Walmes Zeviani
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp",
"reshape", "plyr", "nlme")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra doBy multcomp reshape plyr nlme
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
source("http://dl.dropboxusercontent.com/u/48140237/bandas.R")
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/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] splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] nlme_3.1-110 reshape_0.8.4 plyr_1.8 multcomp_1.2-21
## [5] mvtnorm_0.9-9994 doBy_4.5-10 MASS_7.3-33 survival_2.37-7
## [9] latticeExtra_0.6-24 RColorBrewer_1.0-5 lattice_0.20-29 knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.1 formatR_0.9 grid_3.1.0 lme4_1.0-5 Matrix_1.1-3 methods_3.1.0
## [7] minqa_1.2.2 stringr_0.6.2 tools_3.1.0
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])
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))
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)~.)
## 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)
round(cor(S), 2)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00 -0.03 -0.63 -0.74 -0.54
## [2,] -0.03 1.00 0.53 -0.16 -0.27
## [3,] -0.63 0.53 1.00 0.43 0.26
## [4,] -0.74 -0.16 0.43 1.00 0.76
## [5,] -0.54 -0.27 0.26 0.76 1.00
f <- factor(do.call(rbind, strsplit(rownames(S), "\\."))[,2])
by(S, f, cor)
## INDICES: PC
## V1 V2 V3 V4 V5
## V1 1.00000 -0.05257 -0.7637 -0.8692 -0.6616
## V2 -0.05257 1.00000 0.5347 -0.3156 -0.4437
## V3 -0.76368 0.53471 1.0000 0.5066 0.3067
## V4 -0.86924 -0.31561 0.5066 1.0000 0.6731
## V5 -0.66163 -0.44369 0.3067 0.6731 1.0000
## ---------------------------------------------------------------------------
## INDICES: PD
## V1 V2 V3 V4 V5
## V1 1.00000 0.05694 0.04814 -0.4826 -0.2610
## V2 0.05694 1.00000 0.50713 0.2368 0.2376
## V3 0.04814 0.50713 1.00000 0.2831 0.1594
## V4 -0.48264 0.23681 0.28310 1.0000 0.8992
## V5 -0.26104 0.23763 0.15936 0.8992 1.0000
## 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.07777 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.079
## 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 Prof25
## 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 -0.085
## sistemaPD:Prof10 0.282 0.000 0.000 0.000 -0.417 -0.209 -0.707 -0.432 -0.272 -0.201
## sistemaPD:Prof15 0.454 0.000 0.000 0.000 -0.671 -0.336 -0.432 -0.707 -0.399 -0.202
## sistemaPD:Prof20 0.548 0.000 0.000 0.000 -0.809 -0.405 -0.272 -0.399 -0.707 -0.441
## sistemaPD:Prof25 0.115 0.000 0.000 0.000 -0.169 -0.085 -0.201 -0.202 -0.441 -0.707
## gesso2:Prof10 0.282 0.000 0.000 0.000 -0.208 -0.418 -0.707 -0.432 -0.272 -0.201
## gesso2:Prof15 0.454 0.000 0.000 0.000 -0.335 -0.672 -0.432 -0.707 -0.399 -0.202
## gesso2:Prof20 0.548 0.000 0.000 0.000 -0.404 -0.811 -0.272 -0.399 -0.707 -0.441
## gesso2:Prof25 0.115 0.000 0.000 0.000 -0.085 -0.170 -0.201 -0.202 -0.441 -0.707
## sistemaPD:gesso2:Prof10 -0.200 0.000 0.000 0.000 0.295 0.295 0.500 0.306 0.192 0.142
## sistemaPD:gesso2:Prof15 -0.321 0.000 0.000 0.000 0.474 0.475 0.306 0.500 0.282 0.143
## sistemaPD:gesso2:Prof20 -0.388 0.000 0.000 0.000 0.572 0.573 0.192 0.282 0.500 0.312
## sistemaPD:gesso2:Prof25 -0.081 0.000 0.000 0.000 0.120 0.120 0.142 0.143 0.312 0.500
## ssPD:2 sPD:P10 sPD:P15 sPD:P20 sPD:P25 g2:P10 g2:P15 g2:P20 g2:P25
## bloco2
## bloco3
## bloco4
## sistemaPD
## gesso2
## Prof10
## Prof15
## Prof20
## Prof25
## sistemaPD:gesso2
## sistemaPD:Prof10 0.295
## sistemaPD:Prof15 0.475 0.611
## sistemaPD:Prof20 0.573 0.385 0.564
## sistemaPD:Prof25 0.120 0.285 0.285 0.624
## gesso2:Prof10 0.295 0.500 0.306 0.192 0.142
## gesso2:Prof15 0.475 0.306 0.500 0.282 0.143 0.611
## gesso2:Prof20 0.573 0.192 0.282 0.500 0.312 0.385 0.564
## gesso2:Prof25 0.120 0.142 0.143 0.312 0.500 0.285 0.285 0.624
## sistemaPD:gesso2:Prof10 -0.418 -0.707 -0.432 -0.272 -0.201 -0.707 -0.432 -0.272 -0.201
## sistemaPD:gesso2:Prof15 -0.672 -0.432 -0.707 -0.399 -0.202 -0.432 -0.707 -0.399 -0.202
## sistemaPD:gesso2:Prof20 -0.811 -0.272 -0.399 -0.707 -0.441 -0.272 -0.399 -0.707 -0.441
## sistemaPD:gesso2:Prof25 -0.170 -0.201 -0.202 -0.441 -0.707 -0.201 -0.202 -0.441 -0.707
## 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
## sistemaPD:gesso2:Prof10
## sistemaPD:gesso2:Prof15 0.611
## sistemaPD:gesso2:Prof20 0.385 0.564
## sistemaPD:gesso2:Prof25 0.285 0.285 0.624
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.05309 -0.62628 0.01909 0.49551 2.21989
##
## 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
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 de figura.
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
##-----------------------------------------------------------------------------
## 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 2.873e-15 -0.356080 0.3561
## 19 PC 2 25 3.466e-01 -0.009506 0.7027
## 20 PD 2 25 2.331e-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, by=grid$gesso, f=0.05, layout=c(NA,1),
data=grid, draw=FALSE, panel=segplot.by, key=key,
pch=as.integer(grid$gesso),
ylab="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, by=gri$modelo, f=0.05,
data=gri, draw=FALSE, panel=segplot.by, key=key,
pch=as.integer(gri$modelo),
ylab="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)")
