Exemplo de análise
Experimento em arranjo fatorial com dois tratamentos adicionais: modelo gaussiano com variável transformada
Definições da sessão
##-----------------------------------------------------------------------------
## Pacotes necessários.
pkg <- c("latticeExtra", "plyr", "doBy", "multcomp")
sapply(pkg, require, character.only=TRUE)
##-----------------------------------------------------------------------------
## Definições gráficas.
## Definições de preenchimentos, linhas, etc.
mycol <- brewer.pal(6, "Set1")
ps <- list(
box.rectangle=list(col=1, fill=c("gray70")),
box.umbrella=list(col=1, lty=1),
dot.symbol=list(col=1, pch=19),
dot.line=list(col="gray50", lty=3),
plot.symbol=list(col=1, cex=1.2),
plot.line=list(col=1),
plot.polygon=list(col="gray95"),
superpose.line=list(col=mycol),
superpose.symbol=list(col=mycol, pch=c(1:5)),
superpose.polygon=list(col=mycol),
strip.background=list(col=c("gray80","gray50")))
trellis.par.set(ps)
## show.settings()
## Remove objetos da mémoria para começar vazia.
rm(list=ls())Carregar os dados
##-----------------------------------------------------------------------------
## Endereço.
url <- "http://www.leg.ufpr.br/~walmes/data/zimmermann_fusarium.txt"
if(Sys.info()["user"]=="walmes"){
url <- "~/Dropbox/CursoR/DADOS/zimmermann_fusarium.txt"
}
## Leitura e inspeção.
da <- read.table(file=url, header=TRUE, sep="\t")
names(da) <- substr(names(da), 1, 4)
str(da)## 'data.frame': 55 obs. of 4 variables:
## $ trat: Factor w/ 11 levels "benlate antes",..: 3 1 2 6 4 5 9 7 8 10 ...
## $ fung: Factor w/ 5 levels "benlate","captam",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ apli: Factor w/ 5 levels "antes","misturado",..: 4 1 2 4 1 2 4 1 2 3 ...
## $ infe: int 5 1 3 7 3 6 0 0 0 6 ...## Primeiras linhas.
head(da)## trat fung apli infe
## 1 benlate puro benlate puro 5
## 2 benlate antes benlate antes 1
## 3 benlate misturado benlate misturado 3
## 4 captam puro captam puro 7
## 5 captam antes captam antes 3
## 6 captam misturado captam misturado 6Análise exploratória preliminar
##-----------------------------------------------------------------------------
## Frequência.
cbind(xtabs(~trat, data=da))## [,1]
## benlate antes 5
## benlate misturado 5
## benlate puro 5
## captam antes 5
## captam misturado 5
## captam puro 5
## derosal antes 5
## derosal misturado 5
## derosal puro 5
## polimero 5
## testemunha 5cbind(xtabs(~fung+apli, data=da))## antes misturado polimero puro testemunha
## benlate 5 5 0 5 0
## captam 5 5 0 5 0
## derosal 5 5 0 5 0
## polimero 0 0 5 0 0
## testemunha 0 0 0 0 5##----------------------------------------------------------------------
## IMPORTANTE: não é obrigatório mas conveniente que as testemunhas
## sejam os últimos níveis.
da$apli <- factor(da$apli,
levels=levels(da$apli)[c(1,4,2,3,5)])
levels(da$apli)## [1] "antes" "puro" "misturado" "polimero" "testemunha"levels(da$fung) ## Em ordem (por sorte).## [1] "benlate" "captam" "derosal" "polimero" "testemunha"levels(da$trat) ## Em ordem (por sorte).## [1] "benlate antes" "benlate misturado" "benlate puro"
## [4] "captam antes" "captam misturado" "captam puro"
## [7] "derosal antes" "derosal misturado" "derosal puro"
## [10] "polimero" "testemunha"##----------------------------------------------------------------------
## Gráficos.
xyplot(infe/40~fung, groups=apli, data=da, jitter.x=TRUE,
ylab="Proporção amostral de sementes infectadas (total de 40)",
xlab="Fungicida usado no controle", auto.key=list(columns=2))
Análise dos dados
##-----------------------------------------------------------------------------
## Especificação do modelo.
## Será análisado como modelo gaussiano para comparações com o GLM
## binomial.
## O experimento foi instalado em DIC, assim, motivado pelo desenho e
## termos em estudo, o modelo tem os seguintes efeitos.
##
## \eta_{ij} = \mu+\alpha_{i}+\tau_{j}+\gamma_{ij}, se I
## \eta_{k} = \mu+\delta_{k}, se II
##
## Em que I representa a porção fatorial (com dois efeitos principais e
## interação) e II a porção adicional (as testemunhas).
##----------------------------------------------------------------------
## Contrastes.
K <- C(da$fung, contr=contr.helmert)
K <- C(da$fung, contr=contr.treatment)
K <- attr(K, "contrasts")
## O que o constraste representa?
MASS::fractions(solve(cbind(1, K)))## benlate captam derosal polimero testemunha
## 1 0 0 0 0
## 2 -1 1 0 0 0
## 3 -1 0 1 0 0
## 4 -1 0 0 1 0
## 5 -1 0 0 0 1## Os coeficientes podem ser alterados, caso as hipóteses sejam outras.
## Por exemplo, pode se usar os seguintes contrastes.
U <- rbind(c(1,0,0,0,0), ## Média no primeiro nível.
c(-1,1,0,0,0), ## Contraste entre 1 e 2.
c(-1,0,1,0,0), ## Contraste entre 2 e 3.
c(1,1,1,-3,0), ## {1,2,3} vs 1ra testemunha.
c(1,1,1,0,-3)) ## {1,2,3} vs 2da testemunha
U## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] -1 1 0 0 0
## [3,] -1 0 1 0 0
## [4,] 1 1 1 -3 0
## [5,] 1 1 1 0 -3K <- solve(U)[,-1]
MASS::fractions(K)## [,1] [,2] [,3] [,4]
## [1,] 0 0 0 0
## [2,] 1 0 0 0
## [3,] 0 1 0 0
## [4,] 1/3 1/3 -1/3 0
## [5,] 1/3 1/3 0 -1/3## Atribuindo esses contrastes.
contrasts(da$fung) <- K
contrasts(da$apli) <- K
##----------------------------------------------------------------------
## Ajuste do modelo.
## ctr <- list(fung=contr.helmert, apli=contr.helmert)
## ctr <- list(fung=contr.treatment, apli=contr.treatment)
m0 <- lm(infe/40~fung*apli,
## contrasts=ctr,
data=da)
## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Com transformação justificada para estabilização da variância para
## variáveis de proporção em que se assume binomial.
m0 <- update(m0, asin(sqrt(.))~.)
## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de análise de deviance.
anova(m0, test="F")## Analysis of Variance Table
##
## Response: asin(sqrt(infe/40))
## Df Sum Sq Mean Sq F value Pr(>F)
## fung 4 1.51473 0.37868 30.7953 3.031e-12 ***
## apli 2 0.11973 0.05986 4.8682 0.012304 *
## fung:apli 4 0.19038 0.04760 3.8706 0.008855 **
## Residuals 44 0.54106 0.01230
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## NOTE: os graus de liberdade não são os corretos. Os tratamentros
## adicionais são parâmetros dentro do vetor de parâmetros do primeiro
## fator declarado. No segundo, por confundimento completo, se tornam
## não estimáveis e por isso não alteram os graus de liberdade.O incoveniente do fatorial incompleto
##----------------------------------------------------------------------
## Efeitos não estimáveis.
coef(m0)## (Intercept) fung1 fung2 fung3 fung4
## 0.063512086 0.064242442 -0.031756043 -0.900310438 -1.027930173
## apli1 apli2 apli3 apli4 fung1:apli1
## 0.008761339 0.109345123 NA NA 0.300480902
## fung2:apli1 fung3:apli1 fung4:apli1 fung1:apli2 fung2:apli2
## 0.004585299 NA NA 0.138143533 -0.141101166
## fung3:apli2 fung4:apli2 fung1:apli3 fung2:apli3 fung3:apli3
## NA NA NA NA NA
## fung4:apli3 fung1:apli4 fung2:apli4 fung3:apli4 fung4:apli4
## NA NA NA NA NA## Dimensões diferentes no vetor de estimativas e matrix de
## covariância das estimativas.
c(length(coef(m0)), nrow(vcov(m0)))## [1] 25 11##----------------------------------------------------------------------
## TRUQUE: Obter a matriz do modelo, remover as colunas que correspondem
## aos efeitos não estimáveis e ajudar o modelo passando uma matriz.
X <- model.matrix(m0)
a <- attr(X, "assign") ## Será usado para decompor a deviance.
i <- is.na(coef(m0))
levelplot(t(X[nrow(X):1,]), scales=list(x=list(rot=90)))+
layer(panel.abline(v=which(i), lty=2))
X <- X[,!i]
a <- a[!i]
levelplot(t(X[nrow(X):1,]), scales=list(x=list(rot=90)))
##----------------------------------------------------------------------
## Modelo declarado por matriz, apenas colunas de efeitos estimáveis.
m1 <- update(m0, .~0+X)
## Iguais deviances pois são o mesmo modelo.
c(deviance(m0), deviance(m1))## [1] 0.5410595 0.5410595c(length(coef(m1)), nrow(vcov(m1)))## [1] 11 11Decomposição da soma de quadrados ou deviance
a## [1] 0 1 1 1 1 2 2 3 3 3 3a[4:5] <- max(a)+1:2
L <- diag(length(a))
Ls <- by(data=L, INDICES=a, FUN=as.matrix)
Ls <- Ls[-1] ## Remove intercepto.
names(Ls) <- c("fung", "apli", "fung:apli",
"fung vs test1", "fung vs test2")
Ls <- Ls[c(1,4,5,2,3)]
names(Ls)## [1] "fung" "fung vs test1" "fung vs test2" "apli"
## [5] "fung:apli"lapply(Ls,
function(L){
summary(glht(m1, linfct=L), test=Ftest())
})## $fung
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 2 == 0 0.06424
## 3 == 0 -0.03176
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 0.9726 2 44 0.3861
##
## $`fung vs test1`
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 4 == 0 -0.9003
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 32.96 1 44 8.099e-07
##
## $`fung vs test2`
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 5 == 0 -1.028
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 42.96 1 44 5.154e-08
##
## $apli
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 6 == 0 0.008761
## 7 == 0 0.109345
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 1.501 2 44 0.2341
##
## $`fung:apli`
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 8 == 0 0.300481
## 9 == 0 0.004585
## 10 == 0 0.138144
## 11 == 0 -0.141101
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 3.871 4 44 0.008855## anova(m0, test="F") ## Sequencial.
drop1(m0, test="F", scope=.~.) ## Marginal.## Single term deletions
##
## Model:
## asin(sqrt(infe/40)) ~ fung + apli + fung:apli
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 0.54106 -232.19
## fung 2 0.023919 0.56498 -233.81 0.9726 0.386094
## apli 2 0.036917 0.57798 -232.56 1.5011 0.234081
## fung:apli 4 0.190383 0.73144 -223.60 3.8706 0.008855 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## O p-valor de F do Wald é igual ao p- valor de F da drop1(). São
## iguais porque o modelo é de classe "lm" e nessa classe a LL é uma
## função quadrática, então a estatística F de LL e F de Wald são
## identicas.
## MAS CUIDADO! Essas hipóteses são marginais e algumas delas só fazem
## sentido em casos muito específicos. Por exemplo, não existe sentido
## em testar o efeito de um efeito principal ajustado para interações
## que o contenham. Também não faz sentido testar a média do fatorial
## contra uma testemunha se a interação for signifcativa.Obter as estimativas de proporção de infecção
##----------------------------------------------------------------------
## Criar um modelo de mentira que seja um fatorial completo.
## Atribui indices incrementais as linhas.
rownames(da) <- NULL
dau <- unique(subset(da, select=c(fung, apli)))
Xu <- X[rownames(dau),]
Xu## (Intercept) fung1 fung2 fung3 fung4 apli1
## 1 1 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000
## 2 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 3 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 4 1 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000
## 5 1 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 6 1 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 7 1 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000
## 8 1 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000
## 9 1 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000
## 10 1 0.3333333 0.3333333 -0.3333333 0.0000000 0.3333333
## 11 1 0.3333333 0.3333333 0.0000000 -0.3333333 0.3333333
## apli2 fung1:apli1 fung2:apli1 fung1:apli2 fung2:apli2
## 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 3 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 4 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000
## 5 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 6 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
## 7 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000
## 8 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 9 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000
## 10 0.3333333 0.1111111 0.1111111 0.1111111 0.1111111
## 11 0.3333333 0.1111111 0.1111111 0.1111111 0.1111111## 3x3+2 = 11 celas experimentais. Se fosse um experimento em blocos,
## não haveria linhas únicas ao considerar o termo bloco, então uma rota
## diferente seria necessária para aplicar a média no efeito dos blocos.
ic <- confint(glht(m1, linfct=Xu), calpha=univariate_calpha())$confint
ic## Estimate lwr upr
## 1 7.227342e-02 -0.02767258 0.17221943
## 2 6.351209e-02 -0.03643392 0.16345809
## 3 1.728572e-01 0.07291120 0.27280321
## 4 4.369968e-01 0.33705076 0.53694277
## 5 1.277545e-01 0.02780852 0.22770053
## 6 3.752432e-01 0.27529718 0.47518919
## 7 4.510268e-02 -0.05484332 0.14504869
## 8 3.175604e-02 -0.06818996 0.13170205
## 9 -8.326673e-17 -0.09994601 0.09994601
## 10 4.473808e-01 0.34743480 0.54732681
## 11 4.899207e-01 0.38997471 0.58986672
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.015368pred <- cbind(dau, ic)
str(pred)## 'data.frame': 11 obs. of 5 variables:
## $ fung : Factor w/ 5 levels "benlate","captam",..: 1 1 1 2 2 2 3 3 3 4 ...
## ..- attr(*, "contrasts")= num [1:5, 1:4] 0 1 0 0.333 0.333 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr "benlate" "captam" "derosal" "polimero" ...
## .. .. ..$ : NULL
## $ apli : Factor w/ 5 levels "antes","puro",..: 2 1 3 2 1 3 2 1 3 4 ...
## ..- attr(*, "contrasts")= num [1:5, 1:4] 0 1 0 0.333 0.333 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr "antes" "puro" "misturado" "polimero" ...
## .. .. ..$ : NULL
## $ Estimate: num 0.0723 0.0635 0.1729 0.437 0.1278 ...
## $ lwr : num -0.0277 -0.0364 0.0729 0.3371 0.0278 ...
## $ upr : num 0.172 0.163 0.273 0.537 0.228 ...##----------------------------------------------------------------------
## Gráfico.
segplot(fung:apli~lwr+upr, data=pred, horizontal = FALSE,
ylab="asin(sqrt(Proprorção de infecção))",
xlab="Tratamento", scales=list(x=list(rot=90)),
centers=Estimate, draw.bands=FALSE)
##----------------------------------------------------------------------
## Gráfico aprimorado.
require(wzRfun)## Loading required package: wzRfunpred$desloc <- 0.05*c(0,-1,1,0,0)[match(pred$apli,
levels(da$apli))]
xyplot(asin(sqrt(infe/40))~fung, groups=apli, data=da,
jitter.x=TRUE, amount=0.05,
xlab="Fungicida",
ylab=expression(asin(sqrt("Proporção de infecção"))),
prepanel=prepanel.cbH,
ly=pred$lwr, uy=pred$upr,
auto.key=list(
columns=2, type="o", divide=1,
title="Aplicação", cex.title=1.2,
lines=TRUE, points=FALSE))+
as.layer(xyplot(Estimate~fung, groups=apli,
data=pred, type=c("p","l"),
pch=15,
ly=pred$lwr, uy=pred$upr,
desloc=pred$desloc,
cty="bars",
panel.groups=panel.cbH,
panel=panel.superpose))
Comparações múltiplas
##----------------------------------------------------------------------
## Testar os tratamentos entre si.
Xu## (Intercept) fung1 fung2 fung3 fung4 apli1
## 1 1 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000
## 2 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 3 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 4 1 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000
## 5 1 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 6 1 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 7 1 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000
## 8 1 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000
## 9 1 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000
## 10 1 0.3333333 0.3333333 -0.3333333 0.0000000 0.3333333
## 11 1 0.3333333 0.3333333 0.0000000 -0.3333333 0.3333333
## apli2 fung1:apli1 fung2:apli1 fung1:apli2 fung2:apli2
## 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 3 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 4 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000
## 5 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 6 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
## 7 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000
## 8 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 9 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000
## 10 0.3333333 0.1111111 0.1111111 0.1111111 0.1111111
## 11 0.3333333 0.1111111 0.1111111 0.1111111 0.1111111dau## fung apli
## 1 benlate puro
## 2 benlate antes
## 3 benlate misturado
## 4 captam puro
## 5 captam antes
## 6 captam misturado
## 7 derosal puro
## 8 derosal antes
## 9 derosal misturado
## 10 polimero polimero
## 11 testemunha testemunhaXm <- Xu
rownames(Xm) <- with(dau, paste(fung, apli, sep=":"))
d <- combn(nrow(Xm), 2)
ncol(d)## [1] 55## Matriz com as diferenças entre linhas da anterior, ou seja, os
## contrastes.
Xc <- sapply(as.data.frame(d),
function(l){
K <- list(Xm[l[1],]-Xm[l[2],])
names(K) <- paste(rownames(Xm)[l], collapse="-")
return(K)
})
Xc <- do.call(rbind, Xc)
## Xc
## summary(glht(m1, linfct=Xc), test=adjusted(type="fdr"))
##----------------------------------------------------------------------
## Para facilitar, dado que existe interação, declarar o modelo de um
## fator é equivalente.
m2 <- update(m0, .~trat)
anova(m2)## Analysis of Variance Table
##
## Response: asin(sqrt(infe/40))
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 10 1.82484 0.182484 14.84 4.463e-11 ***
## Residuals 44 0.54106 0.012297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1c(deviance(m2), deviance(m1))## [1] 0.5410595 0.5410595g0 <- summary(glht(m2, linfct=mcp(trat="Tukey")),
test=adjusted(type="fdr"))
## cld(g0)
## str(cld(g0))
cbind(cld(g0, decreasing=TRUE)$mcletters$Letters)## [,1]
## benlate antes "bc"
## benlate misturado "b"
## benlate puro "bc"
## captam antes "bc"
## captam misturado "a"
## captam puro "a"
## derosal antes "bc"
## derosal misturado "c"
## derosal puro "bc"
## polimero "a"
## testemunha "a"Informações da sessão
##-----------------------------------------------------------------------------
## Versões dos pacotes e data do documento.
sessionInfo()## R version 3.2.0 (2015-04-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods stats graphics grDevices utils datasets
## [7] base
##
## other attached packages:
## [1] wzRfun_0.4 multcomp_1.4-0 TH.data_1.0-6
## [4] mvtnorm_1.0-2 doBy_4.5-13 survival_2.38-2
## [7] plyr_1.8.3 latticeExtra_0.6-26 lattice_0.20-31
## [10] RColorBrewer_1.1-2 rmarkdown_0.7 knitr_1.10.5
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.11.6 magrittr_1.5 splines_3.2.0
## [4] MASS_7.3-41 stringr_1.0.0 tools_3.2.0
## [7] grid_3.2.0 htmltools_0.2.6 yaml_2.1.13
## [10] digest_0.6.8 Matrix_1.2-1 formatR_1.2
## [13] codetools_0.2-11 evaluate_0.7 sandwich_2.3-3
## [16] stringi_0.4-1 zoo_1.7-12Sys.time()## [1] "2015-06-23 10:51:04 BRT"