Walmes Zeviani
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "lme4", "plyr",
"aod")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra doBy multcomp lme4 plyr aod
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
source("http://dl.dropboxusercontent.com/u/48140237/apc.R")
source("http://dl.dropboxusercontent.com/u/48140237/equallevels.R")
source("http://dl.dropboxusercontent.com/u/48140237/desdobrglht.R")
## source("http://dl.dropboxusercontent.com/u/48140237/bandas.R")
## source("http://dl.dropboxusercontent.com/u/48140237/segplotby.R")
sessionInfo()
## R version 3.1.0 (2014-04-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=pt_BR.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] aod_1.3 plyr_1.8 lme4_1.0-5 Matrix_1.1-3
## [5] multcomp_1.2-21 mvtnorm_0.9-9994 doBy_4.5-10 MASS_7.3-33
## [9] survival_2.37-7 latticeExtra_0.6-24 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.1 formatR_0.9 grid_3.1.0 minqa_1.2.2 nlme_3.1-110 stringr_0.6.2
## [7] tools_3.1.0
trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Leitura dos dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/farms.txt",
header=TRUE, sep="\t")
str(da)
## 'data.frame': 48 obs. of 4 variables:
## $ farm : int 1 1 1 1 1 1 2 2 2 2 ...
## $ block: int 1 1 2 2 3 3 1 1 2 2 ...
## $ trt : Factor w/ 2 levels "A","B": 1 2 1 2 1 2 1 2 1 2 ...
## $ resp : num 30.9 33.3 30.8 30.9 32.3 ...
##-----------------------------------------------------------------------------
## Editar.
da <- transform(da,
farm=as.factor(farm),
block=as.factor(block))
str(da)
## 'data.frame': 48 obs. of 4 variables:
## $ farm : Factor w/ 8 levels "1","2","3","4",..: 1 1 1 1 1 1 2 2 2 2 ...
## $ block: Factor w/ 3 levels "1","2","3": 1 1 2 2 3 3 1 1 2 2 ...
## $ trt : Factor w/ 2 levels "A","B": 1 2 1 2 1 2 1 2 1 2 ...
## $ resp : num 30.9 33.3 30.8 30.9 32.3 ...
##-----------------------------------------------------------------------------
## Layout.
ftable(farm~block+trt, data=da)
## farm 1 2 3 4 5 6 7 8
## block trt
## 1 A 1 1 1 1 1 1 1 1
## B 1 1 1 1 1 1 1 1
## 2 A 1 1 1 1 1 1 1 1
## B 1 1 1 1 1 1 1 1
## 3 A 1 1 1 1 1 1 1 1
## B 1 1 1 1 1 1 1 1
## ftable(trt~farm+block, data=da)
##-----------------------------------------------------------------------------
## Ver.
xyplot(resp~trt|farm, groups=block, data=da, type="o")
##-----------------------------------------------------------------------------
## Ajuste do modelo que considera apenas o efeito de tratamento como
## fixo, os demais e interações são aleatórios.
## * fixo: trt;
## * aleatório: farm, block dentro de farm e farm interação com trt.
m0 <- lmer(resp~trt+(1|farm)+(1|farm:block)+(1|farm:trt),
data=da, REML=FALSE)
summary(m0)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: resp ~ trt + (1 | farm) + (1 | farm:block) + (1 | farm:trt)
## Data: da
##
## AIC BIC logLik deviance
## 241.4 252.7 -114.7 229.4
##
## Random effects:
## Groups Name Variance Std.Dev.
## farm:block (Intercept) 0.919 0.958
## farm:trt (Intercept) 8.061 2.839
## farm (Intercept) 17.533 4.187
## Residual 1.656 1.287
## Number of obs: 48, groups: farm:block, 24; farm:trt, 16; farm, 8
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 30.211 1.818 16.61
## trtB 0.841 1.467 0.57
##
## Correlation of Fixed Effects:
## (Intr)
## trtB -0.403
## Predição dos efeitos aleatórios.
ranef(m0)
## $`farm:block`
## (Intercept)
## 1:1 -0.06270
## 1:2 -0.69653
## 1:3 0.82623
## 2:1 0.91915
## 2:2 0.02759
## 2:3 -1.06385
## 3:1 0.43810
## 3:2 0.74581
## 3:3 -0.79536
## 4:1 -0.05646
## 4:2 0.07241
## 4:3 0.13553
## 5:1 0.25487
## 5:2 -0.22642
## 5:3 -0.32372
## 6:1 0.72254
## 6:2 -0.09012
## 6:3 -0.78706
## 7:1 0.99247
## 7:2 -0.39615
## 7:3 -0.59603
## 8:1 0.26797
## 8:2 -0.53417
## 8:3 0.22589
##
## $`farm:trt`
## (Intercept)
## 1:A -0.17028
## 1:B 0.75824
## 2:A 2.53631
## 2:B -3.56397
## 3:A 0.09408
## 3:B 3.31561
## 4:A 3.02067
## 4:B -1.69133
## 5:A -1.08595
## 5:B -1.50516
## 6:A -0.14137
## 6:B -1.21573
## 7:A -3.42833
## 7:B 3.43083
## 8:A -0.82514
## 8:B 0.47151
##
## $farm
## (Intercept)
## 1 1.278767
## 2 -2.235080
## 3 7.415796
## 4 2.891205
## 5 -5.635450
## 6 -2.951573
## 7 0.005443
## 8 -0.769108
##
## attr(,"class")
## [1] "ranef.mer"
## Existe interação entre farm e trt? Esse componente de variância é
## diferente de zero?
## Modelo que abandona o termo aleatório farm:trt.
m1 <- lmer(resp~trt+(1|farm)+(1|farm:block), data=da)
## Teste para o termo abandonado.
anova(m0, m1)
## Data: da
## Models:
## m1: resp ~ trt + (1 | farm) + (1 | farm:block)
## m0: resp ~ trt + (1 | farm) + (1 | farm:block) + (1 | farm:trt)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 5 264 274 -127 254
## m0 6 241 253 -115 229 25 1 5.7e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Existe interação farm:trt.
##-----------------------------------------------------------------------------
## Teste para os termos de efeito fixo.
anova(m0)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## trt 1 0.544 0.544 0.33
## Não tem p-valor! É de propósito. Leia:
## https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
## Alternativas (em order de recomendação):
## * Teste de razão de verossimilhanças entre modelos aninhados (um com
## outro sem o termo fixo de interesse, usar REML=FALSE).
## * Teste de Wald (inferência aproximada, pacote aod).
V <- vcov(m0)
b <- fixef(m0)
Terms <- which(attr(m0@pp$X, "assign")==1)
wt <- wald.test(Sigma=V, b=b, Terms=Terms)
wt
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 0.33, df = 1, P(> X2) = 0.57
wt$result
## $chi2
## chi2 df P
## 0.3287 1.0000 0.5664
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
bfb <- unlist(ranef(m0)[1])
bft <- unlist(ranef(m0)[2])
bf <- unlist(ranef(m0)[3])
xyplot(r~f, type=c("p", "smooth"))
xyplot(sqrt(abs(r))~f, type=c("p", "smooth"))
qqmath(r)
qqmath(bfb)
qqmath(bft)
qqmath(bf)
##-----------------------------------------------------------------------------
## Médias ajustadas.
## Formula só da parte de efeito fixo.
f <- nobars(formula(m0))
X <- LSmatrix(lm(f, da), effect=c("trt"))
str(X)
## LSmatrix [1:2, 1:2] 1 1 0 1
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:2] "(Intercept)" "trtB"
## - attr(*, "grid")='data.frame': 2 obs. of 1 variable:
## ..$ trt: chr [1:2] "A" "B"
## - attr(*, "class")= chr [1:2] "LSmatrix" "matrix"
## Estimativas das médias.
summary(glht(m0, linfct=X))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = resp ~ trt + (1 | farm) + (1 | farm:block) + (1 |
## farm:trt), data = da, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 == 0 30.21 1.82 16.6 <1e-10 ***
## 2 == 0 31.05 1.82 17.1 <1e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## Contraste.
summary(glht(m0, linfct=rbind(X[1,]-X[2,])))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = resp ~ trt + (1 | farm) + (1 | farm:block) + (1 |
## farm:trt), data = da, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 == 0 -0.841 1.467 -0.57 0.57
## (Adjusted p values reported -- single-step method)
##-----------------------------------------------------------------------------
## Para desdobrar a interação farm:trt, tratar como efeito fixo.
m0 <- lmer(resp~farm*trt+(1|farm:block),
data=da, REML=FALSE)
summary(m0)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: resp ~ farm * trt + (1 | farm:block)
## Data: da
##
## AIC BIC logLik deviance
## 194.87 228.56 -79.44 158.87
##
## Random effects:
## Groups Name Variance Std.Dev.
## farm:block (Intercept) 0.612 0.783
## Residual 1.104 1.051
## Number of obs: 48, groups: farm:block, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 31.330 0.756 41.4
## farm2 -0.683 1.070 -0.6
## farm3 6.527 1.070 6.1
## farm4 5.050 1.070 4.7
## farm5 -8.013 1.070 -7.5
## farm6 -4.273 1.070 -4.0
## farm7 -4.777 1.070 -4.5
## farm8 -2.783 1.070 -2.6
## trtB 1.833 0.858 2.1
## farm2:trtB -7.510 1.213 -6.2
## farm3:trtB 2.450 1.213 2.0
## farm4:trtB -6.027 1.213 -5.0
## farm5:trtB -1.440 1.213 -1.2
## farm6:trtB -2.140 1.213 -1.8
## farm7:trtB 6.337 1.213 5.2
## farm8:trtB 0.393 1.213 0.3
##
## Correlation of Fixed Effects:
## (Intr) farm2 farm3 farm4 farm5 farm6 farm7 farm8 trtB frm2:B frm3:B frm4:B
## farm2 -0.707
## farm3 -0.707 0.500
## farm4 -0.707 0.500 0.500
## farm5 -0.707 0.500 0.500 0.500
## farm6 -0.707 0.500 0.500 0.500 0.500
## farm7 -0.707 0.500 0.500 0.500 0.500 0.500
## farm8 -0.707 0.500 0.500 0.500 0.500 0.500 0.500
## trtB -0.567 0.401 0.401 0.401 0.401 0.401 0.401 0.401
## farm2:trtB 0.401 -0.567 -0.284 -0.284 -0.284 -0.284 -0.284 -0.284 -0.707
## farm3:trtB 0.401 -0.284 -0.567 -0.284 -0.284 -0.284 -0.284 -0.284 -0.707 0.500
## farm4:trtB 0.401 -0.284 -0.284 -0.567 -0.284 -0.284 -0.284 -0.284 -0.707 0.500 0.500
## farm5:trtB 0.401 -0.284 -0.284 -0.284 -0.567 -0.284 -0.284 -0.284 -0.707 0.500 0.500 0.500
## farm6:trtB 0.401 -0.284 -0.284 -0.284 -0.284 -0.567 -0.284 -0.284 -0.707 0.500 0.500 0.500
## farm7:trtB 0.401 -0.284 -0.284 -0.284 -0.284 -0.284 -0.567 -0.284 -0.707 0.500 0.500 0.500
## farm8:trtB 0.401 -0.284 -0.284 -0.284 -0.284 -0.284 -0.284 -0.567 -0.707 0.500 0.500 0.500
## frm5:B frm6:B frm7:B
## farm2
## farm3
## farm4
## farm5
## farm6
## farm7
## farm8
## trtB
## farm2:trtB
## farm3:trtB
## farm4:trtB
## farm5:trtB
## farm6:trtB 0.500
## farm7:trtB 0.500 0.500
## farm8:trtB 0.500 0.500 0.500
##-----------------------------------------------------------------------------
## Testes de Wald.
V <- vcov(m0)
b <- fixef(m0)
## Termos de efeito fixo.
nobars(formula(m0))
## resp ~ farm * trt
a <- attr(m0@pp$X, "assign"); a
## [1] 0 1 1 1 1 1 1 1 2 3 3 3 3 3 3 3
split(fixef(m0), a)
## $`0`
## (Intercept)
## 31.33
##
## $`1`
## farm2 farm3 farm4 farm5 farm6 farm7 farm8
## -0.6833 6.5267 5.0500 -8.0133 -4.2733 -4.7767 -2.7833
##
## $`2`
## trtB
## 1.833
##
## $`3`
## farm2:trtB farm3:trtB farm4:trtB farm5:trtB farm6:trtB farm7:trtB farm8:trtB
## -7.5100 2.4500 -6.0267 -1.4400 -2.1400 6.3367 0.3933
## Teste para farm.
Terms <- which(a==1)
wald.test(Sigma=V, b=b, Terms=Terms)$result
## $chi2
## chi2 df P
## 299.9 7.0 0.0
## Teste para trt.
Terms <- which(a==2)
wald.test(Sigma=V, b=b, Terms=Terms)$result
## $chi2
## chi2 df P
## 4.56765 1.00000 0.03258
## Teste para farm:trt.
Terms <- which(a==3)
wald.test(Sigma=V, b=b, Terms=Terms)$result
## $chi2
## chi2 df P
## 187.3 7.0 0.0
##-----------------------------------------------------------------------------
## Médias ajustadas.
## Formula só da parte de efeito fixo.
f <- nobars(formula(m0))
X <- LSmatrix(lm(f, da), effect=c("farm", "trt"))
str(X)
## LSmatrix [1:16, 1:16] 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:16] "(Intercept)" "farm2" "farm3" "farm4" ...
## - attr(*, "grid")='data.frame': 16 obs. of 2 variables:
## ..$ farm: chr [1:16] "1" "2" "3" "4" ...
## ..$ trt : chr [1:16] "A" "A" "A" "A" ...
## - attr(*, "class")= chr [1:2] "LSmatrix" "matrix"
grid <- attr(X, "grid")
grid <- equallevels(grid, da)
## Estimativas das médias.
g <- summary(glht(m0, linfct=X), test=adjusted(type="fdr"))
ci <- confint(g, calpha=univariate_calpha())$confint
grid <- cbind(grid, ci)
grid
## farm trt Estimate lwr upr
## 1 1 A 31.33 29.85 32.81
## 2 2 A 30.65 29.16 32.13
## 3 3 A 37.86 36.37 39.34
## 4 4 A 36.38 34.90 37.86
## 5 5 A 23.32 21.83 24.80
## 6 6 A 27.06 25.57 28.54
## 7 7 A 26.55 25.07 28.04
## 8 8 A 28.55 27.06 30.03
## 9 1 B 33.16 31.68 34.65
## 10 2 B 24.97 23.49 26.45
## 11 3 B 42.14 40.66 43.62
## 12 4 B 32.19 30.70 33.67
## 13 5 B 23.71 22.23 25.19
## 14 6 B 26.75 25.27 28.23
## 15 7 B 34.72 33.24 36.21
## 16 8 B 30.77 29.29 32.26
##-----------------------------------------------------------------------------
## Contrastes.
## Parte a matriz de acordo com farm.
L <- by(X, grid$farm, FUN=as.matrix)
L <- lapply(L, "rownames<-", levels(da$trt))
L <- lapply(L, desdobr.glht, m0=m0, focus="trt")
L
## $`1`
## trt estimate lwr upr cld
## 1 A 31.33 29.85 32.81 b
## 2 B 33.16 31.68 34.65 a
##
## $`2`
## trt estimate lwr upr cld
## 1 A 30.65 29.16 32.13 b
## 2 B 24.97 23.49 26.45 a
##
## $`3`
## trt estimate lwr upr cld
## 1 A 37.86 36.37 39.34 b
## 2 B 42.14 40.66 43.62 a
##
## $`4`
## trt estimate lwr upr cld
## 1 A 36.38 34.9 37.86 b
## 2 B 32.19 30.7 33.67 a
##
## $`5`
## trt estimate lwr upr cld
## 1 A 23.32 21.83 24.80 a
## 2 B 23.71 22.23 25.19 a
##
## $`6`
## trt estimate lwr upr cld
## 1 A 27.06 25.57 28.54 a
## 2 B 26.75 25.27 28.23 a
##
## $`7`
## trt estimate lwr upr cld
## 1 A 26.55 25.07 28.04 b
## 2 B 34.72 33.24 36.21 a
##
## $`8`
## trt estimate lwr upr cld
## 1 A 28.55 27.06 30.03 b
## 2 B 30.77 29.29 32.26 a
##-----------------------------------------------------------------------------
## Gráfico.
grid <- ldply(L)
names(grid)[1] <- "farm"
grid <- equallevels(grid, da)
segplot(trt~lwr+upr|farm, centers=estimate,
ylab="Tratamento", xlab="Resposta",
data=grid, draw=FALSE, strip=FALSE, strip.left=TRUE,
layout=c(1,NA), cld=grid$cld,
panel=function(x, y, z, centers, subscripts, cld, ...){
panel.segplot(x, y, z, centers=centers,
subscripts=subscripts, ...)
panel.text(x=centers[subscripts], y=as.numeric(z)[subscripts],
labels=cld[subscripts], pos=4)
})
