Análise de experimento em parcela subdividida
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "gridExtra", "doBy", "multcomp",
"reshape", "plyr", "nlme", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra gridExtra doBy multcomp reshape
## TRUE TRUE TRUE TRUE TRUE TRUE
## plyr nlme wzRfun
## TRUE TRUE TRUE
sessionInfo()
## R version 3.1.1 (2014-07-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.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] splines grid stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.2 nlme_3.1-117 plyr_1.8.1 reshape_0.8.5
## [5] multcomp_1.3-3 TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-10
## [9] MASS_7.3-34 survival_2.37-7 gridExtra_0.9.1 latticeExtra_0.6-26
## [13] RColorBrewer_1.0-5 lattice_0.20-29 knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 lme4_1.1-6 Matrix_1.1-4
## [5] methods_3.1.1 minqa_1.2.3 Rcpp_0.11.2 RcppEigen_0.3.2.0.2
## [9] sandwich_2.3-0 stringr_0.6.2 tools_3.1.1 zoo_1.7-10
trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Leitura dos dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/soja_planta_tcc.txt",
header=TRUE, sep="\t")
str(da)
## 'data.frame': 40 obs. of 9 variables:
## $ sistema : Factor w/ 2 levels "convencional",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ adubpk : int 40 40 40 40 60 60 60 60 90 90 ...
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ alt15.12: num 50.6 48.3 51.6 43.3 49 50.3 51 49 48.6 56.3 ...
## $ alt05.01: num 69.6 59.6 68.3 60 63 69.3 64.3 67 58 61.6 ...
## $ alt26.01: num 76.6 74.6 70.3 68 68.3 74 67 64.6 66.3 67.6 ...
## $ alt1vag : num 12.6 14 14 15.3 12.3 13.6 15.3 12.3 15 14 ...
## $ prod : num 2044 2383 3137 2889 2308 ...
## $ p100 : num 32 30.3 33.7 29 29.4 30.9 30.3 31.3 32 30.8 ...
##-----------------------------------------------------------------------------
## Visualizar.
xyplot(prod~adubpk|bloco, groups=sistema, data=da, type="o",
auto.key=TRUE)
##-----------------------------------------------------------------------------
## Análise considerando o modelo para experimentos em parcela
## subdividida sendo o sistema a parcela e a adubação a subparcela.
## Efeito categórico para adub para começar.
da$adub <- factor(da$adubpk)
## Níveis das parcelas.
da$parcela <- with(da, interaction(bloco, sistema, drop=TRUE))
## Adub categórico.
m0 <- lme(prod~bloco+sistema*adub, random=~1|parcela, data=da,
method="ML")
## Adub linear.
m1 <- update(m0, prod~bloco+sistema*adubpk)
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
b <- unlist(ranef(m0))
grid.arrange(nrow=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r),
qqmath(b))
##-----------------------------------------------------------------------------
## Quadro para o teste dos efeitos fixos (Wald).
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 24 1360.8 <.0001
## bloco 3 3 1.2 0.4383
## sistema 1 3 0.5 0.5250
## adub 4 24 2.9 0.0414
## sistema:adub 4 24 0.9 0.5045
anova(m1, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 9 603.0 618.2 -292.5
## m0 2 15 611.1 636.4 -290.5 1 vs 2 3.963 0.6817
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 30 1663.3 <.0001
## bloco 3 3 1.5 0.3766
## sistema 1 3 0.6 0.4857
## adubpk 1 30 11.9 0.0017
## sistema:adubpk 1 30 1.4 0.2464
## Teste de razão de verossimilhanças para o efeito de sistema.
m2 <- update(m0, prod~bloco+adubpk)
anova(m2, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 7 601.4 613.2 -293.7
## m0 2 15 611.1 636.4 -290.5 1 vs 2 6.341 0.6091
##-----------------------------------------------------------------------------
## Obter os valores preditos com bandas de confiança.
pred <- with(da, expand.grid(bloco=levels(bloco),
sistema=levels(sistema),
adubpk=seq(40,150,by=10)))
str(pred)
## 'data.frame': 96 obs. of 3 variables:
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ sistema: Factor w/ 2 levels "convencional",..: 1 1 1 1 2 2 2 2 1 1 ...
## $ adubpk : num 40 40 40 40 40 40 40 40 50 50 ...
## - attr(*, "out.attrs")=List of 2
## ..$ dim : Named int 4 2 12
## .. ..- attr(*, "names")= chr "bloco" "sistema" "adubpk"
## ..$ dimnames:List of 3
## .. ..$ bloco : chr "bloco=I" "bloco=II" "bloco=III" "bloco=IV"
## .. ..$ sistema: chr "sistema=convencional" "sistema=direto"
## .. ..$ adubpk : chr "adubpk= 40" "adubpk= 50" "adubpk= 60" "adubpk= 70" ...
## Matriz para predição com todos os blocos.
X <- model.matrix(formula(m1)[-2], data=pred)
head(X)
## (Intercept) blocoII blocoIII blocoIV sistemadireto adubpk sistemadireto:adubpk
## 1 1 0 0 0 0 40 0
## 2 1 1 0 0 0 40 0
## 3 1 0 1 0 0 40 0
## 4 1 0 0 1 0 40 0
## 5 1 0 0 0 1 40 40
## 6 1 1 0 0 1 40 40
## Para obter o efeito médio dos blocos
g <- aggregate(X~sistema+adubpk, data=pred, mean)
head(g)
## sistema adubpk (Intercept) blocoII blocoIII blocoIV sistemadireto adubpk
## 1 convencional 40 1 0.25 0.25 0.25 0 40
## 2 direto 40 1 0.25 0.25 0.25 1 40
## 3 convencional 50 1 0.25 0.25 0.25 0 50
## 4 direto 50 1 0.25 0.25 0.25 1 50
## 5 convencional 60 1 0.25 0.25 0.25 0 60
## 6 direto 60 1 0.25 0.25 0.25 1 60
## sistemadireto:adubpk
## 1 0
## 2 40
## 3 0
## 4 50
## 5 0
## 6 60
str(g)
## 'data.frame': 24 obs. of 9 variables:
## $ sistema : Factor w/ 2 levels "convencional",..: 1 2 1 2 1 2 1 2 1 2 ...
## $ adubpk : num 40 40 50 50 60 60 70 70 80 80 ...
## $ (Intercept) : num 1 1 1 1 1 1 1 1 1 1 ...
## $ blocoII : num 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ blocoIII : num 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ blocoIV : num 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ sistemadireto : num 0 1 0 1 0 1 0 1 0 1 ...
## $ adubpk : num 40 40 50 50 60 60 70 70 80 80 ...
## $ sistemadireto:adubpk: num 0 40 0 50 0 60 0 70 0 80 ...
w <- 1:2
pred <- g[,w]
X <- as.matrix(g[,-w])
## Verifica se os nomes correspondem.
colnames(X)==names(fixef(m1))
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## Obtém os erros padrões e limites dos intervalos de confiança.
U <- chol(vcov(m1))
pred$se <- sqrt(apply(X%*%t(U), 1, function(x) sum(x^2)))
zval <- qnorm(p=c(lwr=0.025, fit=0.5, upr=0.975))
me <- outer(pred$se, zval, "*")
fit <- crossprod(t(X), fixef(m1))
pred <- cbind(pred, sweep(me, 1, fit, "+"))
str(pred)
## 'data.frame': 24 obs. of 6 variables:
## $ sistema: Factor w/ 2 levels "convencional",..: 1 2 1 2 1 2 1 2 1 2 ...
## $ adubpk : num 40 40 50 50 60 60 70 70 80 80 ...
## $ se : num 138 138 124 124 111 ...
## $ lwr : num 2400 2323 2462 2422 2522 ...
## $ fit : num 2670 2593 2705 2665 2740 ...
## $ upr : num 2940 2864 2948 2908 2959 ...
##-----------------------------------------------------------------------------
## Representação.
sub <-
"Retas não diferem entre si pelo teste de razão de verossimilhanças (5%)."
xyplot(fit~adubpk, groups=sistema, data=pred, type="l",
auto.key=list(points=FALSE, lines=TRUE,
title="Sistema de cultivo", cex.title=1.1,
corner=c(0.05,0.95)),
xlab=expression("Adubação com PK"~(kg~ha^{-1})),
ylab=expression("Produtividade"~(kg~ha^{-1})),
sub=list(sub, font=3, cex=0.8),
prepanel=prepanel.cbH, cty="bands",
ly=pred$lwr, uy=pred$upr, alpha=0.1, fill=1,
panel.groups=panel.cbH,
panel=panel.superpose)
