Sapling for the statistical design of a trial

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

The statistician designs the statistics before the protocol is final. Sapling simulates the design in R from the assumptions, the statistician adjusts and refines, and the result becomes the statistical section of the protocol and, later, the start of the SAP.

Simulated data are assumptions made visible, not evidence. Every number a simulation gives follows from the assumptions of the design plan. No number of a simulation is a result of a study, and the report says so at its top.

What it is for#

A sample size is a sentence in the protocol: "300 patients give 90% power for a difference of 3 with SD 6.5". The sentence hides what the study will really do: patients stop the drug (an intercurrent event), patients leave (missing data), the estimand says what to do with both, and the analysis may or may not estimate that estimand. A simulation makes these visible before the protocol is final:

• the true value of each estimand under the assumptions (a treatment-policy and a hypothetical estimand have different true values under the same assumptions);
• the power, the type I error, the bias and the coverage of the planned analysis against that true value, each with its Monte Carlo standard error, under the expected assumptions and under pessimistic scenarios;
• the smallest N that reaches the power, under each scenario;
• how other analyses do (ANCOVA with LOCF against MMRM, with or without the strata): --compare;
• the findings: where what the plan states differs from what the simulation found.

The commands#

sapling design new <study> [--kind continuous|tte]   # a study folder with design.qmd (a made-up design to change)
cd <study>
sapling design check                                 # the plan against itself (and against its last simulation)
sapling design simulate [--scenario PESS ...] [--n-sim N] [--seed K] [--jobs 7] [--compare ancova-locf,mmrm-all-data] [--curve | --no-curve]
sapling design simulate --refine [--n-sim 1000]      # add trials to the quick pass (the same seeds' sequence), do not start over
sapling design export [--docx]                       # the statistical section of the protocol
sapling design to-sap <folder> [--from design.qmd]   # a SAP that starts with the blocks of the design
sapling new <study> --from-design <folder>           # the same
node tools/design_proof.ts [--quick]                 # the proof against known answers: DESIGN PROOF PASS

