Survival function

The most common one is the survival function.

Let \(T\) be a non-negative continuous random variable, representing the time until the event of interest. The survival function \(S(t)\) is the probability that a randomly chosen individual is still at risk at time \(t\), where \(0 \le t \le +\infty\). For each \(t\), it is given by

\[ \begin{align*} S(t) &= P(T > t)\\ &= 1 - P(T \le t)\\ &= 1 - F(t)\\ &= 1 - \int^t_0 f(u)\text{d}u, \end{align*} \]

where \(f(\cdot)\) and \(F(\cdot)\) are the density and the cumulative distribution functions of \(T\), respectively.

\(S(t)\) represents, for each time \(t\), the probability that the time until the event is greater than this time \(t\). In other words, it models the probability that the event of interest happens after \(t\).

In the context of our examples mentioned above, it gives the probability that:

• a randomly selected patient will survive beyond time \(t\) or the proportion of patients still alive after time \(t\),
• a cab takes more than \(t\) minutes to arrive, or
• an unemployed person take more than \(t\) months to find a new job.

The survival function \(S(t)\) is:

• a decreasing function,
• taking values in \([0, 1]\) (since it is a probability), and
• equal to 1 at \(t = 0\) (i.e., \(S(0) = 1\)) and 0 at \(t = \infty\) (i.e., \(S(\infty) = 0\)).

Visually we have:

The curve shows the proportion of individuals (or experimental units) who, as time goes on, have not experienced the event of interest. As time progresses, events occur, so the proportion who have not experienced the event decreases.

Cumulative hazard function

The cumulative hazard function, which is the total hazard experienced up to time \(t\), is defined as:

\[H(t) = -log\left(S(t)\right)\]

and has the following properties:

• increasing function,
• taking value in \([0, +\infty]\), and
• \(S(t) = exp(-H(t))\).

Hazard function

The hazard function \(h(t)\), or hazard rate, defines the instantaneous event rate at time \(t\) for an individual still at risk at that time. It can be obtained by

\[ \begin{align*} h(t) &= \lim_{\Delta t \rightarrow 0} \frac{P(t \le T < t + \Delta t | T \ge t)}{\Delta t}\\ &= \frac{d}{dt} H(t)\\ &= \frac{f(t)}{S(t)}. \end{align*} \]

and has the following properties:

• positive function (not necessarily increasing or decreasing)
• the hazard function \(h(t)\) can have many different shapes and is therefore a useful tool to summarize survival data

In the context of cancer research when death is the event of interest, \(h(t)\) measures the instantaneous risk of dying right after time \(t\) given the individual is alive at time \(t\).

To link the hazard rate with the survival function; the survival curve represents the hazard rates. A steeper slope indicates a higher hazard rate because events happen more frequently, reducing the proportion of individuals who have not experienced the event at a faster rate. On the contrary, a gradual and flatter slope indicates a lower hazard rate because events occur less frequently, reducing the proportion of individuals who have not experiences the event at a slower rate. More formally:

\[S(t) = \exp\left(-\int^t_0 h(u) \text{d}u\right).\]

Note that, in contrast to the survival function which focuses on not having an event, the hazard function focuses on the event occurring.

By hand

To understand how it works, let’s first estimate it by hand on the following dataset:¹

where:

• subject is the individual’s identifier
• time is the time to event (in years)²
• event is the event status (0 = censored, 1 = event happened)

Remember that for each subject, we need to know at least 2 pieces of information:

1. the time until the event of interest or the time until the censoring, and
2. whether we have observed the event of interest or if we have observed censoring.

We first need to count the number of distinct event times. Ignoring censored observations, we have 5 distinct event times:

2, 5, 7, 9 and 16

The easiest way to do the calculation by hand is by filling the following table (a table with 5 rows since there are 5 distinct event times):

We fill columns one by one:

• \(y_{(j)}\) = the ordered distinct event times:

2, 5, 7, 9 and 16

So the table becomes:

• \(d_{(j)}\) = the number of observations for each distinct event time. For this, the frequency for each distinct event time is useful:

## time
##  2  5  7  9 16 
##  2  1  1  1  2

The table becomes:

• \(R_{(j)}\) = the remaining number of individuals at risk. For this, the distribution of time (censored and not censored) is useful:

## time
##  2  3  5  7  9 16 18 
##  2  1  2  1  1  2  1

We see that:

• At the beginning there are 10 subjects
• Just before time \(t = 5\), there are 7 subjects left (10 subjects - 2 who had the event - 1 who is censored)
• Just before time \(t = 7\), there are 5 subjects left (= 10 - 2 - 1 - 2)
• Just before time \(t = 9\), there are 4 subjects left (= 10 - 2 - 1 - 2 - 1)
• Just before time \(t = 16\), there are 3 subjects left (= 10 - 2 - 1 - 2 - 1 - 1)

The table becomes:

• \(1 - \frac{d_{(j)}}{R_{(j)}}\) is straightforward, so the table becomes:

The Kaplan-Meier estimator is:

\[\hat{S}_{KM}(t) = \prod_{j:y_{(j)} \le t} \left(1 - \frac{d_{(j)}}{R_{(j)}} \right)\]

