# load package JM
library("JM")

###############
# Practical 1 #
###############

# the linear mixed model
lmeFit.p1 <- lme(log(serBilir) ~ year + drug:year, data = pbc2,
    random = ~ year | id)

# create the indicator for the composite event
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")

# the Cox model
survFit.p1 <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)

# the joint model
jointFit.p1 <- jointModel(lmeFit.p1, survFit.p1, timeVar = "year",
    method = "piecewise-PH-aGH")
    
summary(jointFit.p1)

confint(jointFit.p1, parm = "Longitudinal")
exp(confint(jointFit.p1, parm = "Event"))


# interaction term between serBilir and drug
intFact <- list(value = ~ drug, data = pbc2.id)
jointFit2.p1 <- update(jointFit.p1, interFact = intFact)

summary(jointFit2.p1)

# the joint model under the null H_0: beta_2 = 0
lmeFit2.p1 <- lme(log(serBilir) ~ year, data = pbc2,
    random = ~ year | id)
jointFit3.p1 <- update(jointFit2.p1, lmeObject = lmeFit2.p1)

anova(jointFit3.p1, jointFit2.p1)

# the joint model under the null H_0: gamma = alpha_2 = 0
survFit2.p1 <- coxph(Surv(years, status2) ~ 1, data = pbc2.id, x = TRUE)
jointFit4.p1 <- update(jointFit2.p1, survObject = survFit2.p1, interFact = NULL)

anova(jointFit4.p1, jointFit2.p1)

# the joint model under the null H_0: beta_2 = gamma = alpha_2 = 0
jointFit5.p1 <- update(jointFit2.p1, lmeObject = lmeFit2.p1, 
    survObject = survFit2.p1, interFact = NULL)

anova(jointFit5.p1, jointFit2.p1)


###############
# Practical 2 #
###############

# the linear mixed model
lmeFit.p4 <- lme(pro ~ time + time:treat, data = prothro,
    random = ~ time | id)

# the Cox model
survFit.p4 <- coxph(Surv(Time, death) ~ treat, data = prothros, x = TRUE)

# the joint model
jointFit.p4 <- jointModel(lmeFit.p4, survFit.p4, timeVar = "time",
    method = "piecewise-PH-aGH")

summary(jointFit.p4)


# the data of Patient 155
dataP155 <- prothro[prothro$id == 155, ]

# survival probabilities using only the 1st measurement
sfit <- survfitJM(jointFit.p4, newdata = dataP155[1, ])

sfit
plot(sfit)
plot(sfit, include.y = TRUE)

# we will use a for-loop to update predictions
pdf()
for (i in 1:nrow(dataP155)) {
    data.i <- dataP155[1:i, ]
    sfit.i <- survfitJM(jointFit.p4, newdata = data.i)
    plot(sfit.i, estimator = "mean", include.y = TRUE,
        conf.int = TRUE, fill.area = TRUE, col.area = "lightgrey", 
        ylab2 = "Prothrobin")
}
dev.off()


# predictions of future longitudinal responses, using only the 1st measurement
lfit <- predict(jointFit.p4, newdata = dataP155[1, ], type = "Subject",
    interval = "conf", returnData = TRUE)
    
library(lattice)
xyplot(pred + low + upp ~ time, data = lfit, type = "l",
    lty = c(1,2,2), col = 1, lwd = 2)

# we will use a for-loop to update predictions
pdf()
for (i in 1:nrow(dataP155)) {
    data.i <- dataP155[1:i, ]
    lfit.i <- predict(jointFit.p4, newdata = data.i, type = "Subject",
        interval = "conf", returnData = TRUE)
    print(xyplot(pred + low + upp ~ time, data = lfit.i, type = "l",
        lty = c(1,2,2), col = 1, lwd = 2), 
        xlab = "Time", ylab = "Predicted Prothrobin Level")
}
dev.off()

# AUC at 2 years with 0.5 year window
aucJM(jointFit.p4, newdata = prothro, Tstart = 2, Dt = 0.5)

# prediction error at 2 years with 0.5 year window
prederrJM(jointFit.p4, newdata = prothro, Tstart = 2, Thoriz = 2.5)