Grasiela Bruzamarello Tognon
Walmes Marques Zeviani
Francine Lorena Cuquel
Experimento feito com estacas de B. millefolra 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_milleflora.txt", header=TRUE, sep="\t")
## str(da)
da$trat <- factor(da$trat)
levels(da$trat) <- 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)
## [1] 200
##-----------------------------------------------------------------------------
## 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.51840 1.2322 (Intr)
## dias 0.02606 0.1614 -0.235
## Residual 1.23206 1.1100
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 802 12580 <.0001
## est 3 180 39 <.0001
## trat 4 180 3 0.0343
## dias 1 802 2 0.1254
## I(dias^2) 1 802 1 0.3677
## est:trat 12 180 1 0.1631
## est:dias 3 802 7 0.0001
## est:I(dias^2) 3 802 5 0.0015
## trat:dias 4 802 5 0.0003
## trat:I(dias^2) 4 802 1 0.4887
## est:trat:dias 12 802 4 <.0001
## est:trat:I(dias^2) 12 802 4 <.0001
##-----------------------------------------------------------------------------
## 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': 780 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 -8.83 -8.92 -9 -9.08 -9.15 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num -9.94 -9.96 -9.99 -10.02 -10.06 ...
## $ fit : num -8.83 -8.92 -9 -9.08 -9.15 ...
## $ upr : num -7.72 -7.88 -8.01 -8.13 -8.24 ...
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) 3.11398 1.7646 (Intr)
## dias 0.01211 0.1101 -0.413
## Residual 6.16235 2.4824
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 823 81872 <.0001
## est 3 180 43 <.0001
## trat 4 180 0 0.9064
## dias 1 823 143 <.0001
## est:trat 12 180 1 0.8192
## est:dias 3 823 38 <.0001
## trat:dias 4 823 1 0.6487
## est:trat:dias 12 823 1 0.6012
##-----------------------------------------------------------------------------
## Modelo final.
m1 <- update(m0, y~est*dias)
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 44 5195 5413 -2554
## m1 2 12 5153 5212 -2564 1 vs 2 21.49 0.9204
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 839 80827 <.0001
## est 3 196 42 <.0001
## dias 1 839 144 <.0001
## est:dias 3 839 38 <.0001
##-----------------------------------------------------------------------------
## 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 <- ddply(da, .(est), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
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': 156 obs. of 6 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 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 38.9 38.9 38.9 38.9 38.9 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 38.1 38.1 38.1 38.2 38.2 ...
## $ fit : num 38.9 38.9 38.9 38.9 38.9 ...
## $ upr : num 39.6 39.6 39.6 39.6 39.6 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias, groups=est, data=da, type="p", auto.key=TRUE)+
as.layer(xyplot(fit~dias, groups=est, 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, groups=est, 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) 6.82193 2.6119 (Intr)
## dias 0.03206 0.1791 -0.744
## Residual 4.71118 2.1705
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 827 11868 <.0001
## est 3 180 37 <.0001
## trat 4 180 1 0.3601
## dias 1 827 26 <.0001
## est:trat 12 180 1 0.2421
## est:dias 3 827 36 <.0001
## trat:dias 4 827 1 0.4446
## est:trat:dias 12 827 1 0.3319
##-----------------------------------------------------------------------------
## Modelo final.
m1 <- update(m0, y~est*dias)
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 44 5071 5289 -2491
## m1 2 12 5043 5103 -2510 1 vs 2 36.77 0.2576
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 843 11052 <.0001
## est 3 196 34 <.0001
## dias 1 843 26 <.0001
## est:dias 3 843 36 <.0001
##-----------------------------------------------------------------------------
## 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 <- ddply(da, .(est), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
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': 156 obs. of 6 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 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 16.9 16.9 16.8 16.7 16.7 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 16 16 15.9 15.9 15.9 ...
## $ fit : num 16.9 16.9 16.8 16.7 16.7 ...
## $ upr : num 17.8 17.7 17.7 17.6 17.5 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias, groups=est, data=da, type="p", auto.key=TRUE)+
as.layer(xyplot(fit~dias, groups=est, 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, groups=est, 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) 10.93664 3.3071 (Intr)
## dias 0.05167 0.2273 -0.745
## Residual 8.38498 2.8957
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 816 18827 <.0001
## est 3 180 50 <.0001
## trat 4 180 1 0.2504
## dias 1 816 57 <.0001
## est:trat 12 180 1 0.4848
## est:dias 3 816 37 <.0001
## trat:dias 4 816 1 0.2048
## est:trat:dias 12 816 1 0.2836
##-----------------------------------------------------------------------------
## Modelo final.
m1 <- update(m0, y~est*dias)
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 44 5597 5814 -2754
## m1 2 12 5570 5629 -2773 1 vs 2 37.16 0.2432
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 832 17701 <.0001
## est 3 196 47 <.0001
## dias 1 832 56 <.0001
## est:dias 3 832 36 <.0001
##-----------------------------------------------------------------------------
## 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 <- ddply(da, .(est), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
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': 156 obs. of 6 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 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 26.2 26.2 26.1 26.1 26.1 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 25.1 25.1 25 25 25 ...
## $ fit : num 26.2 26.2 26.1 26.1 26.1 ...
## $ upr : num 27.4 27.3 27.2 27.2 27.1 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias, groups=est, data=da, type="p", auto.key=TRUE)+
as.layer(xyplot(fit~dias, groups=est, 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, groups=est, 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 3449 3753 -1660
## m1 2 44 3691 3900 -1801 1 vs 2 281.9 <.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) 53.2848 7.2996 (Intr)
## dias 0.1282 0.3581 -0.55
## Residual 0.4293 0.6552
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 619 1081.8 <.0001
## est 3 180 3.4 0.0186
## trat 4 180 0.6 0.6833
## dias 1 619 96.0 <.0001
## I(dias^2) 1 619 196.9 <.0001
## est:trat 12 180 0.6 0.8638
## est:dias 3 619 49.9 <.0001
## est:I(dias^2) 3 619 19.9 <.0001
## trat:dias 4 619 2.6 0.0359
## trat:I(dias^2) 4 619 2.9 0.0225
## est:trat:dias 12 619 1.7 0.0563
## est:trat:I(dias^2) 12 619 3.8 <.0001
##-----------------------------------------------------------------------------
## 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': 780 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 12.9 13.1 13.1 13.2 13.3 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 8.43 8.57 8.71 8.83 8.95 ...
## $ fit : num 12.9 13.1 13.1 13.2 13.3 ...
## $ upr : num 17.5 17.5 17.6 17.6 17.7 ...
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.98,0.02),
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': 1050 obs. of 14 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ 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 38.8 40 40.1 40.3 39.9 ...
## $ a : num -11.2 -13.2 -12.9 -13.2 -12.5 ...
## $ b : num 21.2 18.4 17.1 17.3 17.7 ...
## $ AE : num 29.8 29.3 28.3 28.6 28.4 ...
## $ peso : num 9.32 13 21.99 5.63 31.9 ...
## $ escala: int 1 1 1 1 1 1 1 1 1 1 ...
## $ ue : Factor w/ 200 levels "inverno.Controle.1",..: 3 23 43 63 83 103 123 143 163 183 ...
## $ i : int 1 2 3 4 5 6 7 8 9 10 ...
## $ 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 7 11 7 5 9 11 7 7 9 11 ...
## $ dir : num 9 13 7 7 9 13 7 9 9 13 ...
## $ status: num 3 3 0 3 0 3 0 3 0 3 ...
table(dc$status)
##
## 0 3
## 78 122
aggregate(cbind(esq, dir)~trat+est, dc, mean)
## trat est esq dir
## 1 Controle inverno 8.0 8.6
## 2 Ácido cítrico inverno 7.2 8.4
## 3 Metabissulfito de sódio inverno 8.4 8.8
## 4 Ácido cítrico e Sacarose inverno 6.6 7.6
## 5 Metabissulfito de sódio e Sacarose inverno 7.4 8.8
## 6 Controle outono 11.4 12.8
## 7 Ácido cítrico outono 11.8 13.0
## 8 Metabissulfito de sódio outono 11.4 13.0
## 9 Ácido cítrico e Sacarose outono 11.2 12.6
## 10 Metabissulfito de sódio e Sacarose outono 11.0 12.6
## 11 Controle primavera 4.8 5.8
## 12 Ácido cítrico primavera 6.6 7.0
## 13 Metabissulfito de sódio primavera 4.0 5.8
## 14 Ácido cítrico e Sacarose primavera 4.6 5.8
## 15 Metabissulfito de sódio e Sacarose primavera 5.4 6.6
## 16 Controle verão 5.2 6.6
## 17 Ácido cítrico verão 6.8 8.4
## 18 Metabissulfito de sódio verão 6.6 8.2
## 19 Ácido cítrico e Sacarose verão 7.8 8.8
## 20 Metabissulfito de sódio e Sacarose verão 7.4 8.8
##-----------------------------------------------------------------------------
## Ajuste de modelos de sobrevivência.
head(dc)
## est trat rept esq dir status
## inverno.Controle.1 inverno Controle 1 7 9 3
## outono.Controle.1 outono Controle 1 11 13 3
## primavera.Controle.1 primavera Controle 1 7 7 0
## verão.Controle.1 verão Controle 1 5 7 3
## inverno.Ácido cítrico.1 inverno Ácido cítrico 1 9 9 0
## outono.Ácido cítrico.1 outono Ácido cítrico 1 11 13 3
## Resposta com censura intervalar e à direita.
S <- with(dc, Surv(time=esq, time2=dir,
event=status, type="interval"))
dist <- "weibull"
head(S)
## [1] [ 7, 9] [11, 13] 7+ [ 5, 7] 9+ [11, 13]
## 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 496.0 NA NA NA
## 2 est + trat 191 518.8 = -12 -22.8 0.02949
## Efeito de tratamento nulo?
anova(s1, s2)
## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 est + trat 191 518.8 NA NA NA
## 2 est 195 521.3 = -4 -2.449 0.6538
anova(s0, s2)
## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 est * trat 179 496.0 NA NA NA
## 2 est 195 521.3 = -16 -25.25 0.06559
## 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.419 0.12746
## 2 outono 2.582 0.08313
## 3 primavera 1.979 0.09779
## 4 verão 2.097 0.08352
## 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.3044 0.2720 -0.3406 -0.1417
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.3044 0.0488 47.2 <2e-16 ***
## outono == 0 2.5765 0.0384 67.1 <2e-16 ***
## primavera == 0 1.9638 0.0434 45.2 <2e-16 ***
## verão == 0 2.1627 0.0389 55.6 <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.2720 0.0611 -4.46 <0.001 ***
## inverno-primavera == 0 0.3406 0.0644 5.29 <0.001 ***
## inverno-verão == 0 0.1417 0.0620 2.29 0.1007
## outono-primavera == 0 0.6127 0.0575 10.66 <0.001 ***
## outono-verão == 0 0.4137 0.0544 7.61 <0.001 ***
## primavera-verão == 0 -0.1990 0.0581 -3.43 0.0034 **
## ---
## 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 verão. 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.