The Method Lab

This page is the repository's own document docs/method-lab.md, shown as it is, so that it never drifts from the product.

The statistician's harness for each analysis of the plan.

Why#

An AI coding agent is strong because it has a harness: it runs the code, tests it, debugs it and checks the result, and it acts through these tools. A statistician who writes a SAP has the same needs for each analysis: what exactly does the model say, which program computes it, do its assumptions hold, does the test do what the plan says (its type I error, its power), and what does the table look like. The Method Lab gives the statistician, and their AI agent through MCP, these tools for every analysis of the plan. Each result is data (JSON) with a verdict line, and files a reviewer can read.

What it gives for a block#

sapc lab <plan.qmd> <block-id> [--part formula|program|assumptions|simulate|example|all] [--data DIR] [--unblinded] [--nsim N] [--seed S]

A block is an analysis with a model or a test (P_ADAS), a sensitivity or supplementary analysis (SENS_TP), or a table, estimand or method that leads to one. sapc lab <plan.qmd> lists the blocks; sapc lab <plan.qmd> --all runs the lab for all of them and writes the report out/lab/index.html.

The files go to out/lab/<block>/ next to the plan: formula.json, formula.md, formula.html, program.R, program.sas, pseudocode.txt, program.json, assumptions.json and its PNGs, simulation.json and its PNGs, example.json, example/shell.txt, example/data/*.csv (dummy), summary.json (one verdict per part), spec.json (the analysis as the lab read it) and data.R (the program's data steps, which the checks use too).

The last line of sapc lab is SAPLING LAB PASS | CHECK | FAIL: CHECK when a check asks for a look, FAIL when a part could not run.

ICH E9: blinded before unblinding#

ICH E9 (and its blind review) lets the statistician look at the data before unblinding only without the treatment groups. The lab keeps this rule; the mode is in every result and on every table.

An agent never sets unblinded on its own: the tool's description says so.

The pilot (CDISCPILOT01, pilot-reviewed.qmd, the ADaM data of the chain)#

From the repository, R 4.4 with the packages of SAPLING_R_LIBS:

D=~/r-agentic/studies/cdiscpilot01/output/adam      # read only
node src/cli.ts lab examples/cdiscpilot01/pilot-reviewed.qmd --all --part all --data $D      # 11 blocks, about 1 minute
node src/cli.ts lab examples/cdiscpilot01/lab-mmrm.qmd P_ADAS_MMRM --part all --data $D --nsim 200

• Formula of P_ADAS: $Y_i = \mu + \beta\,\mathrm{BASE}i + \gamma{s(i)} + \tau d_i + \varepsilon_i$, $\varepsilon_i \sim N(0,\sigma^2)$ iid; the estimand EST1's difference in means is $\theta(d) = \tau d$; $H_0: \tau = 0$ two-sided at 0.05; $F = \frac{(RSS_0 - RSS_1)/1}{RSS_1/(n-p)} \sim F_{1,n-p}$ (= the Type III test of PROC GLM); CI $\hat\tau \pm t_{0.975,n-p}\,SE$.
• Program: the generated R program on the real ADQSADAS (UNBLINDED, in the test test/lab-r.test.ts only) gives n = 234 (79, 81, 74) and p = 0.2447057 = the chain's P_EFF (ard/p_eff.csv, raw 0.244705673868529, shown 0.245); slope −0.0118 per mg (95% CI −0.0317 to 0.0081). Its SAS is the chain's macro ars_e_ancova: CLASS SITEGR1; MODEL CHG = TRTPN SITEGR1 BASE / SS3; and the Type III row of TRTPN.
• Assumptions of P_ADAS on BLINDED data (TRTA, TRTAN, TRTP, TRTPN, ARM, TRT01P/PN, TRT01A/AN, AVGDD, CUMDOSE removed): residual normality OK (Shapiro-Wilk p = 0.021, skewness 0.38: below the CHECK rule), variance across the SITEGR1 levels OK (p = 0.15), BASE linear OK (p = 0.31), influence OK (max Cook's distance 0.06), independence OK; CHECK: the smallest SITEGR1 level has 8 subjects (fewer than 3 per group); 3 checks DEFERRED to unblinding. P_RACE: CHECK, expected counts down to 0.33 (3 of 9 cells below 5): Fisher's exact test. P_AE: CHECK, 220 of 230 terms have an expected count below 5. P_AGE: CHECK, normality (skewness −0.70).
• Simulation of P_ADAS (1000 per scenario, seed 20261110, 3.6 s): type I error 0.057 (0.044–0.073); power at the plan's assumptions (effect 3.0, SD 6.5, 100 per group, no dropout) 0.915 (0.896–0.931). The block samplesize block SS_ADAS plans 90% for a two-sample t-test (analytic 0.904); the planned test is the slope over 0, 54, 81 mg: with the same 3.0 at the high dose and a linear dose response its information is $n\sum(d-\bar d)^2 = 340200$, the SE of the high-dose effect 0.903 against 0.919, analytic power 0.914: the simulation agrees. With the blinded data's residual SD 5.15 (the baseline and the site explain part of the 6.5), its dropout (34% without a Week 24 value) and the plan's LOCF, the power is 0.950 (0.935–0.962), but the estimate of the high-dose effect is biased by −0.58 against the effect had nobody stopped (coverage 0.876): LOCF carries the earlier, smaller effects forward. Under MNAR (the patients who stop lose δ of the benefit), the bias against the treatment-policy effect goes from −0.52 (δ = 0) to +1.41 (δ = 6). P_ADAS_SLOPE (the same model without a missing-data block: observed cases) has power 0.899.
• MMRM (lab-mmrm.qmd, the change at Weeks 8, 16, 24 without the LOCF rows, covariance=us, fallback={toeph; ar1h; toep; ar1; cs}, Kenward-Roger, REML), BLINDED: us converges (AIC 3114.3, 0.6 s); each fallback converges too (toeph AIC 3112.3, ar1h 3128.8, toep 3132.5, ar1 3146.6, cs 3130.6). Normalized residuals: CHECK (excess kurtosis 2.7). The plan has no MNAR sensitivity analysis for it: CHECK. Simulation (200, 31 s): type I 0.060 (0.035–0.102), power at the plan's assumptions 0.905, with the blinded data 0.890; no bias under random dropout.
• Example runs: dummy data (the tipping point, 100 imputations, 13 deltas; the MMRM with us), and the real data with a dummy randomisation for every block.

The MCP tools#

sapling_lab_formula, sapling_lab_program, sapling_lab_assumptions, sapling_lab_simulate, sapling_lab_example (mcp/server.ts, defined in src/lab/mcp.ts). Input: file, block (omit it for the list of lab blocks), and for the R parts data, unblinded, nsim, seed. Output: {block, analysis, method, dir, verdicts, verdict, <part>} with the part's JSON below. Verdicts:

FORMULA: ANCOVA, test for linear dose response; 5 part(s): model, estimand, hypotheses, statistic, estimator
PROGRAM: R (stats), SAS (PROC MEANS, PROC GLM), pseudocode
ASSUMPTIONS: CHECK (levels_SITEGR1) [BLINDED; 5 ok, 3 deferred to unblinding]
SIMULATION: power 0.950 (0.935–0.962), at the plan's assumptions 0.915; type I 0.057 (0.044–0.073); 1000 simulations, seed 20261110, 3.63 s [BLINDED]
EXAMPLE: the program ran on blinded data (dummy randomisation): n 234, p 0.607, estimate 0.005 (1.1 s) [EXAMPLE]

The assistant's context includes a block's lab results when the lab has run (src/lab/context.ts): "Explain" adds them, and a request that names a block (or has the cursor in it) gets them, for questions such as "Which covariance structure?".

The UI contract (for the editor's Method Lab panel)#

The visual editor shows the lab as tabs on each analysis card: Formula · Program · Assumptions · Simulation · Example. It uses only this API and these JSON files; no editor code is in src/lab/.

API (src/lab/index.ts, Node; the bundle vscode/sapling/src/lab/ is compiled from it):

c is the compiled plan (analyze(path, text) of src/service.ts). Run the R parts off the UI thread (a child process: node src/cli.ts lab <plan> <id> --part <part>, or a worker) and show a progress state; R parts can be cancelled by killing the process.

JSON (each has schema, block, analysis, test, verdict):

Plots are PNG files named in plot/plots, relative to the block's folder. Statuses: OK (green), CHECK (amber: a person looks), DEFERRED (grey: needs unblinding), INFO (grey), ERROR (red). A verdict that starts with ...: ERROR is a part that could not run (no R, no data set): show its message.

Rules for the panel: show the data mode on every tab that has one; never offer "unblinded" as a default or remember it between sessions (ask the statistician each time); the example numbers are never shown without their label; the panel writes nothing into the plan. The bundle (vscode/scripts/bundle.mjs) must also copy r/lab/lab.R for the R parts.

Requirements and settings#

• R ≥ 4.1 (SAPLING_RSCRIPT to choose it) with jsonlite; for the methods: emmeans, mmrm, survival; haven for .xpt and .sas7bdat; ragg for the plots when there is no X11 (else the platform's PNG device). Nothing is installed: SAPLING_R_LIBS names the library (default here: the study-lens renv library of this workspace).
• SAPLING_LAB_DATA: the default data folder (--data). The files are <data set>.rds, .xpt, .sas7bdat or .csv.
• Tests: test/lab.test.ts (no R), test/lab-r.test.ts (R and the pilot data; skipped without them), test/lab-mcp.test.ts.

Limits#

• Multiple imputation without mice or rbmi: the program imputes by sequential Bayesian normal regression (as PROC MI MONOTONE REG) in base R, and combines by Rubin's rules with the Barnard-Rubin df. Jump to reference is approximated by the reference group's regressions on the patient's own observed visits (copy reference for monotone dropout); exact J2R needs rbmi.
• The simulation of a sensitivity block is the simulation of its analysis (with the MNAR shifts); the logistic and CMH simulations use the unstratified table test; binary and time-to-event power needs rates= or hr= on the samplesize block.
• The SAS programs are generated and compared as text with the chain's macros; SAS is not run here.
• The lab's dropout model is random (MCAR, a case of MAR); MNAR enters as shifts of the unobserved values.

See also docs/design.md: the design simulation (several scenarios, intercurrent events and dropout mechanisms, the true value of each estimand, the power curve over N) uses this lab's spec and the R twin's analyses, for a design plan that has no data yet.