#===============================================================================
#. Título: Teste de comparação de médias múltiplas - TCMM
#. Curso : DEX076 - Metodologia e Estatística Experimental - Agronomia
# Autor : José Cláudio Faria
# Data  : 2024-10-09 12:20:11
#===============================================================================
# Objetivos:
# Exemplificar, através do R, as análises e determinações numéricas do curso
# de MEE referentes a análise de variância sob DIC e testes de comparação de
# médias usando pacotes alternativos
# Testes:
#   -Tukey, Dunnet, SNK, Duncan
#===============================================================================

#. Pacotes necessários:
# gplots     (CRAN)
# multcomp   (CRAN)
# TukeyC     (CRAN)
# agricolae  (CRAN)
# ScottKnott (CRAN)

#.. Usando recursos convencionais
# install.packages(c('gplots',
#                    'multcomp',
#                    'TukeyC',
#                    'agricolae',
#                    'ScottKnott'))  # Instala todos os pacotes necessários do CRAN.

#.. Usando o pacote 'pacman' (CRAN)
# install.packages('pacman')
# library(pacman)
# p_load (gplots
#         multcomp
#         TukeyC
#         agricolae
#         ScottKnott)

#. Dados
#.. Exemplo numérico visto em TCMM
dad <- data.frame(
  tra=gl(4,
         5,
         labels=LETTERS[1:4]),
  r=1:5,
  pro=c(25, 26, 20, 23, 21,
        31, 25, 28, 27, 24,
        22, 26, 28, 25, 29,
        33, 29, 31, 34, 28)
)

#.. Exemplo numérico visto em ANOVA e DIC
#dad <- data.frame(
#  tra=gl(4,
#         6,
#         labels=LETTERS[1:4]),
#  r=1:6,
#  pro=c(58, 49, 51, 56, 50, 48,
#        60, 55, 66, 61, 54, 61,
#        59, 47, 44, 49, 62, 60,
#        45, 33, 34, 48, 42, 44)
#)

#. Análise de variância - ANOVA
str(dad)
is.factor(dad$tra)  # Necessário ser fator para fazer ANOVA
is.numeric(dad$pro) # Necessário ser numérico para fazer ANOVA

library(gplots)

# plotmeans (gplots)
plotmeans(pro ~ tra,
          data=dad,
          mean.labels=FALSE,
          digits=3,
          col='blue',
          connect=FALSE,
          ylab='Produção',
          xlab='Tratamentos')

av <- aov(pro ~ tra,
          data=dad)

summary(av)

# cv(av)    cv: Função disponível na página da disciplina e
#               recentente incorporada ao pacote TukeyC >= 1.1-5


# Uso dos pacotes para comparação de médias
#. TukeyHSD - pacote stats
tkh <- TukeyHSD(av, 
                'tra',
                conf.level=95/100)
tkh
plot(tkh)

#. multicomp (pacote)
# Tuckey, Dunnet, ...
library(multcomp)

#.. Tukey
# glht (multcomp): compara
tk_t <- glht(av,
             linfct=mcp(tra='Tukey'))

summary(tk_t)

# Gráficos
plot(tk_t)

old.par <- par(mar=c(5, 5, 6, 2)) # Altera parametros do dispositivo gráfico: plot especial
# cld (multcomp): discrimina com letras - Funcional apenas para teste de Tukey
plot(cld(tk_t))
plot(cld(tk_t, level=12/100))
par(old.par) # Retorna parâmetros do dispositivo gráfico ao normal

#.. Dunnet
#... Opção I - mais automatizada
(tr <- rep(5, 4))
names(tr) <- LETTERS[1:4]
tr

mc_1 <- contrMat(tr,
                 type='Dunnet',
                 base=1)  # Testemunha é o primeiro tratamento
mc_1

t_d1 <- glht(av,
             linfct=mcp(tra=mc_1))
summary(t_d1)
plot(t_d1)

#... Opção II - mais manual
               # A  B  C  D
mc_2 <- rbind(c(-1, 1, 0, 0),
              c(-1, 0, 1, 0),
              c(-1, 0, 0, 1))
rownames(mc_2) <- c('A - B', 'A - C', 'A - D')
colnames(mc_2) <- LETTERS[1:4]
mc_2

t_d2 <- glht(av,
             linfct=mcp(tra=mc_2))

summary(t_d2)
plot(t_d2)


#. agricolae (pacote)
# Tuckey, SNK, Duncan, ...
library(agricolae)

#.. Tukey
# Versão antiga do pacote
#tk <- with(dad,
#           HSD.test(pro, tra,
#                    DFerror=df.residual(av),
#                    MSerror=deviance(av)/df.residual(av),
#                    alpha=5/100))
tk <- HSD.test(av,
               'tra',
               alpha=5/100)

tk

#.. SNK
# Versão antiga do pacote
#snk <- with(dad,
#            SNK.test(pro, tra,
#                     DFerror=df.residual(av),
#                     MSerror=deviance(av)/df.residual(av),
#                     alpha=5/100))
snk <- SNK.test(av,
                'tra',
                alpha=5/100)

