Survival Analysis

Survival Analysis Patterns

Overview

Comprehensive survival analysis methods in R covering Kaplan-Meier estimation, Cox proportional hazards models, parametric survival models, and advanced techniques for time-to-event data.

Basic Survival Objects

Creating Survival Data

library(survival)

# Right-censored data (most common)
surv_obj <- Surv(time = df$time, event = df$status)

# Left-truncated (delayed entry)
surv_obj <- Surv(time = df$entry_time, time2 = df$event_time, event = df$status)

# Interval censoring
surv_obj <- Surv(time = df$left, time2 = df$right, type = "interval2")

# Check structure
head(surv_obj)

Kaplan-Meier Estimation

Basic KM Analysis

# Fit Kaplan-Meier
km_fit <- survfit(Surv(time, status) ~ 1, data = df)

# Summary statistics
summary(km_fit)

# Median survival
km_fit

# Survival at specific times
summary(km_fit, times = c(12, 24, 36, 48, 60))

KM by Groups

# Stratified KM
km_fit <- survfit(Surv(time, status) ~ treatment, data = df)

# Log-rank test
survdiff(Surv(time, status) ~ treatment, data = df)

# Stratified log-rank
survdiff(Surv(time, status) ~ treatment + strata(site), data = df)

# Pairwise comparisons
pairwise_survdiff(Surv(time, status) ~ treatment, data = df)

Publication-Quality KM Plots (survminer)

library(survminer)

# Basic KM plot
ggsurvplot(km_fit, data = df)

# Full featured plot
ggsurvplot(
  km_fit,
  data = df,
  pval = TRUE,               # Add p-value
  conf.int = TRUE,           # Confidence intervals
  risk.table = TRUE,         # Risk table
  risk.table.col = "strata", # Color by group
  surv.median.line = "hv",   # Median line
  ggtheme = theme_bw(),
  palette = c("#E7B800", "#2E9FDF"),
  xlab = "Time (months)",
  ylab = "Survival probability",
  legend.title = "Treatment",
  legend.labs = c("Control", "Treatment")
)

Risk Tables and Annotations

ggsurvplot(
  km_fit,
  data = df,
  risk.table = TRUE,
  risk.table.height = 0.25,
  risk.table.y.text = FALSE,    # Show bars instead of names
  cumevents = TRUE,             # Cumulative events table
  cumcensor = TRUE,             # Cumulative censored table
  tables.height = 0.2,
  break.time.by = 12,           # Time axis breaks
  xlim = c(0, 60),
  surv.scale = "percent"        # Y-axis as percentage
)

Cox Proportional Hazards

Basic Cox Model

# Fit Cox model
cox_fit <- coxph(Surv(time, status) ~ treatment + age + sex, data = df)

# Model summary
summary(cox_fit)

# Hazard ratios with CI
exp(coef(cox_fit))
exp(confint(cox_fit))

# Tidy output
library(broom)
tidy(cox_fit, exponentiate = TRUE, conf.int = TRUE)

Proportional Hazards Assumption

# Test PH assumption
ph_test <- cox.zph(cox_fit)
print(ph_test)

# Plot Schoenfeld residuals
ggcoxzph(ph_test)

# Test for individual variables
ph_test$table

If proportional hazards is not credible, do not rely only on a single Cox hazard ratio or log-rank p-value. Consider prespecified alternatives such as time-varying treatment effects, weighted log-rank tests, milestone survival, or RMST.

Stratified Cox Model

# Stratify by variable (different baseline hazard)
cox_strat <- coxph(
  Surv(time, status) ~ treatment + age + strata(center),
  data = df
)

Time-Varying Covariates

# Start-stop format
cox_tv <- coxph(
  Surv(tstart, tstop, status) ~ treatment + biomarker,
  data = df_long
)

# With time-varying coefficient
cox_tv <- coxph(
  Surv(time, status) ~ treatment + tt(age),
  data = df,
  tt = function(x, t, ...) x * log(t)
)

Forest Plots

library(survminer)

# Forest plot from Cox model
ggforest(cox_fit, data = df)

# Custom forest plot
library(forestplot)
forest_data |>
  forestplot(
    mean = hr,
    lower = hr_lower,
    upper = hr_upper,
    labeltext = c(variable, n, hr_text),
    is.summary = is_summary
  )

Parametric Survival Models

With survreg

# Weibull model
weibull_fit <- survreg(
  Surv(time, status) ~ treatment + age,
  data = df,
  dist = "weibull"
)

