Grasiela Bruzamarello Tognon
Walmes Marques Zeviani
Francine Lorena Cuquel
O experimento avaliou o efeito dos fator calor de campo {com, sem}. As respostas avaliadas foram medidas de características de coloração das estacas (L, a, b, E), o peso e a classificação quanto à viabilidade comercial (nota de 1 à 3). 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: 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.2 rootSolve_1.6.5 multcomp_1.3-3 TH.data_1.0-3
## [5] mvtnorm_0.9-9997 plyr_1.8.1 doBy_4.5-10 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.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 lme4_1.1-6 Matrix_1.1-4
## [5] methods_3.1.1 minqa_1.2.3 Rcpp_0.11.2 RcppEigen_0.3.2.0.2
## [9] sandwich_2.3-0 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("calor_milleflora.txt", header=TRUE, sep="\t")
levels(da$calor) <- c("Com","Sem")
str(da)
## 'data.frame': 100 obs. of 9 variables:
## $ dias : int 1 1 1 1 1 1 1 1 1 1 ...
## $ calor : Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ L : num 41 37.5 42.1 37.6 41.6 ...
## $ a : num -12.72 -10.28 -11.39 -9.39 -12.28 ...
## $ b : num 20.2 14.5 18.9 14.6 18.6 ...
## $ E : num 30.9 23.8 29.9 28.5 29.8 ...
## $ peso : num 49.6 55.3 37.1 32.5 38.5 ...
## $ classif: int 1 2 1 1 1 1 2 1 1 1 ...
##-----------------------------------------------------------------------------
## Ver.
xlab <- "Dia de avaliação"
##--------------------------------------------
## Variável L.
xyplot(L~dias, groups=calor, data=da, jitter.x=TRUE,
auto.key=list(columns=2, title="Calor"), xlab=xlab,
type=c("p","a"),
par.settings=ps)
## xyplot(L~dias|calor, groups=rept, data=da, type=c("o"),
## par.settings=ps, xlab=xlab)
##--------------------------------------------
## Variável a.
xyplot(a~dias, groups=calor, data=da, jitter.x=TRUE,
auto.key=list(columns=2, title="Calor"), xlab=xlab,
type=c("p","a"),
par.settings=ps)
## xyplot(a~dias, groups=calor,
## data=da[-c(60,74,80,94,100),], # sem pontos atípicos
## jitter.x=TRUE,
## auto.key=list(columns=2, title="Calor"), xlab=xlab,
## type=c("p","a"),
## par.settings=ps)
## xyplot(a~dias|calor, groups=rept, data=da, type=c("o"),
## par.settings=ps, xlab=xlab)
##--------------------------------------------
## Variável b.
xyplot(b~dias, groups=calor, data=da, jitter.x=TRUE,
auto.key=list(columns=2, title="Calor"), xlab=xlab,
type=c("p","a"),
par.settings=ps)
## xyplot(b~dias|calor, groups=rept, data=da, type=c("o"),
## par.settings=ps, xlab=xlab)
##--------------------------------------------
## Variável E.
xyplot(E~dias, groups=calor, data=da, jitter.x=TRUE,
auto.key=list(columns=2, title="Calor"), xlab=xlab,
type=c("p","a"),
par.settings=ps)
## xyplot(E~dias|calor, groups=rept, data=da, type=c("o"),
## par.settings=ps, xlab=xlab)
##--------------------------------------------
## Peso.
xyplot(peso~dias, groups=calor, data=da, jitter.x=TRUE,
auto.key=list(columns=2, title="Calor"), xlab=xlab,
type=c("p","a"),
par.settings=ps)
## xyplot(peso~dias|calor, groups=rept, data=da, type=c("o"),
## par.settings=ps, xlab=xlab)
##--------------------------------------------
## Classificação quanto a viabilidade.
xyplot(classif~dias, groups=calor, data=da, jitter.x=TRUE,
auto.key=list(columns=2, title="Calor"), xlab=xlab,
type=c("p","a"),
par.settings=ps)
## xyplot(classif~dias|calor, groups=rept, data=da, type=c("o"),
## par.settings=ps, xlab=xlab)
##-----------------------------------------------------------------------------
## Frequência de estacas viáveis/inviáveis em cada dia.
da$viab <- with(da, cut(classif, c(-Inf,2,Inf), right=TRUE,
labels=c("viável","inviável")))
str(da)
## 'data.frame': 100 obs. of 10 variables:
## $ dias : int 1 1 1 1 1 1 1 1 1 1 ...
## $ calor : Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ L : num 41 37.5 42.1 37.6 41.6 ...
## $ a : num -12.72 -10.28 -11.39 -9.39 -12.28 ...
## $ b : num 20.2 14.5 18.9 14.6 18.6 ...
## $ E : num 30.9 23.8 29.9 28.5 29.8 ...
## $ peso : num 49.6 55.3 37.1 32.5 38.5 ...
## $ classif: int 1 2 1 1 1 1 2 1 1 1 ...
## $ viab : Factor w/ 2 levels "viável","inviável": 1 1 1 1 1 1 1 1 1 1 ...
tb <- xtabs(~calor+dias+viab, data=da)
plot(tb, xlab="Calor", ylab="Dias de avaliação")
##-----------------------------------------------------------------------------
## Análise do peso.
## Testar segmentado linear.
## Espera-se que o sem calor de campo tenha um atraso na queda de peso, o
## que seria representado por um plato linear enquanto que o com calor de
## campo seria apenas linear.
str(da)
## 'data.frame': 100 obs. of 10 variables:
## $ dias : int 1 1 1 1 1 1 1 1 1 1 ...
## $ calor : Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ L : num 41 37.5 42.1 37.6 41.6 ...
## $ a : num -12.72 -10.28 -11.39 -9.39 -12.28 ...
## $ b : num 20.2 14.5 18.9 14.6 18.6 ...
## $ E : num 30.9 23.8 29.9 28.5 29.8 ...
## $ peso : num 49.6 55.3 37.1 32.5 38.5 ...
## $ classif: int 1 2 1 1 1 1 2 1 1 1 ...
## $ viab : Factor w/ 2 levels "viável","inviável": 1 1 1 1 1 1 1 1 1 1 ...
da$y <- da$peso
da$ue <- with(da, interaction(calor, rept, drop=TRUE))
nlevels(da$ue)
## [1] 20
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 536.7 552.4 -262.4
## m2 2 8 523.8 544.7 -253.9 1 vs 2 16.920 0.0002
## m0 3 12 524.7 555.9 -250.3 2 vs 3 7.174 0.1270
m0 <- m2
## Interação de calor com o tempo nos termos quadráticos.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 76 325.2 <.0001
## calor 1 18 4.3 0.053
## poly(dias, 2) 2 76 34.2 <.0001
## calor:poly(dias, 2) 2 76 3.8 0.027
## Estimativas dos termos de efeito fixo.
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 30.720 2.722 76 11.287 6.495e-18
## calorSem 7.971 3.849 18 2.071 5.301e-02
## poly(dias, 2)1 -11.534 2.847 76 -4.052 1.216e-04
## poly(dias, 2)2 -2.931 2.847 76 -1.030 3.065e-01
## calorSem:poly(dias, 2)1 -6.883 4.026 76 -1.710 9.140e-02
## calorSem:poly(dias, 2)2 -8.686 4.026 76 -2.157 3.413e-02
## 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.
## start <- list(b0=40, tx=-5/4, xb=4.5)
## plot(y~dias, data=subset(da, calor=="Sem"))
## with(start, curve(b0+tx*(x-xb)*(x>xb), add=TRUE))
##
## n0 <- nls(y~b0+tx*(dias-xb)*(dias>xb),
## data=subset(db, calor=="Sem"),
## start=start)
## coef(n0)
##
## with(as.list(coef(n0)), curve(b0+tx*(x-xb)*(x>xb), add=TRUE, col=2))
##
## start <- list(b0=30, tx=-5/4, xb=1.5)
## plot(y~dias, data=subset(da, calor=="Com"))
## with(start, curve(b0+tx*(x-xb)*(x>xb), add=TRUE))
##
## n1 <- nls(y~b0+tx*(dias-xb)*(dias>xb),
## data=subset(db, calor=="Com"),
## start=start)
## coef(n1)
##
## with(as.list(coef(n1)), curve(b0+tx*(x-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(ch)
## Valores iniciais obtidos com o ajuste indivídual.
ch <- c(40.61, -8.68900000136728, -1.166,
0.62299999863557, 4.25700400189413, -0.944542355400808)
##-----------------------------------------------------------------------------
## Ajuste dos modelos não lineares mistos.
## Modelo livre, cada valor com seus parâmetros.
n11 <- nlme(y~b0+tx*(dias-xb)*(dias>xb),
data=db,
random=b0~1, fixed=b0+tx+xb~calor,
start=list(fixed=ch),
verbose=TRUE)
##
## **Iteration 1
## LME step: Loglik: -254.5 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.44
##
## PNLS step: RSS = 386.7
## fixed effects:31.92 8.689 -0.543 -0.623 3.312 0.9444
## iterations: 4
##
## Convergence:
## fixed reStruct
## 2.000e+00 1.268e-09
##
## **Iteration 2
## LME step: Loglik: -254.5 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.44
##
## PNLS step: RSS = 386.7
## fixed effects:31.92 8.689 -0.543 -0.623 3.312 0.9444
## iterations: 1
##
## Convergence:
## fixed reStruct
## 0 0
summary(n11)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: y ~ b0 + tx * (dias - xb) * (dias > xb)
## Data: db
## AIC BIC logLik
## 525 545.9 -254.5
##
## Random effects:
## Formula: b0 ~ 1 | ue
## b0.(Intercept) Residual
## StdDev: 8.298 1.966
##
## Fixed effects: b0 + tx + xb ~ calor
## Value Std.Error DF t-value p-value
## b0.(Intercept) 31.92 2.744 75 11.632 0.0000
## b0.calorSem 8.69 3.881 75 2.239 0.0281
## tx.(Intercept) -0.54 0.227 75 -2.395 0.0191
## tx.calorSem -0.62 0.321 75 -1.943 0.0558
## xb.(Intercept) 3.31 1.880 75 1.762 0.0822
## xb.calorSem 0.94 2.018 75 0.468 0.6411
## Correlation:
## b0.(I) b0.clS tx.(I) tx.clS xb.(I)
## b0.calorSem -0.707
## tx.(Intercept) 0.000 0.000
## tx.calorSem 0.000 0.000 -0.707
## xb.(Intercept) -0.073 0.052 -0.819 0.579
## xb.calorSem 0.068 -0.071 0.763 -0.727 -0.932
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.46081 -0.54919 0.04599 0.50328 3.31075
##
## Number of Observations: 100
## Number of Groups: 20
logLik(n11)
## 'log Lik.' -254.5 (df=8)
logLik(m0)
## 'log Lik.' -253.9 (df=8)
## Modelo com a restrição de que xb para com valor é zero.
db$i <- with(db, ifelse(calor=="Com", 0, 1))
n12 <- nlme(y~b0+tx*(dias-xb*i)*(dias>xb*i),
data=db,
random=b0~1, fixed=list(b0+tx~calor, xb~1),
start=list(fixed=ch[-6]),
verbose=TRUE)
##
## **Iteration 1
## LME step: Loglik: -255.2 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.431
##
## PNLS step: RSS = 393.2
## fixed effects:32.76 7.851 -0.4078 -0.7582 4.257
## iterations: 2
##
## Convergence:
## fixed reStruct
## 2.107e+00 2.593e-09
##
## **Iteration 2
## LME step: Loglik: -255.2 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.431
##
## PNLS step: RSS = 393.2
## fixed effects:32.76 7.851 -0.4078 -0.7582 4.257
## iterations: 1
##
## Convergence:
## fixed reStruct
## 0.000e+00 1.551e-16
summary(n12)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: y ~ b0 + tx * (dias - xb * i) * (dias > xb * i)
## Data: db
## AIC BIC logLik
## 524.4 542.6 -255.2
##
## Random effects:
## Formula: b0 ~ 1 | ue
## b0.(Intercept) Residual
## StdDev: 8.297 1.983
##
## Fixed effects: list(b0 + tx ~ calor, xb ~ 1)
## Value Std.Error DF t-value p-value
## b0.(Intercept) 32.76 2.755 77 11.892 0.0000
## b0.calorSem 7.85 3.878 18 2.024 0.0580
## tx.(Intercept) -0.41 0.102 77 -4.009 0.0001
## tx.calorSem -0.76 0.249 77 -3.043 0.0032
## xb 4.26 0.735 77 5.793 0.0000
## Correlation:
## b0.(I) b0.clS tx.(I) tx.clS
## b0.calorSem -0.710
## tx.(Intercept) -0.185 0.131
## tx.calorSem 0.075 -0.054 -0.408
## xb 0.000 -0.062 0.000 -0.665
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.44077 -0.52340 0.06172 0.52120 3.06706
##
## Number of Observations: 100
## Number of Groups: 20
anova(n12, n11)
## Model df AIC BIC logLik Test L.Ratio p-value
## n12 1 7 524.4 542.6 -255.2
## n11 2 8 525.0 545.9 -254.5 1 vs 2 1.328 0.2491
## Modelo com a restrição de que as taxas são iguais.
n13 <- nlme(y~b0+tx*(dias-xb)*(dias>xb),
data=db,
random=b0~1, fixed=list(b0~calor, tx~1, xb~calor),
start=list(fixed=ch[-c(4)]),
verbose=TRUE)
##
## **Iteration 1
## LME step: Loglik: -256.5 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.415
##
## PNLS step: RSS = 406.1
## fixed effects:31.92 8.689 -0.8545 4.657 -1.4
## iterations: 3
##
## Convergence:
## fixed reStruct
## 2.0e+00 9.9e-10
##
## **Iteration 2
## LME step: Loglik: -256.5 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.415
##
## PNLS step: RSS = 406.1
## fixed effects:31.92 8.689 -0.8545 4.657 -1.4
## iterations: 1
##
## Convergence:
## fixed reStruct
## 0 0
summary(n13)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: y ~ b0 + tx * (dias - xb) * (dias > xb)
## Data: db
## AIC BIC logLik
## 527 545.2 -256.5
##
## Random effects:
## Formula: b0 ~ 1 | ue
## b0.(Intercept) Residual
## StdDev: 8.296 2.015
##
## Fixed effects: list(b0 ~ calor, tx ~ 1, xb ~ calor)
## Value Std.Error DF t-value p-value
## b0.(Intercept) 31.92 2.731 76 11.689 0.0000
## b0.calorSem 8.69 3.862 76 2.250 0.0274
## tx -0.85 0.163 76 -5.228 0.0000
## xb.(Intercept) 4.66 0.830 76 5.611 0.0000
## xb.calorSem -1.40 1.023 76 -1.368 0.1755
## Correlation:
## b0.(I) b0.clS tx xb.(I)
## b0.calorSem -0.707
## tx 0.000 0.000
## xb.(Intercept) -0.110 0.078 -0.540
## xb.calorSem 0.089 -0.127 -0.262 -0.433
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.40220 -0.57854 0.00548 0.54349 3.23198
##
## Number of Observations: 100
## Number of Groups: 20
## As taxas são as mesmas?
anova(n11, n13)
## Model df AIC BIC logLik Test L.Ratio p-value
## n11 1 8 525 545.9 -254.5
## n13 2 7 527 545.2 -256.5 1 vs 2 3.917 0.0478
AIC(n11); AIC(m0)
##
## 525
## [1] 523.8
##-----------------------------------------------------------------------------
## Verificar os pressupostos.
## Valores ajustados pelo modelo para cada UE.
plot(augPred(n12), par.settings=ps)
## DIagnóstico.
r <- residuals(n12, type="pearson")
f <- fitted(n12)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p","smooth")),
xyplot(sqrt(abs(r))~f, type=c("p","smooth")),
qqmath(residuals(n12)),
qqmath(unlist(ranef(n12))))
##-----------------------------------------------------------------------------
## 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(n12, 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 32.4 40.6 32.2 40.6 32.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.000" "dias=1.276" "dias=1.552" "dias=1.828" ...
pred$yn <- predict(n12, 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(n12)
## b0.(Intercept) b0.calorSem tx.(Intercept) tx.calorSem xb
## 32.7586 7.8514 -0.4078 -0.7582 4.2570
formula(n12)
## y ~ b0 + tx * (dias - xb * i) * (dias > xb * i)
## <environment: 0x76e0fc0>
## b0.(Intercept) b0.calorSem tx.(Intercept) tx.calorSem xb
## y ~ b0 + tx * (dias - xb * i) * (dias > xb * i)
## 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,])
}
th <- fixef(n12)
## 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 5
colnames(F)==names(th)
## [1] TRUE TRUE TRUE TRUE TRUE
head(F)
## b0.(Intercept) b0.calorSem tx.(Intercept) tx.calorSem xb
## [1,] 1 0 1.000 0 0
## [2,] 1 1 0.000 0 0
## [3,] 1 0 1.276 0 0
## [4,] 1 1 0.000 0 0
## [5,] 1 0 1.552 0 0
## [6,] 1 1 0.000 0 0
## Cálculos do IC.
U <- chol(vcov(n12))
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(n12, 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 32.4 40.6 32.2 40.6 32.1 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ yn : atomic 32.4 40.6 32.2 40.6 32.1 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ ym : atomic 32 39.9 32 40.1 32 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ se : num 2.67 2.66 2.66 2.66 2.66 ...
## $ lwr : num 27.1 35.4 27 35.4 26.9 ...
## $ fit : num 32.4 40.6 32.2 40.6 32.1 ...
## $ upr : num 37.6 45.8 37.5 45.8 37.3 ...
xyplot(y~dias, groups=calor, data=da, key=key, par.settings=ps,
ylab="Peso das estacas (g)", 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))
##-----------------------------------------------------------------------------
## Análise do L.
str(da)
## 'data.frame': 100 obs. of 12 variables:
## $ dias : int 1 1 1 1 1 1 1 1 1 1 ...
## $ calor : Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ L : num 41 37.5 42.1 37.6 41.6 ...
## $ a : num -12.72 -10.28 -11.39 -9.39 -12.28 ...
## $ b : num 20.2 14.5 18.9 14.6 18.6 ...
## $ E : num 30.9 23.8 29.9 28.5 29.8 ...
## $ peso : num 49.6 55.3 37.1 32.5 38.5 ...
## $ classif: int 1 2 1 1 1 1 2 1 1 1 ...
## $ viab : Factor w/ 2 levels "viável","inviável": 1 1 1 1 1 1 1 1 1 1 ...
## $ y : num 49.6 55.3 37.1 32.5 38.5 ...
## $ ue : Factor w/ 20 levels "Com.1","Sem.1",..: 2 4 6 8 10 12 14 16 18 20 ...
da$y <- da$L
da$ue <- with(da, interaction(calor, rept, drop=TRUE))
nlevels(da$ue)
## [1] 20
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")
anova(m1, m2, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 6 527.8 543.5 -257.9
## m2 2 8 524.5 545.4 -254.3 1 vs 2 7.317 0.0258
## m0 3 12 530.6 561.9 -253.3 2 vs 3 1.922 0.7500
m0 <- m2
## Sem interação
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 76 8381 <.0001
## calor 1 18 8 0.0123
## poly(dias, 2) 2 76 16 <.0001
## calor:poly(dias, 2) 2 76 0 0.9176
m0 <- lme(y~calor+poly(dias, 2), random=~1|ue, data=db, method="ML")
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 6 520.7 536.3 -254.3
## m2 2 8 524.5 545.4 -254.3 1 vs 2 0.183 0.9125
## summary(m0)
## VarCorr(m0)
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 43.510 0.6455 78 67.408 6.637e-71
## calorSem -2.567 0.9128 18 -2.813 1.152e-02
## poly(dias, 2)1 14.465 2.8625 78 5.053 2.797e-06
## poly(dias, 2)2 7.746 2.8625 78 2.706 8.358e-03
## plot(augPred(m0), par.settings=ps)
## plot(m0, par.settings=ps)
## plot(residuals(m0))
## qqnorm(residuals(m0))
## qqnorm(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)
## 'data.frame': 60 obs. of 3 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 ...
## $ y : atomic 42.4 39.8 42.3 39.7 42.2 ...
## ..- 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.000" "dias=1.276" "dias=1.552" "dias=1.828" ...
xyplot(y~dias, groups=calor, data=da, key=key, par.settings=ps,
ylab="Medida de coloração L", xlab="Período de observação (dias)")+
as.layer(xyplot(y~dias, groups=calor, data=pred, type="l"))
##-----------------------------------------------------------------------------
## Análise do a.
str(da)
## 'data.frame': 100 obs. of 12 variables:
## $ dias : int 1 1 1 1 1 1 1 1 1 1 ...
## $ calor : Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ L : num 41 37.5 42.1 37.6 41.6 ...
## $ a : num -12.72 -10.28 -11.39 -9.39 -12.28 ...
## $ b : num 20.2 14.5 18.9 14.6 18.6 ...
## $ E : num 30.9 23.8 29.9 28.5 29.8 ...
## $ peso : num 49.6 55.3 37.1 32.5 38.5 ...
## $ classif: int 1 2 1 1 1 1 2 1 1 1 ...
## $ viab : Factor w/ 2 levels "viável","inviável": 1 1 1 1 1 1 1 1 1 1 ...
## $ y : num 41 37.5 42.1 37.6 41.6 ...
## $ ue : Factor w/ 20 levels "Com.1","Sem.1",..: 2 4 6 8 10 12 14 16 18 20 ...
da$y <- da$a # c(60,74,80,94,100)
da$y[c(60,74,80,94,100)] <- NA
da$ue <- with(da, interaction(calor, rept, drop=TRUE))
nlevels(da$ue)
## [1] 20
db <- groupedData(data=da[!is.na(da$y),], 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")
anova(m1, m2, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 6 333.8 349.1 -160.9
## m2 2 8 337.5 357.9 -160.7 1 vs 2 0.3521 0.8386
## m0 3 12 342.5 373.2 -159.3 2 vs 3 2.9431 0.5674
m0 <- m1
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 73 2298.7 <.0001
## calor 1 18 0.2 0.6944
## poly(dias, 1) 1 73 26.3 <.0001
## calor:poly(dias, 1) 1 73 0.1 0.7196
m0 <- lme(y~calor+poly(dias, 1), random=~1|ue, data=db, method="ML")
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 5 332.0 344.7 -161.0
## m1 2 6 333.8 349.1 -160.9 1 vs 2 0.1353 0.713
## summary(m0)
## VarCorr(m0)
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) -10.2520 0.3089 74 -33.1920 3.468e-46
## calorSem -0.2608 0.4314 18 -0.6045 5.530e-01
## poly(dias, 1) 6.1969 1.2024 74 5.1539 2.050e-06
## plot(augPred(m0), par.settings=ps)
## plot(m0, par.settings=ps)
## plot(residuals(m0))
## qqnorm(residuals(m0))
## qqnorm(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)
## 'data.frame': 60 obs. of 3 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 ...
## $ y : atomic -11.1 -11.4 -11.1 -11.3 -11 ...
## ..- 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.000" "dias=1.276" "dias=1.552" "dias=1.828" ...
xyplot(y~dias, groups=calor, data=da, key=key, par.settings=ps,
ylab="Medida de coloração a", xlab="Período de observação (dias)")+
as.layer(xyplot(y~dias, groups=calor, data=pred, type="l"))
##-----------------------------------------------------------------------------
## Análise do b.
str(da)
## 'data.frame': 100 obs. of 12 variables:
## $ dias : int 1 1 1 1 1 1 1 1 1 1 ...
## $ calor : Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ L : num 41 37.5 42.1 37.6 41.6 ...
## $ a : num -12.72 -10.28 -11.39 -9.39 -12.28 ...
## $ b : num 20.2 14.5 18.9 14.6 18.6 ...
## $ E : num 30.9 23.8 29.9 28.5 29.8 ...
## $ peso : num 49.6 55.3 37.1 32.5 38.5 ...
## $ classif: int 1 2 1 1 1 1 2 1 1 1 ...
## $ viab : Factor w/ 2 levels "viável","inviável": 1 1 1 1 1 1 1 1 1 1 ...
## $ y : num -12.72 -10.28 -11.39 -9.39 -12.28 ...
## $ ue : Factor w/ 20 levels "Com.1","Sem.1",..: 2 4 6 8 10 12 14 16 18 20 ...
da$y <- da$b
da$ue <- with(da, interaction(calor, rept, drop=TRUE))
nlevels(da$ue)
## [1] 20
db <- groupedData(data=da[!is.na(da$y),], 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")
anova(m1, m2, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 6 481.3 496.9 -234.6
## m2 2 8 476.0 496.8 -230.0 1 vs 2 9.317 0.0095
## m0 3 12 480.1 511.4 -228.1 2 vs 3 3.852 0.4265
m0 <- m2
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 76 2437.7 <.0001
## calor 1 18 2.3 0.1507
## poly(dias, 2) 2 76 3.3 0.0433
## calor:poly(dias, 2) 2 76 2.6 0.0812
m0 <- lme(y~calor+poly(dias, 2), random=~1|ue, data=db, method="ML")
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 6 477.3 492.9 -232.7
## m2 2 8 476.0 496.8 -230.0 1 vs 2 5.341 0.0692
## summary(m0)
## VarCorr(m0)
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 18.709 0.5146 78 36.354 1.158e-50
## calorSem -1.104 0.7278 18 -1.517 1.467e-01
## poly(dias, 2)1 2.338 2.3097 78 1.012 3.145e-01
## poly(dias, 2)2 5.281 2.3097 78 2.287 2.494e-02
## plot(augPred(m0), par.settings=ps)
## plot(m0, par.settings=ps)
## plot(residuals(m0))
## qqnorm(residuals(m0))
## qqnorm(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)
## 'data.frame': 60 obs. of 3 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 ...
## $ y : atomic 19 17.9 18.9 17.8 18.7 ...
## ..- 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.000" "dias=1.276" "dias=1.552" "dias=1.828" ...
xyplot(y~dias, groups=calor, data=da, key=key, par.settings=ps,
ylab="Medida de coloração b", xlab="Período de observação (dias)")+
as.layer(xyplot(y~dias, groups=calor, data=pred, type="l"))
##-----------------------------------------------------------------------------
## Análise do E.
str(da)
## 'data.frame': 100 obs. of 12 variables:
## $ dias : int 1 1 1 1 1 1 1 1 1 1 ...
## $ calor : Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
## $ rept : int 1 2 3 4 5 6 7 8 9 10 ...
## $ L : num 41 37.5 42.1 37.6 41.6 ...
## $ a : num -12.72 -10.28 -11.39 -9.39 -12.28 ...
## $ b : num 20.2 14.5 18.9 14.6 18.6 ...
## $ E : num 30.9 23.8 29.9 28.5 29.8 ...
## $ peso : num 49.6 55.3 37.1 32.5 38.5 ...
## $ classif: int 1 2 1 1 1 1 2 1 1 1 ...
## $ viab : Factor w/ 2 levels "viável","inviável": 1 1 1 1 1 1 1 1 1 1 ...
## $ y : num 20.2 14.5 18.9 14.6 18.6 ...
## $ ue : Factor w/ 20 levels "Com.1","Sem.1",..: 2 4 6 8 10 12 14 16 18 20 ...
da$y <- da$E
da$ue <- with(da, interaction(calor, rept, drop=TRUE))
nlevels(da$ue)
## [1] 20
db <- groupedData(data=da[!is.na(da$y),], 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")
anova(m1, m2, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 6 526.5 542.1 -257.3
## m2 2 8 521.0 541.8 -252.5 1 vs 2 9.564 0.0084
## m0 3 12 527.8 559.1 -251.9 2 vs 3 1.133 0.8890
m0 <- m2
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 76 4417 <.0001
## calor 1 18 5 0.0361
## poly(dias, 2) 2 76 5 0.0099
## calor:poly(dias, 2) 2 76 1 0.4910
m0 <- lme(y~calor+poly(dias, 2), random=~1|ue, data=db, method="ML")
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 6 518.5 534.1 -253.2
## m2 2 8 521.0 541.8 -252.5 1 vs 2 1.513 0.4692
## summary(m0)
## VarCorr(m0)
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 30.068 0.6123 78 49.108 2.012e-60
## calorSem -1.982 0.8659 18 -2.288 3.442e-02
## poly(dias, 2)1 3.731 2.8602 78 1.305 1.959e-01
## poly(dias, 2)2 8.156 2.8602 78 2.851 5.567e-03
## plot(augPred(m0), par.settings=ps)
## plot(m0, par.settings=ps)
## plot(residuals(m0))
## qqnorm(residuals(m0))
## qqnorm(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)
## 'data.frame': 60 obs. of 3 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 ...
## $ y : atomic 30.5 28.5 30.3 28.3 30.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.000" "dias=1.276" "dias=1.552" "dias=1.828" ...
xyplot(y~dias, groups=calor, data=da, key=key, par.settings=ps,
ylab="Medida de coloração E", xlab="Período de observação (dias)")+
as.layer(xyplot(y~dias, groups=calor, data=pred, type="l"))
##-----------------------------------------------------------------------------
## Análise de sobreviência para a variável classificação.
## Dado com censura intervalar.
da$muda <- as.numeric(da$classif>2)
aggregate(muda~dias+calor, da, mean)
## dias calor muda
## 1 1 Com 0.0
## 2 3 Com 0.0
## 3 5 Com 0.3
## 4 7 Com 0.6
## 5 9 Com 1.0
## 6 1 Sem 0.0
## 7 3 Sem 0.0
## 8 5 Sem 0.0
## 9 7 Sem 0.1
## 10 9 Sem 0.8
myfun <- function(x, y){
Y <- cumsum(y)
wm <- sum(Y==0)
x <- c(x[wm+0:1])
if(length(x)==1) x <- c(0.0001, x)
matrix(x, ncol=2, dimnames=list(NULL, c("esq","dir")))
}
dc <- split(da, f=da$ue)
dc <- do.call(rbind,
lapply(dc,
function(d){
z <- myfun(x=d$dias, y=d$muda)
d <- cbind(d[1,c("calor","rept")])
d <- cbind(d, data.frame(z))
d
}))
str(dc)
## 'data.frame': 20 obs. of 4 variables:
## $ calor: Factor w/ 2 levels "Com","Sem": 1 2 1 2 1 2 1 2 1 2 ...
## $ rept : int 1 1 2 2 3 3 4 4 5 5 ...
## $ esq : int 7 7 7 5 7 7 3 7 7 9 ...
## $ dir : int 9 9 9 7 9 9 5 9 9 NA ...
na <- is.na(dc$dir)
dc$dir[na] <- dc$esq[na]
dc$status <- 3 ## Censura intervalar.
dc$status[na] <- 0 ## Censura à direita.
aggregate(cbind(esq, dir)~calor, dc, mean)
## calor esq dir
## 1 Com 5.2 7.2
## 2 Sem 7.2 8.8
##-----------------------------------------------------------------------------
## Conduzindo uma análise de sobreviência.
S <- with(dc, Surv(time=esq, time2=dir, event=status, type="interval"))
S
## [1] [7, 9] [7, 9] [7, 9] [5, 7] [7, 9] [7, 9] [3, 5] [7, 9] [7, 9] 9+ [5, 7] [7, 9]
## [13] [5, 7] [7, 9] [5, 7] [7, 9] [3, 5] [7, 9] [3, 5] 9+
dist <- "weibull"
s0 <- survreg(S~calor, data=dc, dist=dist)
s1 <- survreg(S~1, data=dc, dist=dist)
anova(s0, s1)
## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 calor 17 43.64 NA NA NA
## 2 1 18 49.53 = -1 -5.891 0.01522
summary(s0)
##
## Call:
## survreg(formula = S ~ calor, data = dc, dist = dist)
## Value Std. Error z p
## (Intercept) 1.930 0.0611 31.57 1.01e-218
## calorSem 0.217 0.0847 2.56 1.04e-02
## Log(scale) -1.842 0.2366 -7.79 6.91e-15
##
## Scale= 0.159
##
## Weibull distribution
## Loglik(model)= -21.8 Loglik(intercept only)= -24.8
## Chisq= 5.89 on 1 degrees of freedom, p= 0.015
## Number of Newton-Raphson Iterations: 5
## n= 20
## exp(1.930)
## exp(1.930+0.217)
## Tempo médio de viabilidade.
predict(s0, newdata=data.frame(calor=levels(dc$calor)), type="response")
## 1 2
## 6.891 8.560
## Tempo para que 50% de estacas da população fiquem inviáveis.
predict(s0, newdata=data.frame(calor=levels(dc$calor)),
type="quantile", p=0.5)
## 1 2
## 6.502 8.077
X <- LSmatrix(s0, effect="calor")
ci <- confint(glht(s0, linfct=X))
ci
##
## Simultaneous Confidence Intervals
##
## Fit: survreg(formula = S ~ calor, data = dc, dist = dist)
##
## Quantile = 2.236
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1 == 0 1.930 1.793 2.067
## 2 == 0 2.147 2.016 2.278
th <- exp(ci$confint[,"Estimate"])
sh <- 1/summary(s0)[[6]]
## sh <- 1/coef(s0, allCoef=TRUE)
## names(s0)
## s0$icoef
curve(pweibull(x, shape=sh, scale=th[1]), 0, 13, ann=FALSE)
curve(pweibull(x, shape=sh, scale=th[2]), add=TRUE, lty=2)
title(xlab="Dias de avaliação", ylab="Probabilidade de passar à inviável")
legend("topleft", bty="n", lty=1:2,
title="Calor", legend=levels(da$calor))
