Grasiela Bruzamarello Tognon
Walmes Marques Zeviani
Francine Lorena Cuquel
O experimento avaliou o efeito dos fator calor de campo {com, sem}. A resposta avaliada foi o total se açúcares. Foram consideradas 10 estacas (unidades experimentais, UE) para cada nível de calor de campo. O registro das respostas foi feito ao longo do tempo (dias, medidas repetidas nas mesmas UE).
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "gridExtra", "nlme", "survival",
"doBy", "plyr", "multcomp", "rootSolve", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra gridExtra nlme survival doBy
## TRUE TRUE TRUE TRUE TRUE TRUE
## plyr multcomp rootSolve wzRfun
## TRUE 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 rootSolve_1.6.5 multcomp_1.3-7 TH.data_1.0-3
## [5] mvtnorm_0.9-9997 plyr_1.8.1 doBy_4.5-11 MASS_7.3-34
## [9] survival_2.37-7 nlme_3.1-117 gridExtra_0.9.1 latticeExtra_0.6-26
## [13] RColorBrewer_1.0-5 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.
## list.files(pattern="*.txt")
da <- read.table("acucares_milleflora.txt", header=TRUE, sep="\t")
levels(da$calor) <- c("Com","Sem")
str(da)
## 'data.frame': 40 obs. of 4 variables:
## $ dias : int 1 1 1 1 1 1 1 1 3 3 ...
## $ calor: Factor w/ 2 levels "Com","Sem": 2 2 2 2 1 1 1 1 2 2 ...
## $ rept : int 1 2 3 4 1 2 3 4 1 2 ...
## $ ast : num 1583 1447 1544 1335 2530 ...
##-----------------------------------------------------------------------------
## Ver.
xlab <- "Dia de avaliação"
xyplot(ast~dias, groups=calor, data=da, jitter.x=TRUE,
auto.key=list(columns=2, title="Calor"), xlab=xlab,
type=c("p","a"),
par.settings=ps)
##-----------------------------------------------------------------------------
## Análise.
str(da)
## 'data.frame': 40 obs. of 4 variables:
## $ dias : int 1 1 1 1 1 1 1 1 3 3 ...
## $ calor: Factor w/ 2 levels "Com","Sem": 2 2 2 2 1 1 1 1 2 2 ...
## $ rept : int 1 2 3 4 1 2 3 4 1 2 ...
## $ ast : num 1583 1447 1544 1335 2530 ...
da$y <- sqrt(da$ast)
da$ue <- with(da, interaction(calor, rept, drop=TRUE))
nlevels(da$ue)
## [1] 8
db <- groupedData(data=da, formula=y~dias|ue, order.groups=FALSE)
m0 <- lme(y~calor*factor(dias), random=~1|ue, data=db, method="ML")
m1 <- lme(y~calor*poly(dias, 1), random=~1|ue, data=db, method="ML")
m2 <- lme(y~calor*poly(dias, 2), random=~1|ue, data=db, method="ML")
## Presença de efeito quadrático.
anova(m1, m2, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 6 229.7412 239.8745 -108.87060
## m2 2 8 205.7236 219.2347 -94.86182 1 vs 2 28.01755 <.0001
## m0 3 12 175.5342 195.8007 -75.76708 2 vs 3 38.18947 <.0001
m0 <- m2
## Interação de calor com o tempo nos termos quadráticos.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 28 5509.310 <.0001
## calor 1 6 5.965 0.0503
## poly(dias, 2) 2 28 60.182 <.0001
## calor:poly(dias, 2) 2 28 13.251 0.0001
## Estimativas dos termos de efeito fixo.
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 31.946007 0.6293803 28 50.757875 4.140646e-29
## calorSem 2.173792 0.8900782 6 2.442248 5.031775e-02
## poly(dias, 2)1 -37.979395 3.9756770 28 -9.552938 2.615813e-10
## poly(dias, 2)2 23.264244 3.9756770 28 5.851643 2.729917e-06
## calorSem:poly(dias, 2)1 19.720209 5.6224564 28 3.507401 1.546412e-03
## calorSem:poly(dias, 2)2 -21.187729 5.6224564 28 -3.768411 7.790929e-04
## Valores ajustados pelo modelo para cada UE.
plot(augPred(m0), par.settings=ps)
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p","smooth")),
xyplot(sqrt(abs(r))~f, type=c("p","smooth")),
qqmath(residuals(m0)),
qqmath(unlist(ranef(m0))))
##-----------------------------------------------------------------------------
## Predição.
pred <- with(db,
expand.grid(calor=levels(calor),
dias=seq(min(dias), max(dias), l=30)))
pred$y <- predict(m0, newdata=pred, level=0)
## str(pred)
xyplot(y~dias, groups=calor, data=da, auto.key=TRUE, par.settings=ps)+
as.layer(xyplot(y~dias, groups=calor, data=pred, type="l"))
## As curvas ajustadas sugerem que o sem calor tem um periodo (de 1 à 5 dias) sem queda
## de peso, depois disso cai a uma taxa igual a do com calor. Para confirmar isso será
## ajustado o modelo platô linear.
##-----------------------------------------------------------------------------
## Ajuste do modelo platô linear para cada nível de calor.
## plot(y~dias, data=subset(da, calor=="Sem"))
## start <- list(b0=38, tx=-1.5, xb=4.5)
## with(start, curve(b0+tx*x*(x<xb)+tx*xb*(x>=xb), add=TRUE))
##
## n0 <- nls(y~b0+tx*dias*(dias<xb)+tx*xb*(dias>=xb),
## data=subset(db, calor=="Sem"),
## start=start)
## coef(n0)
##
## with(as.list(coef(n0)),
## curve(b0+tx*x*(x<xb)+tx*xb*(x>=xb), add=TRUE, col=2))
##
## ##--------------------------------------------
##
## plot(y~dias, data=subset(da, calor=="Com"))
## start <- list(b0=50, tx=-7, xb=3)
## with(start, curve(b0+tx*x*(x<xb)+tx*xb*(x>=xb), add=TRUE))
##
## n1 <- nls(y~b0+tx*dias*(dias<xb)+tx*xb*(dias>=xb),
## data=subset(db, calor=="Com"),
## start=start)
## coef(n1)
##
## with(as.list(coef(n1)),
## curve(b0+tx*x*(x<xb)+tx*xb*(x>=xb), add=TRUE, col=2))
##
## cbind(coef(n0), coef(n1))
## ch <- c(matrix(c(coef(n0), coef(n1)-coef(n0)), ncol=3, by=TRUE))
## dput(round(ch, 5))
##
## Valores iniciais obtidos com o ajuste indivídual.
ch <- c(39.69756, 15.90691, -1.17877, -7.53894, 7.65916, -4.46943)
##-----------------------------------------------------------------------------
## Ajuste dos modelos não lineares mistos.
## Modelo livre, cada valor com seus parâmetros.
n11 <- nlme(y~b0+tx*dias*(dias<xb)+tx*xb*(dias>=xb),
data=db,
random=b0~1, fixed=b0+tx+xb~calor,
start=list(fixed=ch),
verbose=TRUE)
##
## **Iteration 1
## LME step: Loglik: -106.5924 , nlm iterations: 24
## reStruct parameters:
## ue
## 10.13129
##
## PNLS step: RSS = 134.4862
## fixed effects:55.6045 -15.9069 -8.71772 7.53894 3.18973 4.4678
## iterations: 6
##
## Convergence:
## fixed reStruct
## 2.0003647 0.7965753
##
## **Iteration 2
## LME step: Loglik: -81.00918 , nlm iterations: 1
## reStruct parameters:
## ue
## 10.13129
##
## PNLS step: RSS = 134.4862
## fixed effects:55.6045 -15.9069 -8.71772 7.53894 3.18973 4.4678
## iterations: 1
##
## Convergence:
## fixed reStruct
## 0.000000e+00 7.845604e-09
summary(n11)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: y ~ b0 + tx * dias * (dias < xb) + tx * xb * (dias >= xb)
## Data: db
## AIC BIC logLik
## 178.0184 191.5294 -81.00918
##
## Random effects:
## Formula: b0 ~ 1 | ue
## b0.(Intercept) Residual
## StdDev: 7.300402e-05 1.833618
##
## Fixed effects: b0 + tx + xb ~ calor
## Value Std.Error DF t-value p-value
## b0.(Intercept) 55.60447 1.5723145 28 35.36473 0e+00
## b0.calorSem -15.90691 1.8736290 6 -8.48989 1e-04
## tx.(Intercept) -8.71772 0.7031604 28 -12.39791 0e+00
## tx.calorSem 7.53894 0.7374809 28 10.22256 0e+00
## xb.(Intercept) 3.18973 0.1416042 28 22.52565 0e+00
## xb.calorSem 4.46780 1.1771378 28 3.79548 7e-04
## Correlation:
## b0.(I) b0.clS tx.(I) tx.clS xb.(I)
## b0.calorSem -0.839
## tx.(Intercept) -0.894 0.751
## tx.calorSem 0.853 -0.859 -0.953
## xb.(Intercept) -0.351 0.295 0.678 -0.646
## xb.calorSem 0.042 -0.219 -0.082 0.254 -0.120
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.60546219 -0.51331447 0.09503422 0.60219556 2.01647075
##
## Number of Observations: 40
## Number of Groups: 8
logLik(n11)
## 'log Lik.' -81.00918 (df=8)
logLik(m0)
## 'log Lik.' -94.86182 (df=8)
##-----------------------------------------------------------------------------
## Verificar os pressupostos.
## Valores ajustados pelo modelo para cada UE.
plot(augPred(n11), par.settings=ps)
## DIagnóstico.
r <- residuals(n11, type="pearson")
f <- fitted(n11)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p","smooth")),
xyplot(sqrt(abs(r))~f, type=c("p","smooth")),
qqmath(residuals(n11)),
qqmath(unlist(ranef(n11))))
##-----------------------------------------------------------------------------
## Predição.
pred <- with(db,
expand.grid(calor=levels(calor),
dias=seq(min(dias), max(dias), l=30)))
pred$i <- with(pred, ifelse(calor=="Com", 0, 1))
pred$y <- predict(n11, newdata=pred, level=0)
str(pred)
## 'data.frame': 60 obs. of 4 variables:
## $ calor: Factor w/ 2 levels "Com","Sem": 1 2 1 2 1 2 1 2 1 2 ...
## $ dias : num 1 1 1.28 1.28 1.55 ...
## $ i : num 0 1 0 1 0 1 0 1 0 1 ...
## $ y : atomic 46.9 38.5 44.5 38.2 42.1 ...
## ..- attr(*, "label")= chr "Predicted values"
## - attr(*, "out.attrs")=List of 2
## ..$ dim : Named int 2 30
## .. ..- attr(*, "names")= chr "calor" "dias"
## ..$ dimnames:List of 2
## .. ..$ calor: chr "calor=Com" "calor=Sem"
## .. ..$ dias : chr "dias=1.000000" "dias=1.275862" "dias=1.551724" "dias=1.827586" ...
pred$yn <- predict(n11, newdata=pred, level=0)
pred$ym <- predict(m0, newdata=pred, level=0)
## Valores preditos e os observados em cada UE.
xyplot(y~dias, groups=calor, data=da, auto.key=TRUE, par.settings=ps)+
as.layer(xyplot(yn~dias, groups=calor, data=pred, type="l", lwd=2))+
as.layer(xyplot(ym~dias, groups=calor, data=pred, type="l", lwd=2, lty=2))+
as.layer(xyplot(y~dias, groups=ue, data=da, type="l", lty=3, col="gray50"))
## l <- levels(da$calor); nl <- nlevels(da$calor)
l <- c("Com calor","Sem calor"); nl <- length(l)
key <- list(title="Tratamento de calor", cex.title=1.1,
lines=list(col=mycol[1:nl], pch=1, lwd=2),
text=list(l),
columns=nl, divide=1, type="o")
xyplot(y~dias, groups=calor, data=da, key=key, par.settings=ps,
ylab="Peso das estacas (g)", xlab="Dias")+
as.layer(xyplot(yn~dias, groups=calor, data=pred, type="l", lwd=2))
## Então o sem calor a perda de peso tem um padrão com duas fases: até
## aproximadamente 4 dias o peso não cai, após isso cai. No com calor a
## taxa é menor mas existe desde o início.
##-----------------------------------------------------------------------------
## Inclusão das bandas de confiança.
fixef(n11)
## b0.(Intercept) b0.calorSem tx.(Intercept) tx.calorSem xb.(Intercept) xb.calorSem
## 55.604472 -15.906909 -8.717719 7.538944 3.189728 4.467801
formula(n11)
## y ~ b0 + tx * dias * (dias < xb) + tx * xb * (dias >= xb)
## <environment: 0xaee7e50>
## Para obter derivadas numéricas no MLE.
numF <- function(th, xi, Z){
th <- as.matrix(th)
## Z%*%th[c(1,2),]+Z%*%th[c(3,4),]*(xi-Z[,2]*th[5,])*(xi>Z[,2]*th[5,])
Z%*%th[c(1,2),]+Z%*%th[c(3,4),]*xi*(xi<Z%*%th[5:6,])+
Z%*%th[c(3,4),]*Z%*%th[c(5:6),]*(xi>=Z%*%th[5:6,])
}
th <- fixef(n11)
## Verifica se os preditos com numF() correspondem aos do predict().
all(pred$yn==numF(th=th, xi=pred$dias, Z=model.matrix(~calor, pred)))
## [1] TRUE
## Gradiente no MLE.
F <- gradient(numF, x=th, xi=pred$dias,
Z=model.matrix(~calor, pred))
dim(F)
## [1] 60 6
colnames(F)==names(th)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
head(F)
## b0.(Intercept) b0.calorSem tx.(Intercept) tx.calorSem xb.(Intercept) xb.calorSem
## [1,] 1 0 1.000000 0.000000 0 0
## [2,] 1 1 1.000000 1.000000 0 0
## [3,] 1 0 1.275862 0.000000 0 0
## [4,] 1 1 1.275862 1.275862 0 0
## [5,] 1 0 1.551724 0.000000 0 0
## [6,] 1 1 1.551724 1.551724 0 0
## Cálculos do IC.
U <- chol(vcov(n11))
pred$se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
zval <- qnorm(p=c(lwr=0.025, fit=0.5, upr=0.975))
me <- outer(pred$se, zval, "*")
fit <- predict(n11, newdata=pred, level=0)
pred <- cbind(pred, sweep(me, 1, fit, "+"))
str(pred)
## 'data.frame': 60 obs. of 10 variables:
## $ calor: Factor w/ 2 levels "Com","Sem": 1 2 1 2 1 2 1 2 1 2 ...
## $ dias : num 1 1 1.28 1.28 1.55 ...
## $ i : num 0 1 0 1 0 1 0 1 0 1 ...
## $ y : atomic 46.9 38.5 44.5 38.2 42.1 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ yn : atomic 46.9 38.5 44.5 38.2 42.1 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ ym : atomic 44.8 38.6 43.1 38.2 41.4 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ se : num 0.917 0.767 0.8 0.723 0.71 ...
## $ lwr : num 45.1 37 42.9 36.8 40.7 ...
## $ fit : num 46.9 38.5 44.5 38.2 42.1 ...
## $ upr : num 48.7 40 46.1 39.6 43.5 ...
xyplot(y~dias, groups=calor, data=da, key=key, par.settings=ps,
ylab="Raíz da quantidade de açúcares totais",
xlab="Dias")+
as.layer(xyplot(fit~dias, groups=calor, data=pred, type="l",
prepanel=prepanel.cbH, cty="bands", alpha=0.3,
ly=pred$lwr, uy=pred$upr,
panel.groups=panel.cbH,
panel=panel.superpose,
par.settings=ps))
