Standard-Curve Simulation Study — Master Strategy & Operations Guide

Purpose of this document. Single source of truth for the simulation-recovery study that validates the hierarchical Bayesian calibration framework for the Statistics in Medicine manuscript. Consolidates the strategy, design decisions, file-by-file build, and operational learnings so the project can proceed without re-deriving anything from the (very long) development thread.

Status at time of writing. Anchor regeneration complete (Jobs 1–3). Anchor extract stamped submission-ready. Next step: truth generation (01_generate_truth.R) and the staged simulation run.

---

1. The manuscript this serves

Target journal: Statistics in Medicine (methodology motivated by medical application; expects a simulation with known truth).

Framing (decided): Framing A — "the plate is the wrong inferential unit." Three conventions of routine standard-curve practice (per-plate fitting, per-fit form selection, deletion of inconvenient calibrators) are all consequences of treating the plate as the inferential unit; each is individually defensible and collectively expensive. Framing B (the cost of deleting data — the masking audit) is the sharpest evidence inside Framing A. Framing C (model-form uncertainty as a component of measurement uncertainty; the stacking-mixed CDAN estimator) is signposted only, as future work — no equation in the paper.

Why the simulation exists — the one adversarial reader. The referee who says "shrinkage always lowers variance; you've proven nothing about accuracy." Section 4 (the simulation) exists to answer that and nothing else. Design principle throughout: known truth is the arbiter, and the hierarchical estimator must be given every chance to lose.

Primary recovery metric: RMSE on back-calculated log10 concentration. Parameter RMSE is secondary.

---

2. Simulation design (the science)

2.1 Three tiers of "truth" — the inverse-crime firewall

The single most important design choice. If truth is drawn from exactly the model being fit, the hierarchical fit wins by construction and the result is worthless. Three tiers, and the headline claim rests on Tier 2:

Tier	Truth drawn from	Purpose	Expectation
T1 Concordant	the fitted hierarchical model exactly	best case; verifies sampler recovers its own generative model; SBC lives here	Bayes must win, or something is broken
T2 Plausible-perturbed	population params + correlated cross-plate deviations, log-normal (not normal) parameter noise, per-plate γ	the honest test — resembles real assays but violates the fitted model's independence/distributional assumptions	Bayes should still win on RMSE; margin narrower
T3 Misspecified form	a monotone spline outside the fitted 4PL/5PL/Gompertz family	adversarial; neither engine's forms are correct	open — if Bayes loses here, that is a finding and bounds the claim

If compute forces a cut, cut T1 (diagnostic, not evidence). Never cut T3.

2.2 T3 spline specification (decided: "Spline-B")

• Take the best-fitting parametric curve per antigen; evaluate densely.
• Add a smooth monotone (I-spline / monotone-cubic) perturbation; asymptotes preserved.
• Calibrate the perturbation magnitude so the spline's max deviation from its nearest-fitting 5PL equals the anchor's global median per-plate residual RMSE (lack_of_fit_global). This ties "how misspecified" to a measured, defensible quantity rather than an arbitrary knob.
• Convention: the global lack-of-fit median includes pathological antigens (T3 is the adversarial tier; real corpora contain pathological curves).
• Upgrade path: switch to option "C" (real reference curves) when the Sonia Macallister control panel data arrives.

2.3 Scenario grid — 54 cells

dynamic_range (3: low/mid/high) × asymmetry_g (3: symmetric/moderate/danger_zone [0.2,0.6]) × noise (2: assay_typical/elevated) × tier (3: T1/T2/T3) = 54 cells.

• Dynamic range is the spine — §5.2 claims the advantage grows with it; the simulation must confirm-or-contradict against truth.
• N_sim = 1000 replicates per cell. Chosen so a 95% coverage estimate has Monte Carlo SE ≈ 0.69% (distinguishes 95% from ~93.5% cleanly).
• Masking is a paired within-replicate toggle, NOT a grid axis: same truth, refit with and without the 3×-blank rule. Paired differences → tight intervals.

2.4 Metrics (computed against stored truth)