For each \(j\), we thus take the cumulative product:

• \(j_1 = 0.8\)
• \(j_2 = 0.8 \cdot 0.857 = 0.6856\)
• \(j_3 = 0.6856 \cdot 0.8 = 0.54848\)
• \(j_4 = 0.54848 \cdot 0.75 = 0.41136\)
• \(j_5 = 0.41136 \cdot 0.333 = 0.1369829\)

So finally, we have the survival probabilities (rounded to 3 digits):

We can now represent graphically the Kaplan-Meier estimator:

To draw this survival curve, remember that:

• the x-axis corresponds to the time variable in the initial dataset, and
• the y-axis corresponds to the survival probabilities found above.

In R

We now compare our results with the results found in R.

We first create the dataset with the time and event variables:

# create dataset
dat <- data.frame(
  time = c(3, 5, 7, 2, 18, 16, 2, 9, 16, 5),
  event = c(0, 1, 1, 1, 0, 1, 1, 1, 1, 0)
)

We then run the Kaplan-Meier estimator with the survfit() and Surv() functions:

# KM
library(survival)

km <- survfit(Surv(time, event) ~ 1,
  data = dat
)

Notice that the Surv() function accepts two arguments:

1. the time variable, and
2. the event variable.

The ~ 1 in the survfit() function indicates that we estimate the Kaplan-Meier without any grouping. See more on this later in the post.

Finally, we display the results and draw the Kaplan-Meier plot:

# results
summary(km)

## Call: survfit(formula = Surv(time, event) ~ 1, data = dat)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     2     10       2    0.800   0.126       0.5868        1.000
##     5      7       1    0.686   0.151       0.4447        1.000
##     7      5       1    0.549   0.172       0.2963        1.000
##     9      4       1    0.411   0.176       0.1782        0.950
##    16      3       2    0.137   0.126       0.0225        0.834

# plot
plot(km,
  xlab = "Time",
  ylab = "Survival probability",
  conf.int = FALSE
)

The survival probabilities can be found in the survival column. Remark that results by hand and in R are similar (any difference with the results by hand is due to rounding).

Alternatively, we can use the ggsurvplot() function within the {survminer} package:

library(survminer)

# plot
ggsurvplot(km,
  conf.int = FALSE,
  legend = "none"
)

Note that the crosses on the survival curve denote the censored observations.

The advantage with the ggsurvplot() function is that it is easy to draw the median survival directly on the plot:³

ggsurvplot(km,
  conf.int = FALSE,
  surv.median.line = "hv",
  legend = "none"
)

To find the median survival:⁴

summary(km)$table["median"]

## median 
##      9

# or more simply
km

## Call: survfit(formula = Surv(time, event) ~ 1, data = dat)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 10      7      9       5      NA

Suppose that the event of interest is death:

• At time zero, the survival probability is 1 (100% of the subjects are alive).
• The median indicates that the median survival time is 9 years.⁵ This is the time at which the survival \(S(t)\) is 50%. In other words, is the time after which half of the subjects are expected to have died.
• From the plot, we also see that \(S(5) = P(T > 5 \text{ years}) =\) Probability of survival of more than 5 years for these subjects = 75%. This means that 75% of all subjects survive longer than 5 years, and that 25% of all subjects die within the first 5 years.

For the sake of completeness, let’s do another example with a much larger dataset; the tongue dataset within the {KMsurv} package.⁶

# load data
library(KMsurv)
data(tongue)

# preview data
head(tongue)

##   type time delta
## 1    1    1     1
## 2    1    3     1
## 3    1    3     1
## 4    1    4     1
## 5    1   10     1
## 6    1   13     1

• type is the tumor DNA profile (1 = aneuploid tumor, 2 = diploid tumor)
• time is the time to death or on-study time (in weeks)
• delta is the death indicator (0 = alive, 1 = dead)

For this example, we focus on the aneuploid type:

anaploid <- subset(tongue, type == 1)

We can now plot the estimated survival function and estimate the median time to death. Since it is an estimator, we can also construct a confidence interval for the estimated survival at each time \(t\) and for the estimated median survival time.⁷

# results
fit <- survfit(Surv(time, delta) ~ 1,
  data = anaploid,
  conf.type = "log-log"
)

fit

## Call: survfit(formula = Surv(time, delta) ~ 1, data = anaploid, conf.type = "log-log")
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 52     31     93      65     157

# plot
ggsurvplot(fit,
  surv.median.line = "hv",
  legend = "none"
)

The median survival time is estimated to be 93 weeks, with a 95% confidence interval between 65 and 157 weeks.

Kaplan-Meier curves can be seen as descriptive statistics for survival data. We now focus on the second branch of statistics, hypothesis testing which allows to draw conclusions on the population based on a sample (see a quick reminder about the difference between population and sample if you need).

Log-rank test

In this article, we focus on comparing survival between two groups using the log-rank test (also known as Mantel-Cox test). This test is the most common hypothesis test to compare survival between two groups.

