Introduction to Survival Analysis in R

Hello

What are survival data anyway?

Examples from cancer

• Time from diagnosis to death
• Time from surgery to recurrence of disease
• Time from start of treatment to progression of disease
• Time from response to recurrence of disease

Examples from other fields

• Time from HIV infection to development of AIDS
• Time to from diagnosis with heart disease to heart attack
• Time from dicharge from rehabilitation facility to recurrence of substance abuse
• Time from birth to initiation of sexual activity
• Time from production to machine malfunction

A rose by any other name…

Because time-to-event data are common in many fields, it also goes by names besides survival analysis including:

What is censoring?

But isn’t this just binary data??

• Binary data doesn’t have the ability to change depending on the time of analysis, e.g. 5-year survival will have the same value whether it is analyzed at 5 years and 1 day, 5 years and 2 days, 6 years, etc. Either a participant died by 5 years or they didn’t.

• Time-to-event data may have different values depending on the time of analysis, e.g. overall survival will have different values depending on whether it is analyzed at 5 years and 1 day or at 6 years, since additional participants can die between those two time points.

Reasons for censoring

A subject may be censored due to:

• Loss to follow-up
• Withdrawal from study
• No event by end of fixed study period

Other types of censoring

Censoring must be considered in the analysis

Data components for survival analysis

To analyze survival data, we need to know the observed time (\(Y_i\)) and the event indicator (\(\delta_i\)). For a subject (denoted by \(i\)):

• Observed time is the minimum of the event time (\(T_i\)) and censoring time (\(C_i\)) (\(Y_i = \min(T_i, C_i)\))
• Event indicator (\(\delta_i\)) is 1 if the event is observed (i.e. \(T_i \leq C_i\)) and 0 if censored (i.e. \(T_i > C_i\))

Load needed packages

library(dplyr)
library(ggplot2)
library(lubridate)
library(survival)
library(ggsurvfit)
library(gtsummary)

Load example data

To access the example data used throughout this talk, install and load the cancersimdata package from my GitHub repo:

# If needed, install the remotes package first
install.packages("remotes")

# Then install the GitHub repository for the dataset
remotes::install_github("zabore/cancersimdata")

# Finally, load the repository
library(cancersimdata)

Example data background

bc_rt_data is a synthetic dataset based on real breast cancer data. The dataset contains information on 3000 women with T1-2N1M0 breast cancer, who had a mastectomy between 1995-2015.

The original study examined the association between post-mastectomy radiation therapy and disease recurrence.

Example data contents

Relevant variables include:

• rt: PMRT indicator, 1 = yes, 0 = no
• os_event: Death indicator, 1 = dead, 0 = censored
• date_of_mastecomy: Date of mastectomy
• date_last_follow_up_death: Date of last follow-up or death

Example data format

bc_rt_data |> 
  select(rt, os_event, date_of_mastectomy, date_last_follow_up_death) |> 
  print(n = 10)

# A tibble: 3,000 × 4
      rt os_event date_of_mastectomy date_last_follow_up_death
   <dbl>    <dbl> <chr>              <chr>                    
 1     0        0 09/11/2004         07/24/2009               
 2     0        0 08/06/2011         11/07/2016               
 3     0        0 06/10/1998         11/04/2005               
 4     1        0 10/06/2013         04/30/2019               
 5     1        0 11/22/2002         01/17/2008               
 6     0        0 01/31/2013         01/12/2019               
 7     1        0 05/03/2014         08/23/2016               
 8     0        0 02/16/1998         11/12/2004               
 9     1        0 03/04/2014         02/07/2020               
10     1        1 03/06/1998         09/08/2003               
# ℹ 2,990 more rows

Event indicator

It is important to pay attention to the format of the event indicator.

The Surv() function in the survival package accepts by default TRUE/FALSE, where TRUE is event and FALSE is censored; 1/0 where 1 is event and 0 is censored; or 2/1 where 2 is event and 1 is censored. Please take care to ensure the event indicator is properly formatted.

