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_tridentata.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 35.5 42.9 40.8 40.1 38.7 ...
## $ a : num -9.79 -8.67 -9.83 -10.23 -7.67 ...
## $ b : num 13.3 13.7 12.9 14.1 15.2 ...
## $ E : num 21.7 25.8 24.4 25 24.4 ...
## $ peso : num 44.2 20 24.1 29.2 32.9 ...
## $ 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 35.5 42.9 40.8 40.1 38.7 ...
## $ a : num -9.79 -8.67 -9.83 -10.23 -7.67 ...
## $ b : num 13.3 13.7 12.9 14.1 15.2 ...
## $ E : num 21.7 25.8 24.4 25 24.4 ...
## $ peso : num 44.2 20 24.1 29.2 32.9 ...
## $ 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 35.5 42.9 40.8 40.1 38.7 ...
## $ a : num -9.79 -8.67 -9.83 -10.23 -7.67 ...
## $ b : num 13.3 13.7 12.9 14.1 15.2 ...
## $ E : num 21.7 25.8 24.4 25 24.4 ...
## $ peso : num 44.2 20 24.1 29.2 32.9 ...
## $ 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 538.4 554.0 -263.2
## m2 2 8 533.1 553.9 -258.5 1 vs 2 9.316 0.0095
## m0 3 12 540.7 572.0 -258.4 2 vs 3 0.316 0.9887
m0 <- m2
## Interação de calor com o tempo nos termos quadráticos.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 76 224.52 <.0001
## calor 1 18 1.48 0.2389
## poly(dias, 2) 2 76 56.34 <.0001
## calor:poly(dias, 2) 2 76 1.82 0.1688
## Estimativas dos termos de efeito fixo.
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 21.083 2.166 76 9.734 5.309e-15
## calorSem 3.731 3.063 18 1.218 2.389e-01
## poly(dias, 2)1 -19.071 3.193 76 -5.973 7.013e-08
## poly(dias, 2)2 -4.945 3.193 76 -1.549 1.256e-01
## calorSem:poly(dias, 2)1 -7.904 4.516 76 -1.750 8.409e-02
## calorSem:poly(dias, 2)2 -3.436 4.516 76 -0.761 4.490e-01
## 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=30, 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=25, 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(27.3890000050413, -4.37300000734149, -1.47024999513089, 0.534999994359986,
## 4.0812220000996, -0.525933745750815)
## dput(fixef(n12))
ch <- c(24.4542499088563, 2.93475007676411, -0.674249988405072,
-0.796000034655532, 4.08122214761683, 0)
##-----------------------------------------------------------------------------
## 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: -259 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.093
##
## PNLS step: RSS = 484.9
## fixed effects:23.02 4.373 -0.9353 -0.535 3.555 0.5259
## iterations: 3
##
## Convergence:
## fixed reStruct
## 1.000e+00 1.517e-09
##
## **Iteration 2
## LME step: Loglik: -259 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.093
##
## PNLS step: RSS = 484.9
## fixed effects:23.02 4.373 -0.9353 -0.535 3.555 0.5259
## 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
## 534 554.8 -259
##
## Random effects:
## Formula: b0 ~ 1 | ue
## b0.(Intercept) Residual
## StdDev: 6.567 2.202
##
## Fixed effects: b0 + tx + xb ~ calor
## Value Std.Error DF t-value p-value
## b0.(Intercept) 23.016 2.2014 75 10.455 0.0000
## b0.calorSem 4.373 3.1132 75 1.405 0.1643
## tx.(Intercept) -0.935 0.2539 75 -3.683 0.0004
## tx.calorSem -0.535 0.3591 75 -1.490 0.1405
## xb.(Intercept) 3.555 1.1689 75 3.042 0.0032
## xb.calorSem 0.526 1.3488 75 0.390 0.6977
## 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.107 0.076 -0.800 0.566
## xb.calorSem 0.093 -0.107 0.693 -0.755 -0.867
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.1325 -0.4506 -0.0360 0.3514 4.3773
##
## Number of Observations: 100
## Number of Groups: 20
logLik(n11)
## 'log Lik.' -259 (df=8)
logLik(m0)
## 'log Lik.' -258.5 (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: -260.3 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.076
##
## PNLS step: RSS = 500.8
## fixed effects:24.45 2.935 -0.6743 -0.796 4.081
## iterations: 2
##
## Convergence:
## fixed reStruct
## 4.930e-08 1.808e-09
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
## 534.6 552.8 -260.3
##
## Random effects:
## Formula: b0 ~ 1 | ue
## b0.(Intercept) Residual
## StdDev: 6.565 2.238
##
## Fixed effects: list(b0 + tx ~ calor, xb ~ 1)
## Value Std.Error DF t-value p-value
## b0.(Intercept) 24.454 2.2297 76 10.968 0.0000
## b0.calorSem 2.935 3.1259 76 0.939 0.3508
## tx.(Intercept) -0.674 0.1148 76 -5.874 0.0000
## tx.calorSem -0.796 0.2812 76 -2.831 0.0059
## xb 4.081 0.6804 76 5.999 0.0000
## Correlation:
## b0.(I) b0.clS tx.(I) tx.clS
## b0.calorSem -0.713
## tx.(Intercept) -0.257 0.184
## tx.calorSem 0.105 -0.075 -0.408
## xb 0.000 -0.084 0.000 -0.684
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.08144 -0.48480 0.01693 0.48725 4.30613
##
## Number of Observations: 100
## Number of Groups: 20
anova(n12, n11)
## Model df AIC BIC logLik Test L.Ratio p-value
## n12 1 7 534.6 552.8 -260.3
## n11 2 8 534.0 554.8 -259.0 1 vs 2 2.573 0.1087
## 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: -260.2 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.078
##
## PNLS step: RSS = 500.4
## fixed effects:23.02 4.692 -1.113 4.105 -1.355
## iterations: 4
##
## Convergence:
## fixed reStruct
## 1.000000 0.001164
##
## **Iteration 2
## LME step: Loglik: -260.3 , nlm iterations: 1
## reStruct parameters:
## ue
## -1.077
##
## PNLS step: RSS = 500.4
## fixed effects:23.02 4.692 -1.113 4.105 -1.355
## 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
## 534.5 552.8 -260.3
##
## Random effects:
## Formula: b0 ~ 1 | ue
## b0.(Intercept) Residual
## StdDev: 6.565 2.237
##
## Fixed effects: list(b0 ~ calor, tx ~ 1, xb ~ calor)
## Value Std.Error DF t-value p-value
## b0.(Intercept) 23.016 2.1909 76 10.505 0.0000
## b0.calorSem 4.692 3.1406 76 1.494 0.1393
## tx -1.113 0.1372 76 -8.114 0.0000
## xb.(Intercept) 4.105 0.6941 76 5.915 0.0000
## xb.calorSem -1.355 0.9423 76 -1.438 0.1544
## Correlation:
## b0.(I) b0.clS tx xb.(I)
## b0.calorSem -0.698
## tx 0.000 0.000
## xb.(Intercept) -0.156 0.109 -0.514
## xb.calorSem 0.115 -0.240 -0.046 -0.518
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.34836 -0.45810 0.04078 0.45606 4.16491
##
## 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 534.0 554.8 -259.0
## n13 2 7 534.5 552.8 -260.3 1 vs 2 2.523 0.1122
AIC(n11); AIC(m0)
##
## 534
## [1] 533.1
##-----------------------------------------------------------------------------
## 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 23.8 27.4 23.6 27.4 23.4 ...
## ..- 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
## 24.4543 2.9347 -0.6743 -0.7960 4.0812
formula(n12)
## y ~ b0 + tx * (dias - xb * i) * (dias > xb * i)
## <environment: 0xaa1a208>
## 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 23.8 27.4 23.6 27.4 23.4 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ yn : atomic 23.8 27.4 23.6 27.4 23.4 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ ym : atomic 23.2 27.6 23.2 27.6 23.1 ...
## ..- attr(*, "label")= chr "Predicted values"
## $ se : num 2.15 2.14 2.14 2.14 2.14 ...
## $ lwr : num 19.6 23.2 19.4 23.2 19.2 ...
## $ fit : num 23.8 27.4 23.6 27.4 23.4 ...
## $ upr : num 28 31.6 27.8 31.6 27.6 ...
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 35.5 42.9 40.8 40.1 38.7 ...
## $ a : num -9.79 -8.67 -9.83 -10.23 -7.67 ...
## $ b : num 13.3 13.7 12.9 14.1 15.2 ...
## $ E : num 21.7 25.8 24.4 25 24.4 ...
## $ peso : num 44.2 20 24.1 29.2 32.9 ...
## $ 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 44.2 20 24.1 29.2 32.9 ...
## $ 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 526.9 542.5 -257.4
## m2 2 8 530.7 551.6 -257.4 1 vs 2 0.1489 0.9282
## m0 3 12 537.6 568.8 -256.8 2 vs 3 1.1704 0.8829
m0 <- m1
## Sem interação
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 78 7301 <.0001
## calor 1 18 2 0.1915
## poly(dias, 1) 1 78 12 0.0011
## calor:poly(dias, 1) 1 78 0 0.8708
m0 <- lme(y~calor+dias, random=~1|ue, data=db, method="ML")
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 5 524.9 537.9 -257.5
## m2 2 8 530.7 551.6 -257.4 1 vs 2 0.1767 0.9813
## summary(m0)
## VarCorr(m0)
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 41.3771 0.8667 79 47.742 4.982e-60
## calorSem -1.3486 0.9886 18 -1.364 1.893e-01
## dias 0.3503 0.1025 79 3.419 9.967e-04
## 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 41.7 40.4 41.8 40.5 41.9 ...
## ..- 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 35.5 42.9 40.8 40.1 38.7 ...
## $ a : num -9.79 -8.67 -9.83 -10.23 -7.67 ...
## $ b : num 13.3 13.7 12.9 14.1 15.2 ...
## $ E : num 21.7 25.8 24.4 25 24.4 ...
## $ peso : num 44.2 20 24.1 29.2 32.9 ...
## $ 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 35.5 42.9 40.8 40.1 38.7 ...
## $ 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 459.4 475.0 -223.7
## m2 2 8 461.1 482.0 -222.6 1 vs 2 2.211 0.3310
## m0 3 12 467.1 498.4 -221.6 2 vs 3 2.006 0.7347
m0 <- m1
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 78 466.3 <.0001
## calor 1 18 6.5 0.0201
## poly(dias, 1) 1 78 19.3 <.0001
## calor:poly(dias, 1) 1 78 16.0 0.0001
m0 <- lme(y~calor*dias, random=~1|ue, data=db, method="ML")
anova(m0, m1)
## Model df AIC BIC logLik
## m0 1 6 459.4 475 -223.7
## m1 2 6 459.4 475 -223.7
## summary(m0)
## VarCorr(m0)
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) -9.71405 0.7548 78 -12.8705 5.540e-21
## calorSem -0.83600 1.0674 18 -0.7832 4.437e-01
## dias 0.02805 0.1010 78 0.2776 7.820e-01
## calorSem:dias 0.57180 0.1429 78 4.0017 1.421e-04
## 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 -9.69 -9.95 -9.68 -9.78 -9.67 ...
## ..- 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 35.5 42.9 40.8 40.1 38.7 ...
## $ a : num -9.79 -8.67 -9.83 -10.23 -7.67 ...
## $ b : num 13.3 13.7 12.9 14.1 15.2 ...
## $ E : num 21.7 25.8 24.4 25 24.4 ...
## $ peso : num 44.2 20 24.1 29.2 32.9 ...
## $ 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 -9.79 -8.67 -9.83 -10.23 -7.67 ...
## $ 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 494.6 510.2 -241.3
## m2 2 8 498.2 519.1 -241.1 1 vs 2 0.3694 0.8314
## m0 3 12 505.8 537.1 -240.9 2 vs 3 0.3877 0.9835
m0 <- m2
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 76 1071.2 <.0001
## calor 1 18 7.6 0.0128
## poly(dias, 2) 2 76 0.6 0.5294
## calor:poly(dias, 2) 2 76 0.3 0.7678
m0 <- lme(y~calor+dias, random=~1|ue, data=db, method="ML")
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 5 493.1 506.2 -241.6
## m2 2 8 498.2 519.1 -241.1 1 vs 2 0.9212 0.8203
## summary(m0)
## VarCorr(m0)
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 16.42518 0.7861 79 20.8942 6.426e-34
## calorSem -2.62160 0.9343 18 -2.8059 1.169e-02
## dias 0.08352 0.0852 79 0.9803 3.299e-01
## 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 16.5 13.9 16.5 13.9 16.6 ...
## ..- 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 35.5 42.9 40.8 40.1 38.7 ...
## $ a : num -9.79 -8.67 -9.83 -10.23 -7.67 ...
## $ b : num 13.3 13.7 12.9 14.1 15.2 ...
## $ E : num 21.7 25.8 24.4 25 24.4 ...
## $ peso : num 44.2 20 24.1 29.2 32.9 ...
## $ 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 13.3 13.7 12.9 14.1 15.2 ...
## $ 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 565.1 580.7 -276.5
## m2 2 8 569.0 589.9 -276.5 1 vs 2 0.0351 0.9826
## m0 3 12 575.8 607.1 -275.9 2 vs 3 1.2251 0.8739
m0 <- m1
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 78 2028.2 <.0001
## calor 1 18 6.1 0.0238
## poly(dias, 1) 1 78 3.9 0.0528
## calor:poly(dias, 1) 1 78 1.8 0.1880
m0 <- lme(y~calor+dias, random=~1|ue, data=db, method="ML")
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 5 564.9 577.9 -277.4
## m2 2 8 569.0 589.9 -276.5 1 vs 2 1.852 0.6037
## summary(m0)
## VarCorr(m0)
summary(m0)$tTable
## Value Std.Error DF t-value p-value
## (Intercept) 27.0441 1.0471 79 25.828 2.863e-40
## calorSem -2.9386 1.1842 18 -2.482 2.318e-02
## dias 0.2458 0.1257 79 1.955 5.416e-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 27.3 24.4 27.4 24.4 27.4 ...
## ..- 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.5
## 5 9 Com 0.8
## 6 1 Sem 0.0
## 7 3 Sem 0.0
## 8 5 Sem 0.3
## 9 7 Sem 0.7
## 10 9 Sem 0.9
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 5 5 5 9 3 3 9 7 3 ...
## $ dir : int 9 7 7 7 NA 5 5 NA 9 5 ...
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.6 7.2
## 2 Sem 5.2 7.0
##-----------------------------------------------------------------------------
## Conduzindo uma análise de sobreviência.
S <- with(dc, Surv(time=esq, time2=dir, event=status, type="interval"))
S
## [1] [7, 9] [5, 7] [5, 7] [5, 7] 9+ [3, 5] [3, 5] 9+ [7, 9] [3, 5] [5, 7] [7, 9]
## [13] [5, 7] [5, 7] 9+ [7, 9] [3, 5] [3, 5] [3, 5] [5, 7]
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 56.53 NA NA NA
## 2 1 18 56.93 = -1 -0.3962 0.529
summary(s0)
##
## Call:
## survreg(formula = S ~ calor, data = dc, dist = dist)
## Value Std. Error z p
## (Intercept) 2.0182 0.108 18.72 3.61e-78
## calorSem -0.0938 0.149 -0.63 5.28e-01
## Log(scale) -1.2066 0.208 -5.79 6.94e-09
##
## Scale= 0.299
##
## Weibull distribution
## Loglik(model)= -28.3 Loglik(intercept only)= -28.5
## Chisq= 0.4 on 1 degrees of freedom, p= 0.53
## Number of Newton-Raphson Iterations: 4
## 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
## 7.525 6.851
## 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.743 6.140
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 2.018 1.777 2.259
## 2 == 0 1.924 1.692 2.156
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))
