POLLUX trial 2016 - example of survival analysis

Packages today

library(tidyverse)
library(survival)
library(patchwork)
library(broom)
library(survminer)
library(forestmodel)
library(haven)
library(reactable)

base <- "#1f77b4"

Time-to-event data

Image we include 10 patients between 2010 and 2020, we can illustrate them like this:

Time-to-event data

To get the survival function we simply order the observations by their length

The survival function

The survival function is a function that gives the probability that a patient (or other object of interest) will survive past a certain time.

If we assume our data is a sample of \(n\) patients from a population , as \(n \xrightarrow{} \infty\) we approximate the true survivor function.

Censoring

We almost never have fully observed data. Lets say we only have data until 2026

Censoring our example

We can’t draw the full survival curve as not everyone experiences the event (dies). However we still know the follow-up time (survival) for each patient.

Kaplan-Meier estimator

How it works is that censored observations doesn’t “influence” the curve until next event (death), where it goes lower due to fewer at risk. The Kaplan-Meier estimate of the survival function is marked with the dashed line.

Calculating KM by hand

\[\hat{S}(t)=\prod_{i:t_i\leq t}{\left(1-\frac{d_i}{n_i}\right)}\] Data from our case from before:

Example KM

First we need some data, the dataset lung can always be called in R

head(lung)

  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  56   1       0       90        90       NA      15
4    5  210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0

This dataset contains a time (survival) and status (censoring) variable which is what we need to plot a survival curve. We can try our function with this data.

S <- function(data, time, status) {
 
  t_i <- data[["time"]] |> sort() |> unique()
  n_i <- map_int(t_i, ~ nrow(data[data[["time"]] >= .x,]))
  d_i <- map_int(t_i, ~ nrow(data[data[["time"]] == .x & data[["status"]] == 2,]))
  
  tibble(t_i, n_i, d_i, survival = cumprod(1 - d_i / n_i))
}

Example KM

survfit(Surv(time, status) ~ 1, data = lung)

Call: survfit(formula = Surv(time, status) ~ 1, data = lung)

       n events median 0.95LCL 0.95UCL
[1,] 228    165    310     285     363

survfit(Surv(time, status) ~ sex, data = lung)

Call: survfit(formula = Surv(time, status) ~ sex, data = lung)

        n events median 0.95LCL 0.95UCL
sex=1 138    112    270     212     310
sex=2  90     53    426     348     550

fit <- survfit(Surv(time, status) ~ 1, data = lung)
summary(fit, times = c(10, 100, 500))

Call: survfit(formula = Surv(time, status) ~ 1, data = lung)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   10    227       1    0.996 0.00438        0.987        1.000
  100    196      30    0.864 0.02271        0.821        0.910
  500     41     107    0.293 0.03508        0.232        0.371

Example KM

S(time = time, status = status, data = lung) |> 
  head()

# A tibble: 6 × 4
    t_i   n_i   d_i survival
  <dbl> <int> <int>    <dbl>
1     5   228     1    0.996
2    11   227     3    0.982
3    12   224     1    0.978
4    13   223     2    0.969
5    15   221     1    0.965
6    26   220     1    0.961

S(time = time, status = status, data = lung) |> 
  ggplot(aes(t_i, survival)) + 
  geom_step() + theme_minimal(base_size = 26)

survfit(Surv(time, status) ~ 1, data = lung) |> 
  tidy() |> select(1:3,5) |> head()

# A tibble: 6 × 4
   time n.risk n.event estimate
  <dbl>  <dbl>   <dbl>    <dbl>
1     5    228       1    0.996
2    11    227       3    0.982
3    12    224       1    0.978
4    13    223       2    0.969
5    15    221       1    0.965
6    26    220       1    0.961

survfit(Surv(time, status) ~ 1, lung) |> 
  tidy() |> ggplot(aes(time, estimate)) + 
  geom_step() + theme_minimal(base_size = 26)

Confidence intervals on KM

Commonly calculated using Greenwood’s formula:

\[\widehat{\text{Var}}\left(\hat{S}(t)\right) = \left[\hat{S}(t)\right]^2 \sum_{i:\, t_i \le t} \frac{d_i}{n_i(n_i - d_i)}\]

Represent both the uncertainty due to censoring, but also the general sampling bias we have in small samples

survfit(Surv(time, status) ~ 1, data = lung) |> plot()

survfit(Surv(time, status) ~ sex, data = lung) |> plot(conf.int=T)

Publication ready survival curves

However, when plotting survival curves we often want some more info (risk table, censoring, etc.). Therefore i use the package survminer.

lung2 <- lung |> 
  mutate(sex = if_else(sex == 1, "Male","Female"),
         age_group = case_when(between(age, 0,59) ~ "35 - 59 years",
                               between(age, 60,69) ~ "60 - 69 years",
                               between(age, 70,99) ~ "70 - 85 years"))

surv_plot <- function(data, strat, plot_title) {
  
  dat <- data[!is.na(data[[strat]]),]
  dat$time <- dat$time/365
  labs <- sort(unique(dat[[strat]]))
  
  splot <- survminer::ggsurvplot(fit = surv_fit(as.formula(paste0("Surv(time, status) ~", strat)), data = dat), 
                                 data = dat,
                                 palette = c("#0B0405FF", "#40498EFF","#38AAACFF", "#DEF5E5FF"), 
                                 ggtheme = theme_bw(), xlim = c(0,2.5),
                                 break.time.by = 0.5,
                                 pval = TRUE,
                                 pval.method = TRUE,
                                 pval.size = 3.5,
                                 pval.coord = c(0,0.05),
                                 pval.method.coord = c(0,0.12),
                                 risk.table = "nrisk_cumcensor",
                                 risk.table.y.text = FALSE,
                                 risk.table.height = 0.3,
                                 tables.y.text = TRUE,
                                 fontsize = 3, size = 0.2,
                                 conf.int = TRUE,
                                 conf.int.alpha = 0.2,
                                 conf.int.style = "ribbon",
                                 censor = T, censor.shape = 124, censor.size = 2, legend.labs = labs)
  
  splot$plot <- splot$plot + 
    labs(title = paste(plot_title), fill = NULL, color = NULL, x = "Time from diagnosis (years)") +
    scale_y_continuous(labels = scales::percent_format())+
    guides(color = guide_legend(override.aes = list(shape = NA, linewidth=0.6)))+
    theme(panel.grid = element_blank(),
          panel.border = element_rect(color="black", fill=NA, linewidth=0.6),
          legend.justification = c(1,1),
          legend.key.size = unit(4,"mm"),
          legend.spacing.y = unit(0,"mm"),
          legend.background = element_rect(color="black", fill=NA, linewidth = 0.6*0.5),
          legend.position = c(1,1),
          legend.box.margin = margin(t = 0.6*(.pt*72.27/96/4), r = 0.6*(.pt*72.27/96/4)),
          axis.ticks = element_line(color="black"),
          axis.text = element_text(color="black"),
          axis.title = element_text(size = 10, color="black"),
          plot.background = element_rect(fill="white"),
          plot.tag.position = c(0,1))
  
  splot$table <- splot$table +
    labs(title="No. at risk (No. censored)", x = NULL, y = NULL) + 
    coord_cartesian(clip = "off", xlim=c(0,2.5),ylim=c(0,3)) +
    theme(panel.grid = element_blank(),
          plot.title = element_text(size=10),
          panel.border = element_rect(color = NA),
          axis.ticks.x = element_blank(),
          axis.text.x = element_blank())
  splot
}

pl <- list(surv_plot(data = lung2, strat = "age_group", plot_title = "Overall survival by age group"),
           surv_plot(data = lung2, strat = "sex", plot_title = "Overall survival by sex"))

pl[[2]]$table <- pl[[2]]$table + coord_cartesian(clip = "off", xlim=c(0,2.5),ylim=c(-1,2))