Here we see that the documentation stipulates 1 is event (death) and 0 is censored.

Start and end dates

Here, the start and end dates are in the dataset, formatted as character variables.

We need to:

1. Convert them to date format
2. Calculate the follow-up times

Formatting dates

The lubridate package offers a more comprehensive and user-friendly suite of functions for date manipulation. See the lubridate website for details.

The mdy() function converts character values ordered month day year:

bc_rt_data <- 
  bc_rt_data |> 
  mutate(
    date_of_mastectomy = mdy(date_of_mastectomy), 
    date_last_follow_up_death = mdy(date_last_follow_up_death)
  )

Formatted dates

bc_rt_data |> 
  select(rt, os_event, date_of_mastectomy, date_last_follow_up_death) |> 
  print(n = 10)

# A tibble: 3,000 × 4
      rt os_event date_of_mastectomy date_last_follow_up_death
   <dbl>    <dbl> <date>             <date>                   
 1     0        0 2004-09-11         2009-07-24               
 2     0        0 2011-08-06         2016-11-07               
 3     0        0 1998-06-10         2005-11-04               
 4     1        0 2013-10-06         2019-04-30               
 5     1        0 2002-11-22         2008-01-17               
 6     0        0 2013-01-31         2019-01-12               
 7     1        0 2014-05-03         2016-08-23               
 8     0        0 1998-02-16         2004-11-12               
 9     1        0 2014-03-04         2020-02-07               
10     1        1 1998-03-06         2003-09-08               
# ℹ 2,990 more rows

Calculating follow-up times

%--% is a special operator that creates an interval from a specific instant to another instant.

dyears(1) converts the interval to be on the years scale

bc_rt_data <-
  bc_rt_data |> 
  mutate(
    os_years = (date_of_mastectomy %--% date_last_follow_up_death) / dyears(1)
  )

Checking follow-up times

In real-world data it is common to encounter errors in data, such as end dates that come before start dates, etc.

Typically as a quick check, I look at the numeric and visual distribution of follow-up times.

Packages, in summary

Survival object

The Surv() function creates a survival object for use as the response in a model formula. There is one value for each subject that is the survival time, followed by a + if the subject was censored.

 [1] 4.865161+ 5.256674+ 7.403149+ 5.563313+ 5.152635+ 5.946612+ 2.308008+
 [8] 6.737851+ 5.930185+ 5.508556 

We see that that the first 9 subjects were censored at variaous times, and the 10th subject had an event at 5.5 years.

Kaplan-Meier

The Kaplan-Meier method is the most common way to estimate survival times and probabilities. It is a non-parametric approach that results in a step function, where there is a step down each time an event occurs.

Survival curves

The survfit() function creates survival curves using the Kaplan-Meier method based on a formula.

s1 <- survfit(Surv(os_years, os_event) ~ 1, data = bc_rt_data)

Some key components of this survfit object that will be used to create survival curves include:

• time: the timepoints at which the curve has a step, i.e. at least one event occurred
• surv: the survival probability estimate at the corresponding time

Kaplan-Meier curves with ggsurvfit

The ggsurvfit package works best if you create the survfit object using the included survfit2() function, which uses the same syntax to what we saw previously with survfit().

survfit2() tracks the environment from the function call, which allows the plot to have better default values for labeling and p-value reporting.

s2 <- survfit2(Surv(os_years, os_event) ~ 1, data = bc_rt_data)

Survival curves: common mistakes

1. Plotting the y-axis on a scale other than 0, 1
2. Plotting the x-axis beyond the limit of reasonable confidence
3. Adding too much extra information to plot face, e.g. hazard ratios, p-values, median survival time, etc
4. Using default/non-descriptive axis labels
5. Ignoring negative or missing survival times

Estimating x-time survival

summary(survfit(Surv(os_years, os_event) ~ 1, data = bc_rt_data), times = 10)

Call: survfit(formula = Surv(os_years, os_event) ~ 1, data = bc_rt_data)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   10    386     501    0.727  0.0125        0.703        0.752

x-time survival: common mistakes

