survival_analysis.R

# =============================================================================
# Survival Analysis: Kaplan-Meier and Cox Regression
#
# Author: Nikhil Kirtipal
#
# Demonstrates:
#   - Kaplan-Meier estimation and log-rank test
#   - Cox proportional hazards regression
#   - PH assumption testing (Schoenfeld residuals)
#   - Publication-ready summary tables
#
# Data: survival::lung (NCCTG Lung Cancer, 228 patients) - built into the package
#       Loprinzi et al. (1994), J Clin Oncol 12:601-607
#
# Tested: R 4.4.2, survival 3.6-6, survminer 0.5.2, gtsummary 2.5.1,
#         gt 1.3.0, dplyr 1.2.0, tidyr 1.3.2
#
# Replace the data preparation section to use your own clinical dataset.
# =============================================================================

library(survival)
library(survminer)
library(gtsummary)
library(gt)
library(dplyr)
library(tidyr)
#
suppressWarnings(data("lung", package = "survival"))

# ---- 0. Data preparation ----------------------------------------------------
# lung codes status as 1 = censored, 2 = dead. Survival functions expect
# 0 = censored, 1 = event, so subtract 1.
#
# ph.ecog has a single patient at level 3, which gives an unstable univariable
# model, so it is dropped. Complete cases are used throughout so that the
# unadjusted and adjusted models are fitted on the same rows.

dat <- lung %>%
  filter(!is.na(ph.ecog), ph.ecog < 3) %>%
  mutate(
    status  = status - 1,
    sex     = factor(sex, levels = 1:2, labels = c("Male", "Female")),
    ph.ecog = factor(ph.ecog, levels = 0:2,
                     labels = c("Asymptomatic", "Ambulatory", "In bed <50%"))
  ) %>%
  drop_na(time, status, sex, age, ph.ecog, wt.loss)

cat("n =", nrow(dat), " events =", sum(dat$status), "\n")

# Variable labels feed straight through to the gtsummary tables
dat <- dat %>%
  labelled::set_variable_labels(
    time     = "Follow-up time (days)",
    age      = "Age (years)",
    sex      = "Sex",
    ph.ecog  = "ECOG performance status",
    ph.karno = "Karnofsky score (physician)",
    wt.loss  = "Weight loss (lbs, 6 months)"
  )


# ---- 1. Kaplan-Meier --------------------------------------------------------

surv_obj <- Surv(time = dat$time, event = dat$status)
fit <- survfit(surv_obj ~ sex, data = dat)

print(fit)
print(surv_median(fit))

pal <- c("#1B6B70", "#E07A5F")   # teal, coral

km_plot <- ggsurvplot(
  fit,
  data           = dat,
  palette        = pal,
  pval           = TRUE,
  pval.size      = 4.5,
  pval.coord     = c(20, 0.08),
  conf.int       = TRUE,
  conf.int.alpha = 0.15,
  censor.shape   = "|",
  censor.size    = 3,
  size           = 1.1,
  xlab           = "Time (days)",
  ylab           = "Survival probability",
  legend.title   = "",
  legend.labs    = c("Male", "Female"),
  legend         = c(0.85, 0.85),
  risk.table          = TRUE,
  risk.table.height   = 0.24,
  risk.table.title    = "Number at risk",
  risk.table.y.text      = TRUE,
  risk.table.y.text.col  = TRUE,
  risk.table.fontsize = 3.8,
  tables.theme        = theme_cleantable(),
  break.time.by  = 250,
  xlim           = c(0, 1000),
  ggtheme = theme_minimal(base_size = 13) +
    theme(
      panel.grid.minor   = element_blank(),
      panel.grid.major.x = element_blank(),
      panel.grid.major.y = element_line(colour = "grey92", linewidth = 0.4),
      axis.line          = element_line(colour = "grey30", linewidth = 0.4),
      axis.ticks         = element_line(colour = "grey30", linewidth = 0.4),
      legend.text        = element_text(size = 11),
      plot.margin        = margin(10, 14, 6, 10)
    )
)

print(km_plot)

# ---- 2. Baseline characteristics --------------------------------------------

tab1 <- dat %>%
  select(age, sex, ph.ecog, wt.loss) %>%
  tbl_summary(
    by        = sex,
    missing   = "no",
    statistic = list(
      all_categorical() ~ "{n} ({p}%)",
      all_continuous()  ~ "{mean} ± {sd}"
    ),
    digits = list(
      all_categorical() ~ c(0, 1),
      all_continuous()  ~ c(1, 1)
    )
  ) %>%
  bold_labels() %>%
  add_overall(last = TRUE) %>%
  add_p(pvalue_fun = ~style_pvalue(.x, digits = 3))

print(tab1)


# ---- 3. Cox proportional hazards --------------------------------------------

cox <- coxph(Surv(time, status) ~ sex + age + ph.ecog + wt.loss, data = dat)
summary(cox)


# ---- 4. Check the proportional hazards assumption ---------------------------
# The Cox model assumes hazard ratios stay constant over time. This is the
# assumption the entire model rests on, and it is routinely skipped.
#
# cox.zph tests the correlation between scaled Schoenfeld residuals and time.
#   p < 0.05 means the assumption is violated for that term.
#
# If violated: stratify on the offending variable, add a time-varying
# coefficient, or move to an accelerated failure time (AFT) model.

zph <- cox.zph(cox)
print(zph)

