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 vitro via discos de folhas em função da temperatura, com 8 níveis. Foram consideradas 4 repetições para cada temperatura. O experimento foi realizado duas vezes sob o mesmo delineamento (DIC) e níveis de temperatura. 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("disco_pustula.txt", header=TRUE, sep="\t")
str(da)
## 'data.frame': 64 obs. of 20 variables:
## $ bloco: Factor w/ 2 levels "I","II": 1 1 1 1 1 1 1 1 1 1 ...
## $ temp : int 18 18 18 18 20 20 20 20 22 22 ...
## $ rept : int 1 2 3 4 1 2 3 4 1 2 ...
## $ dia3 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ dia4 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ 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 0 0 0 0 0 0 0 0 ...
## $ dia10: int 0 0 0 0 0 0 0 0 21 8 ...
## $ dia11: int 0 0 0 0 0 0 0 4 82 32 ...
## $ dia12: int 0 0 0 0 0 3 5 7 228 123 ...
## $ dia13: int 0 0 0 0 4 15 10 20 451 468 ...
## $ dia14: int 1 0 0 0 47 79 18 29 510 557 ...
## $ dia15: int 3 0 0 1 74 162 27 110 565 563 ...
## $ dia16: int 3 5 0 2 123 238 40 163 570 580 ...
## $ dia17: int 7 8 0 4 273 344 88 247 570 580 ...
## $ dia18: int 7 8 0 4 273 344 88 247 570 580 ...
## $ dia19: int 7 8 7 7 273 344 88 247 570 580 ...
##-----------------------------------------------------------------------------
## Conversões.
## Criar os níveis de unidade experimental.
da <- transform(da,
ue=interaction(bloco,temp,rept),
rept=interaction(bloco,rept),
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)+2
tot <- length(x)+2
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"),
measure.vars=diacol, variable_name="dia")
db$dia <- as.integer(gsub("\\D", "", db$dia))
db <- na.omit(db)
str(db)
## 'data.frame': 1088 obs. of 6 variables:
## $ bloco: Factor w/ 2 levels "I","II": 1 1 1 1 1 1 1 1 1 1 ...
## $ temp : int 18 18 18 18 20 20 20 20 22 22 ...
## $ Temp : Factor w/ 8 levels "18","20","22",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ ue : Factor w/ 64 levels "I.18.1","II.18.1",..: 1 17 33 49 3 19 35 51 5 21 ...
## $ dia : int 3 3 3 3 3 3 3 3 3 3 ...
## $ value: int 0 0 0 0 0 0 0 0 0 0 ...
da <- da[,-diacol]
## str(da)
##-----------------------------------------------------------------------------
cap <-
"Tempo para o aparecimento de pustulas em discos 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 <- "Unidades experimentais"
ylab <- "Tempo para o aparecimento de pustulas nos discos"
dotplot(time2+time~rept|Temp, data=subset(da), layout=c(NA,1),
scales=list(x=list(relation="free", draw=FALSE)),
xlab=xlab, ylab=ylab)+
as.layer(
segplot(rept~time2+time|Temp, data=subset(da), draw=FALSE,
layout=c(NA,1), horizontal=FALSE,
scales=list(x=list(relation="free", draw=FALSE))))
##-----------------------------------------------------------------------------
cap <-
"Número de pustulas 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 pustulas.
ylab <- "Pustulas"
xlab <- "Dias de avaliação"
xyplot(value~dia|Temp, groups=ue, 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 <- "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) 59 281.8 NA NA NA
## 2 bloco + Temp 54 259.9 = 5 21.89 0.0005484
## Estimativas dos parâmetros.
summary(m1)
##
## Call:
## survreg(formula = S ~ bloco + Temp, data = da, dist = dist)
## Value Std. Error z p
## (Intercept) 15.341 0.7257 21.14 3.45e-99
## blocoII 0.410 0.4827 0.85 3.95e-01
## Temp20 -2.546 0.9662 -2.64 8.40e-03
## Temp22 -5.046 0.9662 -5.22 1.76e-07
## Temp24 -6.796 0.9662 -7.03 2.00e-12
## Temp25 -7.296 0.9662 -7.55 4.30e-14
## Temp26 -8.046 0.9662 -8.33 8.23e-17
## Temp28 -7.921 0.9662 -8.20 2.43e-16
## Temp30 -0.692 0.9738 -0.71 4.77e-01
## Log(scale) 0.644 0.0948 6.79 1.10e-11
##
## Scale= 1.9
##
## Gaussian distribution
## Loglik(model)= -129.9 Loglik(intercept only)= -169.9
## Chisq= 79.97 on 8 degrees of freedom, p= 5e-14
## Number of Newton-Raphson Iterations: 5
## n= 64
##-----------------------------------------------------------------------------
## Predição.
pred <- with(da, expand.grid(bloco="I", Temp=levels(Temp)))
pred$time <- predict(m1, newdata=pred, type="quantile", p=0.5)
X <- LSmatrix(lm(time2~bloco+Temp, data=da), effect=c("Temp"))
rownames(X) <- levels(da$Temp)
grid <- equallevels(attr(X, "grid"), da)
## ci <- confint(glht(m1, linfct=X), calpha=univariate_calpha())$confint
## grid <- cbind(grid, ci)
grid <- desdobr.glht(X, m0=m1, focus="Temp", test="fdr")
grid$temp <- as.integer(as.character(grid$Temp))
| 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 discos. | |||
| FV | GL | chi-quadrado | Pr(>chi) |
|---|---|---|---|
| bloco | 1 | 0.72 | 0.3952 |
| Temp | 7 | 161.07 | 0.0000 |
| Tabela 2: Estimativas acompanhadas do intervalo de confiança de 95% para o tempo mediano para aparecimento de pústulas em função da temperatura. Estimativas seguidas de letras iguais indicam contraste entre médias não diferente de zero. | ||||||
| Temp | estimate | lwr | upr | cld | temp | |
|---|---|---|---|---|---|---|
| 1 | 18 | 16 | 14.20 | 16.89 | a | 18 |
| 2 | 20 | 13 | 11.67 | 14.33 | b | 20 |
| 3 | 22 | 11 | 9.17 | 11.83 | c | 22 |
| 4 | 24 | 9 | 7.42 | 10.08 | cd | 24 |
| 5 | 25 | 8 | 6.92 | 9.58 | d | 25 |
| 6 | 26 | 8 | 6.17 | 8.83 | d | 26 |
| 7 | 28 | 8 | 6.29 | 8.96 | d | 28 |
| 8 | 30 | 15 | 13.50 | 16.21 | ab | 30 |
##-----------------------------------------------------------------------------
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. Barras de erro representam o intervalo de confiança de 95%."
##-----------------------------------------------------------------------------
## Número acumulado de pustulas.
ylab <- "Tempo para latência (dias)"
xlab <- expression("Temperatura"~(degree*C))
## xyplot(time~as.integer(as.character(Temp)), data=pred,
## xlab=xlab, ylab=ylab)
xyplot(estimate~temp, data=grid,
xlab=xlab, ylab=ylab,
ly=grid$lwr, uy=grid$upr, cty="bars",
prepanel=prepanel.cbH, panel=panel.cbH)
##-----------------------------------------------------------------------------
## Área abaixo da curva do número de pustulas.
## O bloco dois não será considerado.
dc <- ddply(subset(db, bloco=="I"), .(temp, Temp, ue), summarise,
aac=aac(dia/(19-3), value))
##-----------------------------------------------------------------------------
## Ver.
## xyplot(aac~temp, data=dc)
## xyplot(sqrt(aac)~temp, data=dc)
## xyplot(log(aac+1)~temp, data=dc)
## xyplot(aac^(1/3)~temp, data=dc)
##-----------------------------------------------------------------------------
## Análise.
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))
n0 <- gnls(y~thy*exp(thq*(temp-thx)^2+thc*(temp-thx)^3),
params=thy+thq+thx+thc~1,
data=dc,
start=unlist(start),
weights=varComb(varIdent(form=~1), varPower(value=0.7,
form=~fitted(.))))
summary(n0)
## Generalized nonlinear least squares fit
## Model: y ~ thy * exp(thq * (temp - thx)^2 + thc * (temp - thx)^3)
## Data: dc
## AIC BIC logLik
## 139.4 148.2 -63.7
##
## Combination of variance functions:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 0.4439
##
## Coefficients:
## Value Std.Error t-value p-value
## thy 17.776 0.9098 19.54 0.0000
## thq -0.073 0.0058 -12.53 0.0000
## thx 23.981 0.2243 106.89 0.0000
## thc -0.001 0.0017 -0.47 0.6428
##
## Correlation:
## thy thq thx
## thq -0.540
## thx -0.173 -0.033
## thc 0.114 0.100 -0.807
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.2718 -0.5804 -0.0802 0.7452 1.8118
##
## Residual standard error: 0.8334
## Degrees of freedom: 32 total; 28 residual
intervals(n0)
## Approximate 95% confidence intervals
##
## Coefficients:
## lower est. upper
## thy 15.912167 17.7757201 19.639274
## thq -0.084683 -0.0727806 -0.060878
## thx 23.521461 23.9810113 24.440561
## thc -0.004181 -0.0007786 0.002624
## attr(,"label")
## [1] "Coefficients:"
##
## Variance function:
## lower est. upper
## B.power 0.2028 0.4439 0.6849
## attr(,"label")
## [1] "Variance function:"
##
## Residual standard error:
## lower est. upper
## 0.4688 0.7796 1.2964
##-----------------------------------------------------------------------------
## Predição.
pred <- data.frame(temp=seq(16, 32, by=0.25))
pred$y <- predict(n0, newdata=pred)
## 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(coef(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 3: 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 | t-value | p-value | lower | upper |
|---|---|---|---|---|---|---|
| thy | 17.7757 | 0.9098 | 19.5390 | 0.0000 | 15.9122 | 19.6393 |
| thq | -0.0728 | 0.0058 | -12.5258 | 0.0000 | -0.0847 | -0.0609 |
| thx | 23.9810 | 0.2243 | 106.8935 | 0.0000 | 23.5215 | 24.4406 |
| thc | -0.0008 | 0.0017 | -0.4688 | 0.6428 | -0.0042 | 0.0026 |
##-----------------------------------------------------------------------------
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 pustulas"
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:00 BRT"