Computer exercise 1: Describing and comparing survival data

About this computer exercise.

This exercise covers Chapters 1 and 2 of Collett, Modelling Survival Data in Medical Research (4th ed.). It follows on from the two lectures: Lecture 1 (basic principles, censoring and the Kaplan–Meier estimator) and Lecture 2 (the hazard and cumulative hazard, the Nelson–Aalen estimator and the log-rank test). Where a question builds on a lecture slide, the slide is noted next to the question.

Data and software

library(tidyverse)
library(sta310r)
library(survival)
library(ggsurvfit)
library(gtsummary)
library(broom)

1. From a clinical question to a time-to-event data set

A survival analysis starts with four decisions: the time origin, the time scale, the event and the censoring mechanism. These decisions are part of the estimand, not a technical detail, and are the same ingredients used throughout Lecture 1 (see the censoring and recruitment slides). The help pages (?iud, ?breast, ?prostate, ?melanoma, ?pulmonary, and ?sta310r::myeloma) give the background of each study.

1. For each (or some) of the six data sets, describe:

1. the time origin
2. time scale unit
3. the event of interest
4. what a censored observation means, for example whether the censoring is administrative, due to loss to follow-up, or due to a competing event.

Some help pages state the origin explicitly; where they do not, say what you would assume and why. You may answer in a table.

2. Look at the prostate data set.

1. How many events are there and how many observations in total?

2. Create a survfit object:

(surv_prostate <- survfit(Surv(time, status) ~ 1, data = prostate))
#> Call: survfit(formula = Surv(time, status) ~ 1, data = prostate)
#> 
#>       n events median 0.95LCL 0.95UCL
#> [1,] 38      6     NA      69      NA

What is the median survival? Is it possible to estimate from the data? Interpret the 0.95LCL and 0.95UCL.

3. What can still be estimated from these data?

4. What happens to the Kaplan–Meier curve after the last death at 69 months, and why should the far right-hand tail be interpreted cautiously?

5. Suppose the men who were censored were censored because they became too unwell to attend the clinic. In which direction would the Kaplan–Meier estimate be biased? What is the name of this assumption violation?

3. True or false (no programmiong required)?

1. A patient censored at 10 months is still in the risk set at 12 months.
2. If censoring is informative, the Kaplan–Meier estimate is always biased upwards.
3. The hazard $h(t)$ is the probability of dying at time $t$.
4. Administrative censoring at the end of a study is usually independent censoring, provided the analysis time was specified before the data were examined.
5. Left censoring means that the true event time is smaller than the observed time, and interval censoring means that the event time is known only to lie in an interval.
6. Study time is calendar time in the study, while patient time is measured from each patient’s own time origin. Administrative censoring at a fixed calendar date therefore produces different patient-time censoring times for patients recruited at different dates.

2. The empirical survivor function when there is no censoring

The pulmonary data contain the survival times of 11 men with pulmonary metastasis from osteosarcoma, with no censoring. The empirical survivor function is simply the proportion of men whose survival time is at least $t$: $$\hat S(t) = \frac{\text{number with survival time} \ge t}{n}.$$ The estimated curve is a step function, constant between adjacent death times and decreasing at each death time. This is the “estimation of the survival function” idea from Lecture 1 (slides 15–18), where, with no censoring, $\hat S(t)$ is simply the proportion still event-free at $t$.

4. Explain why this code is useful (interpret what it does, include comments and possibly rewrite it in a more readable way). Then run it and look at the resulting table:

emp_pulm <-
  pulmonary |>
  distinct(time) |>
  arrange(time) |>
  mutate(
    n_at_least = map_int(time, \(t) sum(pulmonary$time >= t)),
    S_hat = n_at_least / nrow(pulmonary)
  )

1. Is S_hat equivalemnt to the survival column from summary(survfit(Surv(pulmonary$time) ~ 1))? If not, why?

2. Plot $\hat S(t)$ above as a step function with ggplot2 (HINT: geom_step(direction = "vh")). Compare your figure to survfit2(Surv(time) ~ 1, pulmonary) |> ggsurvfit() Are they the same?

3. At which times does the curve drop, and by how much? Why is the curve flat between 15 and 17 months?

5. The observed times are 11, 13, 13, 13, 13, 13, 14, 14, 15, 15, 17 months.

Their mean is 13.7272727 months and their median is 13 months. Why is the median usually preferred for survival data? Does the mean or the median describe the “typical” patient better here? What happens to the mean if the largest observation is censored?

6. Five men die at 13 months. With the convention that the curve drops at the event time (the convention used by survfit()), the estimated survival is 0.909 just before 13 months and 0.455 at 13 months. If one of the men had been censored at exactly 13 months instead of dying, should that censored observation be counted before or after the deaths at 13 months when the risk set is formed? What difference would the other convention make to the estimate?

7. (Optional.) When there is no censoring, Greenwood’s formula reduces to the binomial standard error $$\operatorname{se}\{\hat S(t)\} = \sqrt{\frac{\hat S(t)(1-\hat S(t))}{n}}.$$

1. Use this formula to compute the standard error of $\hat S(13)$ and compare it with the std.err reported by survfit.

2. Why is the standard error NaN at the last death time, where $\hat S(t)=0$? Does that mean the estimate is known with certainty?

3. Kaplan–Meier

The iud data contain the times, in weeks, to discontinuation of an intrauterine device because of menstrual bleeding problems, for 18 women. Nine times are events and nine are censored. This is the data set used in Collett’s Examples 2.3–2.10.