p <- pl[[1]]$plot + pl[[2]]$plot + pl[[1]]$table + pl[[2]]$table +
  plot_layout(design = "AB\nCD", height=c(1,0.3))

p

See here for code example on my Linkedin

Alternative endpoints

Currently we have talked about survival, commonly refered as overall survival (OS).

• Overall survival (OS) (most commonly used)
• Progression free survival (PFS/PF$) (common in palliative treatment)
• Disease free survival (DFS) (common in adjuvant treatment)
• Relapse free survival (RFS)
• Event free survival (EFS)
• Time to progression (TTP)

Definitions (what is censoring, what is event, when does time start) of some of the more exotic can differ.

Coding OS and PFS

Genereally to perform a survival analysis 3 dates are needed (to generate the status/event/delta variable and the time/survival variable)

• Date of start (e.g. diagnosis)
• Date of event (e.g. death)
• Date of censoring (e.g. administrative censoring)

In general the status indicator is coded 1 if there is a date of event, and 0 if this is missing. The time variable is coded as time difference between date of start and date of event for all that has the event, and date of start and date of end of study for everyone else.

For public data status and time are often pre-calculated as dates are sensitive. In observational studies made with registry data they have to be coded manually. When changing endpoint we need to redo variables.

Types of censoring

Censoring - You know the subject exists, but exact event time is incomplete.

• Right censoring: Event occurs after observation ends.
  • Administrative censoring: Subjects are censored at a cutoff date, often same as date of as last data update.
  • Loss to follow-up: Subjects status are unknown, e.g. due to moving to another country.
• Left censoring: Event occurred before observation began, but exact time is unknown.
• Interval-censoring: The exact event time unknown, but is known to have occurred between two specific points in time.

Truncation - Systematic sampling bias; subject is entirely missing from data.

• Left truncation: Subjects enter study only if they survive until study start.
• Right truncation: Subjects observed only if event happens before cutoff time.

Illustration of censoring

AMGEN Panitumumab trial

We load the data from a SAS file.

pds <- read_sas("AMG_data/PDS_DSA_20050203/adsl_pds2019.sas7bdat", NULL) |> 
  mutate(DTHMO = round(DTHDY/30.43,1),
         PFSMOCR = round(PFSDYCR/30.43,1))

To include the KRAS subtype we load the dataset with biomarker data.

biomark <- read_sas("AMG_data/PDS_DSA_20050203/biomark_pds2019.sas7bdat", NULL) |> 
  select(SUBJID, BMMTR1)

pds <- pds |> 
  left_join(biomark) |> 
  mutate(KRAS = case_when(BMMTR1 == "Wild-type" ~ "WT",
                          BMMTR1 == "Mutant"    ~ "MT",
                          TRUE                  ~ "Missing"))

pds_WT <- pds |> 
  filter(KRAS == "WT")

AMGEN Panitumumab trial

Here we can see some of the data.

pds |> select(-DTHDY,-PFSDYCR,-BMMTR1,-B_ECOG,-HISSUBTY,-DIAGMONS, -B_METANM, -ATRT) |> 
  reactable(bordered = T, defaultPageSize = 6, style = list(fontSize = "18px"))

OS estimate Panitumumab

survfit(Surv(DTHMO, DTH) ~ TRT, data = pds)

Call: survfit(formula = Surv(DTHMO, DTH) ~ TRT, data = pds)

                           n events median 0.95LCL 0.95UCL
TRT=FOLFOX alone         466    350   19.4    17.4    21.8
TRT=Panitumumab + FOLFOX 469    340   19.6    18.0    21.7

Here we see that for the whole cohort there is no difference in OS.

Wildtype subset:

survfit(Surv(DTHMO, DTH) ~ TRT, data = pds_WT)

Call: survfit(formula = Surv(DTHMO, DTH) ~ TRT, data = pds_WT)

                           n events median 0.95LCL 0.95UCL
TRT=FOLFOX alone         254    176   20.5    17.8    23.6
TRT=Panitumumab + FOLFOX 260    170   23.9    20.3    28.0

