Software 12.4.1 (Events and patients). Schoenfeld’s and Freedman’s formulas for a hazard ratio of \(0.7\), the events needed for four hazard ratios, and the conversion to patients with Simpson’s rule:

> library(survival); options(show.signif.stars = FALSE)
> events <- function(hr, alpha = 0.05, power = 0.8, p = 0.5)
+   (qnorm(1 - alpha / 2) + qnorm(power))^2 / (p * (1 - p) * log(hr)^2)
> freedman <- (qnorm(0.975) + qnorm(0.8))^2 * (1.7 / 0.3)^2
> round(c(schoenfeld = events(0.7), freedman = freedman))
schoenfeld   freedman
       247        252
> round(sapply(c(0.5, 0.6, 0.7, 0.8), events))   # events needed, 80% power
[1]  65 120 247 631
> # probability of an event: accrual a = 24, further follow-up f = 12 months
> S <- function(t, hr = 1) exp(-hr * log(2) / 12 * t)   # control median 12
> pev <- function(hr, a = 24, f = 12)
+   1 - (S(f, hr) + 4 * S(f + a / 2, hr) + S(f + a, hr)) / 6   # Simpson
> pbar <- mean(c(pev(1), pev(0.7)))
> round(c(p_control = pev(1), p_treated = pev(0.7), p_mean = pbar), 4)
p_control p_treated    p_mean
   0.7292    0.6059    0.6675
> ceiling(events(0.7) / pbar)                    # patients needed
[1] 370

Reading the output. \(247\) events (Freedman \(252\)) for \(\Delta = 0.7\); the number rises steeply as \(\Delta \) approaches \(1\). With the given accrual and follow-up, two-thirds of patients will be seen to die, so \(370\) are needed.

Python: lifelines.statistics.sample_size_necessary_under_cph.

Software 12.4.2 (Power by simulation). \(2000\) trials of the design of Example 12.1.2, with exponential survival, uniform accrual over \(24\) months and \(12\) months of further follow-up, each analysed by the log-rank test:

> set.seed(12)
> sim <- function(n, hr = 0.7, a = 24, f = 12) {
+   g <- rep(0:1, length.out = n)
+   t <- rexp(n, log(2) / 12 * hr^g)               # survival times
+   cens <- f + runif(n, 0, a)                     # time from entry to end
+   y <- pmin(t, cens); d <- as.numeric(t <= cens)
+   c(p = 1 - pchisq(survdiff(Surv(y, d) ~ g)$chisq, 1), events = sum(d))
+ }
> r <- replicate(2000, sim(370))
> round(c(power = mean(r["p", ] < 0.05), events = mean(r["events", ])), 3)
  power  events
  0.789 247.018

Reading the output. The simulated power is \(78.9\%\) (with simulation error about \(0.9\%\)), close to the \(80\%\) intended, and the trials see \(247\) deaths on average, the number the formula asked for. Changing the simulated survival times gives power for any pattern of hazards.

Python: the same loop with numpy and lifelines.statistics.logrank_test.

Software 12.4.3 (The PBC trial). The \(312\) randomized patients, deaths by arm, the log-rank test and Kaplan–Meier estimates at \(5\) and \(10\) years (trt is \(1\) for D-penicillamine and \(2\) for placebo):

> pb <- pbc[!is.na(pbc$trt), ]                     # the 312 trial patients
> pb$death <- as.numeric(pb$status == 2)           # transplant = censored
> table(trt = pb$trt, death = pb$death)
   death
trt  0  1
  1 93 65
  2 94 60
> survdiff(Surv(time / 365.25, death) ~ trt, data = pb)
Call:
survdiff(formula = Surv(time/365.25, death) ~ trt, data = pb)

        N Observed Expected (O-E)^2/E (O-E)^2/V
trt=1 158       65     63.2    0.0502     0.102
trt=2 154       60     61.8    0.0513     0.102

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7
> km <- survfit(Surv(time / 365.25, death) ~ trt, data = pb)
> matrix(round(summary(km, times = c(5, 10))$surv, 3), 2,
+        dimnames = list(years = c(5, 10), trt = c("D-pen", "placebo")))
     trt
years D-pen placebo
   5  0.708   0.715
   10 0.425   0.457

Reading the output. Observed and expected deaths are almost equal in both arms (\(65\) against \(63.2\), \(60\) against \(61.8\)), \(p = 0.75\); survival is similar at \(5\) and \(10\) years.

Python: KaplanMeierFitter and logrank_test.

Software 12.4.4 (The Mayo model). The Cox model with age, oedema and three log-transformed laboratory values, its hazard ratios per year and per doubling, the tests of proportional hazards, the concordance index, and the treatment effect adjusted for the model:

> mayo <- coxph(Surv(time, death) ~ age + edema + log(bili) + log(albumin) +
+                 log(protime), data = pb)
> unit <- c(1, 1, log(2), log(2), log(2))         # per year; per doubling
> round(cbind(HR = exp(coef(mayo) * unit), p = summary(mayo)$coef[, 5]), 4)
                 HR      p
age          1.0338 0.0001
edema        2.1943 0.0086
log(bili)    1.8404 0.0000
log(albumin) 0.1199 0.0000
log(protime) 8.0781 0.0032
> round(cox.zph(mayo)$table, 3)
              chisq df     p
age           0.571  1 0.450
edema         4.295  1 0.038
log(bili)     1.512  1 0.219
log(albumin)  0.153  1 0.696
log(protime)  5.352  1 0.021
GLOBAL       10.551  5 0.061
> round(concordance(mayo)$concordance, 3)
[1] 0.844
> m2 <- update(mayo, . ~ . + trt)
> round(summary(m2)$coef["trt", c(2, 5)], 3)       # treatment, adjusted
exp(coef)  Pr(>|z|)
    1.145     0.465

Reading the output. All five predictors are strong (the \(p\)-values shown as \(0\) are below \(0.00005\)). The proportional hazards tests are acceptable overall; \(C = 0.84\); adjusted for prognosis, the drug’s hazard ratio is \(1.15\), \(p = 0.47\).

Python: CoxPHFitter with check_assumptions and concordance_index_.