par(mfrow = c(2, 2))
plot(zph)
par(mfrow = c(1, 1))

# Stratified alternative, if sex were to violate PH:
# cox_strat <- coxph(Surv(time, status) ~ strata(sex) + age + ph.ecog + wt.loss,
#                    data = dat)


# ---- 5. Unadjusted hazard ratios --------------------------------------------

hr_uni <- dat %>%
  select(time, status, sex, age, ph.ecog, wt.loss) %>%
  tbl_uvregression(
    method       = coxph,
    y            = Surv(time, status),
    exponentiate = TRUE,
    pvalue_fun   = ~style_pvalue(.x, digits = 3)
  ) %>%
  modify_column_merge(
    pattern = "{estimate} ({conf.low}, {conf.high})",
    rows    = !is.na(estimate)
  ) %>%
  modify_header(estimate ~ "**HR (95% CI)**") %>%
  bold_labels()

print(hr_uni)


# ---- 6. Adjusted hazard ratios ----------------------------------------------

hr_multi <- cox %>%
  tbl_regression(
    exponentiate = TRUE,
    pvalue_fun   = ~style_pvalue(.x, digits = 3)
  ) %>%
  modify_column_merge(
    pattern = "{estimate} ({conf.low}, {conf.high})",
    rows    = !is.na(estimate)
  ) %>%
  modify_header(estimate ~ "**HR (95% CI)**") %>%
  bold_labels()

print(hr_multi)


# ---- 7. Side-by-side table --------------------------------------------------

hr_table <- tbl_merge(
  tbls        = list(hr_uni, hr_multi),
  tab_spanner = c("**Unadjusted**", "**Adjusted**")
)

print(hr_table)


# ---- 8. Export --------------------------------------------------------------
# ggsave() does not work on a ggsurvplot object (it is a list, not a ggplot),
# so the plot is written through a png device instead.

dir.create("output", showWarnings = FALSE)

tab1_styled <- tab1 %>%
  as_gt() %>%
  tab_header(title = md("**Baseline characteristics**"),
             subtitle = md("NCCTG Lung Cancer cohort, *n* = 212")) %>%
  tab_options(
    table.font.size                   = px(13),
    column_labels.font.weight         = "bold",
    column_labels.border.top.width    = px(2),
    column_labels.border.top.color    = "#1B6B70",
    column_labels.border.bottom.width = px(1.5),
    column_labels.border.bottom.color = "#1B6B70",
    table_body.hlines.color           = "grey92",
    table.border.bottom.color         = "#1B6B70",
    table.border.bottom.width         = px(2),
    data_row.padding                  = px(5)
  ) %>%
  tab_style(style = cell_text(color = "#1B6B70", weight = "bold"),
            locations = cells_column_labels())

hr_styled <- hr_table %>%
  as_gt() %>%
  tab_header(title = md("**Hazard ratios for overall survival**"),
             subtitle = md("Cox proportional hazards, *n* = 212, 150 events")) %>%
  tab_options(
    table.font.size                   = px(13),
    column_labels.font.weight         = "bold",
    column_labels.border.top.width    = px(2),
    column_labels.border.top.color    = "#1B6B70",
    column_labels.border.bottom.width = px(1.5),
    column_labels.border.bottom.color = "#1B6B70",
    table_body.hlines.color           = "grey92",
    table.border.bottom.color         = "#1B6B70",
    table.border.bottom.width         = px(2),
    data_row.padding                  = px(5)
  ) %>%
  tab_style(style = cell_text(color = "#1B6B70", weight = "bold"),
            locations = cells_column_labels())

gtsave(tab1_styled, "output/baseline_table.docx")
gtsave(hr_styled,   "output/hazard_ratios.docx")

png("output/km_curve.png", width = 2400, height = 2100, res = 300)
print(km_plot)
dev.off()

# Save the PH test as a table 

zph_tab <- as.data.frame(zph$table) %>%
  tibble::rownames_to_column("Term") %>%
  mutate(across(c(chisq, p), ~round(.x, 3)))

write.csv(zph_tab, "output/ph_assumption_test.csv", row.names = FALSE)
print(zph_tab)

# Save the Schoenfeld residual plots

png("output/ph_diagnostics.png", width = 2400, height = 1800, res = 300)
par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1), bg = "white")

plot(zph[1], col = "#1B6B70", lwd = 2, main = "Sex")
abline(h = 0, col = "grey60", lty = 3)

plot(zph[2], col = "#1B6B70", lwd = 2, main = "Age")
abline(h = 0, col = "grey60", lty = 3)

plot(zph[3], col = "#1B6B70", lwd = 2, main = "ECOG performance status")
abline(h = 0, col = "grey60", lty = 3)

plot(zph[4], col = "#1B6B70", lwd = 2, main = "Weight loss")
abline(h = 0, col = "grey60", lty = 3)

par(mfrow = c(1, 1))

dev.off()


# =============================================================================
# Notes
#
# Events per variable: aim for at least 10 events per model term. This example
#   has 150 events and 5 terms, which is comfortable. Fitting 9 covariates to
#   50 events is not.
#
# Missing data: coxph drops rows with NA silently. Handled up front here so
#   every model uses the same sample. Otherwise the unadjusted and adjusted
#   estimates come from different n and are not comparable.
#
# Multiple testing: p-values here are unadjusted. Correct them if screening
#   many covariates.
# =============================================================================