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/apc.R")
source("http://dl.dropboxusercontent.com/u/48140237/segplotby.R")
source("http://dl.dropboxusercontent.com/u/48140237/equallevels.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/soja_solo_tcc.txt",
header=TRUE, sep="\t")
str(da)
## 'data.frame': 80 obs. of 16 variables:
## $ sistema: Factor w/ 2 levels "convencional",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ adubpk : int 40 40 40 40 40 40 40 40 60 60 ...
## $ prof : Factor w/ 2 levels "00-20","20-40": 1 1 1 1 2 2 2 2 1 1 ...
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Zn : num 0.8 0.9 0.9 1.2 0.7 0.3 0.7 1.1 1.3 1.1 ...
## $ phcacl2: num 5.1 5.3 4.7 5 4.6 4.3 4.4 4.5 5.8 6.5 ...
## $ phh2o : num 5.9 6.3 5.3 5.6 5.2 5 5 5.1 6.6 7.5 ...
## $ P : int 3 7 2 2 2 2 1 2 4 6 ...
## $ K : num 2.2 1.7 2.8 4.4 1.6 1 1.9 3.8 3 1.5 ...
## $ Al : num 0 0 2.5 0.6 3.1 7.4 6.6 8 0 0 ...
## $ Ca : int 45 58 30 37 27 19 23 22 27 81 ...
## $ Mg : int 34 33 16 18 23 12 13 10 36 33 ...
## $ Hal : int 47 38 72 53 69 89 80 76 29 19 ...
## $ sb : num 81.2 92.7 48.8 59.4 51.6 ...
## $ ctc : num 128 131 121 112 121 ...
## $ v : int 63 70 40 52 42 26 32 32 76 85 ...
##-----------------------------------------------------------------------------
## Visualizar.
xyplot(ctc~adubpk|sistema, groups=prof, data=da, type=c("p","smooth"),
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 e subparcelas.
da$parc <- with(da, interaction(bloco, sistema, drop=TRUE))
da$subp <- with(da, interaction(bloco, sistema, adub, drop=TRUE))
str(da)
## 'data.frame': 80 obs. of 19 variables:
## $ sistema: Factor w/ 2 levels "convencional",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ adubpk : int 40 40 40 40 40 40 40 40 60 60 ...
## $ prof : Factor w/ 2 levels "00-20","20-40": 1 1 1 1 2 2 2 2 1 1 ...
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Zn : num 0.8 0.9 0.9 1.2 0.7 0.3 0.7 1.1 1.3 1.1 ...
## $ phcacl2: num 5.1 5.3 4.7 5 4.6 4.3 4.4 4.5 5.8 6.5 ...
## $ phh2o : num 5.9 6.3 5.3 5.6 5.2 5 5 5.1 6.6 7.5 ...
## $ P : int 3 7 2 2 2 2 1 2 4 6 ...
## $ K : num 2.2 1.7 2.8 4.4 1.6 1 1.9 3.8 3 1.5 ...
## $ Al : num 0 0 2.5 0.6 3.1 7.4 6.6 8 0 0 ...
## $ Ca : int 45 58 30 37 27 19 23 22 27 81 ...
## $ Mg : int 34 33 16 18 23 12 13 10 36 33 ...
## $ Hal : int 47 38 72 53 69 89 80 76 29 19 ...
## $ sb : num 81.2 92.7 48.8 59.4 51.6 ...
## $ ctc : num 128 131 121 112 121 ...
## $ v : int 63 70 40 52 42 26 32 32 76 85 ...
## $ adub : Factor w/ 5 levels "40","60","90",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ parc : Factor w/ 8 levels "I.convencional",..: 5 6 7 8 5 6 7 8 5 6 ...
## $ subp : Factor w/ 40 levels "I.convencional.40",..: 5 6 7 8 5 6 7 8 13 14 ...
## Adub categórico.
m0 <- lme(ctc~bloco+sistema*adub*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)
##-----------------------------------------------------------------------------
## Quadro para o teste dos efeitos fixos (Wald).
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 30 12774 <.0001
## bloco 3 3 2 0.2364
## sistema 1 3 3 0.1747
## adub 4 24 1 0.3130
## prof 1 30 16 0.0004
## sistema:adub 4 24 1 0.4574
## sistema:prof 1 30 5 0.0303
## adub:prof 4 30 1 0.5711
## sistema:adub:prof 4 30 0 0.8630
## Abandono das interações não significativas.
m1 <- update(m0, .~bloco+sistema*prof+adubpk)
anova(m1, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 11 572.1 598.3 -275.1
## m0 2 26 589.9 651.8 -268.9 1 vs 2 12.24 0.6605
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 38 13649 <.0001
## bloco 3 3 3 0.2211
## sistema 1 3 3 0.1645
## prof 1 38 18 0.0002
## adubpk 1 31 4 0.0585
## sistema:prof 1 38 6 0.0222
## Sem efeito da adubação. O que existe é um comportamento distinto da
## ctc na profundidade devido ao sistema.
##-----------------------------------------------------------------------------
## Obter as médias de ctc em cada prof para cada sistema. Pode-se ficar
## a adubação em qualquer nível, bem como bloco.
X <- LSmatrix(lm(formula(m1), data=da), effect=c("sistema","prof"),
at=list(adubpk=120))
head(X)
## (Intercept) blocoII blocoIII blocoIV sistemadireto prof20-40 adubpk sistemadireto:prof20-40
## [1,] 1 0.25 0.25 0.25 0 0 120 0
## [2,] 1 0.25 0.25 0.25 1 0 120 0
## [3,] 1 0.25 0.25 0.25 0 1 120 0
## [4,] 1 0.25 0.25 0.25 1 1 120 1
grid <- attr(X, "grid")
L <- by(X, grid$sistema, apc, lev=c(20,40))
L <- lapply(L, as.matrix); L
## $convencional
## (Intercept) blocoII blocoIII blocoIV sistemadireto prof20-40 adubpk sistemadireto:prof20-40
## 20-40 0 0 0 0 0 -1 0 0
##
## $direto
## (Intercept) blocoII blocoIII blocoIV sistemadireto prof20-40 adubpk sistemadireto:prof20-40
## 20-40 0 0 0 0 0 -1 0 -1
str(L)
## List of 2
## $ convencional: num [1, 1:8] 0 0 0 0 0 -1 0 0
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr "20-40"
## .. ..$ : chr [1:8] "(Intercept)" "blocoII" "blocoIII" "blocoIV" ...
## $ direto : num [1, 1:8] 0 0 0 0 0 -1 0 -1
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr "20-40"
## .. ..$ : chr [1:8] "(Intercept)" "blocoII" "blocoIII" "blocoIV" ...
Xc <- do.call(rbind, L)
rownames(Xc) <- paste0(names(L), "/", rownames(Xc))
Xc
## (Intercept) blocoII blocoIII blocoIV sistemadireto prof20-40 adubpk
## convencional/20-40 0 0 0 0 0 -1 0
## direto/20-40 0 0 0 0 0 -1 0
## sistemadireto:prof20-40
## convencional/20-40 0
## direto/20-40 -1
## Teste para os contrastes.
summary(glht(m1, linfct=Xc))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lme.formula(fixed = ctc ~ bloco + sistema + prof + adubpk + sistema:prof,
## data = da, random = ~1 | parc/subp, method = "ML")
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## convencional/20-40 == 0 2.69 1.99 1.35 0.32
## direto/20-40 == 0 9.77 1.99 4.90 1.9e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## Médias com IC.
ci <- confint(glht(m1, linfct=X), calpha=univariate_calpha())
grid <- cbind(grid, ci$confint)
grid <- equallevels(grid, da)
grid
## sistema prof adubpk Estimate lwr upr
## 1 convencional 00-20 120 124.7 121.1 128.4
## 2 direto 00-20 120 132.2 128.5 135.9
## 3 convencional 20-40 120 122.1 118.4 125.7
## 4 direto 20-40 120 122.4 118.7 126.1
##-----------------------------------------------------------------------------
## Representação.
names(trellis.par.get())
## [1] "grid.pars" "fontsize" "background" "panel.background"
## [5] "clip" "add.line" "add.text" "plot.polygon"
## [9] "box.dot" "box.rectangle" "box.umbrella" "dot.line"
## [13] "dot.symbol" "plot.line" "plot.symbol" "reference.line"
## [17] "strip.background" "strip.shingle" "strip.border" "superpose.line"
## [21] "superpose.symbol" "superpose.polygon" "regions" "shade.colors"
## [25] "axis.line" "axis.text" "axis.components" "layout.heights"
## [29] "layout.widths" "box.3d" "par.xlab.text" "par.ylab.text"
## [33] "par.zlab.text" "par.main.text" "par.sub.text"
l <- levels(grid$prof)
key <- list(type="o", divide=1,
title="Camada do solo (cm)", cex.title=1.1,
lines=list(pch=1:length(l),
col=trellis.par.get("plot.symbol")$col),
text=list(l))
segplot(sistema~lwr+upr, centers=Estimate, by=grid$prof, f=0.05,
data=grid, draw=FALSE, panel=segplot.by, key=key,
pch=as.integer(grid$prof),
xlab="CTC do solo",
ylab="Sistema de cultivo")
##-----------------------------------------------------------------------------
## Os erros padrões de comparações de prof dentro de sistema e de
## sistema dentro de prof.
## Desdobrar prof dentro (fixando) de sistema.
L <- by(X, grid$sistema, apc, lev=c(20,40))
L <- lapply(L, as.matrix)
Xc <- do.call(rbind, L)
rownames(Xc) <- paste0(names(L), "/", rownames(Xc))
summary(glht(m1, linfct=Xc))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lme.formula(fixed = ctc ~ bloco + sistema + prof + adubpk + sistema:prof,
## data = da, random = ~1 | parc/subp, method = "ML")
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## convencional/20-40 == 0 2.69 1.99 1.35 0.32
## direto/20-40 == 0 9.77 1.99 4.90 1.9e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## Desdobrar sistema dentro (fixando) de prof.
L <- by(X, grid$prof, apc, lev=levels(grid$sistema))
L <- lapply(L, as.matrix)
Xc <- do.call(rbind, L)
rownames(Xc) <- paste0(names(L), "/", rownames(Xc))
summary(glht(m1, linfct=Xc))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lme.formula(fixed = ctc ~ bloco + sistema + prof + adubpk + sistema:prof,
## data = da, random = ~1 | parc/subp, method = "ML")
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 00-20/convencional-direto == 0 -7.42 2.46 -3.02 0.005 **
## 20-40/convencional-direto == 0 -0.34 2.46 -0.14 0.987
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summary(glht(m1, linfct=X))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lme.formula(fixed = ctc ~ bloco + sistema + prof + adubpk + sistema:prof,
## data = da, random = ~1 | parc/subp, method = "ML")
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 == 0 124.75 1.88 66.5 <2e-16 ***
## 2 == 0 132.17 1.88 70.4 <2e-16 ***
## 3 == 0 122.06 1.88 65.0 <2e-16 ***
## 4 == 0 122.40 1.88 65.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)