# NOT RUN {
# we construct the composite event indicator (transplantation or death)
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")
# we fit the joint model using splines for the subject-specific
# longitudinal trajectories and a spline-approximated baseline
# risk function
lmeFit <- lme(log(serBilir) ~ ns(year, 2), data = pbc2,
random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModelBayes(lmeFit, survFit, timeVar = "year")
# we will compute survival probabilities for Subject 2 in a dynamic manner,
# i.e., after each longitudinal measurement is recorded
ND <- pbc2[pbc2$id == 2, ] # the data of Subject 2
survPreds <- vector("list", nrow(ND))
for (i in 1:nrow(ND)) {
survPreds[[i]] <- survfitJM(jointFit, newdata = ND[1:i, ])
}
# the default call to the plot method using the first three
# longitudinal measurements
plot(survPreds[[3]])
# we produce the corresponding plot
par(mfrow = c(2, 2), oma = c(0, 2, 0, 2))
for (i in c(1,3,5,7)) {
plot(survPreds[[i]], estimator = "median", conf.int = TRUE,
include.y = TRUE, main = paste("Follow-up time:",
round(survPreds[[i]]$last.time, 1)), ylab = "", ylab2 = "")
}
mtext("log serum bilirubin", side = 2, line = -1, outer = TRUE)
mtext("Survival Probability", side = 4, line = -1, outer = TRUE)
# }
Run the code above in your browser using DataLab