Análise de experimento em parcela subsubdividida
##-----------------------------------------------------------------------------
## 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
## 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.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-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] methods splines grid stats graphics grDevices utils datasets
## [9] base
##
## other attached packages:
## [1] wzRfun_0.3 nlme_3.1-117 plyr_1.8.1 reshape_0.8.5
## [5] multcomp_1.3-7 TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-11
## [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 rmarkdown_0.3.3 knitr_1.7
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 htmltools_0.2.6 Matrix_1.1-4
## [6] Rcpp_0.11.3 sandwich_2.3-0 stringr_0.6.2 tools_3.1.1 yaml_2.1.13
## [11] zoo_1.7-10
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])
grid.arrange(nrow=3,
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 12773.791 <.0001
## bloco 3 3 2.493 0.2364
## sistema 1 3 3.137 0.1747
## adub 4 24 1.260 0.3130
## prof 1 30 15.995 0.0004
## sistema:adub 4 24 0.941 0.4574
## sistema:prof 1 30 5.167 0.0303
## adub:prof 4 30 0.742 0.5711
## sistema:adub:prof 4 30 0.319 0.8630
## Abandono das interações não significativas.
m1 <- update(m0, .~bloco+sistema*prof+adubpk)
## m2 <- update(m0, .~bloco+sistema*prof)
anova(m1, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 11 572.0943 598.2966 -275.0472
## m0 2 26 589.8503 651.7830 -268.9252 1 vs 2 12.244 0.6605
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 38 13649.227 <.0001
## bloco 3 3 2.664 0.2211
## sistema 1 3 3.352 0.1645
## prof 1 38 17.586 0.0002
## adubpk 1 31 3.860 0.0585
## sistema:prof 1 38 5.682 0.0222
## Sem efeito da adubação. O que existe é uma diferença entre camadas
## para a que depende do 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
## [1,] 1 0.25 0.25 0.25 0 0 120
## [2,] 1 0.25 0.25 0.25 1 0 120
## [3,] 1 0.25 0.25 0.25 0 1 120
## [4,] 1 0.25 0.25 0.25 1 1 120
## sistemadireto:prof20-40
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 1
grid <- attr(X, "grid")
rownames(X) <- apply(grid, 1, paste, collapse="|")
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
## 20-40 0 0 0 0 0 -1 0
## sistemadireto:prof20-40
## 20-40 0
##
## $direto
## (Intercept) blocoII blocoIII blocoIV sistemadireto prof20-40 adubpk
## 20-40 0 0 0 0 0 -1 0
## sistemadireto:prof20-40
## 20-40 -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.690 1.994 1.349 0.323
## direto/20-40 == 0 9.775 1.994 4.902 1.89e-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
## convencional|00-20|120 convencional 00-20 120 124.7492 121.0715 128.427
## direto|00-20|120 direto 00-20 120 132.1742 128.4965 135.852
## convencional|20-40|120 convencional 20-40 120 122.0592 118.3815 125.737
## direto|20-40|120 direto 20-40 120 122.3992 118.7215 126.077
##-----------------------------------------------------------------------------
## Representação.
## names(trellis.par.get())
l <- levels(grid$prof)
key <- list(corner=c(x=0.9, y=0.05),
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, groups=grid$prof, f=0.05,
data=grid, draw=FALSE, panel=panel.segplotBy, key=key,
scales=list(y=list(rot=90)),
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 são diferentes pois são fatores com níveis em
## estratos diferentes no modelos misto de parcela subsubdividida.
## 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.690 1.994 1.349 0.323
## direto/20-40 == 0 9.775 1.994 4.902 1.89e-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.425 2.457 -3.022 0.00496 **
## 20-40/convencional-direto == 0 -0.340 2.457 -0.138 0.98712
## ---
## 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|)
## convencional|00-20|120 == 0 124.749 1.876 66.48 <2e-16 ***
## direto|00-20|120 == 0 132.174 1.876 70.44 <2e-16 ***
## convencional|20-40|120 == 0 122.059 1.876 65.05 <2e-16 ***
## direto|20-40|120 == 0 122.399 1.876 65.23 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)