#===============================================================================
# Título: Planejamento de experimento
# Curso : Extensão
# Autor : José Cláudio Faria/UESC/DCET
# Data  : 30/08/2021 21:03:31
# Objetivos:
#===============================================================================
# - Conhecer os recursos básicos do R para auxílio no planejamento e análise
#   de experimentos;
#
# - Exemplo de planejamento de experimento (3x3) no esquema experimental em
#   parcela subdivididas no espaço no delineamento em blocos casualizados:
#   espécie madeirável (parcela) e café conilon (subparcela).
#===============================================================================
# Pacotes necessários:
# - agricolae: planejamento do experimento
# - xlsx     : exportação do planejamento
# - TukeyC   : comparação de média
#===============================================================================

#. Planejamento
library(agricolae)

#... Parcela
EM <- c("Testemunha",
        "Guanandi",
        "Jequitibá")

#... Subparcela
CAFE <- paste('c',
               0:2, 
               sep='')

#... Este planejamento foi feito em 2014. Opção melhor atualmente é usar o pacote "gexp"
res <- design.split(EM,
                    CAFE,
                    r=5,
                    design='rcbd',
                    serie=1,
                    seed=29051956)
str(res)

#. Exportando
plan <- res$book

names(plan) <- c('Parcela',
                 'Sub-parcela',
                 'Bloco',
                 'EM',
                 'CAFE')

#.. Para arquivo texto
write.table(plan,
            "Gilton_cafe.txt",
            row.names=FALSE,
            sep="\t")


#.. Para uma planilha
#install.packages('xlsx')
library(xlsx)
write.xlsx(plan,
           file='Gilton.xlsx',
           sheet='Experimento')


#. Simulação da análise
dad <- plan
dad$y <- rnorm(45, 
               50, 
               5)

dad$CAFE  <- as.factor(dad$CAFE)
str(dad)

#.. ANOVA
av <- aov(y ~ Bloco + EM*CAFE + Error(Bloco/EM),
          data=dad)

summary(av)

#.. TCMM
library(TukeyC)
cv(av)

#... Efeito principal: EM
tk_EM <- TukeyC(av,
                which='EM',
                error='Bloco:EM',
                sig.level=.05)

summary(tk_EM)

#... Efeito principal: CAFE
tk_CAFE <- TukeyC(av,
                  which='CAFE',
                  error='Within',
                  sig.level=0.05)
summary(tk_CAFE)


#.. Testando CAFE (subparcela) dentro de cada nível de EM (parcela)
# Verificando os níveis
attr(av,
     'xlevels')

#... Guanandi/CAFE
tk_Gua_CAFE <- TukeyC(av,
                      which='EM:CAFE',
                      error='Within',
                      fl1=1)  # 1: Guanandi
                              # 2: Jequitibá
                              # 3: Testemunha
summary(tk_Gua_CAFE)

#... Jequitibá/CAFE
tk_Jeq_CAFE <- TukeyC(av,
                      which='EM:CAFE',
                      error='Within',
                      fl1=2)  # 1: Guanandi
                              # 2: Jequitibá
                              # 3: Testemunha
summary(tk_Jeq_CAFE)

#... Testemunha/CAFE
tk_Tes_CAFE <- TukeyC(av,
                      which='EM:CAFE',
                      error='Within',
                      fl1=3)  # 1: Guanandi
                              # 2: Jequitibá
                              # 3: Testemunha
summary(tk_Tes_CAFE)

#.. Testando EM (parcela) dentro de cada nível de CAFE (subparcela)
#... Necessário combinação de erros (feito automaticamente pelas versões mais novas do pacote)
#... c0/EM
tk_co_EM <- TukeyC(av,
                   which='CAFE:EM',
                   error='Within/Bloco:EM',
                   fl1=1)  # c0
                           # c1
                           # c2
summary(tk_co_EM)

#... c1/EM
tk_c1_EM <- TukeyC(av,
                   which='CAFE:EM',
                   error='Within/Bloco:EM',
                   fl1=2)  # c0
                           # c1
                           # c2
summary(tk_c1_EM)

#... c2/EM
tk_c2_EM <- TukeyC(av,
                   which='CAFE:EM',
                   error='Within/Bloco:EM',
                   fl1=3)  # c0
                           # c1
                           # c2
summary(tk_c2_EM)