The intuition behind the test is that if the two groups have different hazard rates, the two survival curves (so their slopes) will differ. More precisely, the log-rank test compares the observed number of events in each group to what would be expected if the survival curves were identical (i.e., if the null hypothesis were true).

Note that, as for the Kaplan-Meier estimator, the log-rank test is a nonparametric test, which makes no assumptions about the survival distributions.

For this example, consider the following dataset:

##    patient group time event
## 1        1     1  4.1     1
## 2        2     1  7.8     0
## 3        3     1 10.0     1
## 4        4     1 10.0     1
## 5        5     1 12.3     0
## 6        6     1 17.2     1
## 7        7     2  9.7     1
## 8        8     2 10.0     1
## 9        9     2 11.1     0
## 10      10     2 13.1     0
## 11      11     2 19.7     1
## 12      12     2 24.1     0

where:

• patient is the patient’s identifier
• group is the group (group 1 or 2)
• time is the time to death (in years)⁹
• event is the event status (0 = censored, 1 = death)

Suppose we are interested in comparing group 1 and 2 in terms of survival, that is, we compare survival curves between the 2 groups:

• \(H_0 : S_1(t) = S_2(t)\) for all \(t\)
• \(H_1 : S_1(t) \ne S_2(t)\) for some \(t\)

It is a statistical test, so if the \(p\)-value < \(\alpha\) (usually 0.05), we reject the null hypothesis and we conclude that survival (or the time to event) is significantly different between the two groups considered.

To perform the log-rank test, the following test statistic will be useful:

\[\begin{eqnarray} U &=& \sum_{j=1}^r w(y_{(j)})\left(O_j - E_j\right) \\ &=& \sum_{j=1}^r w(y_{(j)})\left( d_{(j)1} - \frac{d_{(j)}R_{(j)1}}{R_{(j)}}\right) \end{eqnarray}\]

with \(U^{obs} = \frac{U}{\sqrt{Var(U)}} \sim N(0,1)\) and¹⁰

\[\begin{eqnarray} Var(U) &=& \sum_{j=1}^r w^2(y_{(j)}) \frac{N_{(j)}}{ D_{(j)} }\\ &=& \sum_{j=1}^r w^2(y_{(j)}) \frac{ d_{(j)} \frac{R_{(j)1}}{R_{(j)} } \left( 1 - \frac{R_{(j)1}}{R_{(j)} } \right) \left( R_{(j)} - d_{(j)}\right) }{ R_{(j)} - 1 } \end{eqnarray}\]

## time
##  4.1   10 17.2 
##    1    2    1

## time
##  4.1  7.8   10 12.3 17.2 
##    1    1    2    1    1

## time
##  9.7   10 19.7 
##    1    1    1

## time
##  9.7   10 11.1 13.1 19.7 24.1 
##    1    1    1    1    1    1

## time
##  4.1  9.7   10 17.2 19.7 
##    1    1    3    1    1

## time
##  4.1  7.8  9.7   10 11.1 12.3 13.1 17.2 19.7 24.1 
##    1    1    1    3    1    1    1    1    1    1

dat <- data.frame(
  group = c(rep(1, 6), rep(2, 6)),
  time = c(4.1, 7.8, 10, 10, 12.3, 17.2, 9.7, 10, 11.1, 13.1, 19.7, 24.1),
  event = c(1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0)
)

dat

##    group time event
## 1      1  4.1     1
## 2      1  7.8     0
## 3      1 10.0     1
## 4      1 10.0     1
## 5      1 12.3     0
## 6      1 17.2     1
## 7      2  9.7     1
## 8      2 10.0     1
## 9      2 11.1     0
## 10     2 13.1     0
## 11     2 19.7     1
## 12     2 24.1     0

survdiff(Surv(time, event) ~ group,
  data = dat
)

## Call:
## survdiff(formula = Surv(time, event) ~ group, data = dat)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## group=1 6        4     2.57     0.800      1.62
## group=2 6        3     4.43     0.463      1.62
## 
##  Chisq= 1.6  on 1 degrees of freedom, p= 0.2

fit <- survfit(Surv(time, event) ~ group, data = dat)

ggsurvplot(fit,
  pval = TRUE,
  pval.method = TRUE
)

# load data
data(burn)

# preview data
head(burn)

##   Obs Z1 Z2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 Z11 T1 D1 T2 D2 T3 D3
## 1   1  0  0  0 15  0  0  1  1  0   0   2 12  0 12  0 12  0
## 2   2  0  0  1 20  0  0  1  0  0   0   4  9  0  9  0  9  0
## 3   3  0  0  1 15  0  0  0  1  1   0   2 13  0 13  0  7  1
## 4   4  0  0  0 20  1  0  1  0  0   0   2 11  1 29  0 29  0
## 5   5  0  0  1 70  1  1  1  1  0   0   2 28  1 31  0  4  1
## 6   6  0  0  1 20  1  0  1  0  0   0   4 11  0 11  0  8  1

# fit
fit <- survfit(Surv(T3, D3) ~ Z1, data = burn)

# plot with log-rank test
ggsurvplot(fit,
  pval = TRUE,
  pval.method = TRUE
)