In 1990 the FDA approved levamisole plus fluorouracil (5-FU) for patients just operated on for stage III colon cancer. The approval rested on one randomized trial: 929 patients, assigned to observation, levamisole alone, or the combination, and followed for recurrence and death. Moertel and colleagues reported that the combination cut recurrence by 40% and death by 33%, and that levamisole alone, the drug the investigators expected to work, did nothing measurable. The trial’s data ship with R’s survival package.
ColonAdjuvant rebuilds the headline numbers twice, in Python and in R, proves the verdict in Lean, and asks the questions a regulator would ask that a hazard ratio cannot answer alone.
How it is built
Python, R and Lean steps, run in order by run_all.sh. Each step writes a file the next one reads, so every intermediate result can be opened.
- Pull cancer.rda from a pinned commit of R's survival package, SHA-256 checked python/01
- Fit Cox, log-rank and Kaplan-Meier, in Python and again in R python/02 · python/03
- Compare 74 quantities across languages, then against the abstracts python/04
- Decide approve or not, months of life, the built-in null, subsamples python/05
- Prove the exact log-rank statistic and the verdict (lean/LogRank.lean) python/06 · python/07
Reproducible. The trial file is downloaded from a pinned commit of the survival package and refused if its SHA-256 hash differs, so a maintainer’s later fix cannot change the 929 patients under a rerun. Python, R and Lean are pinned in flake.nix and flake.lock. Python has no reader for R’s .rda format here, so python/rda.py is one, written from R’s own serialization spec. Exporting a CSV from R would have been easier, and would have made the “independent” Python result inherit R’s reading of the file.
Tested. 6 pytest and 2 testthat tests. The .rda reader is tested against a hand-built byte stream rather than its own output; the rest cover the rounding rule and lossless double round-trips. The comparison fails the run if Python and R disagree on any of 74 quantities beyond a relative 10⁻⁶.
The first run did not agree
On the levamisole-alone arm, Python’s lifelines gave a hazard ratio of 0.98412 and R gave 0.98408. Tie handling, the difference everyone expects, turned out to be the smaller effect. From python/03_fit.py:
# lifelines stops its Newton iterations once the Newton decrement is below
# `precision` (default 1e-7). For a coefficient this close to zero (the
# levamisole-alone arm) that leaves the fourth decimal of the hazard ratio
# unconverged; 1e-9 takes it to the same ten decimals as R and statsmodels.
cph = CoxPHFitter().fit(s[["time", "status", "trt"]], "time", "status", fit_options=FIT_OPTIONS)
Against the paper
The full papers are paywalled, so the targets are the abstracts’ numbers, rounded the way the abstracts round.
| Computed | Abstract | ||
|---|---|---|---|
| Patients | 929 | 929 | ✓ |
| Recurrence reduction, combination | 40% | 40% | ✓ |
| Recurrence p | < 0.0001 | < 0.0001 | ✓ |
| Recurrence reduction, levamisole | 2% | 2% | ✓ |
| Death reduction, combination | 31% | 33% | ✗ |
| Death p | 0.0016 | 0.0007 | ✗ |
| Death reduction, levamisole | 3% | 6% | ✗ |
Recurrence matches to the digit, so the pipeline is sound. The death side is weaker in the public data, and adjusted or stratified models (hazard ratios 0.67 to 0.72) do not close it. The package documents its data as “closest to” the final report, not identical. Four of seven, with the reason, is the result I would trust; seven of seven from data that is not the paper’s exact cut would not be.
Approve or not
A paper reports estimates; a regulator says yes or no. The rule, my own stated choice rather than the FDA’s: approve if the log-rank test gives p < 0.001 for recurrence and p < 0.01 for death. Both hold (χ² 19.07 and 9.97, against 10.83 and 6.63).
| Question | Answer |
|---|---|
| Is one hazard ratio a fair summary? | Yes: proportional-hazards tests give p = 0.99 and 0.36 |
| What is the benefit in months? | 8.4 months more life on average over eight years (95% CI 2.8 to 13.8) |
| What does no effect look like here? | Levamisole alone: hazard ratio 0.98, CI 0.79 to 1.21, a placebo built into the trial |
| How much of the trial was needed? | Nearly all of it: with 75% of patients, only 68% of draws approve, and death is the binding endpoint |
Lean computes the log-rank statistic as an exact fraction from the 619 patient rows, a numerator of 505 digits, and its kernel checks the verdict. From lean/LogRank.lean:
theorem recurrence_p_below_0_001 : Frac.lt (Frac.mk' 10828 1000) (chisq recurrence) = true := by decide
theorem death_p_below_0_01 : Frac.lt (Frac.mk' 6635 1000) (chisq death) = true := by decide
Lean does not fit the Cox model: a hazard ratio comes out of an optimiser, and certifying it would need a convergence proof. That is why the rule is stated on the log-rank test, which is a finite sum. For this trial the two rules agree.