However when we only look at WT patients the median survival is 3.4 months longer.

PFS estimate Panitumumab

survfit(Surv(PFSMOCR, PFSCR) ~ TRT, data = pds_WT)

Call: survfit(formula = Surv(PFSMOCR, PFSCR) ~ TRT, data = pds_WT)

                           n events median 0.95LCL 0.95UCL
TRT=FOLFOX alone         254    213    9.2     7.7    10.9
TRT=Panitumumab + FOLFOX 260    212    9.6     9.2    11.3

We find 0.4 months longer median progression free survival for the treated group.

Life tables

When we present the probabilities of survival at different time points we call it life tables.

e.g. here we see that for the FOLFOX alone the 3 year survival rate is 28%.

survfit(Surv(DTHMO, DTH) ~ TRT, data = pds_WT) |> summary(times = c(1,6,12,18,24,36))

Call: survfit(formula = Surv(DTHMO, DTH) ~ TRT, data = pds_WT)

                TRT=FOLFOX alone 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1    252       0    1.000  0.0000        1.000        1.000
    6    225      23    0.908  0.0183        0.873        0.944
   12    182      43    0.734  0.0280        0.681        0.791
   18    137      44    0.556  0.0316        0.497        0.621
   24    101      34    0.416  0.0315        0.359        0.483
   36     33      28    0.282  0.0303        0.228        0.348

                TRT=Panitumumab + FOLFOX 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1    257       3    0.988 0.00662        0.976        1.000
    6    229      21    0.906 0.01819        0.871        0.943
   12    181      44    0.731 0.02796        0.678        0.787
   18    154      26    0.625 0.03067        0.567        0.688
   24    120      31    0.497 0.03184        0.438        0.564
   36     48      40    0.321 0.03067        0.266        0.387

Calculation of the Log Rank test

For each event time \(t_i\), \(O_{i,A}\) and \(O_{i,B}\) are observed number of deaths in group \(A\) (e.g. FOLFOX alone) and \(B\) (e.g. Panitumumab + FOLFOX), respectively. \(N_{i,A}\) and \(N_{i,B}\) are observed number at risk in each group. The expected number deaths in group A assuming no difference between groups:

\[E_{i,A} = N_{i,A} \frac{O_{i,A} + O_{i,B}}{N_{i,A} + N_{i,B}}\]

The Log Rank score statistics:

\[Z = \frac{(O_{A} - E_{A})^2}{\text{var}(O_{A} - E_{A})} \approx \chi^2\]

\[O_{A} - E_{A} = \sum_{i} (O_{i,A} - E_{i,A})\] Choice between \(A\) and \(B\) in calculating the statistic doesn’t matter, difference of \(A\) to \(B\) is the same as \(B\) to \(A\).

Remember Log-rank test says nothing about the magnitude of difference.

Plotting KM-curves

fit1 <- survfit(Surv(DTHMO, DTH) ~ TRT, data = pds_WT)
fit2 <- survfit(Surv(PFSMOCR, PFSCR) ~ TRT, data = pds_WT)

p1 <- ggsurvplot(fit1, risk.table = T, risk.table.y.text = F, break.time.by = 6, xlim = c(0, 48), pval = T, title = "OS", censor = T, 
                 censor.shape = 124, censor.size = 2, palette = c("grey30", base))
p2 <- ggsurvplot(fit2, risk.table = T, risk.table.y.text = F, break.time.by = 6, xlim = c(0, 48), pval = T,title = "PFS", censor = T, 
                 censor.shape = 124, censor.size = 2, palette = c("grey30", base))

p1$plot + p2$plot + p1$table + p2$table + plot_layout(height = c(1, 0.3)) & theme(legend.position = c(0.75, 0.9))

The Hazard Function

The risk of dying before time \(t+dt\) given that you are alive at \(t\) is \(P(t<T\leq t+dt|T>t)\): with \(T\) a random variable denoting the time of death/event.

