Eliane Aparecida Rogovski Czaja
Walmes Marques Zeviani
Louise Larissa May De Mio
O experimento consistiu na avaliação do período de latência e número de pústulas de ferrugem da figueira avaliada in vivo nas folhas de figueira em função da temperatura, com 9 níveis. Foram consideradas de 3 à 5 folhas por planta (unidade experimental) e 3 plantas por temperatura. Para avaliar os 9 níveis, foram feitos 3 experimentos sob o mesmo delineamento (DIC), mudando-se apenas as temperaturas avaliadas. Assim, cada repetição do experimento representa um bloco.
O tempo para aparecimento de pustulas e o valor de número de pústulas (severidade) foram os objetivos de investigação. Como as avaliações foram diárias, o tempo para aparecimento de pústula é uma variável registrada com censura intervalar quando a ocorre entre dois dias e censura à direita quando ocorre após o último dia de avaliação. A partir das avaliações foram obtidos os tempos e tipos de censura. O seguinte modelo foi considerado:
\[ \begin{aligned} Y_{ijk} &\sim Normal(\mu_{ij}, \sigma^2) \newline \mu_{ij} &= \mu+ \alpha_i +\tau_j \newline \end{aligned} \]
em que \(Y_{ijk}\) é o tempo para aparecimento de pústula na unidade experimental \(k\) da temperatura \(j\) no bloco \(i\). O efeito fixo de bloco e manejo são representados por \(\alpha_i\) e \(\tau_j\). A variância residual representada por \(\sigma^2\). O modelo teve seus parâmetros estimados por máxima verossimilhança considerando o tipos de censura intervalar e à direita.
O teste de Wald foi aplicado para avaliar a significância dos termos de efeito fixo. Foram aplicados contrastes entre médias para os níveis de temperatura. O p-valor para os contrastes foram corrigidos pelo método fdr (false discovery rate). Todos os testes de hipótese foram conduzidos sob um nível nominal de significância de 5%.
Para análise da severidade foi considerado à área abaixo da curva do número de pústulas nas folhas ao longo das avaliações dividido pelo número de dias em avaliação o que representa assim a média do número de pústulas. O efeito da temperatura foi representado pelo modelo de regressão não linear:
\[ \begin{aligned} f(x) &= \theta_{y} \exp\{\theta_{q} (x-\theta_{x})^2+\theta_{c}(x-\theta_{x})^3\} \end{aligned} \]
em que \(\theta_y\) representa o valor médio da resposta para a temperatura (\(x\)) ótima \(\theta_x\), \(\theta_q\) está relacionado ao nível de especificidade à temperatura de ótimo (quanto maior \(\theta_q\) menor é a severidade para temperaturas longe do ótimo) e \(\theta_c\) está relacionado à assimetria na resposta para variações de estímulo da temperatura ao redor do ótimo.
O modelo teve seus parâmetros estimados por máxima verossimilhança. Os pressupostos do modelo considerado foram inspecionados por análise gráfica dos resíduos. Transformações da resposta, eliminação de observações ou modelagem da variância foram abordagens considerados para adequar os pressupostos se apresentassem fuga. Todos os testes de hipótese foram conduzidos sob um nível nominal de significância de 5%.
##-----------------------------------------------------------------------------
## Definições da sessão.
## rm(list=ls())
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "plyr",
"reshape", "survival", "nlme")
sapply(pkg, require, character.only=TRUE)
## Funções carregadas da Public do dropbox.
## dburl <- "http://dl.dropboxusercontent.com/u/48140237/"
dburl <- "~/Dropbox/Public/"
extrafun <- c("aac.R","apc.R","desdobrglht.R","equallevels.R","bandas.R")
sapply(paste0(dburl, extrafun), source)
trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Para formatar tabelas.
myformat <- function(table, digits){
tb <- table
for(i in 1:ncol(tb)){
x <- tb[,i]
d <- digits[i]
if(is.numeric(x)){
f <- paste0(" %0.", d, "f")
x <- sprintf(f, x)
} else {
x <- as.character(x)
}
tb[,i] <- x
}
return(tb)
}
## myformat(cars, c(2,4))
##-----------------------------------------------------------------------------
## Leitura.
da <- read.table("folha_pustula.txt", header=TRUE, sep="\t")
str(da)
## 'data.frame': 175 obs. of 17 variables:
## $ bloco: Factor w/ 3 levels "I","II","III": 1 1 1 1 1 1 1 1 1 1 ...
## $ temp : int 20 20 20 20 20 20 20 20 20 20 ...
## $ rept : int 1 1 1 1 1 2 2 2 2 2 ...
## $ folha: Factor w/ 5 levels "A","B","C","D",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ dia5 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ dia6 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ dia7 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ dia8 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ dia9 : int 0 0 1 0 0 0 0 4 2 1 ...
## $ dia10: int 0 0 3 0 3 1 0 23 4 5 ...
## $ dia11: int 2 3 19 3 11 15 5 98 10 19 ...
## $ dia12: int 48 9 25 17 43 35 18 107 41 38 ...
## $ dia13: int 60 11 27 20 57 68 24 110 55 58 ...
## $ dia14: int 60 11 29 20 59 75 25 113 55 60 ...
## $ dia15: int NA NA NA NA NA NA NA NA NA NA ...
## $ dia16: int NA NA NA NA NA NA NA NA NA NA ...
## $ dia17: int NA NA NA NA NA NA NA NA NA NA ...
##-----------------------------------------------------------------------------
## Conversões.
## Criar os níveis de unidade experimental (plantas).
## Criar os níveis de unidade amostral (folhas).
da <- transform(da,
ue=interaction(bloco,temp,rept),
ua=interaction(bloco,temp,rept,folha),
uf=interaction(rept,folha),
Temp=factor(temp))
## Calcular o tempo para a primeira observação não nula de
## pustulas. Esse será o limite superior do intervalo censurado de
## observação.
diacol <- grep("^dia\\d", names(da))
surv <- apply(da[,diacol], 1,
function(x){
x <- na.omit(x)
cs <- cumsum(x)
time <- sum(cs==0)+4
tot <- length(x)+4
isless <- time<tot
c(time2=time+isless, status=ifelse(isless, 3, 0))
})
da <- cbind(da, t(surv))
da$time <- with(da, ifelse(status==0, time2, time2-1))
## str(da)
##-----------------------------------------------------------------------------
## Empilhar nas folhas.
db <- melt(da, id.vars=c("bloco","temp","Temp","ue","ua","uf"),
measure.vars=diacol, variable_name="dia")
db$dia <- as.integer(gsub("\\D", "", db$dia))
db <- na.omit(db)
str(db)
## 'data.frame': 1975 obs. of 8 variables:
## $ bloco: Factor w/ 3 levels "I","II","III": 1 1 1 1 1 1 1 1 1 1 ...
## $ temp : int 20 20 20 20 20 20 20 20 20 20 ...
## $ Temp : Factor w/ 9 levels "15","18","20",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ ue : Factor w/ 135 levels "I.15.1","II.15.1",..: 7 7 7 7 7 34 34 34 34 34 ...
## $ ua : Factor w/ 675 levels "I.15.1.A","II.15.1.A",..: 7 142 277 412 547 34 169 304 439 574 ...
## $ uf : Factor w/ 25 levels "1.A","2.A","3.A",..: 1 6 11 16 21 2 7 12 17 22 ...
## $ dia : int 5 5 5 5 5 5 5 5 5 5 ...
## $ value: int 0 0 0 0 0 0 0 0 0 0 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:300] 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 ...
## .. ..- attr(*, "names")= chr [1:300] "1751" "1752" "1753" "1754" ...
da <- da[,-diacol]
## str(da)
##-----------------------------------------------------------------------------
cap <-
"Tempo para o aparecimento de pústulas nas folhas em função da
temperatura. Pontos ligados por segmentos indicam que o evento foi
observado no intervalo (censura intervalar). Pontos sem segmentos
indicam o evento não ocorreu antes da avaliação (censura à direita)."
##-----------------------------------------------------------------------------
## Enfolhamento.
xlab <- "Folhas"
ylab <- "Tempo para o aparecimento de pústulas nas folhas"
dotplot(time2+time~uf|Temp, data=da, layout=c(NA,1),
scales=list(x=list(relation="free", draw=FALSE)),
xlab=xlab, ylab=ylab)+
as.layer(
segplot(uf~time2+time|Temp, data=da, draw=FALSE,
layout=c(NA,1), horizontal=FALSE,
scales=list(x=list(relation="free", draw=FALSE))))
##-----------------------------------------------------------------------------
cap <-
"Número de pústulas em função dos dias de avaliação. Pontos unidos por linhas
indicam a mesma folha em cada uma das temperaturas."
##-----------------------------------------------------------------------------
## Número acumulado de pústulas.
ylab <- "Pústulas"
xlab <- "Dias de avaliação"
xyplot(value~dia|Temp, groups=uf, data=db, type="o",
xlab=xlab, ylab=ylab)
##-----------------------------------------------------------------------------
## Preparando resposta para análise de sobreviência.
S <- with(da, Surv(time=time, time2=time2, event=status, type="interval"))
##-----------------------------------------------------------------------------
## Ajuste considerando distribuição Weibull.
## dist <- "exponential"
## dist <- "weibull"
dist <- "gaussian"
## Modelo com efeito de temp dado por polinômio.
m0 <- survreg(S~bloco+temp+I(temp^2), data=da, dist=dist)
## Modelo com efeito categórico de temp.
m1 <- survreg(S~bloco+Temp, data=da, dist=dist)
## Existe falta de ajuste?
anova(m0, m1)
## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 bloco + temp + I(temp^2) 169 504.6 NA NA NA
## 2 bloco + Temp 163 504.3 = 6 0.3799 0.999
## Estimativas dos parâmetros.
summary(m0)
##
## Call:
## survreg(formula = S ~ bloco + temp + I(temp^2), data = da, dist = dist)
## Value Std. Error z p
## (Intercept) 96.919 7.4390 13.03 8.43e-39
## blocoII 0.983 0.3443 2.85 4.32e-03
## blocoIII 5.784 0.0000 Inf 0.00e+00
## temp -7.681 0.6286 -12.22 2.46e-34
## I(temp^2) 0.166 0.0131 12.64 1.27e-36
## Log(scale) 0.607 0.0703 8.63 5.91e-18
##
## Scale= 1.83
##
## Gaussian distribution
## Loglik(model)= -252.3 Loglik(intercept only)= -390.3
## Chisq= 276 on 4 degrees of freedom, p= 0
## Number of Newton-Raphson Iterations: 8
## n= 175
##-----------------------------------------------------------------------------
## Estimativa pontual da temperatura ótima para latência.
c0 <- coef(m0)
tmin <- -0.5*c0[4]/c0[5]
tmin
## temp
## 23.18
##-----------------------------------------------------------------------------
## Predição.
pred <- with(da, expand.grid(bloco="I", temp=seq(15, 32, by=0.25)))
pred$time <- predict(m0, newdata=pred, type="quantile", p=0.5)
##-----------------------------------------------------------------------------
## Com bandas.
levels(pred$bloco) <- c("I","II","III")
b <- coef(m0)
b <- b[-3]
F <- model.matrix(~bloco+temp+I(temp^2), pred)
F <- F[,names(b)]
U <- chol(vcov(m0)[names(b),names(b)])
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, "*")
b <- sweep(me, 1, pred$time, "+")
pred <- cbind(pred, b)
| Tabela 1: Quadro de testes de Wald para os termos de efeito fixo do modelo para o tempo até o aparecimento de pústulas em folhas. | |||
| FV | GL | chi-quadrado | Pr(>chi) |
|---|---|---|---|
| bloco | 1 | 8.15 | 0.0043 |
| temp | 1 | 149.30 | 0.0000 |
| I(temp2) | 1 | 159.78 | 0.0000 |
##-----------------------------------------------------------------------------
cap <-
"Tempo mediano estimado para latência em função da temperatura. O tempo
mediano representa o tempo necessário para que 50% das folhas apresentem
latência."
##-----------------------------------------------------------------------------
## Número acumulado de pústulas.
ylab <- "Tempo para latência (dias)"
xlab <- expression("Temperatura"~(degree*C))
xyplot(time~temp, data=pred, type="l",
xlab=xlab, ylab=ylab,
prepanel=prepanel.cbH, cty="bands",
ly=pred$lwr, uy=pred$upr,
panel=panel.cbH)+
layer(panel.abline(v=tmin, lty=2))
##-----------------------------------------------------------------------------
## Área abaixo da curva do número de pústulas.
## Para que não ocorra viés devido o comprimento de observação ter sido
## diferente para os blocos, será selecionado o menor intervalo de
## observação, no caso, 14 dias. Além do mais, as temperaturas extremas
## (15 e 32), não tiveram registro de putulas e não serão consideradas
## nessa análise devido não apresentar variância entre suas observações.
## db <- droplevels(subset(db, dia<=14 & temp>15 & temp<32))
dc <- droplevels(subset(db, dia<=14))
dc <- ddply(dc, .(bloco, temp, Temp, ue, ua, uf), summarise,
aac=aac(dia/14, value))
##-----------------------------------------------------------------------------
## Ver.
## xyplot(aac~temp, groups=bloco, data=dc)
## xyplot(sqrt(aac)~temp, groups=bloco, data=dc)
## xyplot(log(aac+1)~temp, groups=bloco, data=dc)
## xyplot(aac^(1/3)~temp, groups=bloco, data=dc)
##-----------------------------------------------------------------------------
## Análise.
## Beta generalizada.
## dc$y <- log(dc$aac+1)
## plot(y~temp, data=dc)
## start <- list(b0=0.01, tmax=35, tmin=15, thx=1, thn=1)
## with(start, curve(b0*(tmax-x)^thx*(x-tmin)^thn, add=TRUE))
## n0 <- nls(y~b0*(tmax-temp)^thx*(temp-tmin)^thn, data=dc,
## start=start)
## Não dá para ajustar.
dc$y <- dc$aac^(1/2)
start <- list(thy = 3.811, thq = -0.075, thx = 22.862, thc = 0.001)
## plot(y~temp, data=dc)
## with(start, curve(thy*exp(thq*(x-thx)^2+thc*(x-thx)^3), add=TRUE))
dc <- groupedData(y~temp|ue, data=dc, order=FALSE)
n0 <- nlme(y~thy*exp(thq*(temp-thx)^2+thc*(temp-thx)^3),
random=thy~1,
fixed=thy+thq+thx+thc~1,
data=dc,
start=unlist(start),
weights=varComb(varIdent(form=~1), varPower(value=0.7,
form=~fitted(.))))
summary(n0)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: y ~ thy * exp(thq * (temp - thx)^2 + thc * (temp - thx)^3)
## Data: dc
## AIC BIC logLik
## 361.9 384.1 -174
##
## Random effects:
## Formula: thy ~ 1 | ue
## thy Residual
## StdDev: 4.4e-07 0.9189
##
## Combination of variance functions:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 0.6114
## Fixed effects: thy + thq + thx + thc ~ 1
## Value Std.Error DF t-value p-value
## thy 6.282 0.4409 127 14.25 0.0000
## thq -0.086 0.0068 127 -12.54 0.0000
## thx 22.883 0.1873 127 122.17 0.0000
## thc 0.001 0.0013 127 0.92 0.3614
## Correlation:
## thy thq thx
## thq -0.623
## thx -0.198 0.353
## thc 0.297 -0.572 -0.750
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.1952 -0.4996 -0.1969 0.1741 7.2246
##
## Number of Observations: 175
## Number of Groups: 45
intervals(n0)
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## thy 5.41977 6.282235 7.144698
## thq -0.09926 -0.085866 -0.072471
## thx 22.51642 22.882804 23.249190
## thc -0.00135 0.001189 0.003727
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: ue
## lower est. upper
## sd(thy) 7.891e-84 4.4e-07 2.454e+70
##
## Variance function:
## lower est. upper
## B.power 0.5409 0.6114 0.6819
## attr(,"label")
## [1] "Variance function:"
##
## Within-group standard error:
## lower est. upper
## 0.8216 0.9189 1.0277
##-----------------------------------------------------------------------------
## Predição.
pred <- data.frame(temp=seq(15, 32, by=0.25))
pred$y <- predict(n0, newdata=pred, level=0)
## xyplot(y~temp, data=dc)+
## as.layer(xyplot(y~temp, data=pred, type="l"))
##-----------------------------------------------------------------------------
## Inclusão das bandas.
model <- deriv3(~thy*exp(thq*(temp-thx)^2+thc*(temp-thx)^3),
c("thy","thq","thx","thc"),
function(temp, thy,thq,thx,thc){NULL})
pred2 <- pred
m <- do.call(model, args=c(list(temp=pred2$temp),
as.list(fixef(n0))))
F <- attr(m, "gradient")
U <- chol(vcov(n0))
pred2$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(pred2$se, zval, "*")
b <- sweep(me, 1, m, "+")
pred2 <- cbind(pred2, b)
## str(pred2)
| Tabela 2: Quadro de de estimativas dos parâmetros do modelo de regressão não linear para a curva do número médio de pústulas em função da temperatura. | |||||
| Parâmetro | Value | Std.Error | DF | t-value | p-value |
|---|---|---|---|---|---|
| thy | 6.2822 | 0.4409 | 127.0000 | 14.2482 | 0.0000 |
| thq | -0.0859 | 0.0068 | 127.0000 | -12.5392 | 0.0000 |
| thx | 22.8828 | 0.1873 | 127.0000 | 122.1678 | 0.0000 |
| thc | 0.0012 | 0.0013 | 127.0000 | 0.9160 | 0.3614 |
##-----------------------------------------------------------------------------
cap <-
"Raíz do número médio de pústulas ao longo das avaliações em função da
temperatura. Linha representa o valor ajustado pelo modelo não linear
considerado com a respectiva banda de confiança de 95%."
##-----------------------------------------------------------------------------
## Gráfico.
xlab <- "Temperatura"
ylab <- "Raíz do número médio de pústulas"
xyplot(y~temp, data=dc, type="p", xlab=xlab, ylab=ylab)+
as.layer(xyplot(fit~temp, data=pred2, type="l",
prepanel=prepanel.cbH, cty="bands",
ly=pred2$lwr, uy=pred2$upr,
panel=panel.cbH
))
##-----------------------------------------------------------------------------
sessionInfo()
## R version 3.1.0 (2014-04-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 stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] nlme_3.1-117 reshape_0.8.5 plyr_1.8.1 multcomp_1.3-3
## [5] TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-10 MASS_7.3-33
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] aod_1.3 cluster_1.15.2 evaluate_0.5.5 formatR_0.10
## [5] Formula_1.1-1 Gmisc_0.6.2.3 grid_3.1.0 Hmisc_3.14-4
## [9] lme4_1.1-6 Matrix_1.1-4 methods_3.1.0 minqa_1.2.3
## [13] Rcpp_0.11.2 RcppEigen_0.3.2.0.2 sandwich_2.3-0 sp_1.0-15
## [17] stringr_0.6.2 tools_3.1.0 zoo_1.7-10
cbind(Sys.info())
## [,1]
## sysname "Linux"
## release "3.13.0-30-generic"
## version "#54-Ubuntu SMP Mon Jun 9 22:45:01 UTC 2014"
## nodename "brother"
## machine "x86_64"
## login "walmes"
## user "walmes"
## effective_user "walmes"
Sys.time()
## [1] "2014-06-30 17:22:23 BRT"