Grasiela Bruzamarello Tognon
Walmes Marques Zeviani
Francine Lorena Cuquel
Experimento feito com estacas de B. tridentata tendo como fator experimental o tipo de solução para conservação das estacas (4 níveis mais o controle). Cada tratamento foi aplicado em 10 estacas e as variáveis foram observadas ao longo do tempo após aplicação nessas estacas. As medidas são indicadores de qualidade baseados em coloração, peso e nota relacionada à viabilidade comercial. O mesmo experimento foi repetido nas 4 estações do ano (1 rodada de cada estação).
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "gridExtra", "nlme", "survival",
"doBy", "plyr", "multcomp", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra gridExtra nlme survival doBy
## TRUE TRUE TRUE TRUE TRUE TRUE
## plyr multcomp 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.1 multcomp_1.3-6 TH.data_1.0-3 mvtnorm_0.9-9997
## [5] plyr_1.8.1 doBy_4.5-10 MASS_7.3-34 survival_2.37-7
## [9] nlme_3.1-117 gridExtra_0.9.1 latticeExtra_0.6-26 RColorBrewer_1.0-5
## [13] lattice_0.20-29 knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_1.0 lme4_1.1-7 Matrix_1.1-4 methods_3.1.1
## [6] minqa_1.2.3 nloptr_1.0.4 Rcpp_0.11.0 sandwich_2.3-0 stringr_0.6.2
## [11] tools_3.1.1 zoo_1.7-10
citation()
##
## To cite R in publications use:
##
## R Core Team (2014). R: A language and environment for statistical computing. R
## Foundation for Statistical Computing, Vienna, Austria. URL
## http://www.R-project.org/.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {R: A Language and Environment for Statistical Computing},
## author = {{R Core Team}},
## organization = {R Foundation for Statistical Computing},
## address = {Vienna, Austria},
## year = {2014},
## url = {http://www.R-project.org/},
## }
##
## We have invested a lot of time and effort in creating R, please cite it when
## using it for data analysis. See also 'citation("pkgname")' for citing R
## packages.
##-----------------------------------------------------------------------------
## Definições gráficas trellis.
## display.brewer.all()
## mycol <- brewer.pal(6, "Set1")
mycol <- c(brewer.pal(6, "Dark2"), brewer.pal(6, "Set1"))
## names(trellis.par.get())
## dput(trellis.par.get())
ps <- list(box.rectangle=list(col=1, fill=c("gray70")),
box.umbrella=list(col=1, lty=1),
dot.symbol=list(col=1),
dot.line=list(col=1, lty=3),
plot.symbol=list(col=1, cex=0.7),
plot.line=list(col=1),
plot.polygon=list(col="gray80"),
superpose.line=list(col=mycol),
superpose.symbol=list(col=mycol),
superpose.polygon=list(col=mycol),
strip.background=list(col=c("gray90","gray50"))
)
## show.settings()
## show.settings(ps)
##-----------------------------------------------------------------------------
## Ler.
da <- read.table("estacoes_tridentata.txt", header=TRUE, sep="\t")
## str(da)
da$trat <- factor(da$trat)
levels(da$trat) <- c("Controle",
"Metabissulfito de sódio",
"Ácido cítrico",
"Metabissulfito de sódio e Sacarose",
"Ácido cítrico e Sacarose")
da$trat <- factor(da$trat, levels=c("Controle",
"Ácido cítrico",
"Metabissulfito de sódio",
"Ácido cítrico e Sacarose",
"Metabissulfito de sódio e Sacarose"))
## str(da)
da$ue <- with(da, interaction(est, trat, rept, drop=TRUE))
## nlevels(da$ue)
##-----------------------------------------------------------------------------
## a.
da$y <- da$a
useOuterStrips(
xyplot(y~dias|trat+est, groups=rept, data=da,
type=c("p","a")))
m0 <- lme(y~est*trat*(dias+I(dias^2)), random=~1+dias|ue, data=da,
na.action=na.omit, method="ML")
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0); f <- fitted(m0); o <- na.omit(m0$data$y)
par(mfrow=c(2,2))
qqnorm(r)
plot(r~f, type="n")
panel.smooth(f, r)
plot(sqrt(abs(r))~f, type="n")
panel.smooth(f, sqrt(abs(r)))
plot(o~f, type="n")
panel.smooth(f, o)
layout(1)
##-----------------------------------------------------------------------------
## Medidas de ajuste e efeito.
VarCorr(m0)
## ue = pdLogChol(1 + dias)
## Variance StdDev Corr
## (Intercept) 1.881470 1.37167 (Intr)
## dias 0.007195 0.08482 -0.214
## Residual 0.932682 0.96575
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 755 1606.9 <.0001
## est 3 180 1800.8 <.0001
## trat 4 180 1.1 0.3430
## dias 1 755 20.1 <.0001
## I(dias^2) 1 755 5.8 0.0162
## est:trat 12 180 1.6 0.0875
## est:dias 3 755 3.7 0.0123
## est:I(dias^2) 3 755 8.7 <.0001
## trat:dias 4 755 2.8 0.0256
## trat:I(dias^2) 4 755 1.7 0.1506
## est:trat:dias 12 755 2.1 0.0167
## est:trat:I(dias^2) 12 755 1.5 0.1223
##-----------------------------------------------------------------------------
## Modelo final.
m1 <- m0
## anova(m1)
##-----------------------------------------------------------------------------
## Representando o resultado.
## p0 <- with(da,
## expand.grid(est=levels(est),
## trat=levels(trat),
## dias=seq(0,14,by=0.5)))
p0 <- ddply(da, .(est, trat), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
p0$y <- predict(m0, newdata=p0, level=0)
## xyplot(y~dias|est+trat, data=da)+
## as.layer(xyplot(y~dias|est+trat, data=p0,
## type="l"))
##-----------------------------------------------------------------------------
## Predição com bandas de confiança.
## p1 <- with(da,
## expand.grid(est=levels(est),
## trat=levels(trat),
## dias=seq(0,14,by=0.5)))
p1 <- p0
p1$y <- predict(m1, newdata=p1, level=0)
F <- model.matrix(formula(m1), data=p1)
U <- chol(vcov(m1))
se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
tval <- qt(p=c(lwr=0.025, fit=0.5, upr=0.975), df=nrow(da)-ncol(F))
me <- outer(se, tval, "*")
p1 <- cbind(p1, sweep(me, 1, p1$y, "+"))
str(p1)
## 'data.frame': 740 obs. of 7 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat: Factor w/ 5 levels "Controle","Ácido cítrico",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dias: num 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 ...
## $ y : atomic -6.75 -6.87 -6.97 -7.07 -7.17 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num -7.82 -7.9 -7.97 -8.05 -8.12 ...
## $ fit : num -6.75 -6.87 -6.97 -7.07 -7.17 ...
## $ upr : num -5.69 -5.83 -5.97 -6.1 -6.21 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias|est, groups=trat, data=da, type="p", auto.key=TRUE)+
as.layer(xyplot(fit~dias|est, groups=trat, data=p1,
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose))
xyplot(fit~dias|est, groups=trat, data=p1,
auto.key=list(
## corner=c(0.02,0.98),
points=FALSE, lines=TRUE),
xlab="Período de observação (dias)",
ylab="Medida de colocação a",
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose)
##-----------------------------------------------------------------------------
## L.
da$y <- da$L
useOuterStrips(
xyplot(y~dias|trat+est, groups=rept, data=da,
type=c("p","a")))
m0 <- lme(y~est*trat*dias, random=~1+dias|ue, data=da,
na.action=na.omit, method="ML")
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0); f <- fitted(m0); o <- na.omit(m0$data$y)
par(mfrow=c(2,2))
qqnorm(r)
plot(r~f, type="n")
panel.smooth(f, r)
plot(sqrt(abs(r))~f, type="n")
panel.smooth(f, sqrt(abs(r)))
plot(o~f, type="n")
panel.smooth(f, o)
layout(1)
##-----------------------------------------------------------------------------
## Medidas de ajuste e efeito.
VarCorr(m0)
## ue = pdLogChol(1 + dias)
## Variance StdDev Corr
## (Intercept) 4.684e+00 2.1643196 (Intr)
## dias 1.146e-08 0.0001071 0
## Residual 7.905e+00 2.8115282
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 778 47678 <.0001
## est 3 180 4 0.0064
## trat 4 180 3 0.0349
## dias 1 778 38 <.0001
## est:trat 12 180 1 0.5443
## est:dias 3 778 7 0.0002
## trat:dias 4 778 1 0.2601
## est:trat:dias 12 778 2 0.0563
##-----------------------------------------------------------------------------
## Modelo final.
m1 <- update(m0, y~est*dias+trat)
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 44 5257 5473 -2585
## m1 2 16 5239 5318 -2604 1 vs 2 37.67 0.1048
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 794 46444 <.0001
## est 3 192 4 0.0073
## dias 1 794 38 <.0001
## trat 4 192 3 0.0391
## est:dias 3 794 7 0.0002
##-----------------------------------------------------------------------------
## Representando o resultado.
## p0 <- with(da,
## expand.grid(est=levels(est),
## trat=levels(trat),
## dias=seq(0,14,by=0.5)))
p0 <- ddply(da, .(est, trat), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
p0$y <- predict(m0, newdata=p0, level=0)
p1 <- p0
p1$y <- predict(m1, newdata=p1, level=0)
## xyplot(y~dias|est+trat, data=da)+
## as.layer(xyplot(y~dias|est+trat, data=p0,
## type="l"))+
## as.layer(xyplot(y~dias|est+trat, data=p1, col=2,
## type="l"))
## xyplot(y~dias, groups=est, data=da, auto.key=TRUE)+
## as.layer(xyplot(y~dias, groups=est,
## data=subset(p0, trat=="Controle"),
## type="l"))
##-----------------------------------------------------------------------------
## Predição com bandas de confiança.
## p1 <- with(da,
## expand.grid(est=levels(est),
## dias=seq(0,14,by=0.5)))
p1$y <- predict(m1, newdata=p1, level=0)
F <- model.matrix(formula(m1), data=p1)
U <- chol(vcov(m1))
se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
tval <- qt(p=c(lwr=0.025, fit=0.5, upr=0.975), df=nrow(da)-ncol(F))
me <- outer(se, tval, "*")
p1 <- cbind(p1, sweep(me, 1, p1$y, "+"))
str(p1)
## 'data.frame': 740 obs. of 7 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat: Factor w/ 5 levels "Controle","Ácido cítrico",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dias: num 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 ...
## $ y : atomic 39.4 39.4 39.4 39.4 39.5 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 38.3 38.3 38.3 38.3 38.4 ...
## $ fit : num 39.4 39.4 39.4 39.4 39.5 ...
## $ upr : num 40.5 40.5 40.5 40.5 40.5 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias|est, groups=trat, data=da, type="p", auto.key=TRUE)+
as.layer(xyplot(fit~dias|est, groups=trat, data=p1,
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose))
xyplot(fit~dias|est, groups=trat, data=p1,
auto.key=list(
## corner=c(0.02,0.98),
points=FALSE, lines=TRUE),
xlab="Período de observação (dias)",
ylab="Medida de colocação L",
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose)
##-----------------------------------------------------------------------------
## b.
da$y <- da$b
useOuterStrips(
xyplot(y~dias|trat+est, groups=rept, data=da,
type=c("p","a")))
m0 <- lme(y~est*trat*dias, random=~1+dias|ue, data=da,
na.action=na.omit, method="ML")
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0); f <- fitted(m0); o <- na.omit(m0$data$y)
par(mfrow=c(2,2))
qqnorm(r)
plot(r~f, type="n")
panel.smooth(f, r)
plot(sqrt(abs(r))~f, type="n")
panel.smooth(f, sqrt(abs(r)))
plot(o~f, type="n")
panel.smooth(f, o)
layout(1)
##-----------------------------------------------------------------------------
## Medidas de ajuste e efeito.
VarCorr(m0)
## ue = pdLogChol(1 + dias)
## Variance StdDev Corr
## (Intercept) 7.87067 2.805 (Intr)
## dias 0.09062 0.301 -0.326
## Residual 4.91327 2.217
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 780 4908 <.0001
## est 3 180 4 0.0101
## trat 4 180 2 0.2029
## dias 1 780 6 0.0145
## est:trat 12 180 2 0.0128
## est:dias 3 780 7 0.0002
## trat:dias 4 780 3 0.0293
## est:trat:dias 12 780 2 0.1108
##-----------------------------------------------------------------------------
## Modelo final.
## m1 <- update(m0, y~est*dias)
## anova(m0, m1)
m1 <- m0
## anova(m1)
##-----------------------------------------------------------------------------
## Representando o resultado.
## p0 <- with(da,
## expand.grid(est=levels(est),
## trat=levels(trat),
## dias=seq(0,14,by=0.5)))
p0 <- ddply(da, .(est, trat), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
p0$y <- predict(m0, newdata=p0, level=0)
p1 <- p0
p1$y <- predict(m1, newdata=p1, level=0)
## xyplot(y~dias|est+trat, data=da)+
## as.layer(xyplot(y~dias|est+trat, data=p0,
## type="l"))+
## as.layer(xyplot(y~dias|est+trat, data=p1, col=2,
## type="l"))
## xyplot(y~dias, groups=est, data=da, auto.key=TRUE)+
## as.layer(xyplot(y~dias, groups=est,
## data=subset(p0, trat=="Controle"),
## type="l"))
##-----------------------------------------------------------------------------
## Predição com bandas de confiança.
## p1 <- with(da,
## expand.grid(est=levels(est),
## dias=seq(0,14,by=0.5)))
p1$y <- predict(m1, newdata=p1, level=0)
F <- model.matrix(formula(m1), data=p1)
U <- chol(vcov(m1))
se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
tval <- qt(p=c(lwr=0.025, fit=0.5, upr=0.975), df=nrow(da)-ncol(F))
me <- outer(se, tval, "*")
p1 <- cbind(p1, sweep(me, 1, p1$y, "+"))
str(p1)
## 'data.frame': 740 obs. of 7 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat: Factor w/ 5 levels "Controle","Ácido cítrico",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dias: num 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 ...
## $ y : atomic 11.9 11.9 12 12 12 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 9.89 9.95 10 10.05 10.1 ...
## $ fit : num 11.9 11.9 12 12 12 ...
## $ upr : num 13.9 13.9 13.9 13.9 13.9 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias|est, groups=trat, data=da, type="p", auto.key=TRUE)+
as.layer(xyplot(fit~dias|est, groups=trat, data=p1,
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose))
xyplot(fit~dias|est, groups=trat, data=p1,
auto.key=list(
## corner=c(0.02,0.98),
points=FALSE, lines=TRUE),
xlab="Período de observação (dias)",
ylab="Medida de colocação b",
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose)
##-----------------------------------------------------------------------------
## AE.
da$y <- da$AE
useOuterStrips(
xyplot(y~dias|trat+est, groups=rept, data=da,
type=c("p","a")))
m0 <- lme(y~est*trat*dias, random=~1+dias|ue, data=da,
na.action=na.omit, method="ML")
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0); f <- fitted(m0); o <- na.omit(m0$data$y)
par(mfrow=c(2,2))
qqnorm(r)
plot(r~f, type="n")
panel.smooth(f, r)
plot(sqrt(abs(r))~f, type="n")
panel.smooth(f, sqrt(abs(r)))
plot(o~f, type="n")
panel.smooth(f, o)
layout(1)
##-----------------------------------------------------------------------------
## Medidas de ajuste e efeito.
VarCorr(m0)
## ue = pdLogChol(1 + dias)
## Variance StdDev Corr
## (Intercept) 11.3110 3.3632 (Intr)
## dias 0.0459 0.2142 -0.141
## Residual 9.9258 3.1505
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 780 8467 <.0001
## est 3 180 1 0.3135
## trat 4 180 2 0.1082
## dias 1 780 15 0.0001
## est:trat 12 180 2 0.0905
## est:dias 3 780 6 0.0005
## trat:dias 4 780 2 0.1019
## est:trat:dias 12 780 2 0.0862
##-----------------------------------------------------------------------------
## Modelo final.
## m1 <- update(m0, y~(est+trat+dias)^2)
## anova(m0, m1)
m1 <- m0
## anova(m1)
##-----------------------------------------------------------------------------
## Representando o resultado.
## p0 <- with(da,
## expand.grid(est=levels(est),
## trat=levels(trat),
## dias=seq(0,14,by=0.5)))
p0 <- ddply(da, .(est, trat), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
p0$y <- predict(m0, newdata=p0, level=0)
p1 <- p0
p1$y <- predict(m1, newdata=p1, level=0)
## xyplot(y~dias|est+trat, data=da)+
## as.layer(xyplot(y~dias|est+trat, data=p0,
## type="l"))+
## as.layer(xyplot(y~dias|est+trat, data=p1, col=2,
## type="l"))
## xyplot(y~dias, groups=est, data=da, auto.key=TRUE)+
## as.layer(xyplot(y~dias, groups=est,
## data=subset(p0, trat=="Controle"),
## type="l"))
##-----------------------------------------------------------------------------
## Predição com bandas de confiança.
## p1 <- with(da,
## expand.grid(est=levels(est),
## dias=seq(0,14,by=0.5)))
p1$y <- predict(m1, newdata=p1, level=0)
F <- model.matrix(formula(m1), data=p1)
U <- chol(vcov(m1))
se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
tval <- qt(p=c(lwr=0.025, fit=0.5, upr=0.975), df=nrow(da)-ncol(F))
me <- outer(se, tval, "*")
p1 <- cbind(p1, sweep(me, 1, p1$y, "+"))
str(p1)
## 'data.frame': 740 obs. of 7 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat: Factor w/ 5 levels "Controle","Ácido cítrico",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dias: num 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 ...
## $ y : atomic 21.2 21.3 21.3 21.3 21.4 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 18.7 18.7 18.8 18.9 18.9 ...
## $ fit : num 21.2 21.3 21.3 21.3 21.4 ...
## $ upr : num 23.8 23.8 23.8 23.8 23.8 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias|est, groups=trat, data=da, type="p", auto.key=TRUE)+
as.layer(xyplot(fit~dias|est, groups=trat, data=p1,
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose))
xyplot(fit~dias|est, groups=trat, data=p1,
auto.key=list(
## corner=c(0.02,0.98),
points=FALSE, lines=TRUE),
xlab="Período de observação (dias)",
ylab="Medida de colocação E",
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose)
##-----------------------------------------------------------------------------
## Peso.
da$y <- (da$peso)
useOuterStrips(
xyplot(y~dias|trat+est, groups=rept, data=subset(da, !is.na(peso)),
type=c("o")))
m0 <- lme(y~est*trat*(dias+I(dias^2)), random=~1+dias|ue, data=da,
na.action=na.omit, method="ML")
m1 <- lme(y~est*trat*(dias), random=~1+dias|ue, data=da,
na.action=na.omit, method="ML")
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 64 2498 2799 -1185
## m1 2 44 2750 2957 -1331 1 vs 2 292.2 <.0001
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0); f <- fitted(m0); o <- na.omit(m0$data$y)
par(mfrow=c(2,2))
qqnorm(r)
plot(r~f, type="n")
panel.smooth(f, r)
plot(sqrt(abs(r))~f, type="n")
panel.smooth(f, sqrt(abs(r)))
plot(o~f, type="n")
panel.smooth(f, o)
layout(1)
##-----------------------------------------------------------------------------
## Medidas de ajuste e efeito.
VarCorr(m0)
## ue = pdLogChol(1 + dias)
## Variance StdDev Corr
## (Intercept) 37.53845 6.1269 (Intr)
## dias 0.03797 0.1948 -0.332
## Residual 0.11104 0.3332
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 582 1613.5 <.0001
## est 3 180 19.4 <.0001
## trat 4 180 0.8 0.5445
## dias 1 582 63.1 <.0001
## I(dias^2) 1 582 257.4 <.0001
## est:trat 12 180 0.8 0.6719
## est:dias 3 582 6.1 0.0004
## est:I(dias^2) 3 582 10.8 <.0001
## trat:dias 4 582 2.4 0.0453
## trat:I(dias^2) 4 582 2.6 0.0354
## est:trat:dias 12 582 1.2 0.2951
## est:trat:I(dias^2) 12 582 3.1 0.0003
##-----------------------------------------------------------------------------
## Modelo final.
m1 <- m0
## anova(m1)
##-----------------------------------------------------------------------------
## Representando o resultado.
## p0 <- with(da,
## expand.grid(est=levels(est),
## trat=levels(trat),
## dias=seq(0,14,by=0.5)))
p0 <- ddply(da, .(est, trat), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
p0$y <- predict(m0, newdata=p0, level=0)
p1 <- p0
p1$y <- predict(m1, newdata=p1, level=0)
## xyplot(y~dias|est+trat, data=da)+
## as.layer(xyplot(y~dias|est+trat, data=p0,
## type="l"))
##-----------------------------------------------------------------------------
## Predição com bandas de confiança.
p1 <- p0
F <- model.matrix(formula(m1), data=p1)
U <- chol(vcov(m1))
se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
tval <- qt(p=c(lwr=0.025, fit=0.5, upr=0.975), df=nrow(da)-ncol(F))
me <- outer(se, tval, "*")
p1 <- cbind(p1, sweep(me, 1, p1$y, "+"))
str(p1)
## 'data.frame': 740 obs. of 7 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat: Factor w/ 5 levels "Controle","Ácido cítrico",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dias: num 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 ...
## $ y : atomic 14.1 14.2 14.2 14.2 14.2 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 10.3 10.4 10.4 10.5 10.5 ...
## $ fit : num 14.1 14.2 14.2 14.2 14.2 ...
## $ upr : num 17.9 17.9 18 18 18 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias|est, groups=trat, data=da, type="p", auto.key=TRUE)+
as.layer(xyplot(fit~dias|est, groups=trat, data=p1,
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose))
xyplot(fit~dias|est, groups=trat, data=p1,
auto.key=list(
## corner=c(0.02,0.98),
points=FALSE, lines=TRUE),
xlab="Período de observação (dias)",
ylab="Peso das estacas (g)",
type="l", ## col=1,
cty="bands", alpha=0.3,
ly=p1$lwr, uy=p1$upr,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose)
##-----------------------------------------------------------------------------
## Escala.
da$y <- da$escala
da$y <- with(da, ifelse(y<=2, 0, 1))
useOuterStrips(
xyplot(y~dias|trat+est, groups=rept, data=subset(da, !is.na(y)),
type=c("o")))
useOuterStrips(
xyplot(y~dias|trat+est, data=subset(da, !is.na(y)),
type=c("a","p")))
## Fazer por meio de análise de sobreviência.
## Dado com censura intervalar.
da$muda <- as.numeric(da$escala>2)
## aggregate(muda~dias+est+trat, da, mean)
myfun <- function(x, y){
Y <- cumsum(y)
wm <- sum(Y==0)
x <- c(x[wm+0:1])
status <- 3
if(length(x)==1){
x[2] <- x[1]
status <- 3
}
if(is.na(x[2])){
x[2] <- x[1]
status <- 0
}
matrix(c(x, status), ncol=3,
dimnames=list(NULL, c("esq","dir","status")))
}
## myfun(x=1:5, y=c(0,0,0,0,0))
## myfun(x=1:5, y=c(0,0,0,1,1))
## myfun(x=1:5, y=c(1,1,1,1,1))
str(da)
## 'data.frame': 1000 obs. of 13 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ dias : int 1 1 1 1 1 1 1 1 1 1 ...
## $ trat : Factor w/ 5 levels "Controle","Ácido cítrico",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ L : num 32.1 32.4 45.7 42.9 39.5 ...
## $ a : num -7.55 -11.56 -9.89 -9.79 -8.43 ...
## $ b : num 12 17.9 20.4 18 13.3 ...
## $ AE : num 18.1 25 32.5 29 23.3 ...
## $ peso : num 41.2 15.8 17.1 16.6 14.4 ...
## $ escala: int 1 1 1 1 1 1 1 1 1 1 ...
## $ ue : Factor w/ 200 levels "inverno.Controle.1",..: 2 22 42 62 82 102 122 142 162 182 ...
## $ y : num 0 0 0 0 0 0 0 0 0 0 ...
## $ muda : num 0 0 0 0 0 0 0 0 0 0 ...
dc <- split(da, f=da$ue)
dc <- do.call(rbind,
lapply(dc,
function(d){
z <- myfun(x=d$dias, y=d$muda)
d <- cbind(d[1,c("est","trat","rept")])
d <- cbind(d, data.frame(z))
d
}))
str(dc)
## 'data.frame': 200 obs. of 6 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ trat : Factor w/ 5 levels "Controle","Ácido cítrico",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ rept : int 1 1 1 1 1 1 1 1 1 1 ...
## $ esq : num 11 7 5 5 5 9 7 5 5 7 ...
## $ dir : num 11 9 7 7 7 9 9 7 7 9 ...
## $ status: num 0 3 3 3 3 0 3 3 3 3 ...
table(dc$status)
##
## 0 3
## 75 125
aggregate(cbind(esq, dir)~trat+est, dc, mean)
## trat est esq dir
## 1 Controle inverno 6.4 8.0
## 2 Ácido cítrico inverno 8.0 9.6
## 3 Metabissulfito de sódio inverno 8.2 9.2
## 4 Ácido cítrico e Sacarose inverno 9.6 10.6
## 5 Metabissulfito de sódio e Sacarose inverno 9.0 10.2
## 6 Controle outono 8.4 9.0
## 7 Ácido cítrico outono 7.8 8.6
## 8 Metabissulfito de sódio outono 8.4 9.0
## 9 Ácido cítrico e Sacarose outono 7.8 9.0
## 10 Metabissulfito de sódio e Sacarose outono 8.8 9.0
## 11 Controle primavera 6.0 7.4
## 12 Ácido cítrico primavera 8.2 9.0
## 13 Metabissulfito de sódio primavera 6.4 8.0
## 14 Ácido cítrico e Sacarose primavera 6.8 8.0
## 15 Metabissulfito de sódio e Sacarose primavera 7.2 8.4
## 16 Controle verão 4.2 6.0
## 17 Ácido cítrico verão 4.6 6.6
## 18 Metabissulfito de sódio verão 4.8 6.4
## 19 Ácido cítrico e Sacarose verão 4.4 6.0
## 20 Metabissulfito de sódio e Sacarose verão 4.2 6.2
##-----------------------------------------------------------------------------
## Ajuste de modelos de sobrevivência.
head(dc)
## est trat rept esq dir status
## inverno.Controle.1 inverno Controle 1 11 11 0
## outono.Controle.1 outono Controle 1 7 9 3
## primavera.Controle.1 primavera Controle 1 5 7 3
## verão.Controle.1 verão Controle 1 5 7 3
## inverno.Ácido cítrico.1 inverno Ácido cítrico 1 5 7 3
## outono.Ácido cítrico.1 outono Ácido cítrico 1 9 9 0
## Resposta com censura intervalar e à direita.
S <- with(dc, Surv(time=esq, time2=dir,
event=status, type="interval"))
dist <- "weibull"
head(S)
## [1] 11+ [ 7, 9] [ 5, 7] [ 5, 7] [ 5, 7] 9+
## Sequência de modelos.
s0 <- survreg(S~est*trat, data=dc, dist=dist)
s1 <- survreg(S~est+trat, data=dc, dist=dist)
s2 <- survreg(S~est, data=dc, dist=dist)
s3 <- survreg(S~1, data=dc, dist=dist)
## Interação nula?
anova(s0, s1)
## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 est * trat 179 448.6 NA NA NA
## 2 est + trat 191 465.6 = -12 -17.06 0.1475
## Efeito de tratamento nulo?
anova(s1, s2)
## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 est + trat 191 465.6 NA NA NA
## 2 est 195 468.7 = -4 -3.05 0.5495
anova(s0, s2)
## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 est * trat 179 448.6 NA NA NA
## 2 est 195 468.7 = -16 -20.11 0.2154
## summary(s0)
## coef(s0)
##-----------------------------------------------------------------------------
## Tempo médio de viabilidade (tempo para passar de 2 para 3 na escala)
p0 <- with(da,
expand.grid(est=levels(est),
trat=levels(trat)))
p0 <- cbind(p0,
as.data.frame(predict(s0, newdata=p0,
se.fit=TRUE, type="link")))
## xyplot(fit~trat|est, data=p0)
## str(p0)
subset(p0, trat=="Controle", select=c("est","fit","se.fit"))
## est fit se.fit
## 1 inverno 2.219 0.07955
## 2 outono 2.437 0.12882
## 3 primavera 2.082 0.08357
## 4 verão 1.738 0.07896
## Passando para escala real.
p0 <- transform(p0,
Fit=exp(fit),
lwr=exp(fit-1.96*se.fit),
upr=exp(fit+1.96*se.fit))
segplot(trat~lwr+upr|est, data=p0, centers=Fit,
draw=FALSE)
col <- trellis.par.get()$superpose.line$col[1:nlevels(p0$est)]
key <- list(divide=1, type="o",
title="Estação", cex.title=1.1, columns=2,
lines=list(col=col, pch=19, cex=0.8),
text=list(levels(p0$est)))
col <- col[as.integer(p0$est)]
segplot(trat~lwr+upr, groups=p0$est, data=p0, centers=Fit,
panel=panel.segplotBy, f=0.1, col=col, key=key,
draw=FALSE,
xlab="Tempo para tornar inviável (dias)",
ylab="Solução")
##-----------------------------------------------------------------------------
## Teste de médias.
coef(s2)
## (Intercept) estoutono estprimavera estverão
## 2.37469 0.03406 -0.18313 -0.60750
X <- diag(4); X[,1] <- 1
rownames(X) <- levels(da$est)
Xc <- apc(X)
summary(glht(s2, linfct=X))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: survreg(formula = S ~ est, data = dc, dist = dist)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## inverno == 0 2.3747 0.0405 58.6 <2e-16 ***
## outono == 0 2.4088 0.0584 41.2 <2e-16 ***
## primavera == 0 2.1916 0.0414 53.0 <2e-16 ***
## verão == 0 1.7672 0.0365 48.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summary(glht(s2, linfct=Xc))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: survreg(formula = S ~ est, data = dc, dist = dist)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## inverno-outono == 0 -0.0341 0.0705 -0.48 0.962
## inverno-primavera == 0 0.1831 0.0577 3.17 0.008 **
## inverno-verão == 0 0.6075 0.0547 11.11 <0.001 ***
## outono-primavera == 0 0.2172 0.0705 3.08 0.011 *
## outono-verão == 0 0.6416 0.0700 9.17 <0.001 ***
## primavera-verão == 0 0.4244 0.0554 7.65 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## Não tem diferença entre inverno e outono. No entanto, cuidado, não tem
## repetição de níveis. Não pode ser atribuída essa diferença ao efeito
## de estação pois pode estar confundido com outros fatores.