Análise de experimento em parcela subsubdividida na profundidade
##-----------------------------------------------------------------------------
## 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/sistema_gesso_solo.txt",
header=TRUE, sep="\t")
## str(da)
##-----------------------------------------------------------------------------
## Editar.
da <- transform(da,
bloco=as.factor(bloco),
gesso=as.factor(gesso),
Prof=as.factor(prof),
parc=interaction(bloco, sistema),
subp=interaction(bloco, sistema, gesso))
str(da)
## 'data.frame': 80 obs. of 21 variables:
## $ sistema: Factor w/ 2 levels "PC","PD": 2 2 2 2 2 2 2 2 2 2 ...
## $ gesso : Factor w/ 2 levels "0","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ prof : int 5 5 5 5 10 10 10 10 15 15 ...
## $ bloco : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
## $ Al : num 0 0 0 0 0 0 0 0 0 0 ...
## $ MO : num 60.4 42.4 51.7 52.4 31.7 39.9 32.3 40.7 28.7 31 ...
## $ pHKCl : num 5.2 5.6 6.4 4.2 6.7 6.2 6.6 3.9 6.2 5.3 ...
## $ pHCaCl : num 5.5 5.7 6.1 4.3 6.2 6.3 6.5 4 6.2 5.8 ...
## $ pHH2O : num 6 6.4 7 5 7.2 7.1 7.4 4.6 7.2 6.4 ...
## $ Ca : num 69.9 77.3 91.4 79.4 89.2 86.9 95.5 90.5 76.6 57.7 ...
## $ Mg : num 32.3 30.3 28.6 30.4 36 32.6 26 31 28.3 19.6 ...
## $ Hal : int 45 29 18 30 16 19 16 17 20 34 ...
## $ S : num 8.2 2.8 3.8 4.9 7.2 2.7 7.2 2.7 6.6 7.8 ...
## $ P : int 19 17 22 16 9 4 8 8 4 2 ...
## $ K : num 4.8 4.7 8 6.9 1.7 2.3 3.3 2 0.9 1.2 ...
## $ SB : num 107 112 128 117 127 ...
## $ CTC : num 152 141 146 147 143 ...
## $ SAT : int 70 79 88 80 89 87 89 88 84 70 ...
## $ Prof : Factor w/ 5 levels "5","10","15",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ parc : Factor w/ 8 levels "1.PC","2.PC",..: 5 6 7 8 5 6 7 8 5 6 ...
## $ subp : Factor w/ 16 levels "1.PC.0","2.PC.0",..: 5 6 7 8 5 6 7 8 5 6 ...
##-----------------------------------------------------------------------------
## Ver.
da$y <- da$P
xyplot(y~prof|sistema, groups=gesso, data=da, type=c("p","a"),
auto.key=TRUE, jitter.x=TRUE)
##-----------------------------------------------------------------------------
## Ajuste do modelo de parcela subdivididas (assume independência entre
## observações de uma mesma parcela na profundidade).
m0 <- lme(y~bloco+sistema*gesso*Prof, random=~1|parc/subp, data=da,
method="ML")
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
bp <- unlist(ranef(m0)[1])
bs <- unlist(ranef(m0)[2])
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r), qqmath(bp), qqmath(bs))
## Embora esteja considerando só a parte de efeito fixo serve de indicação.
MASS::boxcox(lm(formula(m0), data=da))
## Embora esteja considerando só a parte de efeito aleatório um nível
## anterior ao de resíduo, serve de indicação.
MASS::boxcox(lm(y~subp, data=da))
##-----------------------------------------------------------------------------
## Modelo com a variável transformada. Refazer a análise dos resíduos.
m0 <- update(m0, log(y)~.)
r <- residuals(m0, type="pearson")
f <- fitted(m0)
bp <- unlist(ranef(m0)[1])
bs <- unlist(ranef(m0)[2])
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r), qqmath(bp), qqmath(bs))
## Teste de Wald para os termos de efeito fixo.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 48 1103.4 <.0001
## bloco 3 3 4.3 0.1312
## sistema 1 3 32.3 0.0108
## gesso 1 6 0.8 0.4006
## Prof 4 48 197.2 <.0001
## sistema:gesso 1 6 0.0 0.9324
## sistema:Prof 4 48 6.6 0.0003
## gesso:Prof 4 48 4.4 0.0042
## sistema:gesso:Prof 4 48 1.4 0.2381
##-----------------------------------------------------------------------------
## Mas e a correlação serial ao longo do perfil?
plot(ACF(m0))
S <- split(data.frame(r=r, prof=da$prof), da$subp)
S <- lapply(S,
function(i){
i$r[order(i$prof)]
})
S <- do.call(rbind, S)
splom(S, type=c("p","smooth"))
round(cor(S), 2)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00 0.13 -0.38 -0.42 0.69
## [2,] 0.13 1.00 0.30 -0.10 0.11
## [3,] -0.38 0.30 1.00 -0.05 -0.33
## [4,] -0.42 -0.10 -0.05 1.00 0.16
## [5,] 0.69 0.11 -0.33 0.16 1.00
## Como se comporta dentro de cada sistema de cultivo?
f <- factor(do.call(rbind, strsplit(rownames(S), "\\."))[,2])
by(S, f, cor)
## INDICES: PC
## V1 V2 V3 V4 V5
## V1 1.0000 0.0315 -0.6331 -0.68371 0.68730
## V2 0.0315 1.0000 0.5395 -0.24830 0.17403
## V3 -0.6331 0.5395 1.0000 0.47558 -0.12570
## V4 -0.6837 -0.2483 0.4756 1.00000 -0.03008
## V5 0.6873 0.1740 -0.1257 -0.03008 1.00000
## -------------------------------------------------------------------
## INDICES: PD
## V1 V2 V3 V4 V5
## V1 1.0000 0.23038 -0.1248 -0.11472 0.70587
## V2 0.2304 1.00000 0.1789 -0.02529 0.08042
## V3 -0.1248 0.17887 1.0000 -0.43709 -0.60038
## V4 -0.1147 -0.02529 -0.4371 1.00000 0.41162
## V5 0.7059 0.08042 -0.6004 0.41162 1.00000
## Estrutura de correlação não estruturada (tem k*(k-1)/2 parâmetros).
m1 <- update(m0, correlation=corSymm(form=~as.integer(Prof)|parc/subp))
summary(m1)
## Linear mixed-effects model fit by maximum likelihood
## Data: da
## AIC BIC logLik
## 64.52 150.3 3.74
##
## Random effects:
## Formula: ~1 | parc
## (Intercept)
## StdDev: 0.02511
##
## Formula: ~1 | subp %in% parc
## (Intercept) Residual
## StdDev: 0.07775 0.354
##
## Correlation Structure: General
## Formula: ~as.integer(Prof) | parc/subp
## Parameter estimate(s):
## Correlation:
## 1 2 3 4
## 2 0.634
## 3 0.053 0.407
## 4 -0.378 -0.198 -0.036
## 5 0.940 0.659 0.129 -0.078
## Fixed effects: log(y) ~ bloco + sistema + gesso + Prof + sistema:gesso + sistema:Prof + gesso:Prof + sistema:gesso:Prof
## Value Std.Error DF t-value p-value
## (Intercept) 3.397 0.2246 48 15.12 0.0000
## bloco2 -0.390 0.1048 3 -3.72 0.0338
## bloco3 -0.152 0.1048 3 -1.45 0.2426
## bloco4 -0.374 0.1048 3 -3.57 0.0375
## sistemaPD -0.257 0.3044 3 -0.85 0.4600
## gesso2 -0.119 0.3037 6 -0.39 0.7084
## Prof10 -0.115 0.1793 48 -0.64 0.5232
## Prof15 -1.552 0.2886 48 -5.38 0.0000
## Prof20 -2.821 0.3482 48 -8.10 0.0000
## Prof25 -2.821 0.0729 48 -38.70 0.0000
## sistemaPD:gesso2 0.620 0.4294 6 1.44 0.1990
## sistemaPD:Prof10 -0.859 0.2536 48 -3.39 0.0014
## sistemaPD:Prof15 -0.492 0.4082 48 -1.20 0.2344
## sistemaPD:Prof20 0.257 0.4924 48 0.52 0.6037
## sistemaPD:Prof25 -0.089 0.1031 48 -0.87 0.3910
## gesso2:Prof10 -0.494 0.2536 48 -1.95 0.0574
## gesso2:Prof15 -0.142 0.4082 48 -0.35 0.7291
## gesso2:Prof20 0.695 0.4924 48 1.41 0.1647
## gesso2:Prof25 0.119 0.1031 48 1.16 0.2536
## sistemaPD:gesso2:Prof10 -0.432 0.3587 48 -1.21 0.2339
## sistemaPD:gesso2:Prof15 -0.878 0.5773 48 -1.52 0.1347
## sistemaPD:gesso2:Prof20 -1.094 0.6964 48 -1.57 0.1227
## sistemaPD:gesso2:Prof25 -0.620 0.1458 48 -4.25 0.0001
## Correlation:
## (Intr) bloco2 bloco3 bloco4 sstmPD gesso2 Prof10 Prof15 Prof20
## bloco2 -0.233
## bloco3 -0.233 0.500
## bloco4 -0.233 0.500 0.500
## sistemaPD -0.678 0.000 0.000 0.000
## gesso2 -0.676 0.000 0.000 0.000 0.499
## Prof10 -0.399 0.000 0.000 0.000 0.295 0.295
## Prof15 -0.643 0.000 0.000 0.000 0.474 0.475 0.611
## Prof20 -0.775 0.000 0.000 0.000 0.572 0.573 0.385 0.564
## Prof25 -0.162 0.000 0.000 0.000 0.120 0.120 0.285 0.285 0.624
## sistemaPD:gesso2 0.478 0.000 0.000 0.000 -0.705 -0.707 -0.209 -0.336 -0.405
## sistemaPD:Prof10 0.282 0.000 0.000 0.000 -0.417 -0.209 -0.707 -0.432 -0.272
## sistemaPD:Prof15 0.454 0.000 0.000 0.000 -0.671 -0.336 -0.432 -0.707 -0.399
## sistemaPD:Prof20 0.548 0.000 0.000 0.000 -0.809 -0.405 -0.272 -0.399 -0.707
## sistemaPD:Prof25 0.115 0.000 0.000 0.000 -0.169 -0.085 -0.201 -0.202 -0.441
## gesso2:Prof10 0.282 0.000 0.000 0.000 -0.208 -0.418 -0.707 -0.432 -0.272
## gesso2:Prof15 0.454 0.000 0.000 0.000 -0.335 -0.672 -0.432 -0.707 -0.399
## gesso2:Prof20 0.548 0.000 0.000 0.000 -0.404 -0.811 -0.272 -0.399 -0.707
## gesso2:Prof25 0.115 0.000 0.000 0.000 -0.085 -0.170 -0.201 -0.202 -0.441
## sistemaPD:gesso2:Prof10 -0.200 0.000 0.000 0.000 0.295 0.295 0.500 0.306 0.192
## sistemaPD:gesso2:Prof15 -0.321 0.000 0.000 0.000 0.474 0.475 0.306 0.500 0.282
## sistemaPD:gesso2:Prof20 -0.388 0.000 0.000 0.000 0.572 0.573 0.192 0.282 0.500
## sistemaPD:gesso2:Prof25 -0.081 0.000 0.000 0.000 0.120 0.120 0.142 0.143 0.312
## Prof25 ssPD:2 sPD:P10 sPD:P15 sPD:P20 sPD:P25 g2:P10 g2:P15
## bloco2
## bloco3
## bloco4
## sistemaPD
## gesso2
## Prof10
## Prof15
## Prof20
## Prof25
## sistemaPD:gesso2 -0.085
## sistemaPD:Prof10 -0.201 0.295
## sistemaPD:Prof15 -0.202 0.475 0.611
## sistemaPD:Prof20 -0.441 0.573 0.385 0.564
## sistemaPD:Prof25 -0.707 0.120 0.285 0.285 0.624
## gesso2:Prof10 -0.201 0.295 0.500 0.306 0.192 0.142
## gesso2:Prof15 -0.202 0.475 0.306 0.500 0.282 0.143 0.611
## gesso2:Prof20 -0.441 0.573 0.192 0.282 0.500 0.312 0.385 0.564
## gesso2:Prof25 -0.707 0.120 0.142 0.143 0.312 0.500 0.285 0.285
## sistemaPD:gesso2:Prof10 0.142 -0.418 -0.707 -0.432 -0.272 -0.201 -0.707 -0.432
## sistemaPD:gesso2:Prof15 0.143 -0.672 -0.432 -0.707 -0.399 -0.202 -0.432 -0.707
## sistemaPD:gesso2:Prof20 0.312 -0.811 -0.272 -0.399 -0.707 -0.441 -0.272 -0.399
## sistemaPD:gesso2:Prof25 0.500 -0.170 -0.201 -0.202 -0.441 -0.707 -0.201 -0.202
## g2:P20 g2:P25 sPD:2:P10 sPD:2:P15 sPD:2:P20
## bloco2
## bloco3
## bloco4
## sistemaPD
## gesso2
## Prof10
## Prof15
## Prof20
## Prof25
## sistemaPD:gesso2
## sistemaPD:Prof10
## sistemaPD:Prof15
## sistemaPD:Prof20
## sistemaPD:Prof25
## gesso2:Prof10
## gesso2:Prof15
## gesso2:Prof20
## gesso2:Prof25 0.624
## sistemaPD:gesso2:Prof10 -0.272 -0.201
## sistemaPD:gesso2:Prof15 -0.399 -0.202 0.611
## sistemaPD:gesso2:Prof20 -0.707 -0.441 0.385 0.564
## sistemaPD:gesso2:Prof25 -0.441 -0.707 0.285 0.285 0.624
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.05313 -0.62629 0.01904 0.49545 2.21993
##
## Number of Observations: 80
## Number of Groups:
## parc subp %in% parc
## 8 16
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 26 87.32 149.2 -17.66
## m1 2 36 64.52 150.3 3.74 1 vs 2 42.8 <.0001
## Estrutura de continous AR(1).
m2 <- update(m0, correlation=corCAR1(form=~prof|parc/subp))
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 26 87.32 149.2 -17.66
## m2 2 27 88.13 152.4 -17.06 1 vs 2 1.192 0.2749
AIC(m1)
## [1] 64.52
AIC(m2)
## [1] 88.13
##-----------------------------------------------------------------------------
## Veja como os resultados mudam completamente.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 48 1103.4 <.0001
## bloco 3 3 4.3 0.1312
## sistema 1 3 32.3 0.0108
## gesso 1 6 0.8 0.4006
## Prof 4 48 197.2 <.0001
## sistema:gesso 1 6 0.0 0.9324
## sistema:Prof 4 48 6.6 0.0003
## gesso:Prof 4 48 4.4 0.0042
## sistema:gesso:Prof 4 48 1.4 0.2381
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 48 23832 <.0001
## bloco 3 3 6 0.0810
## sistema 1 3 14 0.0333
## gesso 1 6 38 0.0008
## Prof 4 48 2276 <.0001
## sistema:gesso 1 6 27 0.0020
## sistema:Prof 4 48 19 <.0001
## gesso:Prof 4 48 8 <.0001
## sistema:gesso:Prof 4 48 5 0.0012
##-----------------------------------------------------------------------------
## Qual foi a correlação estimada?
## Embora seja aceitável a correlação de 63% entre 1-2 e de 40% entre
## 2-3, as correlações negativas com a camada 4 as correlações
## superiores à 60% de entre 1-5 e 2-5 parecem bem estranhas. Será que o
## fato dos sistemas serem diferentes quando ao revolvimento do solo
## (PD=não mexe, PC=mexe) não justifica uma especificação de covariância
## ao longo do perfil que seja separada por sistema?
summary(m1)$modelStruct$corStruct
## Correlation structure of class corSymm representing
## Correlation:
## 1 2 3 4
## 2 0.634
## 3 0.053 0.407
## 4 -0.378 -0.198 -0.036
## 5 0.940 0.659 0.129 -0.078
##-----------------------------------------------------------------------------
## O interesse é saber como o gesso altera o nível dos nutrientes em
## cada sistema de manejo. Então obter valores preditos para
## profundidade em cada nível de gesso e sistema.
X <- LSmatrix(lm(formula(m1), data=da),
effect=c("sistema","gesso","Prof"))
grid <- attr(X, "grid")
## Médias com IC.
ci <- confint(glht(m1, linfct=X), calpha=univariate_calpha())
grid <- cbind(grid, ci$confint)
grid <- equallevels(grid, da)
grid
## sistema gesso Prof Estimate lwr upr
## 1 PC 0 5 3.168e+00 2.811560 3.5237
## 2 PD 0 5 2.910e+00 2.554241 3.2664
## 3 PC 2 5 3.049e+00 2.692420 3.4046
## 4 PD 2 5 3.411e+00 3.054923 3.7671
## 5 PC 0 10 3.052e+00 2.696237 3.4084
## 6 PD 0 10 1.936e+00 1.579520 2.2917
## 7 PC 2 10 2.439e+00 2.083246 2.7954
## 8 PD 2 10 1.510e+00 1.153984 1.8661
## 9 PC 0 15 1.615e+00 1.259287 1.9714
## 10 PD 0 15 8.664e-01 0.510354 1.2225
## 11 PC 2 15 1.354e+00 0.997945 1.7101
## 12 PD 2 15 3.466e-01 -0.009506 0.7027
## 13 PC 0 20 3.466e-01 -0.009506 0.7027
## 14 PD 0 20 3.466e-01 -0.009506 0.7027
## 15 PC 2 20 9.222e-01 0.566140 1.2783
## 16 PD 2 20 4.479e-01 0.091860 0.8040
## 17 PC 0 25 3.466e-01 -0.009506 0.7027
## 18 PD 0 25 -7.216e-16 -0.356080 0.3561
## 19 PC 2 25 3.466e-01 -0.009506 0.7027
## 20 PD 2 25 -1.776e-15 -0.356080 0.3561
##-----------------------------------------------------------------------------
## Representação.
## l <- sapply(levels(grid$gesso),
## function(l){
## eval(substitute(expression(gesso~(ton~ha^{-1})),
## list(gesso=l)))
## })
l <- levels(grid$gesso)
key <- list(type="o", divide=1, columns=length(l),
title=expression("Gesso agrícola"~(ton~ha^{-1})),
cex.title=1.1,
lines=list(pch=1:length(l),
col=trellis.par.get("plot.symbol")$col),
text=list(l))
## Para que os pontos sejam corretamente representados, deve-se
## reordenar os dados.
grid <- grid[with(grid, order(sistema, Prof, gesso)),]
segplot(Prof~lwr+upr|sistema, horizontal=FALSE,
centers=Estimate, groups=grid$gesso, f=0.05, layout=c(NA,1),
data=grid, draw=FALSE, panel=panel.segplotBy, key=key,
pch=as.integer(grid$gesso),
ylab="log Fósforo",
xlab="Camada do solo (cm)")
##-----------------------------------------------------------------------------
## Como ficariam os resultados se fosse assumido independência?
grie <- attr(X, "grid")
ci <- confint(glht(m0, linfct=X), calpha=univariate_calpha())
grie <- cbind(grie, ci$confint)
grie <- equallevels(grie, da)
grie <- grie[with(grie, order(sistema, Prof, gesso)),]
grid$modelo <- "correlação"
grie$modelo <- "independente"
gri <- rbind(grid, grie)
gri$modelo <- factor(gri$modelo)
levels(gri$modelo)
## [1] "correlação" "independente"
l <- levels(gri$modelo)
key <- list(type="o", divide=1, columns=length(l),
title="Modelo",
cex.title=1.1,
lines=list(pch=1:length(l),
col=trellis.par.get("plot.symbol")$col),
text=list(l))
gri <- gri[with(gri, order(sistema, gesso, Prof, modelo)),]
segplot(Prof~lwr+upr|sistema*gesso, horizontal=FALSE,
centers=Estimate, groups=gri$modelo, f=0.05,
data=gri, draw=FALSE, panel=panel.segplotBy, key=key,
pch=as.integer(gri$modelo),
ylab="log Fósforo",
xlab="Camada do solo (cm)")
##-----------------------------------------------------------------------------
## Ligando as médias para dar inteção de continuidade.
d <- (as.integer(gri$modelo)-median(seq_along(l)))/5
key <- list(type="o", divide=1, columns=length(l),
title="Modelo",
cex.title=1.1,
lines=list(pch=1:length(l),
col=trellis.par.get("superpose.symbol")$col[1:length(l)]),
text=list(l))
xyplot(Estimate~Prof|sistema*gesso, type="o",
data=gri,
groups=modelo,
ly=gri$lwr, uy=gri$upr,
cty="bars", desloc=d, length=0.02,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose,
key=key, pch=as.integer(gri$modelo),
ylab="Fósforo",
xlab="Camada do solo (cm)")
