# NOT RUN {
# linear mixed model fit
fitLME <- lme(log(serBilir) ~ drug * year, data = pbc2,
random = ~ year | id)
# survival regression fit
fitSURV <- coxph(Surv(years, status2) ~ drug, data = pbc2.id,
x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModelBayes(fitLME, fitSURV, timeVar = "year")
DF <- with(pbc2, expand.grid(drug = levels(drug),
year = seq(min(year), max(year), len = 100)))
Ps <- predict(fitJOINT, DF, interval = "confidence", return = TRUE)
require(lattice)
xyplot(pred + low + upp ~ year | drug, data = Ps,
type = "l", col = c(2,1,1), lty = c(1,2,2), lwd = 2,
ylab = "Average log serum Bilirubin")
# Subject-specific predictions
ND <- pbc2[pbc2$id == 2, ]
Ps.ss <- predict(fitJOINT, ND, type = "Subject",
interval = "confidence", return = TRUE)
xyplot(pred + low + upp ~ year | id, data = Ps.ss,
type = "l", col = c(2,1,1), lty = c(1,2,2), lwd = 2,
ylab = "Average log serum Bilirubin")
# }
Run the code above in your browser using DataLab