4.4 Kaplan–Meier estimation in software

Software 4.4.1 (The estimator by hand). Example 4.1.1: a loop over the distinct event times computing the number at risk, the deaths and the running product; then survfit for comparison:

> library(survival)
> tm <- c(3, 5, 6, 8, 10, 12, 15, 18, 22, 24)       # the ten patients of
> ev <- c(1, 0, 1, 1, 0, 1, 0, 1, 0, 0)             # Section 3, 0 = censored
> S <- 1; out <- NULL
> for (tj in sort(unique(tm[ev == 1]))) {           # event times only
+   n <- sum(tm >= tj); dj <- sum(tm == tj & ev == 1)
+   S <- S * (1 - dj / n)
+   out <- rbind(out, c(t = tj, n_risk = n, deaths = dj, S = round(S, 4)))
+ }
> out
      t n_risk deaths      S
[1,]  3     10      1 0.9000
[2,]  6      8      1 0.7875
[3,]  8      7      1 0.6750
[4,] 12      5      1 0.5400
[5,] 18      3      1 0.3600
> km <- survfit(Surv(tm, ev) ~ 1)
> summary(km)$surv                                  # the same steps
[1] 0.9000 0.7875 0.6750 0.5400 0.3600

Reading the output. The table reproduces the hand calculation, and survfit gives the same five values: \(0.9\), \(0.7875\), \(0.675\), \(0.54\), \(0.36\).

Python: lifelines.KaplanMeierFitter().fit(tm, ev).survival_function_.

Software 4.4.2 (Survival, intervals and the median). The lung data: the estimated survival at \(180\), \(365\) and \(730\) days with the number at risk and the \(95\%\) interval; the median with its interval; and Greenwood’s standard error at one year by hand, against survfit’s:

> fit <- survfit(Surv(time, status) ~ 1, data = lung)
> s <- summary(fit, times = c(180, 365, 730))
> round(cbind(time = s$time, n_risk = s$n.risk, surv = s$surv,
+             lower = s$lower, upper = s$upper), 3)
     time n_risk  surv lower upper
[1,]  180    160 0.722 0.666 0.783
[2,]  365     65 0.409 0.345 0.486
[3,]  730     13 0.116 0.072 0.187
> fit                                               # median and its 95% CI
Call: survfit(formula = Surv(time, status) ~ 1, data = lung)

       n events median 0.95LCL 0.95UCL
[1,] 228    165    310     285     363
> use <- fit$time <= 365                      # Greenwood by hand, at 1 year
> terms <- fit$n.event / (fit$n.risk * (fit$n.risk - fit$n.event))
> v <- sum(terms[use])
> round(c(greenwood_se = s$surv[2] * sqrt(v), survfit_se = s$std.err[2]), 5)
greenwood_se   survfit_se
     0.03582      0.03582

Reading the output. One-year survival is \(40.9\%\) (\(34.5\%\) to \(48.6\%\)). The median is \(310\) days (\(285\) to \(363\)). Greenwood’s sum gives a standard error of \(0.03582\) at one year, exactly the value survfit uses.

Python: kmf.median_survival_time_ and kmf.confidence_interval_.

Software 4.4.3 (Curves by group and median follow-up). Kaplan–Meier estimates for men and women, with the numbers of subjects and events, the medians and their intervals; then the reverse Kaplan–Meier estimate of median follow-up:

> bysex <- survfit(Surv(time, status) ~ sex, data = lung)
> tab <- summary(bysex)$table
> tab[, c("records", "events", "median", "0.95LCL", "0.95UCL")]
      records events median 0.95LCL 0.95UCL
sex=1     138    112    270     212     310
sex=2      90     53    426     348     550
> rev_km <- survfit(Surv(time, status == 1) ~ 1, data = lung)  # reverse KM
> summary(rev_km)$table["median"]                   # median follow-up, days
median
   588

Reading the output. Men: \(138\) patients, \(112\) deaths, median \(270\) days (\(212\) to \(310\)). Women: \(90\) patients, \(53\) deaths, median \(426\) days (\(348\) to \(550\)). The intervals barely overlap, which suggests a real difference. The median follow-up is \(588\) days.

Python: fit a KaplanMeierFitter to each group and plot both on one axis.

Questions on this section

Stuck on something here? Ask below and it stays attached to this topic.