\[h(t) = \frac{P(t<T\leq t+dt|T>t)}{dt}\] The hazard function represents the probability that a subject will experience an event in the next small interval of time, divided by the length of this interval, given that they have survived up to time \(t\). Instantaneous rate of death, given that you are alive at \(t\)

Cumulative hazard function and survival function are given by:

\[H(t) = \int_{0}^{t} h(u) du \qquad\qquad S(t) = exp(-H(t))\]

Sometimes the hazard function is written as: \(\lambda(t)\) and the cumulative hazard function as: \(\Lambda(t)\)

Hazard example from book

What do you think the hazard looks like for this population?

The Hazard Function

Examples of 4 hazard functions and their corresponding survival curves:

The Hazard Function for our example

As the hazard is the instantaneous rate of failure it is 0 in between events, and infinite when an event hits. This is a problem going from discrete to continuous.

For this reason the hazard makes sense fitting some distribution, but this brings additional assumptions.

Calculating hazard rates

If the hazard rate is assumed to be constant over time, you can calculate the crude hazard rate (incidence rate) for each group individually.

\[\text{Hazard rate} = \frac{N_{event}}{\text{person-time at risk}}\]

epiR::epi.conf(cbind(ev = sum(pds_WT[pds_WT$TRT=="FOLFOX alone",]$DTH), 
                     pt = sum(pds_WT[pds_WT$TRT=="FOLFOX alone",]$DTHMO)), 
               ctype = "inc.rate", method = "exact")

          est      lower      upper
ev 0.03270221 0.02804917 0.03790659

epiR::epi.conf(cbind(ev = sum(pds_WT[pds_WT$TRT=="Panitumumab + FOLFOX",]$DTH), 
                     pt = sum(pds_WT[pds_WT$TRT=="Panitumumab + FOLFOX",]$DTHMO)), 
               ctype = "inc.rate", method = "exact")

          est     lower      upper
ev 0.02927451 0.0250392 0.03402102

Cox Proportional Hazards model

Hazard rate ratio (HR) Can be derived from the Proportional Hazards model, where the hazards can change over time, but only proportionally to each other.

The proportional hazard assumption is that all individuals have the same hazard function, but a unique scaling factor. The shape of the hazard function is the same for all individuals, and only a scalar multiple changes per individual.

\[H_i(t) = a_ih(t)\]

At the core of the assumption is that \(a_i\) is not time varying, that is, \(a_i(t)=a_i\). Furthermore, the ratio of this with another subject (called the hazard rate ratio (HR)) is constant for all \(t\):

\[\frac{h_i(t)}{h_j(t)} = \frac{a_i h(t)}{a_j h(t)} = \frac{a_i}{a_j}\]

Cox Proportional Hazards model

The Proportional Hazards model, a regression model:

\[h(t, X) = \exp(X\beta)h_b(t)\] Where \(\beta\) and the baseline hazard \(h_b(t)\) are the unknown parameters estimated from data (non-parametric).

Cox Proportional Hazards model

What makes the CoxPH model special is that it does not assumes a specific shape or mathematical distribution for baseline hazard. It is estimated using partial likelihood (as opposed to maximum likelihood) where we don’t need to estimate the baseline hazard.

coxph(Surv(DTHMO,DTH) ~ TRT, data = pds_WT) |> summary()

Call:
coxph(formula = Surv(DTHMO, DTH) ~ TRT, data = pds_WT)

  n= 514, number of events= 346 

                           coef exp(coef) se(coef)      z Pr(>|z|)
TRTPanitumumab + FOLFOX -0.1149    0.8915   0.1077 -1.066    0.286

                        exp(coef) exp(-coef) lower .95 upper .95
TRTPanitumumab + FOLFOX    0.8915      1.122    0.7219     1.101

Concordance= 0.515  (se = 0.014 )
Likelihood ratio test= 1.14  on 1 df,   p=0.3
Wald test            = 1.14  on 1 df,   p=0.3
Score (logrank) test = 1.14  on 1 df,   p=0.3