• Primary: back-calc RMSE on log10 concentration, stratified by curve zone (near_lower / mid / near_upper — the U-shape).
• Parameter-level (secondary): bias, RMSE, variance/bias decomposition (so a referee can see any RMSE win is not variance-only — the exact objection the section answers).
• Coverage: interval coverage at 50/80/95%; concentration interval coverage.
• Recovery: recovery rate and RMSE-on-recovered, reported together (recovering a sample wrongly is not a win — the Feng 2011 phenomenon under known truth).
• LOQ calibration: of samples declared within the quantifiable window, what fraction have true concentration recovered within the CV budget.
• Masking (paired): Δ(RMSE), Δ(recovery), Δ(fit-failure), masked − unmasked.
• SBC (T1 only): rank uniformity per parameter.

2.5 Staged execution — fail cheaply

Stage	Cells	N_sim	Gate
0 smoke	1 cell	5	plumbing only
1 sanity	mid-DR, all g, T1	50	Bayes MUST win on concordant truth + SBC ranks ~ uniform (incl. b histogram → sigma_log_b decision)
2 honest	DR × g, T2	100	inspect T2 recovery RMSE → decides whether the full grid runs. If Bayes wins only marginally/loses, replan §5 as precision-not-accuracy paper
3 full	all 54	1000	includes T3; days of runtime; resumable

Stage 2 is the decision gate that can reframe the paper.

2.6 Anticipation table (where results will land)

Result	If claims hold	If contradicted
T2 recovery RMSE	Bayes < NLS, margin grows with DR	flat/reversed → §5.2 is a shrinkage artifact; paper becomes precision+limits
Variance/bias split	Bayes win survives removing variance component	win is variance-only → referee is right; reframe
Coverage (SBC + intervals)	≈ nominal	undercoverage → priors need work before submission
T3 misspecified	Bayes degrades gracefully, still ≤ NLS	Bayes worse → bound claim to "when a reasonable form set is available"
Masking Δ(RMSE), paired	masking inflates RMSE, drops correct recoveries	harmless under truth → §5.5 loses its known-truth backbone
g recovery, danger zone	g recovered, wide but honest intervals	unrecoverable → C3 stays parked

---

3. The anchor (foundation of everything)

3.1 What the anchor is