Table of x-time survival probability

We can produce nice tables of x-time survival probability estimates using the tbl_survfit() function from the gtsummary package:

survfit(Surv(os_years, os_event) ~ 1, 
        data = bc_rt_data) |> 
  tbl_survfit(
    times = 10,
    label_header = 
      "**10-year survival (95% CI)**"
  )

Estimating median survival time

Call: survfit(formula = Surv(os_years, os_event) ~ 1, data = bc_rt_data)

        n events median 0.95LCL 0.95UCL
[1,] 3000    526     NA      NA      NA

Median survival: common mistakes

Using a “naive” estimate.

Summarize the median survival time among the 526 patients who died:

# A tibble: 1 × 1
  median_surv
        <dbl>
1        4.79

You get an incorrect estimate of median survival time of 4.8 years when you ignore the fact that censored patients also contribute follow-up time.

Recall the correct estimate of median survival was not reached.

Table of median survival

We can produce nice tables of median survival time estimates using the tbl_survfit() function from the gtsummary package:

survfit(Surv(os_years, os_event) ~ 1, 
        data = bc_rt_data) |> 
  tbl_survfit(
    probs = 0.5,
    label_header = 
      "**Median survival (95% CI)**"
  )

Comparing survival times between groups

Conducting the log-rank test

survdiff(Surv(os_years, os_event) ~ rt, data = bc_rt_data)

Call:
survdiff(formula = Surv(os_years, os_event) ~ rt, data = bc_rt_data)

n=2823, 177 observations deleted due to missingness.

        N Observed Expected (O-E)^2/E (O-E)^2/V
rt=0 1147      233      195      7.55      12.7
rt=1 1676      246      284      5.17      12.7

 Chisq= 12.7  on 1 degrees of freedom, p= 4e-04 

The Cox regression model

We may want to quantify an effect size for a single variable, or include more than one variable into a regression model to account for the effects of multiple variables.

The Cox regression model is a semi-parametric model that can be used to fit univariate and multivariable regression models that have survival outcomes.

\[h(t|X_i) = h_0(t) \exp(\beta_1 X_{i1} + \cdots + \beta_p X_{ip})\]

\(h(t)\): hazard, or the instantaneous rate at which events occur \(h_0(t)\): underlying baseline hazard

Cox regression assumptions

Some key assumptions of the model:

• non-informative censoring
• proportional hazards

Note that parametric regression models for survival outcomes are also available, but they won’t be addressed in this tutorial.

How to interpret a hazard ratio

The quantity of interest from a Cox regression model is a hazard ratio (HR), which represents the instantaneous rate of occurrence of the event of interest in those who are still at risk for the event.

If you have a regression parameter \(\beta\), then HR = \(\exp(\beta)\).

A HR < 1 indicates reduced hazard of event whereas a HR > 1 indicates an increased hazard of event.

Fitting Cox models

Fit regression models using the coxph() function from the survival package, which takes a Surv() object on the left hand side and has standard syntax for regression formulas in R on the right hand side.

mod1 <- coxph(Surv(os_years, os_event) ~ rt, data = bc_rt_data)

Table of Cox model results

We can creates tables of results using the tbl_regression() function from the gtsummary package, with the option to exponentiate set to TRUE to return the hazard ratio rather than the log hazard ratio:

mod1 |> 
  tbl_regression(
    exp = TRUE,
    label = list(rt ~ "PMRT")
    ) 

HR = 0.72 implies that receipt of PMRT is associated with 0.72 times the hazard of death as compared to no PMRT.

Cox regression: common mistakes

1. Overfitting. This occurs when there are too few events to support the number of included variables. Rule of thumb is 10-15 events per degree of freedom.
2. Interpreting a hazard as a risk - they are related, but they are not the same.
3. Overlooking the proportional hazards assumption

Connect with me

zabore2@ccf.org

https://www.emilyzabor.com/

https://github.com/zabore

https://www.linkedin.com/in/emily-zabor-59b902b7/

https://bsky.app/profile/zabore.bsky.social/

Further reading