Adjusted CoxPH

We can add covariates to the model

\[h(t) = exp(x_1\beta_1+x_2\beta_2+...+x_j\beta_j)h_b(t)\]

coxph(Surv(DTHMO,DTH) ~ TRT + AGE + SEX, data = pds) |> summary()

Call:
coxph(formula = Surv(DTHMO, DTH) ~ TRT + AGE + SEX, data = pds)

  n= 935, number of events= 690 

                             coef exp(coef)  se(coef)      z Pr(>|z|)  
TRTPanitumumab + FOLFOX -0.034986  0.965619  0.076324 -0.458   0.6467  
AGE                      0.006954  1.006978  0.003880  1.792   0.0731 .
SEXMale                 -0.058426  0.943248  0.079945 -0.731   0.4649  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                        exp(coef) exp(-coef) lower .95 upper .95
TRTPanitumumab + FOLFOX    0.9656     1.0356    0.8315     1.121
AGE                        1.0070     0.9931    0.9993     1.015
SEXMale                    0.9432     1.0602    0.8064     1.103

Concordance= 0.522  (se = 0.012 )
Likelihood ratio test= 3.57  on 3 df,   p=0.3
Wald test            = 3.54  on 3 df,   p=0.3
Score (logrank) test = 3.54  on 3 df,   p=0.3

Stratified CoxPH

If we know a covariate might void the proportional hazards assumption we can estimate separate baseline hazards across the strata, in this case gender. This is known as stratified cox.

coxph(Surv(DTHMO,DTH) ~ TRT + SEX, data = pds_WT) |> broom::tidy(exp=T, conf.int=T)

# A tibble: 2 × 7
  term                   estimate std.error statistic p.value conf.low conf.high
  <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 TRTPanitumumab + FOLF…    0.895     0.108    -1.03    0.302    0.724      1.11
2 SEXMale                   0.912     0.111    -0.830   0.407    0.733      1.13

coxph(Surv(DTHMO,DTH) ~ TRT + strata(SEX), data = pds_WT) |> broom::tidy(exp=T, conf.int=T)

# A tibble: 1 × 7
  term                   estimate std.error statistic p.value conf.low conf.high
  <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 TRTPanitumumab + FOLF…    0.898     0.108    -0.998   0.318    0.727      1.11

Plotting CoxPH

We can also plot the Cox model as a forest plot (remember x-axis should be log)

forest_model(coxph(Surv(DTHMO,DTH) ~ TRT + AGE + SEX + RACE + B_ECOG + KRAS + LIVERMET, data = pds), 
             format_options = forest_model_format_options(banded = F), factor_separate_line = T)

Examples of use

Examples of use

Examples of use

Extra - Partial Maximum Likelihood for COXPH

\[L(\beta) = \prod_{I\in D}\frac{exp(Z_I^T \beta)}{\sum_{k\in R(t_1)}exp(Z_k^T \beta)}\] A cool trick is that this partial maximum likelihood can be easily estimated using an optimizer in R, like GLM models that rely on MLE.

lung <- lung |> filter(!duplicated(time)) # classic Cox does not handle ties, Efron/Breslow

partial_lik <- function(beta, data) {
  da <- arrange(data, time)
  db <- da[da$status == 2,]
  aa <- map_dbl(db$time, ~ if_else(db[db$time == .x, "sex"] == 1, 1, exp(beta)))
  bb <- map_dbl(db$time, ~ nrow(da[da$time >= .x & da$sex == 1,]))
  cc <- map_dbl(db$time, ~ nrow(da[da$time >= .x & da$sex == 2,]) * exp(beta))
  prod(aa / (bb + cc))
}

lung |> optimize(partial_lik, c(-9, 9), maximum = TRUE, tol = 0.000001, data = _) |> pluck("maximum")

[1] -0.4754432

lung |> coxph(Surv(time, status) ~ sex, data = _) |> pluck("coefficients", 1)

[1] -0.4754431