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. 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: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods splines grid stats graphics grDevices utils datasets
## [9] base
##
## other attached packages:
## [1] wzRfun_0.3 multcomp_1.3-7 TH.data_1.0-3 mvtnorm_0.9-9997
## [5] plyr_1.8.1 doBy_4.5-11 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.7
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_1.0 Matrix_1.1-4 Rcpp_0.11.3 sandwich_2.3-0
## [6] stringr_0.6.2 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("proteina.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] 80
##-----------------------------------------------------------------------------
## Ver.
db <- droplevels(subset(da, especie=="B. milleflora"))
str(db)
## 'data.frame': 420 obs. of 7 variables:
## $ especie: Factor w/ 1 level "B. milleflora": 1 1 1 1 1 1 1 1 1 1 ...
## $ 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 2 2 2 2 3 3 ...
## $ rept : int 1 2 3 4 1 2 3 4 1 2 ...
## $ ua : num 3.05 2.41 4.06 3.39 2.46 ...
## $ ue : Factor w/ 80 levels "inverno.Controle.1",..: 3 23 43 63 7 27 47 67 11 31 ...
db$y <- sqrt(db$ua)
useOuterStrips(
xyplot(y~dias|trat+est, groups=rept, data=db,
type=c("p","a")))
##-----------------------------------------------------------------------------
## m0 <- lme(y~est*trat*(dias+I(dias^2)), random=~1+dias|ue, data=db,
## na.action=na.omit, method="ML")
m0 <- lme(y~est*trat*dias, random=~1+dias|ue, data=db,
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.929162e-12 2.632330e-06 (Intr)
## dias 3.811442e-12 1.952292e-06 0
## Residual 3.930117e-02 1.982452e-01
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 320 18593.603 <.0001
## est 3 60 123.900 <.0001
## trat 4 60 21.588 <.0001
## dias 1 320 211.491 <.0001
## est:trat 12 60 5.658 <.0001
## est:dias 3 320 43.014 <.0001
## trat:dias 4 320 2.166 0.0726
## est:trat:dias 12 320 5.810 <.0001
##-----------------------------------------------------------------------------
## Modelo final.
m1 <- m0
## anova(m1)
##-----------------------------------------------------------------------------
## Representando o resultado.
p0 <- ddply(db, .(est, trat), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
p0$y <- predict(m0, newdata=p0, level=0)
##-----------------------------------------------------------------------------
## Predição com bandas de confiança.
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(db)-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 1.52 1.53 1.54 1.54 1.55 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 1.36 1.37 1.38 1.4 1.41 ...
## $ fit : num 1.52 1.53 1.54 1.54 1.55 ...
## $ upr : num 1.68 1.69 1.69 1.69 1.69 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias|est, groups=trat, data=db, 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="Atividade da Polifenoloxidase (min/mg de proteína)",
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)
##-----------------------------------------------------------------------------
## Ver.
db <- droplevels(subset(da, especie=="B. tridentata"))
str(db)
## 'data.frame': 400 obs. of 7 variables:
## $ especie: Factor w/ 1 level "B. tridentata": 1 1 1 1 1 1 1 1 1 1 ...
## $ 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 2 2 2 2 3 3 ...
## $ rept : int 1 2 3 4 1 2 3 4 1 2 ...
## $ ua : num 1.147 1.212 1.089 1.136 0.507 ...
## $ ue : Factor w/ 80 levels "inverno.Controle.1",..: 3 23 43 63 7 27 47 67 11 31 ...
db$y <- sqrt(db$ua)
useOuterStrips(
xyplot(y~dias|trat+est, groups=rept, data=db,
type=c("p","a")))
##-----------------------------------------------------------------------------
## m0 <- lme(y~est*trat*(dias+I(dias^2)), random=~1+dias|ue, data=db,
## na.action=na.omit, method="ML")
m0 <- lme(y~est*trat*dias, random=~1+dias|ue, data=db,
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.624224e-12 2.573757e-06 (Intr)
## dias 8.610672e-14 2.934395e-07 0
## Residual 3.059528e-02 1.749151e-01
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 300 13224.492 <.0001
## est 3 60 10.283 <.0001
## trat 4 60 17.298 <.0001
## dias 1 300 61.291 <.0001
## est:trat 12 60 12.639 <.0001
## est:dias 3 300 8.206 <.0001
## trat:dias 4 300 5.666 2e-04
## est:trat:dias 12 300 7.667 <.0001
##-----------------------------------------------------------------------------
## Modelo final.
m1 <- m0
## anova(m1)
##-----------------------------------------------------------------------------
## Representando o resultado.
p0 <- ddply(db, .(est, trat), summarise,
dias=seq(min(dias)-0.5, max(dias)+0.5, by=0.25))
p0$y <- predict(m0, newdata=p0, level=0)
##-----------------------------------------------------------------------------
## Predição com bandas de confiança.
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(db)-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 1.37 1.36 1.34 1.33 1.32 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ lwr : num 1.23 1.23 1.22 1.21 1.2 ...
## $ fit : num 1.37 1.36 1.34 1.33 1.32 ...
## $ upr : num 1.5 1.48 1.47 1.45 1.43 ...
col <- trellis.par.get()$superpose.line$col[1:4]
xyplot(y~dias|est, groups=trat, data=db, 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="Atividade da Polifenoloxidase (min/mg de proteína)",
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)
