# ========================================================
# SURVIVAL ANALYSIS IN R
# The lung cancer data from the survival package
# ========================================================
# Question: how long did people with advanced lung cancer
# survive after their follow-up began, and did survival
# differ between men and women?
# ========================================================


# ========================================================
# SECTION 1: PACKAGE AND DATA
# ========================================================

# install.packages("survival")   # run once
library(survival)   # Surv(), survfit(), survdiff(), coxph()

# The lung data come with the survival package: 228 people
# with advanced lung cancer, from the North Central Cancer
# Treatment Group
dim(lung)
head(lung)


# ========================================================
# SECTION 2: PREPARE THE VARIABLES
# ========================================================
# time:   days from the start of follow-up to death, or to
#         the last contact for people still alive
# status: 1 = censored (alive at last contact), 2 = died
# sex:    1 = male, 2 = female
# age:    age in years
surv_data <- data.frame(
  time   = lung$time,
  status = lung$status,
  sex    = lung$sex,
  age    = lung$age)

# The event indicator: 1 = died, 0 = censored
surv_data$event <- NA
surv_data$event[surv_data$status == 2] <- 1
surv_data$event[surv_data$status == 1] <- 0
table(surv_data$event)

# Sex as a factor, with men as the reference group
surv_data$sex <- factor(surv_data$sex, levels = c(1, 2),
                        labels = c("Male", "Female"))
table(surv_data$sex)

# An ordinary summary treats every censored time as if it
# were a time of death
summary(surv_data$time)


# ========================================================
# SECTION 3: THE SURVIVAL OBJECT
# ========================================================
# Surv() joins each person's time and event status
surv_object <- Surv(surv_data$time, surv_data$event)
head(surv_object, 10)


# ========================================================
# SECTION 4: THE KAPLAN-MEIER ESTIMATE FOR EVERYONE
# ========================================================
# "~ 1" asks for one curve for the whole sample
km_all <- survfit(Surv(time, event) ~ 1, data = surv_data)
km_all

# Estimated survival at about 6 months, 1 year and 2 years
summary(km_all, times = c(180, 365, 730))

# The Kaplan-Meier curve: mark.time = TRUE draws a small
# tick at each censored time
plot(km_all, mark.time = TRUE,
     xlab = "Days since the start of follow-up",
     ylab = "Proportion surviving",
     main = "Kaplan-Meier estimate, all patients")
abline(h = 0.5, lty = 2)   # half of the group surviving


# ========================================================
# SECTION 5: KAPLAN-MEIER ESTIMATES BY SEX
# ========================================================
km_sex <- survfit(Surv(time, event) ~ sex, data = surv_data)
km_sex

plot(km_sex, mark.time = TRUE, lwd = 2,
     col = c("steelblue", "firebrick"),
     xlab = "Days since the start of follow-up",
     ylab = "Proportion surviving",
     main = "Kaplan-Meier estimates by sex")
legend("topright", legend = c("Male", "Female"),
       col = c("steelblue", "firebrick"), lwd = 2)


# ========================================================
# SECTION 6: THE LOG-RANK TEST
# ========================================================
# Compares the observed number of deaths in each group with
# the number expected if both groups had the same survival
survdiff(Surv(time, event) ~ sex, data = surv_data)


# ========================================================
# SECTION 7: COX PROPORTIONAL HAZARDS REGRESSION
# ========================================================
cox_model <- coxph(Surv(time, event) ~ sex + age,
                   data = surv_data)
summary(cox_model)


# ========================================================
# SECTION 8: CHECK THE PROPORTIONAL HAZARDS ASSUMPTION
# ========================================================
# A small p-value suggests that a hazard ratio changes
# over follow-up time
ph_check <- cox.zph(cox_model)
ph_check

# The estimated log hazard ratio for sex across follow-up,
# with the single value from the Cox model as a dashed line
plot(ph_check[1], resid = FALSE,
     main = "Proportional hazards check: sex")
abline(h = coef(cox_model)["sexFemale"], lty = 2,
       col = "firebrick")