# Other distributions
exponential_fit <- survreg(..., dist = "exponential")
lognormal_fit <- survreg(..., dist = "lognormal")
loglogistic_fit <- survreg(..., dist = "loglogistic")

# AFT interpretation
summary(weibull_fit)

With flexsurv (more distributions)

library(flexsurv)

# Generalized gamma (very flexible)
ggamma_fit <- flexsurvreg(
  Surv(time, status) ~ treatment + age,
  data = df,
  dist = "gengamma"
)

# Compare distributions
compare_distributions <- list(
  exponential = flexsurvreg(Surv(time, status) ~ treatment, data = df, dist = "exp"),
  weibull = flexsurvreg(Surv(time, status) ~ treatment, data = df, dist = "weibull"),
  lognormal = flexsurvreg(Surv(time, status) ~ treatment, data = df, dist = "lognormal"),
  gamma = flexsurvreg(Surv(time, status) ~ treatment, data = df, dist = "gamma"),
  gengamma = flexsurvreg(Surv(time, status) ~ treatment, data = df, dist = "gengamma")
)

# AIC comparison
sapply(compare_distributions, AIC)

Flexible Parametric Models (rstpm2)

library(rstpm2)

# Royston-Parmar spline model
stpm_fit <- stpm2(
  Surv(time, status) ~ treatment + age,
  data = df,
  df = 4  # degrees of freedom for baseline
)

# Time-varying effects
stpm_tv <- stpm2(
  Surv(time, status) ~ treatment + age,
  data = df,
  df = 4,
  tvc = list(treatment = 2)  # time-varying treatment effect
)

Competing Risks

Cause-Specific Hazards

# Fit separate Cox models for each cause
cs_cause1 <- coxph(
  Surv(time, status == 1) ~ treatment + age,
  data = df
)

cs_cause2 <- coxph(
  Surv(time, status == 2) ~ treatment + age,
  data = df
)

Fine-Gray Subdistribution Hazards

library(cmprsk)

# Cumulative incidence function
cif <- cuminc(
  ftime = df$time,
  fstatus = df$status,
  group = df$treatment
)

# Plot CIF
ggcompetingrisks(cif)

# Fine-Gray regression
fg_fit <- crr(
  ftime = df$time,
  fstatus = df$status,
  cov1 = model.matrix(~ treatment + age, df)[, -1],
  failcode = 1  # event of interest
)

summary(fg_fit)

With tidycmprsk

library(tidycmprsk)

# Cumulative incidence
cuminc(Surv(time, status) ~ treatment, data = df) |>
  ggcuminc()

# Fine-Gray model (tidy interface)
crr(Surv(time, status) ~ treatment + age, data = df, failcode = 1)

Restricted Mean Survival Time

library(survRM2)

# RMST comparison
rmst_result <- rmst2(
  time = df$time,
  status = df$status,
  arm = df$treatment,
  tau = 60  # restriction time
)

print(rmst_result)

# RMST regression
library(survRM2)
rmst_reg <- rmst2(
  time = df$time,
  status = df$status,
  arm = df$treatment,
  covariates = df[, c("age", "sex")],
  tau = 60
)

Multi-State Models

library(mstate)

# Define transition matrix
tmat <- transMat(
  x = list(
    c(2, 3),  # from state 1
    c(3),     # from state 2
    c()       # state 3 absorbing
  ),
  names = c("Healthy", "Illness", "Death")
)

# Prepare data
msdata <- msprep(
  time = c(NA, "illness_time", "death_time"),
  status = c(NA, "illness_status", "death_status"),
  data = df,
  trans = tmat
)

# Fit Cox model
ms_cox <- coxph(
  Surv(Tstart, Tstop, status) ~ treatment + strata(trans),
  data = msdata
)

# Transition probabilities
pt <- probtrans(msfit(ms_cox, newdata = new_patient), predt = 0)

Survival with tidymodels (censored)

library(censored)

# Model specification
cox_spec <- proportional_hazards() |>
  set_engine("survival") |>
  set_mode("censored regression")

# Workflow
surv_wf <- workflow() |>
  add_formula(Surv(time, status) ~ treatment + age + sex) |>
  add_model(cox_spec)

# Fit
surv_fit <- fit(surv_wf, data = train_data)

# Predict survival probability
predict(surv_fit, new_data, type = "survival", eval_time = c(12, 24, 36))

# Predict hazard
predict(surv_fit, new_data, type = "hazard", eval_time = c(12, 24, 36))

Key Packages Summary