R - Survival Analysis

Modern Survival Analysis in R

Survival analysis is a statistical method for analyzing the expected duration of time until a specific event occurs. It’s often called time-to-event analysis. For example, it can be used to predict the time to patient recovery after a treatment, customer churn in a subscription service, or failure time for a mechanical part.

The foundational package for this in R is survival. For modern, publication-quality visualizations, we’ll use survminer, which extends ggplot2. We’ll also use dplyr from the tidyverse for data manipulation, which is a modern best practice.

Setup: Installing and Loading Packages

First, ensure you have a modern R development environment like RStudio. Then, install the necessary packages. You only need to do this once.

## Make sure you have the tidyverse for data handling and plotting
# install.packages("tidyverse")

## Install the core survival analysis packages
# install.packages("survival")
# install.packages("survminer")

Core Concepts: The Kaplan-Meier Estimator

The most common starting point is the Kaplan-Meier (KM) curve. It visualizes the probability of an individual surviving past a certain time point. To create it, we need two key functions:

• Surv(time, event): From the survival package. This function creates a special survival object. time is the follow-up duration, and event is a binary indicator (e.g., 1 for event occurred, 0 for censored). Censoring is crucial: it means the event of interest has not been observed for that subject by the end of the study (e.g., the patient was still alive).
• survfit(formula, data): This function fits the Kaplan-Meier model. The formula typically looks like Surv(time, event) ~ group, where group is a categorical variable to compare survival curves (e.g., treatment vs. placebo).

Example: Analyzing the ‘pbc’ Dataset

We will use the built-in pbc dataset from the survival package, which contains data on patients with Primary Biliary Cholangitis. Our goal is to estimate the survival probability over time.

The key columns are time (days to death, transplant, or study end) and status (0=censored, 1=transplant, 2=dead). For our analysis, the ‘event’ is death (status == 2).

# Load the modern libraries we'll use
library(tidyverse)
library(survival)
library(survminer)

# Inspect the data using modern tools
data(pbc)
# Use glimpse() for a compact summary
glimpse(pbc)

# The 'status' column needs cleaning. It has values 0, 1, 2.
# Let's define the event as death (status == 2).
# We create a new column 'event' which is TRUE if status is 2, FALSE otherwise.
pbc <- pbc %>%
  mutate(event = (status == 2))

# Let's check the first few rows of our modified data
head(pbc)

Executing the code above prepares our data. The glimpse() output would look something like this:

Rows: 418
Columns: 21
$ id      <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ...
$ time    <int> 400, 4500, 1012, 1925, 1504, 2503, 1832, 2466, 2400, 51, ...
$ status  <int> 2, 0, 2, 2, 1, 2, 0, 2, 0, 2, ...
$ trt     <int> 1, 1, 1, 1, 2, 2, 2, 2, 1, 2, ...
$ age     <dbl> 58.77, 56.45, 70.07, 54.74, 38.11, 66.26, 55.53, 53.06, ...
$ sex     <fct> f, f, m, f, f, f, f, f, f, f, ...
# ... other columns
$ event   <lgl> TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, ...

Creating and Visualizing the Kaplan-Meier Curve

Now we fit the model and use the modern ggsurvplot() for a clean and informative visualization.

# 1. Fit the Kaplan-Meier model
# The formula `Surv(time, event) ~ 1` means we are calculating the overall survival for all patients combined.
fit <- survfit(Surv(time, event) ~ 1, data = pbc)

# Print the model summary
print(fit)

# 2. Visualize the result with ggsurvplot
# This provides a much richer plot than the base R plot() function.
ggsurvplot(
  fit,
  data = pbc,
  risk.table = TRUE,          # Add a table showing number at risk
  conf.int = TRUE,            # Show confidence intervals
  pval = TRUE,                # Show p-value of the log-rank test
  ggtheme = theme_minimal(),    # Use a clean ggplot2 theme
  legend.title = "Overall Survival",
  surv.median.line = "hv"     # Add a line for the median survival time
)

The print(fit) command provides a summary:

Call: survfit(formula = Surv(time, event) ~ 1, data = pbc)

      n  events  median 0.95LCL 0.95UCL
    418     161    3395    3090    3853

This tells us we have 418 subjects, 161 events (deaths), and the median survival time is 3395 days. The ggsurvplot will generate a high-quality chart showing the survival probability decreasing over time, complete with confidence bands and a risk table below the plot.

Next Steps & Best Practices

• Comparing Groups: You can compare survival between groups, e.g., Surv(time, event) ~ sex. ggsurvplot will automatically create separate curves for males and females and perform a log-rank test to see if the difference is statistically significant.
• Cox Proportional Hazards Model: To understand how continuous predictors (like age or bilirubin levels) affect survival, the next step is to use a Cox Proportional Hazards model (coxph() function). This is a regression technique for survival data.
• Assumption Checking: For Cox models, it’s crucial to check the proportional hazards assumption. The cox.zph() function helps with this.
• Further Learning: Explore the vignettes for the survival and survminer packages for more advanced examples and detailed explanations.