Software 5.4.1 (Nelson–Aalen by hand). Example 5.1.2: the cumulative hazard and its standard error at each event time, the survival estimate \(\exp (-\hat H)\), and the Kaplan–Meier estimate for comparison:

> library(survival)
> tm <- c(3, 5, 6, 8, 10, 12, 15, 18, 22, 24)
> ev <- c(1, 0, 1, 1, 0, 1, 0, 1, 0, 0)
> km <- survfit(Surv(tm, ev) ~ 1)
> k <- km$n.event > 0                     # rows at event times
> H <- cumsum(km$n.event[k] / km$n.risk[k])   # Nelson-Aalen
> v <- cumsum(km$n.event[k] / km$n.risk[k]^2) # its variance
> round(cbind(t = km$time[k], H_NA = H, se = sqrt(v), exp_minus_H = exp(-H),
+             KM = km$surv[k]), 4)
      t   H_NA     se exp_minus_H     KM
[1,]  3 0.1000 0.1000      0.9048 0.9000
[2,]  6 0.2250 0.1601      0.7985 0.7875
[3,]  8 0.3679 0.2146      0.6922 0.6750
[4,] 12 0.5679 0.2933      0.5667 0.5400
[5,] 18 0.9012 0.4440      0.4061 0.3600

Reading the output. \(\hat H\) rises by \(d_j/n_j\) at each death, reaching \(0.901\) at \(18\) months with standard error \(0.444\). Every value of \(\exp (-\hat H)\) exceeds the Kaplan–Meier value, as Proposition 5.1.1 requires, by up to \(0.046\) at the last death.

Python: lifelines.NelsonAalenFitter().fit(tm, ev).cumulative_hazard_.

Software 5.4.2 (A Weibull check). The Nelson–Aalen estimate for lung, a least-squares line of \(\log \hat H\) on \(\log t\) through the death times, and a comparison of \(\hat H(365)\) with \(-\log \hat S_{\text {KM}}(365)\):

> fit <- survfit(Surv(time, status) ~ 1, data = lung)
> ok <- fit$n.event > 0 & fit$surv > 0
> Hl <- cumsum(fit$n.event / fit$n.risk)[ok]
> x <- log(fit$time[ok]); y <- log(Hl)
> round(coef(lm(y ~ x)), 4)               # slope estimates the Weibull shape
(Intercept)           x
    -7.3392      1.2096
> round(c(H_at_365 = Hl[max(which(fit$time[ok] <= 365))],
+         minus_log_KM = -log(summary(fit, times = 365)$surv)), 4)
    H_at_365 minus_log_KM
      0.8883       0.8934

Reading the output. The slope \(1.21\) confirms an increasing hazard, as the Weibull fit found (\(1.32\)). At one year \(\hat H = 0.888\) and \(-\log \hat S_{\text {KM}} = 0.893\): with \(65\) or more at risk throughout the first year, the two estimators nearly agree.

Python: numpy.polyfit(numpy.log(t), numpy.log(H), 1).

Software 5.4.3 (An actuarial life table). Example 5.3.1: deaths and censorings by quarter-year, the numbers entering each interval, the conditional probabilities of death with the half-interval correction, and the survival at the end of each interval; then the Kaplan–Meier estimates at \(91\), \(182\), \(364\) and \(728\) days:

> br <- seq(0, 1092, by = 91)             # quarter-year intervals
> grp <- cut(lung$time, br, right = FALSE)
> d <- unname(tapply(lung$status == 2, grp, sum))       # deaths in each
> cen <- unname(tapply(lung$status == 1, grp, sum))     # and censored
> n <- nrow(lung) - c(0, cumsum(d + cen)[-length(d)])   # entering each
> nprime <- n - cen / 2                   # effective number at risk
> q <- d / nprime; S <- cumprod(1 - q)
> tab <- data.frame(start = head(br, -1), n, deaths = d, censored = cen,
+                   q = round(q, 4), S_end = round(S, 4))
> tab[1:8, ]
  start   n deaths censored      q  S_end
1     0 228     27        0 0.1184 0.8816
2    91 201     38        6 0.1919 0.7124
3   182 157     28       23 0.1924 0.5753
4   273 106     27       12 0.2700 0.4200
5   364  67     13        7 0.2047 0.3340
6   455  47     10        5 0.2247 0.2589
7   546  32      7        3 0.2295 0.1995
8   637  22      8        0 0.3636 0.1270
> round(summary(fit, times = c(91, 182, 364, 728))$surv, 4)   # KM, to compare
[1] 0.8816 0.7081 0.4092 0.1157

Reading the output. The quarterly probability of death rises from \(0.12\) in the first quarter to about \(0.19\)–\(0.27\) thereafter, the increasing hazard again. The life table’s survival at the ends of the intervals, \(0.882\), \(0.712\), …, \(0.420\) at one year and \(0.127\) at two, tracks the Kaplan–Meier values \(0.882\), \(0.708\), \(0.409\), \(0.116\) within a little over a percentage point.

Python: lifelines.utils.survival_table_from_events gives the counts, from which the table follows.