Software 1.5.1 (A survival object). The first patients of the lung data, the coding of the status variable, the survival object, and the numbers of events and censored times:

> library(survival)
> head(lung[, c("time", "status", "age", "sex")], 5)
  time status age sex
1  306      2  74   1
2  455      2  68   1
3 1010      1  56   1
4  210      2  57   1
5  883      2  60   1
> table(lung$status)                    # 1 = censored, 2 = died

  1   2
 63 165
> y <- Surv(lung$time, lung$status == 2)
> head(y, 10)                           # a + marks a censored time
 [1]  306   455  1010+  210   883  1022+  310   361   218   166
> c(events = sum(lung$status == 2), censored = sum(lung$status == 1))
  events censored
     165       63
> round(mean(lung$status == 1), 3)      # proportion censored
[1] 0.276

Reading the output. The third patient was followed for \(1010\) days and was still alive, printed \(1010+\). There are \(165\) deaths and \(63\) censored times, so \(27.6\%\) of the observations are censored: too many to ignore.

Python: with lifelines, the pair is kept as two columns, durations and event_observed.

Software 1.5.2 (The cost of ignoring censoring). Simulated exponential lifetimes with mean \(10\), each censored at an independent uniform time on \((0, 20)\); the fraction of deaths observed, and four estimates of the mean:

> set.seed(1)
> n <- 10000
> t_true <- rexp(n, rate = 1/10)        # true lifetimes, mean 10
> cens <- runif(n, 0, 20)               # censoring time, uniform on (0, 20)
> time <- pmin(t_true, cens); event <- as.numeric(t_true <= cens)
> round(mean(event), 3)                 # fraction of deaths observed
[1] 0.571
> round(c(true = mean(t_true), all_times = mean(time),
+         deaths_only = mean(time[event == 1]),
+         rate_method = sum(time) / sum(event)), 3)
       true   all_times deaths_only rate_method
      9.984       5.616       4.674       9.844

Reading the output. Only \(57\%\) of the deaths are observed. The mean of all recorded times is \(5.6\) and the mean of the uncensored times \(4.7\), both far below the true mean of about \(10\). Dividing the total time at risk by the number of deaths gives \(9.84\), close to the truth.

Python: numpy.random.exponential and numpy.minimum give the same simulation.

Software 1.5.3 (Follow-up by group). The median recorded time for the deceased and the censored patients in lung, and the numbers and proportions of deaths by sex (\(1\) = male, \(2\) = female):

> with(lung, tapply(time, status, median))   # median time by status
  1   2
284 226
> with(lung, table(sex, status))
   status
sex   1   2
  1  26 112
  2  37  53
> round(prop.table(with(lung, table(sex, status)), 1), 3)
   status
sex     1     2
  1 0.188 0.812
  2 0.411 0.589

Reading the output. \(81\%\) of the men died during the study but only \(59\%\) of the women. This suggests that women survived longer, but the comparison is not yet fair: the groups may have been followed for different lengths of time. The Kaplan–Meier curves and the log-rank test of later sections make such comparisons properly.

Python: pandas.crosstab(df.sex, df.status, normalize=’index’).