The CDISC pilot study tested a skin-patch drug, xanomeline, in mild-to-moderate Alzheimer’s disease: placebo, a low dose and a high dose, scored on the ADAS-Cog, where higher is worse. Its 2006 SAS analysis went into a report of the kind the FDA receives, and the report is public. Table 14-3.01 is the trial’s primary result: an ANCOVA of the change in score at week 24, 70 printed numbers in all.
AncovaThreeWays rebuilds that table from the raw files in Python and, separately, in R, and matches all 70. It then proves every one of them in Lean, in exact arithmetic, and asks whether the 2006 answer holds when a third of the patients are no longer carried forward.
How it is built
Python, R and Lean steps, run in order by run_all.sh. Python and R each build their own patient list from the raw files, so a data-handling bug cannot hide in both.
- Pull CDISC pilot ADaM transport files, SHA-256 checked python/01
- Fit the week-24 ANCOVA, built from raw files in Python and again in R python/02 · python/03
- Compare all 70 values of Table 14-3.01, rounded as SAS prints python/04
- Refit MMRM on observed visits, from scratch in Python against R's mmrm python/05 – python/07
- Prove the table in exact fractions, every value a theorem (lean/Table.lean) python/08 · python/09
Reproducible. Python, R and Lean are pinned in flake.nix and flake.lock, and each run writes a session-info file per language. The transport files are checked against a pinned SHA-256 fingerprint. Every number is saved to 17 significant digits and rounded once, at the comparison, by SAS’s rule: half away from zero.
Tested. 7 pytest and 3 testthat tests cover the data selection, the rounding rule and the MMRM against cases with closed-form answers. The comparison against SAS is itself a test: one value off and the run fails.
Three ways to the same 70 numbers
| Placebo | Low dose | High dose | |
|---|---|---|---|
| Mean change at week 24 (SD) | 2.5 (5.80) | 2.0 (5.55) | 1.5 (4.26) |
| Difference vs placebo (SE) | −0.5 (0.82) | −1.0 (0.84) | |
| 95% interval | (−2.1, 1.1) | (−2.7, 0.7) | |
| p vs placebo | 0.569 | 0.233 |
Python and R agree to 12 decimals, and the reason is known. SAS, R and Python name three conventions differently: the type of sums of squares, how group means are weighted, the degrees of freedom. The README shows why none of them can move the result for this model. On the way, an exact comparison caught three bugs that would never have changed the verdict and would each have broken the match: R’s write.csv keeping 15 digits, pandas’ fast CSV reader reading 2 of 702 values one bit off, and pandas reading a SAS zero as 5.4 × 10⁻⁷⁹, so placebo’s 0 mg stopped equalling 0.
The “Range” row hides fractions. It prints whole numbers, but 14 of the 468 underlying scores are prorated totals such as 56.72, which SAS printed as 57. The match needed them found, not assumed away.
Lean proves the table. Every score is an exact fraction, so the whole ANCOVA runs without rounding; the high-dose difference has a 59-digit numerator. Lean’s kernel then checks that each exact value rounds to the number SAS printed in 2006. From lean/Table.lean:
/-- All 70 published values of Table 14-3.01, checked by the kernel. -/
theorem table_14_3_01 : checks.all id = true := by decide +kernel
theorem table_count : checks.length = 70 := by decide
Change one published value from −1.0 to −1.1 and the proof fails. What is not proved is stated: intervals and p-values need a t quantile, which is not a fraction, so that one input comes from scipy with a ±10⁻⁹ margin, and every check must hold at both ends.
The 2006 method against today’s
By week 24 only 155 of 234 patients were still measured. For the other 79 the 2006 analysis carried each patient’s last score forward (LOCF), which assumes a patient who left at week 8 stopped declining at week 8. A mixed model for repeated measures (MMRM) uses only the 539 scores actually taken.
| Week 24 | LOCF (2006) | MMRM |
|---|---|---|
| High dose vs placebo | −1.01 (SE 0.84), p = 0.23 | −0.83 (SE 1.07), p = 0.44 |
| Low dose vs placebo | −0.47 (SE 0.82), p = 0.57 | −0.59 (SE 1.01), p = 0.56 |
Same verdict, very different margin: the old method made the high dose look about twice as convincing as the observed visits support. No Python package fits an MMRM the way SAS and R’s mmrm do, so python/utils_mmrm.py writes one, with REML and Satterthwaite degrees of freedom, in about 150 lines. Python and R compute the same objective to 10 decimals; their defaults part in the fifth because R’s optimiser stops at a gradient of 0.007, and told to keep going it lands where Python did.
The MMRM still assumes that patients who left resembled those who stayed. A treatment-policy estimand would need reference-based imputation (rbmi), and the README names that as the next step rather than claiming it.