# Load packages JM
library("JM")

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


#################
# Section 3.2.1 #
#################

S.KM <- survfit(Surv(years, status2) ~ 1, data = pbc2.id)
H.NA <- survfit(Surv(years, status2) ~ 1, data = pbc2.id,
    type = "fleming-harrington")

par(mfrow = c(1, 2))
plot(S.KM, xlab = "Time (years)", ylab = "Survival Probability",
    main = "Kaplan-Meier Estimator", mark.time = FALSE)
plot(H.NA, xlab = "Time (years)", ylab = "Cumulative Hazard",
    main = "Nelson-Aalen Estimator", fun = "cumhaz", mark.time = FALSE)


#################
# Section 3.3.1 #
#################

coxFit <- coxph(Surv(years, status2) ~ drug + age + sex,
    data = pbc2.id)

summary(coxFit)


###############
# Section 3.5 #
###############

tdCox.pro <- coxph(Surv(start, stop, event) ~ pro + treat,
    data = prothro)

summary(tdCox.pro)