Software 6.4.1 (The log-rank test by hand). Example 6.1.1: a loop over the death times accumulating observed and expected deaths in group \(1\) and the hypergeometric variance; then survdiff:

> library(survival)
> o <- ovarian                          # 26 patients, two treatments
> tt <- sort(unique(o$futime[o$fustat == 1]))   # distinct death times
> O1 <- E1 <- V <- 0
> for (t in tt) {
+   n <- sum(o$futime >= t); n1 <- sum(o$futime >= t & o$rx == 1)
+   d <- sum(o$futime == t & o$fustat == 1)
+   d1 <- sum(o$futime == t & o$fustat == 1 & o$rx == 1)
+   O1 <- O1 + d1; E1 <- E1 + d * n1 / n
+   if (n > 1) V <- V + d * (n1 / n) * (1 - n1 / n) * (n - d) / (n - 1)
+ }
> round(c(O1 = O1, E1 = E1, V = V, chisq = (O1 - E1)^2 / V), 4)
    O1     E1      V  chisq
7.0000 5.2335 2.9362 1.0627
> survdiff(Surv(futime, fustat) ~ rx, data = o)
Call:
survdiff(formula = Surv(futime, fustat) ~ rx, data = o)

      N Observed Expected (O-E)^2/E (O-E)^2/V
rx=1 13        7     5.23     0.596      1.06
rx=2 13        5     6.77     0.461      1.06

 Chisq= 1.1  on 1 degrees of freedom, p= 0.3

Reading the output. The hand calculation gives \(O_1 = 7\), \(E_1 = 5.23\), \(V = 2.94\) and \(\chi ^2 = 1.06\), the value in the last column of survdiff’s table, with \(p = 0.3\). The column \((O - E)^2/E\) is an older approximation that is smaller and should not be used.

Python: lifelines.statistics.logrank_test(t1, t2, e1, e2).

Software 6.4.2 (Effects, weights and strata). The log-rank test for sex in lung with its \(p\)-value and the observed-to-expected hazard ratio; the Peto–Peto version; then the veteran regimens unstratified, stratified by cell type, and by the Peto–Peto test:

> sd <- survdiff(Surv(time, status) ~ sex, data = lung)
> round(c(chisq = sd$chisq, p = pchisq(sd$chisq, 1, lower.tail = FALSE)), 5)
   chisq        p
10.32674  0.00131
> hr <- (sd$obs[2] / sd$exp[2]) / (sd$obs[1] / sd$exp[1]); round(hr, 3)
[1] 0.59
> pp <- survdiff(Surv(time, status) ~ sex, data = lung, rho = 1)  # Peto-Peto
> round(pp$chisq, 3)
[1] 12.714
> v1 <- survdiff(Surv(time, status) ~ trt, data = veteran)
> v2 <- survdiff(Surv(time, status) ~ trt + strata(celltype), data = veteran)
> round(c(unstratified = v1$chisq, stratified = v2$chisq), 3)
unstratified   stratified
       0.008        0.702
> v3 <- survdiff(Surv(time, status) ~ trt, data = veteran, rho = 1)
> round(c(peto_peto = v3$chisq,
+         p = pchisq(v3$chisq, 1, lower.tail = FALSE)), 3)
peto_peto         p
    0.871     0.351

Reading the output. For sex, \(\chi ^2 = 10.3\) (\(p = 0.0013\)) and \(\widehat {\text {HR}} = 0.59\); the Peto–Peto statistic is larger, \(12.7\), because the curves separate early. For the veteran regimens, none of the three tests finds a difference (\(p = 0.93\), \(0.40\), \(0.35\)).

Python: logrank_test(..., weightings=’peto’).

Software 6.4.3 (Three groups and a trend). Example 6.2.1: the numbers, observed and expected deaths for ECOG scores \(0\), \(1\), \(2\); the \(k\)-group test; and the trend test with scores \(0, 1, 2\) using the covariance matrix that survdiff returns:

> lg <- subset(lung, ph.ecog %in% 0:2)
> e3 <- survdiff(Surv(time, status) ~ ph.ecog, data = lg)
> round(cbind(n = e3$n, O = e3$obs, E = e3$exp), 2)
            n  O     E
ph.ecog=0  63 37 53.90
ph.ecog=1 113 82 83.09
ph.ecog=2  50 44 26.00
> round(c(chisq = e3$chisq, p = pchisq(e3$chisq, 2, lower.tail = FALSE)), 6)
    chisq         p
18.012097  0.000123
> w <- 0:2; U <- sum(w * (e3$obs - e3$exp))    # trend test, scores 0, 1, 2
> trend <- U^2 / c(t(w) %*% e3$var %*% w)
> round(c(trend = trend, p = pchisq(trend, 1, lower.tail = FALSE)), 6)
   trend        p
16.44654  0.00005

Reading the output. Patients with score \(2\) had \(44\) deaths against \(26\) expected, those with score \(0\) had \(37\) against \(54\). The overall test gives \(\chi ^2_2 = 18.0\); the trend test \(\chi ^2_1 = 16.4\) captures almost all of it with one degree of freedom.

Python: lifelines.statistics.multivariate_logrank_test.