snk

#.. Duncan
# Versão antiga do pacote
#dc <- with(dad,
#           duncan.test(pro, tra,
#                       DFerror=df.residual(av),
#                       MSerror=deviance(av)/df.residual(av),
#                       alpha=5/100))

dc <- duncan.test(av,
                  'tra',
                  alpha=5/100)
dc

#. TukeyC (pacote)
#.. Tukey

library(TukeyC)
cv(av)   # A função TukeyC foi incorporada ao pacote TukeyC >= 1.1-5

tk <- TukeyC(av,
             sig.level=5/100)

summary(tk)

summary(tk,
        complete=FALSE)

plot(tk,
     dispersion='mm')

# Informações adicionais
str(tk)
names(tk)
tk$out$Result
tk$out$Diff_Prob
tk$out$MSD

# Dificuldade de interpretação quando o número de tratamentos é grande
data(CRD2,
     package='TukeyC')

str(CRD2)

#... Análise detalhada: ANOVA seguida de TCMM
av2 <- aov(y ~ x,
           data=CRD2$dfm)

tk2 <- TukeyC(av2)

summary(tk2)

par(mar=c(8, 5, 6, 2))
plot(tk2,
     id.las=2,
     disp='cip',
     xlab='Tratamentos',
     ylab='Médias')

#... Análise detalhada mais rápida
# Usando os vetores (x e y) da lista (CRD2)
par(mar=c(8, 5, 6, 2))
with(CRD2,
     plot(TukeyC(aov(y ~ x)),
          yl=FALSE,
          disp='cip',
          id.las=2,
          xlab='Tratamentos',
          ylab='Médias'))

# Usando o data.frame (dfm) da lista (CRD2)
par(mar=c(8, 5, 6, 2))
with(CRD2$dfm,
     plot(TukeyC(aov(y ~ x)),
          yl=FALSE,
          disp='cip',
          id.las=2,
          xlab='Tratamentos',
          ylab='Médias'))

#. ScottKnott (pacote)
#library(pacman)
#p_load_gh('ivanalaman/ScottKnott')

#.. ScottKnott
library(ScottKnott)

sk <- SK(av)

summary(sk)

plot(sk,
     disp='cip')

# Informações adicionais
names(sk)
sk

# Um contexto mais adequado para o teste SK
data(CRD2,
     package='ScottKnott')

str(CRD2)

#... Análise detalhada: ANOVA seguida de TCMM
av2 <- aov(y ~ x,
           data=CRD2$dfm)

sk2 <- SK(av2)

summary(sk2)

par(mar=c(8, 5, 6, 2))
plot(sk2,
     id.las=2,
     disp='cip',
     xlab='Tratamentos',
     ylab='Médias')

#... Análise detalhada mais rápida
# Usando os vetores (x e y) da lista (CRD2)
par(mar=c(8, 5, 6, 2))
with(CRD2,
     plot(SK(aov(y ~ x)),
          yl=FALSE,
          disp='cip',
          id.las=2,
          xlab='Tratamentos',
          ylab='Médias'))

# Usando o data.frame (dfm) da lista (CRD2)
par(mar=c(8, 5, 6, 2))
with(CRD2$dfm,
     plot(SK(aov(y ~ x)),
          yl=FALSE,
          disp='cip',
          id.las=2,
          xlab='Tratamentos',
          ylab='Médias'))

#. Alternativas e recursos adicionais
# Recuperar algumas informações úteis contidas no objeto 'av'
# Soma de quadrados: total
(SQDtot <- anova(av)$Sum[1] + anova(av)$Sum[2])

# Graus de liberdade: total
(GLtot <- anova(av)$Df[1] + anova(av)$Df[2])

# Soma de quadrados: tra
(SQDtra <- anova(av)$Sum[1])

# Graus de liberdade: tra
(GLtra <- anova(av)$Df[1])

# Soma de quadrados: resíduo
(SQDRes <- anova(av)$Sum[2])

# Graus de liberdade: resíduo
(GLres <- anova(av)$Df[2])

# Quadrado médio: tra
(QMtra <- anova(av)$Mean[1])

# Quadrado médio: resíduo
(QMres <- anova(av)$Mean[2])

# Cálculos alternativos de algumas medidas básicas
# Calcular médias dos tratamentos e armazenar na variável m
(m <- with(dad,
           tapply(pro,
                  tra,
                  mean)))
sort(m,
     d=TRUE)

# Calcular variâncias dos tratamentos e armazenar na variável vs
(vs <- with(dad,
            tapply(pro,
                   tra,
                   var)))
sort(vs,
     d=TRUE)

# SQDtot em relação á média geral
(SQDtot <- with(dad,
                sum((pro - mean(pro))^2)))

# Calcular os graus de liberdade total
(GLtot <- length(dad$pro) - 1)

# Quadrado médio residuo
(QMres <- sum(resid(av)^2) / av$df.residual)

# Coeficiente de variação
(cv.man <- 100 * sqrt(QMres) / mean(dad$pro))

round(cv.man, 2)