drizopoulos

library("JMbayes")
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")

lmeFit <- lme(log(serBilir) ~ ns(year, 2) + age + sex, data = pbc2,
              random = ~ ns(year, 2) | id)

coxFit <- coxph(Surv(years, status2) ~ sex * age, 
                data = pbc2.id, x = TRUE)