Rafael Henrique de Tonissi e Buschinelli de Goes
Walmes Marques Zeviani
Os dados são de um experimento que avaliou o comportamento ingestivo em função da adição de copaíba ao concentrado oferecido para os animais. O experimentos foi instalado segundo um delineamento de quadrado latino com 4 animais (linhas), 4 períodos (colunas) e 4 níveis de óleo de copaíba: 0, 2.9, 5.8 e 8.7 g/dia. Cada uma das 16 celas experimentais recebeu uma quantidade conhecida de concetrado à 9h30. A cada 15 minutos o conteúdo de substrato restante no comedouro foi medido para se ter conhecimento da curva de ingestão (peso de concentrado ingerido em função do tempo).
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "gridExtra", "nlme",
"doBy", "multcomp", "plyr", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra gridExtra nlme doBy multcomp
## TRUE TRUE TRUE TRUE TRUE TRUE
## plyr wzRfun
## TRUE TRUE
## Instruções para instalação do pacote wzRfun em:
## https://github.com/walmes/wzRfun
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 plyr_1.8.1 multcomp_1.3-6 TH.data_1.0-3
## [5] mvtnorm_0.9-9997 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
##-----------------------------------------------------------------------------
## Ler.
## list.files(pattern="*.txt")
da <- read.table("http://www.leg.ufpr.br/~walmes/data/compingest.txt",
header=TRUE, sep="\t",
stringsAsFactors=FALSE)
## Conversão para fator.
da <- transform(da, animal=factor(animal),
periodo=factor(periodo))
## Passando a representação de horas para minutos após cocho cheio.
da$h <- with(da, as.POSIXct(hora, format="%H:%M"))
da$h <- with(da, h-as.POSIXct("09:30", format="%H:%M"))
units(da$h) <- "mins"
da$h <- as.numeric(da$h)
## Criando nível da unidade experimental.
da$ue <- with(da, interaction(periodo, animal, drop=TRUE))
str(da)
## 'data.frame': 144 obs. of 8 variables:
## $ animal : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 2 2 ...
## $ piquete : int 1 1 1 1 1 1 1 1 2 2 ...
## $ periodo : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ copaiba : num 0 0 0 0 0 0 0 0 2.9 2.9 ...
## $ hora : chr "9:30" "9:45" "10:00" "10:15" ...
## $ pesoconc: int 230 130 62 62 20 20 0 0 470 329 ...
## $ h : num 0 15 30 45 60 75 90 105 0 15 ...
## $ ue : Factor w/ 16 levels "1.1","2.1","3.1",..: 1 1 1 1 1 1 1 1 5 5 ...
## Criando uma cópia que será fator.
da$copa <- as.factor(da$copaiba)
levels(da$copa)
## [1] "0" "2.9" "5.8" "8.7"
##-----------------------------------------------------------------------------
## Tabela de frequência dos registros.
xtabs(~animal+periodo, da)
## periodo
## animal 1 2 3 4
## 1 8 7 9 12
## 2 8 8 9 11
## 3 8 8 9 11
## 4 8 9 9 10
xtabs(~animal+copa, da)
## copa
## animal 0 2.9 5.8 8.7
## 1 8 7 9 12
## 2 11 8 8 9
## 3 9 11 8 8
## 4 9 9 10 8
##-----------------------------------------------------------------------------
## Ver.
xyplot(pesoconc~h|copa, groups=animal, data=da,
type="b", as.table=TRUE,
xlab="Minutos", ylab="Conteúdo de concentrado no comedouro")
A análise exploratória indica que nas unidades experimentais o número de registros no tempo não foi o mesmo. Ainda, pelo gráficos, pode-se perceber que existe considerável variábilidade entre unidades e bem menor variabilidade entre observações dentro das unidades. A tendência apresentada, conforme se esperava, é de decaimento do conteúdo de concentrado em função do tempo.
Sem manisfestar preocupação sobre usar um modelo de decaímento, por enquanto, será ajustado um polinômio de grau 2. Pela análise exploratória, considera-se este o suficiente para representar o efeito do tempo sobre o conteúdo do cocho. Além disso, será considerado o efeito de unidade experimental (16 níveis, 16-1 graus de liberdade) ao invés de efeito de linhas e colunas (2*(4-1) graus de liberdade), o que na realidade equivale à considerar o efeito da interação entre linha e coluna, que são as unidades experimentais. Fica então um experimento longitudinal com efeito fixo do nível de copaíba e tempo e o efeito de unidade experimental sobre os parâmetros do modelo polinomial de grau 2.
##-----------------------------------------------------------------------------
## Especificando o modelo.
da <- groupedData(data=da, pesoconc~h|ue)
str(da)
## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 144 obs. of 9 variables:
## $ animal : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 2 2 ...
## $ piquete : int 1 1 1 1 1 1 1 1 2 2 ...
## $ periodo : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ copaiba : num 0 0 0 0 0 0 0 0 2.9 2.9 ...
## $ hora : chr "9:30" "9:45" "10:00" "10:15" ...
## $ pesoconc: int 230 130 62 62 20 20 0 0 470 329 ...
## $ h : num 0 15 30 45 60 75 90 105 0 15 ...
## $ ue : Ord.factor w/ 16 levels "1.1"<"2.1"<"1.4"<..: 1 1 1 1 1 1 1 1 11 11 ...
## $ copa : Factor w/ 4 levels "0","2.9","5.8",..: 1 1 1 1 1 1 1 1 2 2 ...
## - attr(*, "formula")=Class 'formula' length 3 pesoconc ~ h | ue
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "FUN")=function (x)
## - attr(*, "order.groups")= logi TRUE
## Considerou-se o efeito de copaíba como categórico porque são apenas 4
## níveis. O ajuste de modelos contínuos, como polinomiais, para
## representar o efeito baseado em poucos níveis pode não ser
## justificável.
m0 <- lme(pesoconc~copa*(h+I(h^2)),
random=~1+h+I(h^2)|ue,
data=da, method="ML")
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(nrow=1,
xyplot(r~f, type=c("p","smooth")),
xyplot(sqrt(abs(r))~f, type=c("p","smooth")))
grid.arrange(nrow=1,
xyplot(r~f|da$ue)+layer(panel.abline(h=0, lty=2)),
qqmath(~r|da$ue))
## Estimativas dos efeitos fixos (na parametrização considera).
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 403.422980 52.814525 120 7.63849 5.916e-12
## copa2.9 -5.371338 74.693014 12 -0.07191 9.439e-01
## copa5.8 11.853222 74.738952 12 0.15859 8.766e-01
## copa8.7 -1.684750 74.648051 12 -0.02257 9.824e-01
## h -4.865141 1.059849 120 -4.59041 1.098e-05
## I(h^2) 0.018156 0.005202 120 3.49044 6.754e-04
## copa2.9:h 0.727599 1.500113 120 0.48503 6.285e-01
## copa5.8:h 1.385002 1.512774 120 0.91554 3.617e-01
## copa8.7:h 0.820790 1.489405 120 0.55109 5.826e-01
## copa2.9:I(h^2) -0.003395 0.007414 120 -0.45787 6.479e-01
## copa5.8:I(h^2) -0.009180 0.007641 120 -1.20135 2.320e-01
## copa8.7:I(h^2) -0.006994 0.007213 120 -0.96966 3.342e-01
## Estimativas dos parâmetros de covariância, parte aleatória.
VarCorr(m0)
## ue = pdLogChol(1 + h + I(h^2))
## Variance StdDev Corr
## (Intercept) 9.713e+03 98.556935 (Intr) h
## h 3.472e+00 1.863288 -0.325
## I(h^2) 6.212e-05 0.007881 0.231 -0.984
## Residual 8.298e+02 28.806477
## Quadro de testes de Wald para termos de efeito fixo.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 120 145.74 <.0001
## copa 3 12 0.15 0.9287
## h 1 120 82.21 <.0001
## I(h^2) 1 120 25.62 <.0001
## copa:h 3 120 0.70 0.5565
## copa:I(h^2) 3 120 0.58 0.6320
## Modelo reduzido com o abandono do efeito de copaíba. A diferença em
## grau de ajuste entre esses modelos representa a magnitude do efeito
## de copaíba no comportamento ingestivo.
m1 <- update(m0, .~h+I(h^2))
## Razão de verossimilhanças para testar a hipótese nula de efeito nulo
## de copaíba.
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 19 1534 1590 -747.8
## m1 2 10 1519 1548 -749.3 1 vs 2 3.035 0.9629
## Hipótese não rejeitada.
##-----------------------------------------------------------------------------
## Predição para ver a curva média de comportamento ingestivo ao nível
## das unidades experimentais e dos efeito fixos.
## Combinações presentes.
## with(da, unique(cbind(ue, copa)))
## range(da$h)
## Para a parte de efeito aleatório.
pred1 <- ddply(da, .(ue), summarise,
h=seq(min(h), max(h), by=5))
pred1 <- merge(pred1, unique(da[,c("ue","copa")]), by="ue")
pred1$y1 <- predict(m0, newdata=pred1, level=1)
str(pred1)
## 'data.frame': 400 obs. of 4 variables:
## $ ue : Ord.factor w/ 16 levels "1.1"<"2.1"<"1.4"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ h : num 0 5 10 15 20 25 30 35 40 45 ...
## $ copa: Factor w/ 4 levels "0","2.9","5.8",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ y1 : atomic 203 185 167 150 134 ...
## ..- attr(*, "label")= chr "Predicted values"
## Para a parte de efeito fixo.
m <- median(with(da, tapply(h, ue, max)))
pred0 <- with(da, expand.grid(h=seq(0, m, 5), copa=levels(copa)))
pred0$y0 <- predict(m0, newdata=pred0, level=0)
## Legenda.
key <- list(title="Nível de predição", cex.title=1.1,
lines=list(lwd=1:2, col=c("gray50","black"), lty=2:1),
text=list(c("Unidade experimental", "Copaíba")))
## Texto para as tarjas.
fl <- sapply(levels(da$copa),
function(x){
a <- substitute(Copaíba~x~g~dia^{-1}, list(x=x))
a
})
fl <- do.call(expression, fl)
## Gráfico.
xyplot(pesoconc~h|copa, groups=animal, data=da,
key=key, as.table=TRUE,
strip=strip.custom(factor.levels=fl),
ylab="Conteúdo de concentrado no comedouro (g)",
xlab="Período após abastecimento do comedouro (min)")+
as.layer(xyplot(y1~h|copa, groups=ue, data=pred1,
type="l", lty=2, col="gray50"))+
as.layer(xyplot(y0~h|copa, data=pred0, type="l", col=1, lwd=2))
##-----------------------------------------------------------------------------
## Bandas de confiança para as curvas populacionais (nível de efeito
## fixo).
U <- chol(vcov(m0))
X <- model.matrix(formula(m0)[-2], data=pred0)
pred0$se <- sqrt(apply(X%*%t(U), 1, function(x) sum(x^2)))
zval <- qnorm(p=c(lwr=0.025, fit=0.5, upr=0.975))
me <- outer(pred0$se, zval, "*")
fit <- crossprod(t(X), fixef(m0))
pred0 <- cbind(pred0, sweep(me, 1, fit, "+"))
str(pred0)
## 'data.frame': 100 obs. of 7 variables:
## $ h : num 0 5 10 15 20 25 30 35 40 45 ...
## $ copa: Factor w/ 4 levels "0","2.9","5.8",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ y0 : atomic 403 380 357 335 313 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ se : num 50.6 49 47.9 47.2 47 ...
## $ lwr : num 304 284 263 242 221 ...
## $ fit : num 403 380 357 335 313 ...
## $ upr : num 503 476 450 427 405 ...
xyplot(y0~h, groups=copa, data=pred0, type="l", #col=1, lwd=2,
ylab="Conteúdo de concentrado no comedouro (g)",
xlab="Período após abastecimento do comedouro (min)",
auto.key=list(columns=4, points=FALSE, lines=TRUE,
title=expression(Copaíba~(g~dia^{-1})),
cex.title=1.1),
prepanel=prepanel.cbH, cty="bands",
ly=pred0$lwr, uy=pred0$upr, alpha=0.1, fill=1,
panel.groups=panel.cbH,
panel=panel.superpose)## +
## layer_(panel.abline(v=seq(0, 120, 20),
## h=seq(0, 500, 100),
## col="gray50", lty=2),
## under=TRUE)
Ao conderar tanto a análise (teste de hipótese) quanto os gráficos, verifica-se que não existe efeito da adição de copaíba, considerando os níveis estudados, no comportamento ingestivo de concetrado.