Teste de comparação de médias múltiplas

## Objetivos:
## Exemplificar, através do R, as análises e determinações numéricas 
## referentes a análise de variância sob DIC e testes de comparação 
## de médias usando pacotes alternativos

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

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

##===============================
## Dados
##===============================
tra <- gl(4,
          6,
          labels=LETTERS[1:4])

r   <- rep(1:6, 4)

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)

dad <- data.frame(tra,
                  r,
                  pro)

rm(tra, r, pro)   # Removendo objetos não mais necessários


##===============================
## Análise de variância
##===============================
is.data.frame(dad)
is.factor(dad$tra)  # Necessário ser fator para fazer ANOVA
is.numeric(dad$pro) # Necessário ser numerico 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)


##===============================
## Teste de comparação de médias
##===============================
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Usando a função TukeyHSD do pacote stats
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
tkh <- TukeyHSD(av, 
               'tra')
tkh
plot(tkh)

##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Usando o pacote multicomp
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
id.tra <- with(dad,
               tapply(pro,
                      tra,
                      mean))
## Também pode ser como abaixo
# id.tra <- rep(6, 4);  names(id.tra) <- LETTERS[1:4]

library(multcomp)

## contrMat (multcomp): Gera matriz de comparações
mc <- contrMat(id.tra,
               type='Tukey')
mc

## glht (multcomp): Compara
tkm <- summary(glht(av,
                    linfct=mcp(tra=mc)))

tkm
plot(tkm)

##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Usando o pacote agricolae
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
library(agricolae)

## HSD.test (agricolae)
with(dad,
     HSD.test(pro, tra,
              DFerror=df.residual(av),
              MSerror=deviance(av)/df.residual(av),
              alpha=0.05))

##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Usando o pacote TukeyC
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
library(TukeyC)

## TukeyC (TukeyC)
tkc <- TukeyC(av)

summary(tkc)
summary(tkc, complete=FALSE)
plot(tkc)

## Informações adicionais
names(tkc)
tkc
tkc$Result
tkc$Diff_Prob

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

(GLtot <- anova(av)$Df[1] + anova(av)$Df[2])

(SQDtra <- anova(av)$Sum[1])

(GLtra <- anova(av)$Df[1])

(SQDRes <- anova(av)$Sum[2])

(GLres <- anova(av)$Df[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)

## Calcular a 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)

## Coeficiente de variação
QMres <- sum(resid(av)^2) / av$df.residual
cv <- 100 * sqrt(QMres) / mean(dad$pro)
round(cv, 2)