Survival analysis consists of statistical methods that help us understand and predict how long it takes for an event to occur. These methods help researchers analyze time-to-event data, revealing not just if something happens, but when and how events unfold over time. In this article, we’ll walk through how to perform survival analysis in R using a lung cancer dataset.

Setting Up Your Environment

Before we begin, you’ll need to install and load the necessary R packages:

# Install required packages if you haven't already
install.packages("survival")
install.packages("ggplot2")
install.packages("survminer")  # Will automatically install ggpubr

# Load required libraries
library(survival)    # Core survival analysis functions
library(ggplot2)    # For data visualization
library(survminer)  # For enhanced survival plots

The survival package provides the core functions for survival analysis, while ggplot2 and survminer offer comprehensive visualization capabilities. The survminer package, which automatically loads ggpubr, provides specialized functions for creating publication-ready survival plots with minimal code.

About the Dataset

For this tutorial, we’ll use the lung dataset that comes built into the survival package. This dataset is automatically loaded when we attach the survival package. You might see a warning message when explicitly loading the data with data(“lung”) – this is normal and simply indicates that the dataset is already available.

Let’s verify our data is loaded:

# Check if the dataset exists
exists("lung")  # Should return TRUE

# Look at the first few rows
head(lung)

Image By Author

The dataset contains information about lung cancer patients, including their survival time, status (alive/dead), and various clinical characteristics. Here’s what the key variables represent:

• time: Survival time in days
• status: Patient status (1=censored/alive, 2=dead)
• ph.ecog: ECOG performance score (0=fully active, 1=restricted but ambulatory, 2=ambulatory but unable to work, 3-4=increasing disability)
• sex: Patient’s sex (1=male, 2=female)
• age: Patient’s age in years
• ph.karno: Karnofsky performance score rated by physician (0-100)
• pat.karno: Karnofsky performance score rated by patient (0-100)
• meal.cal: Calories consumed at meals
• wt.loss: Weight loss in last 6 months

This rich dataset will allow us to demonstrate key survival analysis techniques while working with real-world medical data.

Creating Basic Survival Curves

Now, let’s create our first visualization – the Kaplan-Meier survival curve:

# Create survival object
surv_obj <- Surv(time = lung$time, event = lung$status == 2)

# Fit Kaplan-Meier model
km_fit <- survfit(surv_obj ~ 1, data = lung)

# Create survival curve with all titles centered
ggsurvplot(km_fit,
           title = "Overall Survival Curve",
           xlab = "Time (days)",
           ylab = "Survival Probability",
           conf.int = TRUE,
           risk.table = TRUE,
           font.title = c(16, "bold"),
           font.x = c(14),
           font.y = c(14),
           font.tickslab = c(12),
           palette = "#2E9FDF",
           ggtheme = theme_minimal() + 
                     theme(plot.title = element_text(hjust = 0.5),
                           axis.title.x.bottom = element_text(hjust = 0.5)),  
           legend = "none",
           
           # Adjust spacing and proportions
           surv.plot.height = 0.7,
           risk.table.height = 0.3,
           risk.table.y.text.col = TRUE,
           risk.table.y.text = FALSE,
           tables.theme = theme_minimal() +
                         theme(plot.title = element_text(hjust = 0.5)),  
           break.time.by = 250,
           risk.table.pos = "out")

Image By Author

Analyzing Cumulative Risk Over Time

While survival curves show us the probability of surviving past a certain time point, examining cumulative hazard provides a different perspective on risk patterns. Let’s create a Nelson-Aalen estimate of the cumulative hazard, stratified by sex:

# Create cumulative hazard plot by sex
# Fit Nelson-Aalen models for different sex groups
na_fit <- survfit(surv_obj ~ sex, data = lung, type = "fleming-harrington")

# Create enhanced cumulative hazard plot
ggsurvplot(na_fit,
           fun = "cumhaz",
           title = "Cumulative Hazard by Sex",
           xlab = "Time (days)",
           ylab = "Cumulative Hazard",
           conf.int = TRUE,
           palette = c("#E7B800", "#2E9FDF"),
           ggtheme = theme_minimal() + 
                     theme(plot.title = element_text(hjust = 0.5),
                           axis.title.x = element_text(hjust = 0.5)),
           font.title = c(16, "bold"),
           font.x = c(14),
           font.y = c(14),
           font.tickslab = c(12),
           break.time.by = 250,
           legend.title = "Sex",
           legend.labs = c("Male", "Female"))