One real-data hierarchical fit that serves double duty: (1) the §5 application fit (paper's headline numbers), and (2) the generative anchor supplying population parameters, cross-plate SDs, the 5×5 parameter correlation matrix, the variance function, the lack-of-fit yardstick, and dynamic-range coverage for all 54 simulation cells. Making these the same fit gives internal consistency by construction.

3.2 Anchor data scope

• Study: cloneB_GaPs (a clone of MADI P3 GAPS), project_id 90.
• 10 experiments / readout classes: ADCD, ADCP, ADNP, FcgR2aH, FcgR3bNA1, IgG_AW, IgG1, IgG2, IgG3, IgG4.
• 96 (experiment × antigen) groups, 15 plates each, 1,440 curve_ids total. (ADCP/ADNP have only 4 pertussis antigens each → 60 curve_ids; all others 11 antigens → 165.)

3.3 Source anchoring map (decided)

Each (experiment, antigen) is anchored to ONE standard source:

Readout class	Anchor source
IgG_AW, IgG1–4	NIBSC06_140
FcgR2aH, FcgR3bNA1, ADCD	NIBSC06_140
ADCP, ADNP	Sando (Sandoglobulin)

Source strings are normalized (dirty → canonical): anything starting NIBSC* → NIBSC06_140; SANDO* → Sando; SD/STD/STANDARD* → SD. Unrecognized tokens halt (never silently bucketed). In cloneB the column was already clean (one raw spelling per bucket), but the normalizer is retained for future clones.

3.4 Full-Bayes prior decision

Moved from the deployed empirical-Bayes priors to full-Bayes, design-standardized priors. These are R-only changes (data-block hyperparameters), not Stan-model edits — confirmed by reading the .stan files (priors are all passed as data; no hard d ≥ max(y) constraint in the model). Specifically:

• Prior locations anchored to the known calibrator design, not observed y (removes empirical-Bayes peeking).
• Hard d ≥ max(y) pull dropped (prior_d_mu set freely).
• γ (variance-function exponent) pooled at the antigen level.

What the deployed Stan models actually implement (not the idealized spec). Cross-plate SD priors are partially hardcoded:

• sigma_a/d/c ~ normal(0, prior_*_sigma * 0.5) — derived, not free.
• sigma_log_b ~ normal(0, 0.5) — fixed constant, ignores passed data.
• sigma_log_g ~ normal(0, prior_log_g_plate_sd) — data-driven.
• sigma_obs, sigma_blank ~ normal(0, prior_a_sigma) — coupled to a's width.

The sigma_log_b fixed-constant decision (DEFERRED). The Hill slope is the paper's headline parameter (17× reproducibility). Its cross-plate pooling being a fixed constant is the one scientifically questionable constraint. Decision deferred pending the T1 SBC rank histogram for b — if b's ranks are non-uniform, implement a free_crossplate_sd switch in the Stan models (a recompile + worker-image rebuild). If uniform, document as a known property. free_crossplate_sd = 0 (deployed default) for now.

3.5 Three gate-zero fixes baked into the anchor

The anchor manifest stamps fixes = c("stacking_weights_api", "theta_base_floor", "full_bayes_priors"); truth generation refuses to proceed (assert_anchor_complete()) unless all three are present:

1. loo::stacking_weights() API fix (genuine weights, not collapsed to hard selection).
2. theta_base variance-floor fix (σ > 0 at background).
3. Full-Bayes priors (only stamped AFTER SBC passes).

---

4. The three-job anchor pipeline

Hard-separated, with a gating stop between Job 2 and Job 3.

JOB 1 — refit_anchor() — fit real data via the SIM API

• Submits one job per experiment (10 sequential) to the i-spi-compute-sim stack; sampling = n_draws_predict = 3500 (LOQ-damping; both must match or the local-vs-production LOQ comparison is confounded); cdan_cv_threshold = 20 (integer %); blank_option = "ignored"; full-Bayes prior overrides.
• Verifies 1,440 curve_ids present in calib_fit; writes anchor_manifest.rds (NOT stamped).

JOB 2 — sbc_on_anchor() — LOCAL curveRbayes SBC (n_sbc = 300)

• Does not touch the API or DB. Draws true params from the deployed model's actual priors (sigma_log_b from half_normal(0, 0.5) per Stan hardcode), simulates, fits via fit_bayes_single, records rank of truth in posterior.
• Gates Job 3 via evaluate_sbc() / PASS_RULE (not advisory).
• Reports the sigma_log_b finding for the deferred Stan-switch decision.

JOB 3 — extract_anchor() — build the generative anchor

• Only runs on a PASS verdict. Reads calib_* for population means, cross-plate SDs, variance function, lack-of-fit yardstick, form selection, dynamic range.
• Captures joint draws for the 5×5 correlation matrix (Convention 1: from the stacked posterior so g is always present) via one local fit_bayes_single on a representative antigen (IgG1/pt). Falls back to identity + warning if unavailable (→ T2 perturbation would be independent; re-run extract to fix).
• Calls assert_anchor_complete(), writes anchor_extract.rds, stamps full_bayes_priors = TRUE.
• Prints the anchor-inspection banner — the human gate before truth generation.

The SBC PASS_RULE (the gate that protects the whole compute budget)

• Primary test: ECDF-difference band at 95% (two-sided KS statistic, handles ties; epsilon self-adjusts to n_sbc). KS scalar stored, not gating.
• Per-parameter, not global (failure modes differ).
• Carve-out (narrow, deliberate): d_upper may fail in the tails and still pass (expected after dropping the d ≥ max(y) pull), but an interior violation (15–85% quantile range) still HALTS. Empty-interior = FAIL (a distribution with no interior is severely non-uniform).
• HALT params: a_lower, b_slope, c_ec50, g_asymmetry, gamma_variance, sigma_log_b, theta_base, theta_prop.
• Fixed-constant note attached to b_slope and sigma_log_b failures → evidence for the free_crossplate_sd switch.
• Do NOT retighten priors to suppress a failure — that reintroduces the empirical-Bayes undercoverage. Loosen the offending scale and re-run.
• evaluate_sbc() is unit-tested (test_evaluate_sbc.R, 22 checks, 6 synthetic scenarios) independently of any Stan fit.

---

5. Infrastructure

5.1 The dedicated SIM stack

• Separate clone: i-spi-compute-sim (own API + Redis + 16-core worker, replicas: 1, WORKER_CORES=16, 16Gi). Isolated queue — no collision with production i-spi-compute or the batch calculator.
• Base URL: https://madi-preprod.dartmouth.edu/i-spi-compute-sim (ROOT_PATH /i-spi-compute-sim; Traefik i-spi-compute-sim-stripprefix@file middleware confirmed present).
• Reuses the existing i-spi-compute SealedSecret (same API_KEY/REDIS_AUTH) and madi-lumi-reader DB secret. Same physical database as production.

5.2 Environment variables (session-side)