Exemplo de análise
Experimento em arranjo fatorial com dois tratamentos adicionais: modelo linear generalizado
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.
## Haja visto que a resposta consiste no número de sementes infectatas
## num total de 40 unidades, será considerado inicialmente um model
## para resposta do tipo binomial. Esse modelo supõe independencia entre
## os ensaios bernoulli e probablidade de infecção constante na mesma
## unidade experimental.
## 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)
## Como função dos tratamentos primeiro (opcional).
m0 <- glm(cbind(infe, 40-infe)~fung*apli,
## contrasts=ctr,
family=binomial,
data=da)
## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Existem celas com 0 de proporção amostral.
## Avaliação da suposição de equidispersão.
cbind(deviance(m0), df.residual(m0))## [,1] [,2]
## [1,] 75.68084 44pchisq(deviance(m0), df=df.residual(m0))## [1] 0.9979064m0 <- update(m0, family=quasibinomial)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de análise de deviance.
anova(m0, test="F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(infe, 40 - infe)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 54 289.012
## fung 4 181.175 50 107.837 25.7929 4.778e-11 ***
## apli 2 25.276 48 82.562 7.1967 0.001977 **
## fung:apli 4 6.881 44 75.681 0.9796 0.428474
## ---
## 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
## -4.5951199 1.4170660 -0.6981850 -11.9912242 -12.6309668
## apli1 apli2 apli3 apli4 fung1:apli1
## 0.9315582 1.6506809 NA NA 0.7636633
## fung2:apli1 fung3:apli1 fung4:apli1 fung1:apli2 fung2:apli2
## -0.2333732 NA NA -0.3300818 -17.7712541
## 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] 75.68084 75.68084c(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] ## 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 1.4171
## 3 == 0 -0.6982
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 1.755 2 44 0.1847
##
## $`fung vs test1`
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 4 == 0 -11.99
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 0.0002008 1 44 0.9888
##
## $`fung vs test2`
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 5 == 0 -12.63
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 0.0002228 1 44 0.9882
##
## $apli
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 6 == 0 0.9316
## 7 == 0 1.6507
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 1.455 2 44 0.2444
##
## $`fung:apli`
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 8 == 0 0.7637
## 9 == 0 -0.2334
## 10 == 0 -0.3301
## 11 == 0 -17.7713
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 0.5103 4 44 0.7284## anova(m0, test="F") ## Sequencial.
drop1(m0, test="F", scope=.~.) ## Marginal.## Single term deletions
##
## Model:
## cbind(infe, 40 - infe) ~ fung * apli
## Df Deviance F value Pr(>F)
## <none> 75.681
## fung 2 83.287 2.2110 0.1216
## apli 2 81.791 1.7763 0.1812
## fung:apli 4 82.562 1.0001 0.4177## O p-valor de F do Wald está próximo do valor de F da drop1(). Seriam
## iguais se o modelo fosse de classe "lm" porque 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.Reduzindo o modelo
##----------------------------------------------------------------------
## Os testes apresentados apontam não haver interação entre fungicida,
## além do mais, indica pouco efeito de fungicidas e tipos de
## aplicação. Só foram rejeitadas hipóteses considerando as
## testemunhas. Sendo assim, um modelo reduzido que considera toda
## porção fatorial como um único nível será ajustado.
## da$man <- da$fung
## levels(da$man) <- c(rep("ben:cap:der", 3), "pol", "tes")
m2 <- update(m0, .~fung+apli)
anova(m2, m0, test="F")## Analysis of Deviance Table
##
## Model 1: cbind(infe, 40 - infe) ~ fung + apli
## Model 2: cbind(infe, 40 - infe) ~ fung * apli
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 48 82.562
## 2 44 75.681 4 6.8807 0.9796 0.4285anova(m2, test="F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(infe, 40 - infe)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 54 289.012
## fung 4 181.175 50 107.837 24.0842 5.697e-11 ***
## apli 2 25.276 48 82.562 6.7199 0.002673 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(m2, test="F", scope=.~.)## Single term deletions
##
## Model:
## cbind(infe, 40 - infe) ~ fung + apli
## Df Deviance F value Pr(>F)
## <none> 82.562
## fung 2 177.295 27.5384 1.081e-08 ***
## apli 2 107.837 7.3474 0.001645 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##----------------------------------------------------------------------
## O modelo m2 ainda contém NA's. Proceder da mesma forma que
## anteriormente.
coef(m2)## (Intercept) fung1 fung2 fung3 fung4
## -4.648190 1.566605 -1.764020 -6.980987 -7.620730
## apli1 apli2 apli3 apli4
## 1.505112 1.305856 NA NAX <- model.matrix(m2)
a <- attr(X, "assign")
i <- is.na(coef(m2))
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)))
m3 <- update(m2, .~0+X)
c(deviance(m2), deviance(m3))## [1] 82.56155 82.56155Obter 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
## 1 0.0000000
## 2 0.0000000
## 3 1.0000000
## 4 0.0000000
## 5 0.0000000
## 6 1.0000000
## 7 0.0000000
## 8 0.0000000
## 9 1.0000000
## 10 0.3333333
## 11 0.3333333## 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(m3, linfct=Xu), calpha=univariate_calpha())$confint
ic ## Na escala do preditor linear.## Estimate lwr upr
## 1 -3.143078 -3.875339 -2.4108176
## 2 -4.648190 -5.675736 -3.6206447
## 3 -3.342334 -4.099487 -2.5851813
## 4 -1.576474 -2.047126 -1.1058208
## 5 -3.081585 -3.920045 -2.2431259
## 6 -1.775729 -2.275667 -1.2757919
## 7 -4.907098 -6.493339 -3.3208566
## 8 -6.412210 -8.157039 -4.6673808
## 9 -5.106354 -6.704954 -3.5077534
## 10 -1.450010 -1.934480 -0.9655406
## 11 -1.236763 -1.691901 -0.7816239
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964##----------------------------------------------------------------------
## Passar o inverso da função de ligação.
ic <- m2$family$linkinv(ic)
pred <- 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.04136 0.00949 0.03415 0.1713 0.04387 ...
## $ lwr : num 0.02033 0.00342 0.01631 0.11434 0.01945 ...
## $ upr : num 0.0824 0.0261 0.0701 0.2487 0.0959 ...##----------------------------------------------------------------------
## Gráfico.
segplot(fung:apli~lwr+upr, data=pred, horizontal = FALSE,
ylab="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(infe/40~fung, groups=apli, data=da,
jitter.x=TRUE, amount=0.05,
xlab="Fungicida",
ylab="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
##----------------------------------------------------------------------
## Qual o melhor fungicida?
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
## 1 0.0000000
## 2 0.0000000
## 3 1.0000000
## 4 0.0000000
## 5 0.0000000
## 6 1.0000000
## 7 0.0000000
## 8 0.0000000
## 9 1.0000000
## 10 0.3333333
## 11 0.3333333dau## 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 <- by(Xu, INDICES=dau$fung, FUN=colMeans)
Xm <- do.call(rbind, Xm)
Xm## (Intercept) fung1 fung2 fung3 fung4
## benlate 1 0.0000000 0.0000000 0.0000000 0.0000000
## captam 1 1.0000000 0.0000000 0.0000000 0.0000000
## derosal 1 0.0000000 1.0000000 0.0000000 0.0000000
## polimero 1 0.3333333 0.3333333 -0.3333333 0.0000000
## testemunha 1 0.3333333 0.3333333 0.0000000 -0.3333333
## apli1 apli2
## benlate 0.3333333 0.3333333
## captam 0.3333333 0.3333333
## derosal 0.3333333 0.3333333
## polimero 0.3333333 0.3333333
## testemunha 0.3333333 0.3333333## Só as linhas que interessam para a hipótese.
Xm <- Xm[1:3,]
d <- combn(nrow(Xm), 2)
## Matriz com as diferenças entre linhas da anterior, ou seja, os
## contrastes.
Xc <- sapply(as.data.frame(d), #MARGIN=2, #SIMPLIFY=FALSE,
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## (Intercept) fung1 fung2 fung3 fung4 apli1 apli2
## V1.benlate-captam 0 -1 0 0 0 0 0
## V2.benlate-derosal 0 0 -1 0 0 0 0
## V3.captam-derosal 0 1 -1 0 0 0 0summary(glht(m3, linfct=Xc), test=adjusted(type="single-step"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = cbind(infe, 40 - infe) ~ X - 1, family = quasibinomial,
## data = da)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## V1.benlate-captam == 0 -1.5666 0.3812 -4.110 <0.001 ***
## V2.benlate-derosal == 0 1.7640 0.8632 2.044 0.0937 .
## V3.captam-derosal == 0 3.3306 0.8132 4.096 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)##----------------------------------------------------------------------
## Qual melhor método de aplicação?
Xm <- by(Xu, INDICES=dau$apli, FUN=colMeans)
Xm <- do.call(rbind, Xm)
Xm## (Intercept) fung1 fung2 fung3 fung4
## antes 1 0.3333333 0.3333333 0.0000000 0.0000000
## puro 1 0.3333333 0.3333333 0.0000000 0.0000000
## misturado 1 0.3333333 0.3333333 0.0000000 0.0000000
## polimero 1 0.3333333 0.3333333 -0.3333333 0.0000000
## testemunha 1 0.3333333 0.3333333 0.0000000 -0.3333333
## apli1 apli2
## antes 0.0000000 0.0000000
## puro 1.0000000 0.0000000
## misturado 0.0000000 1.0000000
## polimero 0.3333333 0.3333333
## testemunha 0.3333333 0.3333333## Só as linhas que interessam para a hipótese.
Xm <- Xm[1:3,]
d <- combn(nrow(Xm), 2)
## Matriz com as diferenças entre linhas da anterior, ou seja, os
## contrastes.
Xc <- sapply(as.data.frame(d), #MARGIN=2, #SIMPLIFY=FALSE,
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## (Intercept) fung1 fung2 fung3 fung4 apli1 apli2
## V1.antes-puro 0 0 0 0 0 -1 0
## V2.antes-misturado 0 0 0 0 0 0 -1
## V3.puro-misturado 0 0 0 0 0 1 -1summary(glht(m3, linfct=Xc), test=adjusted(type="single-step"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = cbind(infe, 40 - infe) ~ X - 1, family = quasibinomial,
## data = da)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## V1.antes-puro == 0 -1.5051 0.4769 -3.156 0.00436 **
## V2.antes-misturado == 0 -1.3059 0.4851 -2.692 0.01881 *
## V3.puro-misturado == 0 0.1993 0.3279 0.608 0.81291
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)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:34:42 BRT"