Speed. sapling design simulate without --n-sim is a first look: 200 trials per scenario, printed with its Monte Carlo errors, and it says how to refine. --refine continues from it (trial i of scenario k always has the seed seed + 1e6 k + i, and the trials of a scenario are kept in sim/<scenario>.rds, so the 200 are not wasted; a refined run is the same as a direct run of the same N). The power curve is computed when asked (--curve) or when the stated power is not met (then only for the expected case, around the N it needs); --no-curve never. An MMRM fit takes 0.13 s (80% of it inside nlme::gls, the C likelihood; the twin's own overhead is small), so the MMRM example costs about 1850 CPU-s at 1000 trials with the curve of the made-up design, 4 to 5 minutes on 7 free cores. A compiled MMRM likelihood is a later slice. The ANCOVA, time-to-event and multiplicity examples take seconds.

The design plan#

design.qmd is a plan of the language of sap.qmd (docs/markdown.md) with kind: design under sapling:: it holds the statistical design only (the design block, the objectives and endpoints, the estimands with their intercurrent events and strategies, the analysis populations, the primary and key secondary analyses with their models, the sample size, the multiplicity, the interim analyses, the missing data and the sensitivity analyses). With kind: design a plan has no protocol to cite, and that is all it changes: it is the smallest change (a plan without protocol: already had no citation rules; the compiler warned that it could not check citations, and now says so as information for a design plan). Nothing is forked: sapling check, sapling review, the Method Lab and the editor read it as they read a SAP, and its blocks are the SAP's blocks.

The one new block is the assumptions block, with its parts effect, correlation, accrual, dropout, ice and scenario (reference: docs/markdown.md, "The design plan and the assumptions block"). Every number has a basis=, because the protocol must justify it. The examples: examples/design/htn-sbp/design.qmd (a repeated measure, MMRM, a treatment-policy primary estimand and a hypothetical supplementary one) and examples/design/onc-pfs/design.qmd (time to event, Cox).

What is simulated#

The operating characteristics of an analysis are computed against the true value of every estimand, its own first: the bias of ANCOVA with LOCF against the treatment-policy and against the hypothetical estimand are both in the report.

Multiplicity and group-sequential looks are simulated (below). Not simulated yet (said by design check as a note where the plan has it): non-inferiority and equivalence margins, multiple imputation (the planned sensitivity analyses: the comparison covers MMRM, ANCOVA, LOCF and on-treatment analyses), more than two groups, covariates other than the baseline and class factors, informative censoring of a time to event, strategies other than treatment policy and hypothetical for a continuous endpoint (composite for a response only), non-proportional hazards, and the planned Kenward-Roger degrees of freedom (see below).

Multiplicity#

When the plan has a multiplicity block with two or more hypotheses and the assumptions give each hypothesis its effect ([...]{.effect hypothesis="H1" difference="-5" sd="12" basis="..."}) and the correlation of the endpoints ([...]{.correlation between="endpoints" rho="0.6" basis="..."}), sapling design simulate also runs the family (programs/simulate_multiplicity.R, r/design/multiplicity.R): K endpoints on the same subjects (normal, equal correlation), each tested with its planned analysis (the ANCOVA or t-test of its analysis block, through lsmeans()/ttest() of r/sapling/), then the plan's procedure decides (hierarchy = fixed sequence in the order of the family, bonferroni with the weights, holm, hochberg, graphical with the plan's weights and transitions (Bretz et al. 2009), co-primary: each at alpha, all must be rejected). Per scenario: the probability of rejecting each hypothesis, all, at least one, and the family-wise error (a true null rejected). The scenarios: the expected effects, the global null, each hypothesis alone null (a partial null), and the plan's own scenarios with zero="H2" (the hypotheses whose effect is zero in them). The family runs 5 x --n-sim trials (1000 to 5000) because a t-test is cheap. The main simulation is then of the primary analysis only; the other hypotheses' analyses are in the family program (fns). Direction and sidedness: an estimate in the wrong direction never rejects (the direction of benefit is the sign of the hypothesis' expected effect); with sided=2 the two-sided p-value counts in the right direction only, so under the global null the chance of a rejection in the right direction is alpha/2 (a rejection in the wrong direction, a type III error, is not counted as a success nor as a false rejection); with sided=1 the one-sided p-value is p/2 in the right direction. Only single-visit ANCOVA/t-test endpoints are in a family so far (an MMRM hypothesis is a problem with its fix). design.fwer-above-alpha is an error. The example: examples/design/htn-two/design.qmd.

Interim analyses (group sequential)#

When the plan has interim blocks for the primary hypothesis (one block per look: information=, events=, purpose="efficacy and futility", spending= an alpha spending function (obrien-fleming, pocock, hwang-shih-decani with gamma=), futility_z= a non-binding futility boundary, z= or alpha= a stated boundary, the language of docs/markdown.md), sapling design simulate also runs the looks (programs/simulate_interim.R, r/design/interim.R, r/design/gsd.R). The efficacy boundaries are computed in base R from the spending function by the recursive numerical integration of Armitage, McPherson and Rowe, at the one-sided level (a two-sided 5% test: one-sided 2.5%). At each look the planned analysis of r/sapling/ runs on the data so far and its z (the sign of the benefit times qnorm(1 - p/2)) is compared with the boundaries. For an endpoint measured once, a look is at the first ceiling(t N) subjects in the order of entry (the outcome is available at once: no delay); for a time to event, a look is at the calendar time of the ceiling(t D)th event, the final analysis too (D: the events of the first interim block divided by its information, or Schoenfeld's). Per scenario: the probability of stopping for efficacy and for futility at each look, of ending without a rejection, the overall power (the type I error under the null), and the expected number of subjects, events and months. Not simulated: an MMRM over several visits (a problem with its fix), binding futility, a look with a delay of the outcome, sample size re-estimation.

The program: adjust it by hand#

out/design/programs/simulate_<analysis>.R is plain R. The assumptions are the list A (the same numbers as design.qmd), the analyses are functions of r/sapling/ written out as the generated programs write them (mmrm(), lsmeans(), cox()), the run is one call of design_run(). Open it, change a number or a model, run it again:

Rscript --vanilla out/design/programs/simulate_P_SBP.R

The generator (what a trial looks like) is r/design/generate.R, the runner r/design/run.R. A change in the program is a change of the design: make it in design.qmd too (the program carries the hash of the assumptions it was made from). Each trial of scenario k has the seed seed + 1e6 k + i: the results are the same for any number of workers. One example trial is written as data files (example/adsl.csv, advs.csv or adtte.csv, adtte_ontrt.csv, adrs.csv) to look at.

What differs from the study's program: the simulation uses ddfm = "residual" for the MMRM (Kenward-Roger in the R twin takes 24 s a fit; --mmrm-df kr asks for it) and the structures of the R twin (un, csh, cs, arh1, ar1; a planned toeph fallback is dropped and noted). Both change the p-value by a little at N of 100 or more.

The check and the proposals#

sapling design check reads the plan against itself and, when out/design/results.json is of the same assumptions (same hash), against the simulation:

Each finding is data (src/design_findings.ts): {id, severity, block, title, found, stated, fix, change?, evidence?}. change is the suggested change as the new text of the block (set an attribute, replace a block, insert, or a decision in words). The findings never write to the plan. The proposal workflow (src/proposals.ts: the AI proposes, a person accepts) needs one small hook, described in the report of the build: a function that turns a finding with a change into a Proposal of kind replace (target: the block with change.block, text: change.text, expect: change.was, reason: the finding's found and fix), compiled in place by checkBlock as every proposal is, so that "Accept" is the statistician's click.

The export: the statistical section of the protocol#

sapling design export [--docx] writes protocol-statistics.md (and .docx): plain protocol prose made deterministically from the blocks, with no AI. After each paragraph an HTML comment names the blocks it comes from (<!-- from: A_BASE, SS_SBP -->), and protocol-statistics.trace.json lists them. The sections:

1. Objectives, endpoints and estimands (a table with the five attributes of ICH E9(R1))
2. Statistical hypotheses
3. Sample size determination (the numbers, the assumptions with their basis, the method, the result, and a paragraph on the simulations with the operating characteristics under the scenarios)
4. Analysis sets
5. Statistical analyses: general considerations, primary, secondary, sensitivity, handling of intercurrent events and missing data
6. Multiplicity
7. Interim analyses

The headings follow the statistical considerations of the ICH M11 template as far as Sapling knows it. Check the headings and their numbers against your company's template before the text goes into a protocol; the order and the wording of a template are the sponsor's.

From design to SAP#

sapling design to-sap <folder> (or sapling new <folder> --from-design <design folder>) writes sap.qmd with the blocks of the design (kind: design removed, a note on top) and .sapling/design-origin.json (the hash of the design and its blocks). When the protocol exists: add protocol: to the front matter; sapling check then lists each block with no source, and each gets its citation as usual ([]{.cite source="protocol:<section>" words="..."}); a block the protocol does not say is a change of the statistician. The comparison of the design with the protocol's statistical section (a block whose number or wording differs from the sentence it cites) is the rule no-source/cite of the compiler plus a difference report from design-origin.json; the report of differences is not made yet.

The proof#

The rule of the proof. A check compares a simulated number with theory and passes when the difference is within 3.5 Monte Carlo standard errors plus a stated slack for the theory's own approximation (printed with each check as its tolerance, and its reason when there is one). The proof makes about 150 comparisons, so a check at |z| of 2 to 3 is expected a few times by chance and one at |z| above 3.5 is not; a passing check at |z| above 2 is a reason to run it again with more trials and seeds of its own, not to accept it: the type I checks of the group-sequential designs therefore run 30000 trials of the null with seeds that no other check uses (the 0.0323 at |z| = 2.3 that a 3000-trial run showed was noise on its first seeds).

node tools/design_proof.ts runs the whole path (the language, the spec, the R program, the r/sapling/ analyses, the operating characteristics) on plans of known answers and compares what was simulated with what theory says, in Monte Carlo standard errors. A simulation that disagrees with theory is a bug to fix, not a finding. The verdict is DESIGN PROOF PASS; test/design.test.ts and test/design-r.test.ts are the unit tests.

The full run (23 checks, 690 s on 7 cores; --quick runs the same checks with fewer trials, about 3 minutes) on 2026-11-14:

Multiplicity (--only mult; 2 independent endpoints, n = 40 per group, the t-test as the known-answer analysis, 8000 trials per scenario, 15 checks, 414 s). A rejection counts in the right direction only, so every null probability below is halved: under a two-sided 5% test the chance of a rejection in the right direction is alpha/2. The exact values come from power.t.test.

No published graphical example with a known answer is in the proof: the two equivalences above are the answers I can state with certainty.

Group-sequential looks (--only gsd, 54 checks, 161 s). The boundaries are computed in base R by numerical integration (r/design/gsd.R); gsDesign and rpact are not installed here, so the published values are the proof:

The looks themselves are checked against the exact probabilities of the same boundaries (the recursion of gsd_probs() with the drift of the design: delta sqrt(n1 n2 / N) / sd for means, sqrt(events) |log HR| / 2 for a time to event): the probability of stopping for efficacy and for futility at each look under the null and under the stated effect, the overall type I error and the power, for a means design (with and without futility) and a time-to-event design with two interim looks (events-driven). All 54 are within 3.5 Monte Carlo standard errors plus a stated slack. Two questions the first proof run raised, and their answers (checked with 20000 to 40000 trials; the engine's own path and a plain-R replica agree):

1. Is the overall one-sided type I error of the means design (85 to 170 subjects, three looks) above 0.025? No. The 0.0323 at 3000 trials and 0.0298 at 6000 were the same first seeds: Monte Carlo noise. With 24000 further trials the engine gives 0.0248 (sequential seeds) and 0.0255 (independent streams, 24000), a plain-R replica 0.0246 (40000), SE 0.0008 to 0.001. The boundaries are applied to the p-value scale (the z of a look is qnorm(1 - p/2) of the two-sided p-value of the planned t-test), so the small-sample distribution of the statistic is respected and each look's marginal level is exact; boundaries of normal theory applied to the t statistic itself would be the version that inflates the type I error, and this is not what the engine does.
2. Is the early look's stopping probability (0.232 against 0.255) a fault? It is real, not noise (20000 trials: 0.226 against 0.255), and it is the other side of the same fact. At the first look (85 subjects, 83 degrees of freedom) a t statistic reaches the very small p-value of the O'Brien-Fleming boundary (p = 0.003) less often than a normal statistic at the same effect does: in a plain replica 0.234 for the t-test's p-value, 0.252 for a normal statistic on the same random subsets (arm imbalance in the subset costs almost nothing), theory 0.255. The overall power is hardly touched (0.891 against 0.898). A design with a larger n per look is exact (below: 1000 subjects, 500 and 750 at the looks). The simulation says it when it happens (design.interim-small-sample: "the boundaries are computed for a normal statistic ... the early looks stop less often than the normal theory says") and its numbers, not the normal ones, are the ones to size the study with.

The information actually observed (the default). The boundary of a look of a spending-function design is computed from the information observed (Lan and DeMets 1983): for a time to event the events at the look over the planned events (the last look too, when follow-up ends with fewer events than planned: it spends what is left of alpha); for an endpoint measured once, the subjects with a value over the planned analysable number (N x (1 - dropout)). --information planned uses the nominal fractions. Proof: a time to event with follow-up that ends at 237 of the 330 planned events: with the information observed the type I error is 0.0250 (SE 0.0035, alpha is 0.025); with the planned fractions it is 0.0240, also at or below alpha (the difference is too small to resolve with this many trials: reported, not asserted).

The example onc-pfs is a two-sided 5% design, so its efficacy boundaries are one-sided 2.5%; with non-binding futility boundaries the exact probability of rejecting under the null is 2.47%. The simulation gave 3.1% at 1000 trials (SE 0.5%) and 2.1% at 12000 trials (SE 0.14%). That is Monte Carlo variation around 2.47%, not a shortfall of the Cox Wald test or of the loss to follow-up: the same design with no loss to follow-up gives 0.0239 with the Wald test and 0.0243 with the log-rank test, and with 5% loss and the log-rank test 0.0267 (10000 trials each, SE 0.0016).

(z: the difference in Monte Carlo standard errors; a check passes within 3.5 standard errors plus the stated slack of the theory's own approximation. An earlier quick run showed MNAR bias -0.373 against -0.6 at 1.55 standard errors with 200 trials; with 600 it is -0.609: Monte Carlo noise, not a bug.)

What a statistician should know before trusting it#

• The simulation is as good as its assumptions: it shows what the design implies if they are true. It does not check the assumptions against data.
• The SD of effect is the SD of the change after the baseline (what ANCOVA and MMRM power use); with a baseline_coef the marginal SD of the change is larger.
• Intercurrent events and dropout are independent of the outcome except where the assumptions say (MAR depends on the last value). An event that depends on the outcome, or a gradual loss of the effect after it, is not in this version.
• The strategy of an estimand decides its true value; the records of the analysis decide what the analysis does. When they disagree the report shows the bias: that is the finding, not a bug.
• 1000 trials give a Monte Carlo standard error of 1.0 to 1.6 percentage points on a power: look at it before reading a power of 88% against 90%.
• The power curve's smallest N is a probit fit of the simulated powers: its error is a few subjects.

The shape of out/design/results.json#

Stable, because the editor reads it (src/design.ts: DesignResults; the exports of src/design.ts and src/design_findings.ts keep their names and shapes). schema: "sapling.design.results/1"; new fields are added, none is renamed or removed.

{ schema, dir, program,
  spec:      DesignSpec    { plan, study, assumptions, hash, kind ("continuous"|"tte"|"binary"), arms[2], alloc[2], n, alpha, sided, fav,
                             stated {n, power, effect, sd, dropout, block}, paramcd, visits[{label, n}], A (the assumptions as R reads them),
                             analyses[{id, estimand, title, test, planned, where, onTreatment, formula, mmrm?, covars, classes, response}],
                             estimands[{id, title, level, analysis, strategies{event: strategy}}], scenarios[{id, label, over, curve, basis}],
                             assumptionBlock, multiplicity?, interim?, notes[], problems[] },
  settings:  DesignSettings { nSim, seed, jobs, compare[], curve, nSimCurve, mmrmDf, ... },
  run:       { seconds, n_sim, seed, jobs, r, packages{}, assumptions_hash },
  scenarios: [ { scenario, label, null, params, n_sim, seed, seconds, truths{estimand: value}, facts{name: [mean...]},
                 analyses{ analysisId: { own: estimandId, vs{ estimandId: DesignOc } } } } ],
  curve:     [ { scenario, label, points[{n, power, se, failed}], target, n_fit, n_grid, n_sim } ],
  multiplicity?: { scenarioId: { scenario, label, n_sim, true_null[], effects{}, reject{h: p}, reject_se{}, all, all_se, any, any_se, fwer, fwer_se } },
  interim?:  { bounds {t[], upper[], lower[] (null: none), D}, scenarios{ scenarioId: { scenario, label, null, n_sim, p_eff[], p_fut[], p_eff_se[], p_fut_se[],
                                                                                    p_none, power, power_se, exp_n, exp_events, exp_time, theory } } },
  methods, files[] }
DesignOc = { n_sim, n_ok, failed, fallback, rejection, rejection_se, truth, mean_est, sd_est, q[5], bias, bias_se, rmse, rmse_se,
             coverage, coverage_se, mean_se, mean_ci_width }

out/design/findings.json is the list of DesignFinding (id, severity, block, title, found, stated, fix, change?, evidence?).

One harness: the Method Lab and the design simulation#

The Method Lab (docs/method-lab.md) tests one analysis of a plan: its type I error and power at the plan's single assumption. The design simulation extends it and does not sit beside it: it reads the analysis through the lab's spec (labSpec() of src/lab/spec.ts: the test, the model, the records, the visits, the covariance, the estimand), and its analyses are the functions of the R twin (r/sapling/) that the generated programs call. What it adds is the assumptions block (scenarios, dropout mechanisms, intercurrent events with their effect), the true value of each estimand, the comparison of analyses against it, the power curve over N and the report. The lab's own simulation (r/lab/lab.R) keeps its own code for a plan that has data; moving it onto r/design/ is a next step.