Walmes Zeviani
##-----------------------------------------------------------------------------
## Definições da sessão.
## rm(list=ls())
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "plyr", "aod",
"survival")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra doBy multcomp plyr aod survival
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## Funções extras carregadas do dropbox.
fun <- c("apc.R", "desdobrglht.R", "equallevels.R")
dropbox <- "http://dl.dropboxusercontent.com/u/48140237/"
sapply(paste0(dropbox, fun), source)
## http://dl.dropboxusercontent.com/u/48140237/apc.R
## value ?
## visible FALSE
## http://dl.dropboxusercontent.com/u/48140237/desdobrglht.R
## value ?
## visible FALSE
## http://dl.dropboxusercontent.com/u/48140237/equallevels.R
## value ?
## visible FALSE
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] methods splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] aod_1.3 plyr_1.8.1 multcomp_1.3-3 TH.data_1.0-3
## [5] mvtnorm_0.9-9997 doBy_4.5-10 MASS_7.3-33 survival_2.37-7
## [9] latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29 knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.0 lme4_1.1-6
## [5] Matrix_1.1-3 minqa_1.2.3 nlme_3.1-117 Rcpp_0.11.2
## [9] RcppEigen_0.3.2.0.2 sandwich_2.3-0 stringr_0.6.2 tools_3.1.0
## [13] zoo_1.7-10
trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Leitura dos dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/estaca_viabilidade.txt",
header=TRUE, sep="\t")
da$trat <- factor(da$trat)
str(da)
## 'data.frame': 200 obs. of 6 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ esq : int 7 9 9 7 9 9 3 9 9 9 ...
## $ dir : int 9 9 9 9 9 9 5 9 9 9 ...
## $ status: int 3 0 0 3 0 0 3 0 0 0 ...
##-----------------------------------------------------------------------------
## Ver.
dotplot(rept~esq+dir|est+trat, data=da,
ylab="Estaca", xlab="Dias")+
as.layer(segplot(rept~esq+dir|est+trat, data=da))
##-----------------------------------------------------------------------------
## Criando vetor que contém informação dos tempos para os eventos e as
## censuras.
S <- with(da,
Surv(time=esq,
time2=dir,
event=status,
type="interval"))
head(S, 10)
## [1] [7, 9] 9+ 9+ [7, 9] 9+ 9+ [3, 5] 9+ 9+ 9+
## Opções são várias: "exponential", "gamma", etc.
dist <- "weibull"
## survreg's scale = 1/(rweibull shape)
## survreg's intercept = log(rweibull scale)
## Modelo com interação.
s0 <- survreg(S~est*trat, data=da, dist=dist)
## Só efeitos aditivos (interação abandonada).
s1 <- survreg(S~est+trat, data=da, dist=dist)
## Só efeito de estação (trat abandonado).
s2 <- survreg(S~est, data=da, dist=dist)
## Sem efeitos de termos experimentais.
s3 <- survreg(S~1, data=da, dist=dist)
## Comparando os modelos o que implica em testar os termos abandonados
## de um para o outro.
anova(s3, s2, s1, s0)
## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 1 198 594.5 NA NA NA
## 2 est 195 521.3 = 3 73.183 8.882e-16
## 3 est + trat 191 518.8 = 4 2.449 6.538e-01
## 4 est * trat 179 496.0 = 12 22.798 2.949e-02
## Outra forma e aplicar a função dropterm() que testa o abandono de
## termos de efeito fixo mas sem fazer reajustes, via aproximação de Wald.
dropterm(s0, test="Chisq")
## Warning: no 'nobs' method is available
## Warning: no 'nobs' method is available
## Single term deletions
##
## Model:
## S ~ est * trat
## Df AIC LRT Pr(Chi)
## <none> 538
## est:trat 12 537 22.8 0.029 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dropterm(s1, test="Chisq")
## Warning: no 'nobs' method is available
## Warning: no 'nobs' method is available
## Warning: no 'nobs' method is available
## Single term deletions
##
## Model:
## S ~ est + trat
## Df AIC LRT Pr(Chi)
## <none> 537
## est 3 605 74.1 5.7e-16 ***
## trat 4 531 2.4 0.65
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Estimativas dos parâmetros.
summary(s0)
##
## Call:
## survreg(formula = S ~ est * trat, data = da, dist = dist)
## Value Std. Error z p
## (Intercept) 2.4191 0.1275 18.979 2.52e-80
## estoutono 0.1631 0.1514 1.077 2.81e-01
## estprimavera -0.4405 0.1600 -2.754 5.89e-03
## estverão -0.3218 0.1530 -2.103 3.55e-02
## trat2 -0.1975 0.1552 -1.272 2.03e-01
## trat3 0.1012 0.1989 0.509 6.11e-01
## trat4 -0.1851 0.1599 -1.158 2.47e-01
## trat5 -0.2252 0.1519 -1.482 1.38e-01
## estoutono:trat2 0.2505 0.1976 1.267 2.05e-01
## estprimavera:trat2 0.4964 0.2423 2.049 4.05e-02
## estverão:trat2 0.2183 0.1935 1.128 2.59e-01
## estoutono:trat3 -0.1280 0.2292 -0.558 5.77e-01
## estprimavera:trat3 -0.3487 0.2358 -1.479 1.39e-01
## estverão:trat3 -0.0983 0.2295 -0.428 6.68e-01
## estoutono:trat4 0.1734 0.1982 0.875 3.82e-01
## estprimavera:trat4 0.1131 0.2073 0.546 5.85e-01
## estverão:trat4 0.3794 0.2069 1.833 6.67e-02
## estoutono:trat5 0.1745 0.1896 0.920 3.57e-01
## estprimavera:trat5 0.2087 0.2018 1.034 3.01e-01
## estverão:trat5 0.3218 0.1934 1.664 9.62e-02
## Log(scale) -1.5258 0.0884 -17.263 8.92e-67
##
## Scale= 0.217
##
## Weibull distribution
## Loglik(model)= -248 Loglik(intercept only)= -297.2
## Chisq= 98.43 on 19 degrees of freedom, p= 1e-12
## Number of Newton-Raphson Iterations: 6
## n= 200
##-----------------------------------------------------------------------------
## Obter as estimativas para cada tratamento em cada estação. A lm()
## aqui é só para obter a matriz.
## Médias para a estação primavera.
X <- LSmatrix(lm(dir~est*trat, data=da), effect=c("trat","est"))
grid <- attr(X, "grid")
L <- by(X, INDICES=grid$est, FUN=as.matrix)
L <- lapply(L, "rownames<-", levels(da$trat))
g <- lapply(L, desdobr.glht, m0=s0, focus="trat", test="fdr")
grid <- ldply(g); names(grid)[1] <- "est"
grid <- equallevels(grid, da)
grid
## est trat estimate lwr upr cld
## 1 inverno 1 2.419 2.169 2.669 a
## 2 inverno 2 2.222 2.045 2.398 a
## 3 inverno 3 2.520 2.213 2.827 a
## 4 inverno 4 2.234 2.042 2.426 a
## 5 inverno 5 2.194 2.030 2.358 a
## 6 outono 1 2.582 2.419 2.745 a
## 7 outono 2 2.635 2.459 2.811 a
## 8 outono 3 2.555 2.403 2.708 a
## 9 outono 4 2.570 2.408 2.733 a
## 10 outono 5 2.531 2.379 2.684 a
## 11 primavera 1 1.979 1.787 2.170 ab
## 12 primavera 2 2.278 1.970 2.585 a
## 13 primavera 3 1.731 1.577 1.885 b
## 14 primavera 4 1.907 1.731 2.082 ab
## 15 primavera 5 1.962 1.785 2.139 ab
## 16 verão 1 2.097 1.934 2.261 a
## 17 verão 2 2.118 1.963 2.273 a
## 18 verão 3 2.100 1.946 2.255 a
## 19 verão 4 2.292 2.098 2.485 a
## 20 verão 5 2.194 2.030 2.358 a
## Cuidado! Essas estimativas ainda não são as médias ou medianas de
## tempo para tornar inviável. Estão na escala do preditor linear de uma
## particular parametrização da distribuição Weibull. Operações tem que
## ser feitas para chegar aos valores desejados ou interpretáveis.
## Com essas igualdades podemos chegar às curvas de "sobrevivência".
## survreg's scale = 1/(rweibull shape)
## survreg's intercept = log(rweibull scale)
grid$wscale <- exp(grid$estimate)
grid$wshape <- 1/s0$scale
##-----------------------------------------------------------------------------
## As curvas de sobrevivência.
time <- seq(0, 15, by=0.5)
cur <- ddply(grid, .(est, trat), summarise,
time=time, prob=pweibull(time, shape=wshape, scale=wscale))
str(cur)
## 'data.frame': 620 obs. of 4 variables:
## $ est : Factor w/ 4 levels "inverno","outono",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ time: num 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 ...
## $ prob: num 0.00 6.08e-07 1.47e-05 9.51e-05 3.57e-04 ...
## Curvas de "sobrevivência" que expressam a probabilidade de ficar
## inviável como função do tempo.
xyplot(prob~time|est, groups=trat, data=cur, type="l",
auto.key=list(columns=3, lines=TRUE, points=FALSE),
xlab="Dias", ylab=expression(Pr(Inviável)))
