Walmes Zeviani
Milson Evaldo Serafim
Material complementar ao artigo com título ??? submetido à revista ???, volume ???, número ???. A citação será incluída assim que for determinado o volume e paginação.
%% Citação LaTex.
@article{SerafimMaracuja2014,
author = {},
title = {Cinza de cana em substrato de "terra de barranco" para
produção de mudas de maracujazeiro},
doi = {},
journal = {},
month = {},
pages = {},
publisher = {},
url = {},
year = {2014}
}
Donwnload:
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
## Instruções para instalação do wzRfun em:
## browseURL("https://github.com/walmes/wzRfun")
pkg <- c("lattice", "latticeExtra", "grid", "gridExtra",
"reshape", "plyr", "doBy", "multcomp", "MASS", "nlme",
"wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
## lattice latticeExtra grid gridExtra reshape plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## doBy multcomp MASS nlme wzRfun
## TRUE TRUE TRUE TRUE TRUE
##-----------------------------------------------------------------------------
## Funções para essa sessão.
strip.combined <- function(which.given, which.panel,
factor.levels, ...){
if(which.given==1){
panel.rect(0, 0, 1, 1, col="grey90", border=1)
panel.text(x=1, y=0.5, pos=2,
lab=paste("família:",
factor.levels[which.panel[which.given]]))
}
if(which.given==2){
panel.text(x=0, y=0.5, pos=4,
lab=paste("agregado:",
factor.levels[which.panel[which.given]]))
}
}
##-----------------------------------------------------------------------------
## Configurações de cores da lattice.
colp <- brewer.pal(11, "Spectral")
colp <- colorRampPalette(colp, space="rgb")
mycol <- brewer.pal(8, "Dark2")
ps <- list(
box.dot=list(col="black", pch="|"),
box.rectangle=list(col="black", fill=c("gray70")),
box.umbrella=list(col="black", lty=1),
dot.symbol=list(col="black"),
dot.line=list(col="gray90", lty=1),
plot.symbol=list(col="black", cex=0.8),
plot.line=list(col="black", lty=1.2),
plot.polygon=list(col="gray80"),
superpose.line=list(col=mycol),
superpose.symbol=list(col=mycol),
superpose.polygon=list(col=mycol),
regions=list(col=colp(100)),
strip.background=list(col=c("gray90","gray50")),
strip.shingle=list(col=c("gray50","gray90"))
)
## trellis.par.set(ps)
rm(list=c("colp", "mycol", "ps"))
lattice.options(default.theme=modifyList(
standard.theme(color=FALSE),
list(strip.background=list(col="gray90"))))
##-----------------------------------------------------------------------------
## Ler.
url <- "http://www.leg.ufpr.br/~walmes/data/maracuja_altura.txt"
mara <- read.table(url, header=TRUE, sep="\t")
mara <- transform(mara, fam=factor(familia), blc=factor(bloco),
agr=factor(agregado, labels=c("4-10 mm","0-2 mm")),
dta=as.Date(data))
mara$dias <- as.numeric(mara$dta-min(mara$dta))
str(mara)
## 'data.frame': 1344 obs. of 12 variables:
## $ agregado: Factor w/ 2 levels "<2","4 a 10": 2 2 2 2 2 2 2 2 2 2 ...
## $ familia : Factor w/ 2 levels "F29","F48": 1 1 1 1 1 1 1 1 1 1 ...
## $ cinza : num 0 0 0 0 0 0 1.5 1.5 1.5 1.5 ...
## $ bloco : int 1 1 2 2 3 3 1 1 2 2 ...
## $ rept : int 1 2 1 2 1 2 1 2 1 2 ...
## $ data : Factor w/ 8 levels "2012/02/24","2012/03/01",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ altura : num 6 6.5 4.5 4 4.5 4.5 8 7.2 6 7.4 ...
## $ fam : Factor w/ 2 levels "F29","F48": 1 1 1 1 1 1 1 1 1 1 ...
## $ blc : Factor w/ 3 levels "1","2","3": 1 1 2 2 3 3 1 1 2 2 ...
## $ agr : Factor w/ 2 levels "4-10 mm","0-2 mm": 2 2 2 2 2 2 2 2 2 2 ...
## $ dta : Date, format: "2012-02-24" "2012-02-24" ...
## $ dias : num 0 0 0 0 0 0 0 0 0 0 ...
mara$id <- with(mara, interaction(blc, fam, agr, cinza, dias))
nlevels(mara$id) # 3*2*2*7*
## [1] 672
mara <- ddply(mara, .(blc, fam, agr, cinza, dias),
summarise, altura=mean(altura))
str(mara)
## 'data.frame': 672 obs. of 6 variables:
## $ blc : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ fam : Factor w/ 2 levels "F29","F48": 1 1 1 1 1 1 1 1 1 1 ...
## $ agr : Factor w/ 2 levels "4-10 mm","0-2 mm": 1 1 1 1 1 1 1 1 1 1 ...
## $ cinza : num 0 0 0 0 0 0 0 0 1.5 1.5 ...
## $ dias : num 0 6 14 21 28 35 42 47 0 6 ...
## $ altura: num 6.65 7.3 12.25 40.25 62.75 ...
##-----------------------------------------------------------------------------
## Dados de experimento com 4 fatores:
## * fam: categórica a família do maracujá;
## * agr: categórica a classe de agregado usado para encher os vasos de
## cultivo das plantas;
## * cinza: métrica a dose de cinza aplicado ao solo;
## * dias: métrica o dia de avaliação da planta;
## Deseja-se verificar o efeito de fam*agr*cinza no crescimento das
## plantas.
##-----------------------------------------------------------------------------
## Ver.
## xyplot(altura~cinza|fam*agr, groups=dias, data=mara, type=c("p","a"))
xyplot(altura~dias|fam*agr, groups=cinza, data=mara, type=c("p","a"))
## xyplot(log(altura)~cinza|fam*agr, groups=dias, data=mara, type=c("p","a"))
## xyplot(log(altura)~dias|fam*agr, groups=cinza, data=mara, type=c("p","a"))
##-----------------------------------------------------------------------------
## Análise só com a última medida.
## xyplot(altura~cinza|fam, groups=agr, data=subset(mara, dias==47),
## type=c("p","a"))
## xyplot(altura~cinza|fam, groups=agr,
## data=subset(mara[-c(144,456,112,656),], dias==47),
## type=c("p","a"))
##-----------------------------------------------------------------------------
## Modelo.
m0 <- lm(altura~blc+fam*agr*cinza,
data=subset(mara[-c(144,456,112,656),], dias==47))
## Resíduos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de anova.
anova(m0)
## Analysis of Variance Table
##
## Response: altura
## Df Sum Sq Mean Sq F value Pr(>F)
## blc 2 1009 505 5.93 0.0042 **
## fam 1 66 66 0.77 0.3824
## agr 1 230 230 2.70 0.1049
## cinza 1 82 82 0.97 0.3280
## fam:agr 1 9 9 0.10 0.7527
## fam:cinza 1 634 634 7.46 0.0080 **
## agr:cinza 1 670 670 7.88 0.0065 **
## fam:agr:cinza 1 756 756 8.89 0.0039 **
## Residuals 70 5953 85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Usando modelos não lineares mistos.
## Medidas repetidas nos indivíduos que são as parcelas.
mara$parcela <- with(mara, interaction(blc, fam, agr, cinza, sep="_"))
## nlevels(mara$parcela) # 3*2*2*7
marag <- groupedData(formula=altura~dias|parcela, data=mara, order=FALSE)
str(marag)
## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 672 obs. of 7 variables:
## $ blc : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ fam : Factor w/ 2 levels "F29","F48": 1 1 1 1 1 1 1 1 1 1 ...
## $ agr : Factor w/ 2 levels "4-10 mm","0-2 mm": 1 1 1 1 1 1 1 1 1 1 ...
## $ cinza : num 0 0 0 0 0 0 0 0 1.5 1.5 ...
## $ dias : num 0 6 14 21 28 35 42 47 0 6 ...
## $ altura : num 6.65 7.3 12.25 40.25 62.75 ...
## $ parcela: Factor w/ 84 levels "1_F29_4-10 mm_0",..: 1 1 1 1 1 1 1 1 13 13 ...
## - attr(*, "formula")=Class 'formula' length 3 altura ~ dias | parcela
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "FUN")=function (x)
## - attr(*, "order.groups")= logi FALSE
##-----------------------------------------------------------------------------
## Modelo sem efeito dos fatores experimentais.
nn0 <- nlme(altura~SSlogis(dias, A, B, C),
fixed=A+B+C~1,
random=A+B+C~1,
data=marag,
start=list(fixed=c(148, 31, 6.8)))
## Quadro com as estimativas dos parâmetros.
summary(nn0)$tTable
## Value Std.Error DF t-value p-value
## A 148.856 0.9895 586 150.43 0.000e+00
## B 31.428 0.4684 586 67.09 3.399e-277
## C 6.813 0.1429 586 47.67 8.126e-204
VarCorr(nn0)
## parcela = pdLogChol(list(A ~ 1,B ~ 1,C ~ 1))
## Variance StdDev Corr
## A 11.811 3.437 A B
## B 16.928 4.114 0.586
## C 1.084 1.041 -0.097 0.750
## Residual 24.159 4.915
## intervals(nn0)
## plot(nn0)
## pairs(nn0)
## par(mfrow=c(1,3)); apply(ranef(nn0), 2, qqnorm); layout(1)
outl <- which(abs(residuals(nn0, type="pearson"))>3); outl
## 1_F29_4-10 mm_3 1_F29_4-10 mm_3 1_F29_4-10 mm_48 1_F29_4-10 mm_48
## 23 24 54 55
marag <- marag[-outl,]
##-----------------------------------------------------------------------------
## O intercepto não é zero, não sei porque, mas incorporar no modelo.
## modelo: Int+A/(1+exp(-(x-x0)/S)).
## modelo de efeitos aditivos, aleatório em A e B e modelagem da
## variância.
nn1 <- nlme(altura~int+A/(1+exp(-(dias-d50)/S)),
fixed=list(
int~1,
A~blc+fam+agr+cinza,
d50~blc+fam+agr+cinza,
S~blc+fam+agr+cinza),
random=A+d50+S~1,
data=marag,
start=list(
fixed=c(5,
120,0,0,0,0,0,
30,0,0,0,0,0,
4,0,0,0,0,0)))
## Quadro de testes de Wald para termos de efeito fixo.
anova(nn1, type="marginal")
## numDF denDF F-value p-value
## int 1 566 407.8 <.0001
## A.(Intercept) 1 566 2222.8 <.0001
## A.blc 2 566 2.1 0.1190
## A.fam 1 566 0.2 0.6833
## A.agr 1 566 0.8 0.3845
## A.cinza 1 566 0.0 0.8276
## d50.(Intercept) 1 566 1339.3 <.0001
## d50.blc 2 566 4.0 0.0186
## d50.fam 1 566 0.7 0.4012
## d50.agr 1 566 0.1 0.7521
## d50.cinza 1 566 0.0 0.8301
## S.(Intercept) 1 566 605.9 <.0001
## S.blc 2 566 0.9 0.4009
## S.fam 1 566 1.0 0.3130
## S.agr 1 566 0.5 0.4638
## S.cinza 1 566 0.6 0.4325
## summary(nn1)$tTable
VarCorr(nn1)
## parcela = pdLogChol(list(A ~ 1,d50 ~ 1,S ~ 1))
## Variance StdDev Corr
## A.(Intercept) 95.926 9.7942 A.(In) d50.(I
## d50.(Intercept) 8.767 2.9609 -0.308
## S.(Intercept) 0.480 0.6929 -0.436 0.338
## Residual 9.203 3.0337
## plot(nn1)
## par(mfrow=c(1,3)); apply(ranef(nn1), 2, qqnorm); layout(1)
## intervals(nn1)
## Pelo modelo acima ninguém interfere no crescimento da planta.
##-----------------------------------------------------------------------------
## Modelo com interações.
nn2 <- nlme(altura~int+A/(1+exp(-(dias-d50)/S)),
fixed=list(
int~1,
A~blc+fam*agr*cinza,
d50~blc+fam*agr*cinza,
S~blc+fam*agr*cinza),
random=A+d50+S~1,
data=marag,
start=list(
fixed=c(5.2,
120,0,0,0,0,0,0,0,0,0,
30,0,0,0,0,0,0,0,0,0,
4,0,0,0,0,0,0,0,0,0)))
anova(nn2, type="marginal")
## numDF denDF F-value p-value
## int 1 554 402.4 <.0001
## A.(Intercept) 1 554 1535.3 <.0001
## A.blc 2 554 2.6 0.0777
## A.fam 1 554 0.3 0.5882
## A.agr 1 554 4.5 0.0340
## A.cinza 1 554 2.7 0.0984
## A.fam:agr 1 554 3.9 0.0489
## A.fam:cinza 1 554 0.5 0.4877
## A.agr:cinza 1 554 13.2 0.0003
## A.fam:agr:cinza 1 554 6.1 0.0141
## d50.(Intercept) 1 554 856.3 <.0001
## d50.blc 2 554 3.9 0.0201
## d50.fam 1 554 1.3 0.2503
## d50.agr 1 554 0.0 0.9590
## d50.cinza 1 554 0.4 0.5494
## d50.fam:agr 1 554 0.0 0.9844
## d50.fam:cinza 1 554 0.9 0.3517
## d50.agr:cinza 1 554 0.1 0.8176
## d50.fam:agr:cinza 1 554 0.0 0.9192
## S.(Intercept) 1 554 412.4 <.0001
## S.blc 2 554 0.9 0.3989
## S.fam 1 554 0.0 0.8967
## S.agr 1 554 0.0 0.8572
## S.cinza 1 554 0.4 0.5098
## S.fam:agr 1 554 0.4 0.5220
## S.fam:cinza 1 554 0.0 0.8462
## S.agr:cinza 1 554 0.1 0.7825
## S.fam:agr:cinza 1 554 0.1 0.7967
## summary(nn2)$tTable
VarCorr(nn2)
## parcela = pdLogChol(list(A ~ 1,d50 ~ 1,S ~ 1))
## Variance StdDev Corr
## A.(Intercept) 79.7499 8.930 A.(In) d50.(I
## d50.(Intercept) 8.5178 2.919 -0.282
## S.(Intercept) 0.4706 0.686 -0.427 0.342
## Residual 9.1425 3.024
## plot(nn2)
## par(mfrow=c(1,3)); apply(ranef(nn2), 2, qqnorm); layout(1)
## intervals(nn2)
## pairs(nn2)
## Esse parece ser um bom modelo, existe interação tripla apenas para A,
## o resto não muda.
##-----------------------------------------------------------------------------
## Modelo que é uma redução do nn2.
nn3 <- nlme(altura~int+A/(1+exp(-(dias-d50)/S)),
fixed=list(
int~1,
A~blc+fam*agr*cinza,
d50~blc,
S~blc),
random=A+d50+S~1,
data=marag,
start=list(
fixed=c(5.2,
120,0,0,0,0,0,0,0,0,0,
30,0,0,
4,0,0)))
anova(nn3, type="marginal")
## numDF denDF F-value p-value
## int 1 568 412.3 <.0001
## A.(Intercept) 1 568 1560.5 <.0001
## A.blc 2 568 2.8 0.0627
## A.fam 1 568 0.5 0.4791
## A.agr 1 568 4.6 0.0318
## A.cinza 1 568 3.0 0.0845
## A.fam:agr 1 568 4.0 0.0456
## A.fam:cinza 1 568 0.7 0.3922
## A.agr:cinza 1 568 13.1 0.0003
## A.fam:agr:cinza 1 568 6.3 0.0125
## d50.(Intercept) 1 568 2629.2 <.0001
## d50.blc 2 568 3.7 0.0244
## S.(Intercept) 1 568 1189.7 <.0001
## S.blc 2 568 0.8 0.4355
anova(nn3, nn2) ## Ótimo!!
## Model df AIC BIC logLik Test L.Ratio p-value
## nn3 1 24 4077 4185 -2015
## nn2 2 38 4099 4270 -2012 1 vs 2 6.153 0.9625
summary(nn3)$tTable
## Value Std.Error DF t-value p-value
## int 5.0501 0.2487 568 20.3059 2.494e-69
## A.(Intercept) 134.6841 3.4094 568 39.5034 4.492e-165
## A.blc2 -1.3152 2.7302 568 -0.4817 6.302e-01
## A.blc3 -6.2798 2.7980 568 -2.2444 2.519e-02
## A.famF48 2.9420 4.1540 568 0.7082 4.791e-01
## A.agr0-2 mm 8.8981 4.1339 568 2.1525 3.178e-02
## A.cinza 0.2514 0.1455 568 1.7280 8.453e-02
## A.famF48:agr0-2 mm -11.7084 5.8452 568 -2.0031 4.564e-02
## A.famF48:cinza -0.1692 0.1976 568 -0.8564 3.922e-01
## A.agr0-2 mm:cinza -0.7410 0.2045 568 -3.6241 3.160e-04
## A.famF48:agr0-2 mm:cinza 0.6993 0.2791 568 2.5057 1.250e-02
## d50.(Intercept) 30.0030 0.5851 568 51.2756 2.803e-215
## d50.blc2 0.4303 0.8261 568 0.5209 6.026e-01
## d50.blc3 2.1504 0.8315 568 2.5862 9.953e-03
## S.(Intercept) 5.6847 0.1648 568 34.4923 1.897e-141
## S.blc2 -0.1759 0.2262 568 -0.7777 4.371e-01
## S.blc3 0.1189 0.2316 568 0.5131 6.081e-01
VarCorr(nn3)
## parcela = pdLogChol(list(A ~ 1,d50 ~ 1,S ~ 1))
## Variance StdDev Corr
## A.(Intercept) 80.9886 8.9994 A.(In) d50.(I
## d50.(Intercept) 8.8802 2.9800 -0.273
## S.(Intercept) 0.4895 0.6996 -0.413 0.343
## Residual 9.1305 3.0217
## plot(nn3)
## par(mfrow=c(1,3)); apply(ranef(nn3), 2, qqnorm); layout(1)
## intervals(nn3)
## pairs(nn3)
##-----------------------------------------------------------------------------
## Fazendo a predição.
pred <- expand.grid(blc="1",
fam=levels(mara$fam),
agr=levels(mara$agr),
cinza=seq(0,50, l=2),
dias=seq(0,60,l=30))
pred$y <- predict(nn3, newdata=pred, level=0)
## xyplot(y~dias|fam,
## groups=interaction(cinza, agr, sep=": "),
## data=pred, type=c("l"),
## auto.key=list(columns=4, points=FALSE, lines=TRUE))
xyplot(y~dias|fam*agr, groups=cinza, data=pred, type=c("l"))
## xyplot(altura~dias|fam*agr, groups=cinza, data=marag, xlim=c(0,70))+
## as.layer(xyplot(y~dias|fam*agr, groups=cinza, data=pred,
## type=c("l")))
## xyplot(altura~dias|fam*agr, data=subset(marag, blc=="1"), xlim=c(0,70))+
## as.layer(xyplot(y~dias|fam*agr, groups=cinza, data=pred,
## type=c("l")))
##-----------------------------------------------------------------------------
## Predição.
pred <- expand.grid(blc="1",
fam=levels(mara$fam),
agr=levels(mara$agr),
cinza=seq(0,50,l=20),
dias=seq(0,60,l=30))
pred$y <- predict(nn3, newdata=pred, level=0)
## colr <- brewer.pal(11, "Spectral")
## colr <- c("white","black")
## colr <- colorRampPalette(colr, space="rgb")
## Largura da página da coluna da revista (8 cm = 3.150 pol) e da folha
## (17 cm = 6.690 pol).
colr <- c("white","black")
colr <- colorRampPalette(colr, space="rgb")
wireframe(y~dias+cinza|fam*agr, data=pred,
scales=list(arrows=FALSE), drape=TRUE, zlim=c(0, 155),
zlab=list(expression(Altura~da~planta~(cm)), rot=90,
hjust=0.5, vjust=-0.2),
xlab=list(expression(Dias), vjust=1, rot=30),
ylab=list(expression(Cinza~(t~ha^{-1})), vjust=1, rot=-40),
strip=strip.combined, par.strip.text=list(lines=0.6),
par.settings=list(
## regions=list(col=gray.colors(100)),
regions=list(col=colr(100)),
layout.widths=list(xlab.axis.padding=3)))
##-----------------------------------------------------------------------------
## Estimativas pontuais das assíntotas.
summ <- summary(nn3)$tTable
summ[grep("cinza", rownames(summ)),]
## Value Std.Error DF t-value p-value
## A.cinza 0.2514 0.1455 568 1.7280 0.084530
## A.famF48:cinza -0.1692 0.1976 568 -0.8564 0.392154
## A.agr0-2 mm:cinza -0.7410 0.2045 568 -3.6241 0.000316
## A.famF48:agr0-2 mm:cinza 0.6993 0.2791 568 2.5057 0.012500
coefs <- fixef(nn3)
vcovs <- vcov(nn3)
## names(coefs)
idA <- grep("^A\\.", names(coefs))
coefsA <- coefs[idA]
vcovsA <- vcovs[idA, idA]
## length(coefsA); length(coef(m0))
X0 <- LSmatrix(m0, effect=c("agr","fam"), at=list(cinza=0))
b0 <- X0%*%coefsA
v0 <- sqrt(diag(X0%*%vcovsA%*%t(X0)))
ic0 <- sapply(1:4, function(i) b0[i]+c(-1,0,1)*1.96*v0[i])
ic0 <- t(ic0)
colnames(ic0) <- c("Lower","Estimate","Upper")
rownames(ic0) <- paste(c(outer(levels(mara$fam), levels(mara$agr),
paste)), "cinza 0")
X48 <- LSmatrix(m0, effect=c("agr","fam"), at=list(cinza=48))
b48 <- X48%*%coefsA
v48 <- sqrt(diag(X48%*%vcovsA%*%t(X48)))
ic48 <- sapply(1:4, function(i) b48[i]+c(-1,0,1)*1.96*v48[i])
ic48 <- t(ic48)
colnames(ic48) <- c("Lower","Estimate","Upper")
rownames(ic48) <- paste(c(outer(levels(mara$fam), levels(mara$agr),
paste)), "cinza 48")
ics <- cbind(trat=c(rownames(ic0), rownames(ic48)),
data.frame(rbind(ic0, ic48)))
p1 <-
segplot(reorder(trat, Estimate)~Lower+Upper, data=ics,
draw.bands=FALSE, centers=Estimate,
segments.fun=panel.arrows,
ends="both", angle=90, length=.05,
par.settings=simpleTheme(pch=19, col=1),
ylab=NULL,
xlab=expression(Ganho~em~altura~(cm)),
panel=function(x, y, z, ...) {
panel.abline(h=z, col="grey", lty="dashed")
panel.abline(v=0, col=1, lty="dashed")
panel.segplot(x, y, z, ...)})
## print(p1)
##-----------------------------------------------------------------------------
## Mas o contraste é uma matriz menos a outra, então.
X <- X48-X0
b <- X%*%coefsA
v <- X%*%vcovsA%*%t(X)
t <- b/sqrt(diag(v))
Z <- model.matrix(~fam*agr,
expand.grid(fam=levels(mara$fam),
agr=levels(mara$agr)))
d <- coefs[grep("cinza", names(coefs))]
u <- vcovs[grep("cinza", names(coefs)), grep("cinza", names(coefs))]
d <- c(Z%*%d)
u <- sqrt(diag(Z%*%u%*%t(Z)))
ic <- sapply(1:4, function(i) d[i]+c(-1,0,1)*1.96*u[i])
ic <- t(ic)
colnames(ic) <- c("Lower","Estimate","Upper")
rownames(ic) <- c(outer(levels(mara$fam), levels(mara$agr), paste))
ic <- cbind(trat=rownames(ic), as.data.frame(ic))
p2 <-
segplot(reorder(trat, Estimate)~Lower+Upper, data=ic,
draw.bands=FALSE, centers=Estimate, segments.fun=panel.arrows,
ends="both", angle=90, length=.05,
par.settings=simpleTheme(pch=19, col=1),
ylab=NULL,
xlab=expression(Coeficiente~de~inclinação),
panel=function(x, y, z, ...) {
panel.abline(h=z, col="grey", lty="dashed")
panel.abline(v=0, col=1, lty="dashed")
panel.segplot(x, y, z, ...)
})
## print(p2)
plot(p1, split=c(1,1,1,2), more=TRUE)
grid.text(label="A", 0.025, 0.975)
plot(p2, split=c(1,2,1,2), more=FALSE)
grid.text(label="B", 0.025, 0.475)
##-----------------------------------------------------------------------------
## Ler tabela de dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/maracuja_planta.txt",
header=TRUE, sep="\t")
da$bloc <- factor(da$bloc)
str(da)
## 'data.frame': 168 obs. of 9 variables:
## $ agreg: Factor w/ 2 levels "[0,2]","[4,10]": 2 2 2 2 2 2 2 2 2 2 ...
## $ fam : Factor w/ 2 levels "F29","F48": 1 1 1 1 1 1 1 1 1 1 ...
## $ cinz : num 0 0 0 0 0 0 1.5 1.5 1.5 1.5 ...
## $ bloc : Factor w/ 3 levels "1","2","3": 1 1 2 2 3 3 1 1 2 2 ...
## $ rept : int 1 2 1 2 1 2 1 2 1 2 ...
## $ mfpa : num 43.5 43.7 31.8 23.4 36.4 ...
## $ mspa : num 11.1 11.13 8.24 5.39 9.04 ...
## $ ds : num 1.24 1.25 1.23 1.23 1.27 1.18 1.27 1.28 1.01 1.19 ...
## $ cav : num 6470 7260 9446 9046 8468 ...
u <- unique(da$cinz)
cbind(u, u/1.5, sqrt(u/1.5), log(u/1.5, base=2))
## u
## [1,] 0.0 0 0.000 -Inf
## [2,] 1.5 1 1.000 0
## [3,] 3.0 2 1.414 1
## [4,] 6.0 4 2.000 2
## [5,] 12.0 8 2.828 3
## [6,] 24.0 16 4.000 4
## [7,] 48.0 32 5.657 5
## Doses de cinza em uma escala com distância mais regular entre níveis.
da$cin <- sqrt(da$cinz/1.5)
##-----------------------------------------------------------------------------
## Ver.
p1 <- xyplot(mfpa~cin|fam, groups=agreg, data=da, type=c("p","smooth"),
auto.key=list(columns=2))
p2 <- xyplot(mspa~cin|fam, groups=agreg, data=da, type=c("p","smooth"),
auto.key=list(columns=2))
p3 <- xyplot(ds~cin|fam, groups=agreg, data=da, type=c("p","smooth"),
auto.key=list(columns=2))
p4 <- xyplot(cav~cin|fam, groups=agreg, data=da, type=c("p","smooth"),
auto.key=list(columns=2))
grid.arrange(p1, p2, nrow=2)
grid.arrange(p3, p4, nrow=2)
##-----------------------------------------------------------------------------
## Especificação e ajuste dos modelos em batelada.
resp <- c("mfpa", "mspa", "ds", "cav")
form0 <- lapply(paste(resp, "~bloc+fam*agreg*(cinz+I(cinz^2))"),
as.formula)
names(form0) <- resp
## Ajustes.
m0 <- lapply(form0, lm, data=da)
##-----------------------------------------------------------------------------
## Quadros de anova.
lapply(m0, anova)
## $mfpa
## Analysis of Variance Table
##
## Response: mfpa
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 2104 1052 6.55 0.00185 **
## fam 1 89 89 0.55 0.45842
## agreg 1 4130 4130 25.74 1.1e-06 ***
## cinz 1 1848 1848 11.52 0.00088 ***
## I(cinz^2) 1 426 426 2.66 0.10515
## fam:agreg 1 258 258 1.61 0.20690
## fam:cinz 1 431 431 2.69 0.10320
## fam:I(cinz^2) 1 559 559 3.48 0.06384 .
## agreg:cinz 1 840 840 5.23 0.02350 *
## agreg:I(cinz^2) 1 27 27 0.17 0.68484
## fam:agreg:cinz 1 252 252 1.57 0.21173
## fam:agreg:I(cinz^2) 1 123 123 0.77 0.38307
## Residuals 154 24710 160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $mspa
## Analysis of Variance Table
##
## Response: mspa
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 86 43.1 5.27 0.0061 **
## fam 1 0 0.3 0.04 0.8446
## agreg 1 80 80.3 9.80 0.0021 **
## cinz 1 62 61.8 7.54 0.0067 **
## I(cinz^2) 1 5 5.0 0.61 0.4371
## fam:agreg 1 1 1.2 0.15 0.6972
## fam:cinz 1 21 20.5 2.51 0.1155
## fam:I(cinz^2) 1 21 21.3 2.60 0.1091
## agreg:cinz 1 63 63.0 7.70 0.0062 **
## agreg:I(cinz^2) 1 14 13.6 1.66 0.1993
## fam:agreg:cinz 1 11 11.2 1.37 0.2431
## fam:agreg:I(cinz^2) 1 7 7.4 0.90 0.3437
## Residuals 154 1261 8.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $ds
## Analysis of Variance Table
##
## Response: ds
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 0.044 0.022 5.06 0.00742 **
## fam 1 0.053 0.053 12.20 0.00063 ***
## agreg 1 0.024 0.024 5.50 0.02025 *
## cinz 1 0.815 0.815 187.38 < 2e-16 ***
## I(cinz^2) 1 0.090 0.090 20.58 1.1e-05 ***
## fam:agreg 1 0.013 0.013 2.89 0.09093 .
## fam:cinz 1 0.036 0.036 8.22 0.00473 **
## fam:I(cinz^2) 1 0.032 0.032 7.44 0.00713 **
## agreg:cinz 1 0.026 0.026 5.93 0.01602 *
## agreg:I(cinz^2) 1 0.017 0.017 3.91 0.04990 *
## fam:agreg:cinz 1 0.007 0.007 1.52 0.21968
## fam:agreg:I(cinz^2) 1 0.000 0.000 0.02 0.89143
## Residuals 154 0.670 0.004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $cav
## Analysis of Variance Table
##
## Response: cav
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 4.70e+06 2.35e+06 1.54 0.21665
## fam 1 1.30e+07 1.30e+07 8.52 0.00405 **
## agreg 1 2.71e+08 2.71e+08 177.72 < 2e-16 ***
## cinz 1 5.51e+05 5.51e+05 0.36 0.54852
## I(cinz^2) 1 2.36e+06 2.36e+06 1.55 0.21500
## fam:agreg 1 4.24e+06 4.24e+06 2.78 0.09720 .
## fam:cinz 1 5.82e+05 5.82e+05 0.38 0.53749
## fam:I(cinz^2) 1 1.40e+07 1.40e+07 9.20 0.00284 **
## agreg:cinz 1 2.82e+06 2.82e+06 1.85 0.17541
## agreg:I(cinz^2) 1 1.36e+04 1.36e+04 0.01 0.92471
## fam:agreg:cinz 1 1.91e+07 1.91e+07 12.55 0.00053 ***
## fam:agreg:I(cinz^2) 1 2.92e+05 2.92e+05 0.19 0.66213
## Residuals 154 2.34e+08 1.52e+06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Quadro de estimativas dos efeitos.
lapply(m0, summary)
## $mfpa
##
## Call:
## FUN(formula = X[[1L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.07 -8.11 1.07 7.51 32.86
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.94758 3.46877 11.52 < 2e-16 ***
## bloc2 -2.62321 2.39385 -1.10 0.27487
## bloc3 -8.46589 2.39385 -3.54 0.00054 ***
## famF48 8.48798 4.49937 1.89 0.06111 .
## agreg[4,10] -0.71655 4.49937 -0.16 0.87368
## cinz 1.21207 0.46468 2.61 0.00999 **
## I(cinz^2) -0.01879 0.00948 -1.98 0.04921 *
## famF48:agreg[4,10] -12.40053 6.36307 -1.95 0.05313 .
## famF48:cinz -1.18219 0.65716 -1.80 0.07399 .
## famF48:I(cinz^2) 0.02599 0.01341 1.94 0.05439 .
## agreg[4,10]:cinz -0.64216 0.65716 -0.98 0.33002
## agreg[4,10]:I(cinz^2) 0.00444 0.01341 0.33 0.74109
## famF48:agreg[4,10]:cinz 1.09051 0.92937 1.17 0.24245
## famF48:agreg[4,10]:I(cinz^2) -0.01658 0.01896 -0.87 0.38307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.7 on 154 degrees of freedom
## Multiple R-squared: 0.31, Adjusted R-squared: 0.251
## F-statistic: 5.31 on 13 and 154 DF, p-value: 7.8e-08
##
##
## $mspa
##
## Call:
## FUN(formula = X[[2L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.316 -1.897 0.214 2.066 6.930
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.312995 0.783599 11.88 <2e-16 ***
## bloc2 -0.245000 0.540773 -0.45 0.6511
## bloc3 -1.627500 0.540773 -3.01 0.0031 **
## famF48 1.028492 1.016413 1.01 0.3132
## agreg[4,10] 0.113221 1.016413 0.11 0.9114
## cinz 0.174329 0.104973 1.66 0.0988 .
## I(cinz^2) -0.002196 0.002141 -1.03 0.3066
## famF48:agreg[4,10] -2.029210 1.437426 -1.41 0.1601
## famF48:cinz -0.248126 0.148454 -1.67 0.0967 .
## famF48:I(cinz^2) 0.005485 0.003028 1.81 0.0720 .
## agreg[4,10]:cinz -0.074416 0.148454 -0.50 0.6169
## agreg[4,10]:I(cinz^2) -0.000727 0.003028 -0.24 0.8107
## famF48:agreg[4,10]:cinz 0.257012 0.209946 1.22 0.2228
## famF48:agreg[4,10]:I(cinz^2) -0.004068 0.004283 -0.95 0.3437
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.86 on 154 degrees of freedom
## Multiple R-squared: 0.228, Adjusted R-squared: 0.163
## F-statistic: 3.49 on 13 and 154 DF, p-value: 8.92e-05
##
##
## $ds
##
## Call:
## FUN(formula = X[[3L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27187 -0.03054 0.00381 0.03321 0.17906
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.26e+00 1.81e-02 69.96 <2e-16 ***
## bloc2 -2.10e-02 1.25e-02 -1.68 0.0942 .
## bloc3 1.86e-02 1.25e-02 1.50 0.1367
## famF48 4.52e-02 2.34e-02 1.93 0.0556 .
## agreg[4,10] -3.50e-02 2.34e-02 -1.50 0.1369
## cinz -8.08e-03 2.42e-03 -3.34 0.0011 **
## I(cinz^2) 9.00e-05 4.94e-05 1.82 0.0701 .
## famF48:agreg[4,10] -5.85e-02 3.31e-02 -1.76 0.0796 .
## famF48:cinz -9.29e-03 3.42e-03 -2.72 0.0074 **
## famF48:I(cinz^2) 1.41e-04 6.98e-05 2.03 0.0446 *
## agreg[4,10]:cinz 5.06e-03 3.42e-03 1.48 0.1413
## agreg[4,10]:I(cinz^2) -9.08e-05 6.98e-05 -1.30 0.1953
## famF48:agreg[4,10]:cinz 2.20e-03 4.84e-03 0.46 0.6493
## famF48:agreg[4,10]:I(cinz^2) -1.35e-05 9.87e-05 -0.14 0.8914
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.066 on 154 degrees of freedom
## Multiple R-squared: 0.633, Adjusted R-squared: 0.602
## F-statistic: 20.4 on 13 and 154 DF, p-value: <2e-16
##
##
## $cav
##
## Call:
## FUN(formula = X[[4L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2658 -767 -123 801 3246
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6018.438 337.907 17.81 <2e-16 ***
## bloc2 -10.196 233.195 -0.04 0.9652
## bloc3 349.768 233.195 1.50 0.1357
## famF48 -1070.239 438.303 -2.44 0.0157 *
## agreg[4,10] 1370.280 438.303 3.13 0.0021 **
## cinz -130.181 45.267 -2.88 0.0046 **
## I(cinz^2) 2.134 0.923 2.31 0.0222 *
## famF48:agreg[4,10] 1934.023 619.854 3.12 0.0022 **
## famF48:cinz 201.080 64.017 3.14 0.0020 **
## famF48:I(cinz^2) -3.205 1.306 -2.45 0.0152 *
## agreg[4,10]:cinz 73.283 64.017 1.14 0.2541
## agreg[4,10]:I(cinz^2) -0.317 1.306 -0.24 0.8086
## famF48:agreg[4,10]:cinz -122.468 90.534 -1.35 0.1781
## famF48:agreg[4,10]:I(cinz^2) 0.809 1.847 0.44 0.6621
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1230 on 154 degrees of freedom
## Multiple R-squared: 0.586, Adjusted R-squared: 0.551
## F-statistic: 16.8 on 13 and 154 DF, p-value: <2e-16
##-----------------------------------------------------------------------------
## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m0[["mfpa"]]); layout(1) ## ok!
par(mfrow=c(2,2)); plot(m0[["mspa"]]); layout(1) ## ok!
par(mfrow=c(2,2)); plot(m0[["ds"]]); layout(1) ## bom.
par(mfrow=c(2,2)); plot(m0[["cav"]]); layout(1) ## ok!
## par(mfrow=c(2,2))
## lapply(m0, plot, which=1)
## layout(1)
## par(mfrow=c(2,2))
## lapply(m0, plot, which=2)
## layout(1)
## par(mfrow=c(2,2))
## lapply(m0, plot, which=3)
## layout(1)
##-----------------------------------------------------------------------------
## Modelos após abandono de termos não significativos.
form1 <- list(mfpa=mfpa~bloc+agreg*cinz,
mspa=mspa~bloc+agreg*cinz,
ds=ds~bloc+(fam+agreg+cinz)^3+(fam+agreg)*I(cinz^2),
cav=cav~bloc+fam*agreg*(cinz+I(cinz^2)))
m1 <- lapply(form1, lm, data=da, contrast=list(bloc=contr.sum))
lapply(m1, anova)
## $mfpa
## Analysis of Variance Table
##
## Response: mfpa
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 2104 1052 6.34 0.0022 **
## agreg 1 4130 4130 24.89 1.5e-06 ***
## cinz 1 1848 1848 11.14 0.0010 **
## agreg:cinz 1 840 840 5.06 0.0258 *
## Residuals 162 26875 166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $mspa
## Analysis of Variance Table
##
## Response: mspa
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 86 43.1 5.21 0.0064 **
## agreg 1 80 80.3 9.69 0.0022 **
## cinz 1 62 61.8 7.46 0.0070 **
## agreg:cinz 1 63 63.0 7.61 0.0065 **
## Residuals 162 1342 8.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $ds
## Analysis of Variance Table
##
## Response: ds
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 0.044 0.022 5.10 0.0072 **
## fam 1 0.053 0.053 12.27 0.0006 ***
## agreg 1 0.024 0.024 5.54 0.0199 *
## cinz 1 0.815 0.815 188.57 < 2e-16 ***
## I(cinz^2) 1 0.090 0.090 20.71 1.1e-05 ***
## fam:agreg 1 0.013 0.013 2.91 0.0899 .
## fam:cinz 1 0.036 0.036 8.27 0.0046 **
## agreg:cinz 1 0.026 0.026 5.97 0.0157 *
## fam:I(cinz^2) 1 0.032 0.032 7.49 0.0069 **
## agreg:I(cinz^2) 1 0.017 0.017 3.93 0.0492 *
## fam:agreg:cinz 1 0.007 0.007 1.53 0.2182
## Residuals 155 0.670 0.004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $cav
## Analysis of Variance Table
##
## Response: cav
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 4.70e+06 2.35e+06 1.54 0.21665
## fam 1 1.30e+07 1.30e+07 8.52 0.00405 **
## agreg 1 2.71e+08 2.71e+08 177.72 < 2e-16 ***
## cinz 1 5.51e+05 5.51e+05 0.36 0.54852
## I(cinz^2) 1 2.36e+06 2.36e+06 1.55 0.21500
## fam:agreg 1 4.24e+06 4.24e+06 2.78 0.09720 .
## fam:cinz 1 5.82e+05 5.82e+05 0.38 0.53749
## fam:I(cinz^2) 1 1.40e+07 1.40e+07 9.20 0.00284 **
## agreg:cinz 1 2.82e+06 2.82e+06 1.85 0.17541
## agreg:I(cinz^2) 1 1.36e+04 1.36e+04 0.01 0.92471
## fam:agreg:cinz 1 1.91e+07 1.91e+07 12.55 0.00053 ***
## fam:agreg:I(cinz^2) 1 2.92e+05 2.92e+05 0.19 0.66213
## Residuals 154 2.34e+08 1.52e+06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova entre modelos sequenciais.
sapply(names(m1), simplify=FALSE,
function(i){
anova(m0[[i]],m1[[i]])
})
## $mfpa
## Analysis of Variance Table
##
## Model 1: mfpa ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Model 2: mfpa ~ bloc + agreg * cinz
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 154 24710
## 2 162 26875 -8 -2165 1.69 0.11
##
## $mspa
## Analysis of Variance Table
##
## Model 1: mspa ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Model 2: mspa ~ bloc + agreg * cinz
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 154 1261
## 2 162 1342 -8 -80.6 1.23 0.29
##
## $ds
## Analysis of Variance Table
##
## Model 1: ds ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Model 2: ds ~ bloc + (fam + agreg + cinz)^3 + (fam + agreg) * I(cinz^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 154 0.67
## 2 155 0.67 -1 -8.13e-05 0.02 0.89
##
## $cav
## Analysis of Variance Table
##
## Model 1: cav ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Model 2: cav ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 154 2.34e+08
## 2 154 2.34e+08 0 1.79e-07
##-----------------------------------------------------------------------------
## Valores preditos.
gridlist <- list(bloc="1",
fam=levels(da$fam),
agreg=levels(da$agreg),
cin=seq(0,6,l=100),
cinz=seq(0,50,l=100))
m1 <- lapply(m1,
function(i){
predvars <- attr(terms(i), "term.labels")
pred <- do.call(
expand.grid,
gridlist[predvars[!sapply(gridlist[predvars],
is.null)]])
i$newdata <- pred
i
})
mypredict <- function(m){
cbind(m$newdata,
predict(m, newdata=m$newdata,
interval="confidence")-
coef(m)["bloc1"])
}
all.pred <- lapply(m1, mypredict)
## str(all.pred)
## lapply(all.pred, names)
##-----------------------------------------------------------------------------
## Gráficos.
xlab <- expression(Cinza~(t~ha^{-1}))
ylab <- list(expression(Massa~fresca~de~parte~aérea~(g)),
expression(Massa~seca~de~parte~aérea~(g)),
expression(Densidade~do~solo~(Mg~t^{-1})),
expression(Água~consumida~no~ciclo~(mL)))
names(ylab) <- names(m1)
##-----------------------------------------------------------------------------
## Mfpa.
p1 <-
xyplot(mfpa~cinz, groups=agreg, data=da,
ylab=ylab[["mfpa"]], xlab=xlab,
strip=strip.custom(bg="gray90"),
auto.key=list(columns=2, points=FALSE, lines=TRUE,
title="Classe de agregado (mm)", cex.title=1.2))+
as.layer(with(all.pred[["mfpa"]],
xyplot(fit~cinz, groups=agreg, type="l",
ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)))
##-----------------------------------------------------------------------------
## Mspa.
p2 <-
xyplot(mspa~cinz, groups=agreg, data=da,
ylab=ylab[["mspa"]], xlab=xlab,
strip=strip.custom(bg="gray90"),
auto.key=list(columns=2, points=FALSE, lines=TRUE,
title="Classe de agregado (mm)", cex.title=1.2))+
as.layer(with(all.pred[["mspa"]],
xyplot(fit~cinz, groups=agreg, type="l",
ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)))
##-----------------------------------------------------------------------------
## Ds.
p3 <-
xyplot(ds~cinz|fam, groups=agreg, data=da, layout=c(NA,1),
ylab=ylab[["ds"]], xlab=xlab,
strip=strip.custom(bg="gray90"),
auto.key=list(columns=2, points=FALSE, lines=TRUE,
title="Classe de agregado (mm)", cex.title=1.2))+
as.layer(with(all.pred[["ds"]],
xyplot(fit~cinz|fam, groups=agreg, type="l",
ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)))
##-----------------------------------------------------------------------------
## Cav.
p4 <-
xyplot(cav~cinz|fam, groups=agreg, data=da, layout=c(NA,1),
ylab=ylab[["cav"]], xlab=xlab,
strip=strip.custom(bg="gray90"),
auto.key=list(columns=2, points=FALSE, lines=TRUE,
title="Classe de agregado (mm)", cex.title=1.2))+
as.layer(with(all.pred[["cav"]],
xyplot(fit~cinz|fam, groups=agreg, type="l",
ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)))
##-----------------------------------------------------------------------------
## Arranjo 1.
grid.arrange(p1, p2, ncol=2)
grid.text(label="A", x=0.025, y=0.975)
grid.text(label="B", x=0.5+0.025, y=0.975)
##-----------------------------------------------------------------------------
## Arranjo 2.
grid.arrange(p3, p4, nrow=2)
grid.text(label="A", x=0.025, y=0.975)
grid.text(label="B", x=0.025, y=0.975-0.5)