Image By Author

This visualization reveals several important patterns in our lung cancer data:

1. Males (shown in yellow) consistently show a higher cumulative hazard compared to females (in blue), suggesting gender may be an important prognostic factor.

2. The increasing slope of both curves indicates that risk accumulates more rapidly during certain time periods, particularly after 500 days.

3. The widening confidence intervals (shaded areas) towards the later time points reflect increased uncertainty in our estimates, typically due to fewer patients remaining in the study.

The survfit() function with type = “fleming-harrington” provides a robust estimate of cumulative hazard, while ggsurvplot() from the survminer package helps us create publication-quality visualizations with minimal code.

Quantifying Risk Factors

To better understand which factors significantly influence survival time, we’ll use a Cox Proportional Hazards model. This powerful statistical technique allows us to quantify the impact of multiple variables simultaneously while accounting for censored observations:

# Prepare data for Cox model
# Scale continuous variables for better visualization
lung_prep <- lung %>%
  mutate(age_scaled = scale(age),
         wt_scaled = scale(wt.loss))

# Fit Cox model
cox_fit <- coxph(Surv(time, status == 2) ~ age_scaled + sex + ph.ecog + wt_scaled, 
                 data = lung_prep)

# Extract hazard ratios
hazard_ratios <- exp(coef(cox_fit))

# Create data frame for plotting
hr_df <- data.frame( 
  factor = names(hazard_ratios), 
  hr = hazard_ratios 
) 

# Create the plot using ggplot2
ggplot(hr_df, aes(x = hr, y = factor)) +
  geom_bar(stat = "identity", 
          aes(fill = hr > 1),
          width = 0.7) +
  scale_fill_manual(values = c("green3", "red2")) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "black", alpha = 0.3) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    legend.position = "none",
    # Add some padding to the right margin
    plot.margin = margin(5, 30, 5, 5, "pt")
  ) +
  labs(
    title = "Risk Factors in Lung Cancer Survival",
    x = "Hazard Ratio",
    y = ""
  ) +
  # Adjust text positioning
  geom_text(aes(label = sprintf("%.2f", hr)), 
            hjust = ifelse(hr_df$hr > 1, -0.1, 1.2),
            size = 4) +
  # Extend x-axis limits to accommodate labels
  scale_x_continuous(limits = c(0, max(hr_df$hr) * 1.2))

Image By Author

The hazard ratios plot above reveals several key insights:

1. ECOG Performance Score (ph.ecog) has the strongest effect, with a hazard ratio of 1.67, indicating that higher scores are associated with a 67% increase in mortality risk.

2. Age (age_scaled) shows a moderate effect, with a hazard ratio of 1.13, suggesting slightly higher risk with increasing age.

3. Both sex and weight (wt_scaled) appear protective (hazard ratios < 1), with female gender associated with a 45% reduction in risk (hazard ratio = 0.55) and higher weight showing a modest protective effect (hazard ratio = 0.89).

The dashed vertical line at 1.0 serves as a reference – factors with bars extending to the right of this line increase risk, while those to the left decrease risk. Green bars indicate protective factors (hazard ratio < 1) while red bars indicate risk factors (hazard ratio > 1).

Conclusion: From Analysis to Clinical Insights

Throughout this tutorial, we’ve demonstrated how R provides robust tools for survival analysis, allowing us to extract meaningful insights from time-to-event data. Using the lung cancer dataset, we’ve showcased three key aspects of survival analysis:

1. Basic survival curves that visualize overall patient outcomes
2. Cumulative hazard analysis that reveals how risk patterns differ between groups
3. Cox Proportional Hazards modeling that quantifies the impact of various risk factors

The survival package, combined with survminer‘s visualization capabilities, makes it straightforward to perform sophisticated analyses with relatively simple code. To extend this analysis further, you might consider:

• Including interaction terms in your Cox model
• Performing time-dependent covariate analysis
• Validating the proportional hazards assumption
• Conducting stratified analyses for specific patient subgroups

Remember that while R provides powerful tools for survival analysis, the interpretation of results should always be guided by domain expertise and careful consideration of the underlying data structure.