The Kaplan–Meier estimator multiplies the survival factors at each event time (Lecture 1, slides 19–22): $$\hat S(t)=\prod_{j:\,t_{(j)}\le t}\left(1-\frac{d_j}{n_j}\right),$$ where $d_j$ is the number of events at the event time $t_{(j)}$ and $n_j$ the number at risk just before it.

8.

1. Reproduce Collett’s Table 2.2. At each event time $t_j$, compute the number at risk $n_j$, the number of events $d_j$, the survival factor $(n_j-d_j)/n_j$ and the product-limit estimate.

2. Report $\hat S(12), \hat S(24), \hat S(52), \hat S(104)$ and $\hat S(107)$.

9.

1. The largest observed time (max(iud$time)) is 107 weeks, but there is one event and two censorings at that time. Is $\hat S(t)$ defined for $t>107$? What does the survfit object report?
2. What would the estimate be beyond 107 if all three observations at 107 had been censored? What if the only observation at 107 had been an event?
3. Why does the curve stop at the largest observed time rather than continue to, say, 120 weeks?

10. The Kaplan–Meier-type hazard estimate in the interval from $t_{(j)}$ to $t_{(j+1)}$ is $$\hat h(t)=\frac{d_j}{n_j\tau_j},\qquad \tau_j=t_{(j+1)}-t_{(j)}.$$ The estimate is not defined after the final death time, because the last interval is open-ended.

1. Compute $\hat h(t)$ for the event times in the iud data by:

event_times <- sort(unique(iud$time[iud$status == 1]))
(haz <- tibble::tibble(
  t = event_times,
  tau = c(diff(event_times), NA_real_),
  n = map_int(event_times, \(t) sum(iud$time >= t)),
  d = map_int(event_times, \(t) sum(iud$time == t & iud$status == 1)),
  h = d / (n * tau)
))
#> # A tibble: 9 × 5
#>       t   tau     n     d        h
#>   <int> <dbl> <int> <int>    <dbl>
#> 1    10     9    18     1  0.00617
#> 2    19    11    15     1  0.00606
#> 3    30     6    13     1  0.0128 
#> 4    36    23    12     1  0.00362
#> 5    59    16     8     1  0.00781
#> 6    75    18     7     1  0.00794
#> 7    93     4     6     1  0.0417 
#> 8    97    10     5     1  0.02   
#> 9   107    NA     3     1 NA

Make sure you understand the formula and the code. Then interpret the results (you may use a graph). Why is the ungrouped estimate so much more irregular than a cumulative summary such as the survivor function?

2. Why is it usually better to report the cumulative hazard or the survivor function than a series of raw hazard estimates?

11. A colleague ignores censoring and estimates 2-year survival in the iud data as “the proportion of women who did not discontinue within 104 weeks”.

1. Compare this with the Kaplan–Meier estimate $\hat S(104)$.
2. In which direction is the naive estimate biased, and why?
3. Would the problem be worse or less serious in a study with 80 % censoring?
4. Write one sentence interpreting $\hat S(104)$ correctly.

4. Hazard and cumulative hazard

12. The cumulative hazard is $H(t)=-\log S(t)$, and the Nelson–Aalen estimate is the cumulative sum of the event fractions $d_j/n_j$.

1. Compare the estimated survival by Kaplan-Meier and Nelson-Aalen for the sta310r::myeloma data. (HINT: use argument stype, see ?survfit.formula())?

Note that sta310r::myeloma and survival::myeloma are two different data sets so use the ::-notation to access the correct version!

1. How does the estimated survival differ comparing Kaplan–Meier and Nelson–Aalen
2. What does the slope of the cumulative hazard curve represent?
3. The two estimates of $S(t)$ differs only slightly. Is the difference practically important compared with the width of the confidence intervals?

5. Comparing two groups

The breast data contain the survival times of 45 women with breast cancer, classified by whether the tumour stained positively or negatively with the lectin HPA. There are 13 women in the negative group (5 deaths) and 32 in the positive group (21 deaths). Collett uses these data in Examples 2.11–2.14.

13.

1. Plot the Kaplan–Meier curves by staining group, with a risk table and a line at the median (HINT: look at the example: https://www.danieldsjoberg.com/ggsurvfit/index.html).

2. Report $\hat S(24)$, $\hat S(60)$ and $\hat S(120)$ (You may use gtsummary::tbl_survival()).

3. Interpret the curves in clinical language.

4. Is the median estimable in both groups? If not, what should be reported instead?

5. Do the estimated curves cross? What does the answer suggest about the relative suitability of the log-rank, Wilcoxon and Peto–Peto tests? (Collett, Section 2.6.5; Lecture 2, the log-rank test and the Wilcoxon test)

6. Perform a log-rank test to compare the two groups. Report the test statistic and p-value.

6. Three or more groups

The melanoma data contain survival times of 30 patients from two immunotherapy groups, classified into three age groups. Collett uses the data in Examples 2.15 and 2.16.

14. Using only the BCG treatment group, compare survival across the three age groups with a log-rank test.

1. What is the null hypothesis? How many degrees of freedom does the test have?
2. Report the test statistic and p-value. Which age group contributes most to the test statistic?
3. Interpret the test results

15. Now compare the two treatments (treatment variable), first without and then with stratification by age group (HINT: strata(age) in the formula).

1. Report the test statistic, degrees of freedom and p-value for both tests.
2. What does the stratified test estimate that the unstratified test does not? What assumption does the stratified test make about the treatment effect across age groups?
3. Why are the two results so similar here? Check the distribution of age group within each treatment group.
4. Give an example of a situation in which stratification would change the conclusion.