Decisions 45-50: scenario tuning, NC truths, oracle own estimand, separate PVI phases, validation

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Validation, NCE analysis, per-replicate oracle for landmark/naive, RCT truths script, resumable experiment runner
2026-10-27 09:20:22 -04:00 · 2026-10-27 09:19:59 -04:00 · 2026-10-27 09:01:13 -04:00 · 2026-10-27 08:15:30 -04:00 · 2026-10-27 08:15:18 -04:00 · 2026-10-27 07:54:19 -04:00
28 changed files with 16737 additions and 0 deletions
@@ -0,0 +1,6 @@
+.Rproj.user
+.Rhistory
+.RData
+results/
+*.parquet
+.DS_Store
@@ -0,0 +1,346 @@
+**Why surgical trials fail, and what each target trial emulation answers instead: a known-truth simulation of early peripheral vascular intervention for claudication**
+
+*Use case: PVI within 6 months of a new claudication diagnosis, compared with initial conservative management; Medicare-like claims; 5-year major adverse limb events (MALE)*
+
+Analysis plan, draft v1.0, 2026-10-27. Every decision cited as "D#" is recorded, with its reason, in `docs/decisions.md`. The readings are in `docs/literature/`.
+
+# 1. The paper in one paragraph
+
+Whether to perform peripheral vascular intervention (PVI) early for claudication is one of the most common decisions in vascular surgery. It has never been tested in a randomized trial on outcomes that matter to the limb. The trials that do exist measured walking distance and quality of life in 56–240 patients, many recruited poorly, and in their conservative arms patients crossed over to intervention. Observational studies in Medicare suggest early PVI leads to more interventions and more chronic limb-threatening ischemia (CLTI). They use different designs, though (landmark, clone–censor–weight), and those designs answer different questions. We build one simulated Medicare-like population of claudicants and their physicians in which the true effect of early PVI is known by construction. In it we switch on, one at a time, the features that make surgical trials hard:
+
+- variation in surgeon skill;
+- surgeon practice style;
+- surgeons who also decide the outcome (reintervention, CLTI coding);
+- symptom-driven deferral and crossover;
+- unmeasured symptom severity.
+
+We then run naive, landmark, sequential-trial and clone–censor–weight (CCW) analyses, with four ways of handling the physician, alongside the RCT that could realistically be run. Each method is scored against the question it actually answers and against the question a surgeon is asking. We also test whether negative controls can tell when the answer is wrong.
+
+# 2. The gap
+
+## 2.1 What is already known
+
+- **Trials of claudication are small and short, and none measures hard outcomes.** CLEVER (n = 111; planned 252), ERASE (212), IRONIC (158), MIMIC (127), SUPER (240, stopped early, 29% crossover in year 1) and OBACT (56 randomized of 434 screened) all used walking or QoL endpoints. Where follow-up was long, the benefit faded by 2–5 years (ERASE, IRONIC). The Cochrane review pooled no hard limb outcomes (Fakhry 2018).
+
+- **Surgical RCTs fail for structural reasons.** 21% of surgical RCTs stop early, mostly for poor recruitment (Chapman 2014). Surgical trials are discontinued more often than medical trials (43% vs 27%; Rosenthal 2015). Comparisons of surgery with non-surgery recruit about twice as slowly (Cook 2008). Crossover can make ITT uninterpretable (SPORT; Weinstein 2006). The intervention's effect depends on the operator (Birkmeyer 2013), and trial surgeons are selected (Rothwell 2005; Cutteridge 2025).
+
+- **In timing trials, the "conservative" arm becomes delayed surgery.** Examples: hernia watchful waiting reached 68% crossover in the long term (Fitzgibbons 2013); 62% of the ADAM surveillance arm were repaired (Lederle 2002); 85% of CAESAR by 54 months (Cao 2011).
+
+- **Medicare evidence on early PVI.** 3.2% of new claudicants have early PVI. The median physician rate is 0%, and 5.6% of physicians are outliers (Hicks 2020). Early PVI is associated with late PVI (22.5% vs 3.6%; aHR 6.89) and with CLTI (16.4% vs 7.8%) (Sorber 2023). A CCW emulation found more CLTI and MALE (Dun, VESS 2026). The Sorber design is in effect a 6-month landmark analysis.
+
+- **Methods.** Immortal time (Suissa 2008), landmark (Anderson 1983; Dafni 2011), CCW (Hernán 2018; Maringe 2020; grace periods: Wanis 2024) and sequential trials (Hernán 2008; Keogh 2023) are each well described. Duchesneau 2022 compared landmark and CCW against a plasmode truth in cancer surgery. Unmeasured confounding at the cluster level has been studied outside TTE (Li 2020; Cafri 2020). So has bias amplification from conditioning on an instrument (Myers 2011).
+
+## 2.2 The gap this paper fills
+
+**No study has put landmark, sequential-trial and CCW analyses of one surgical timing question side by side, each scored against its own true estimand, in a world where the surgeon drives treatment, modifies its effect, and partly decides the outcome. The study should also run alongside the RCT that could realistically be done.** So surgeons reading a Medicare TTE cannot tell:
+
+1. which question the analysis answered;
+2. how far that is from the question they care about;
+3. which surgery-specific threats could still bias it;
+4. whether a clean negative control means anything.
+
+## 2.3 Why it matters
+
+- Claudication PVI is high-volume, increasingly done in office-based labs, and varies widely between physicians (Hicks 2020; Dun 2024). The question will be answered by claims analyses, not trials.
+
+- Journals now expect TTE (TARGET 2025). Readers in surgery need to know what each design estimates and where it breaks, in clinical language.
+
+# 3. Aims and hypotheses
+
+| **Aim** | **Question** | **Hypotheses** | **Primary output** |
+|---|---|---|---|
+| **Aim 1. What each method answers** | In a clean world (all confounders measured, one physician type), how far is each method's own estimand from the clinically relevant one (question gap), and is each method unbiased for its own estimand? | H1a: Each TTE method is unbiased for its own estimand. H1b: Landmark and sequential-trial ITT have large question gaps, while CCW-B1 has none by definition. H1c: The ITT of the realistic RCT is diluted towards the null by crossover, and is shifted by surgeon selection. | Fig 1 (base panel), Table 2 |
+| **Aim 2. Surgery-specific threats** | How does each threat change the estimation bias of each method × physician variant? | H2a: Cluster-robust SEs fix coverage, not bias. H2b: A physician random intercept in the weight models removes only part of the confounding by physician; within-physician analysis (FE) removes it when style affects the outcome directly, but amplifies residual patient confounding when style acts only on treatment (Z-bias). H2c: Unmeasured severity biases every observational method; no physician variant fixes it. H2d: Crossover driven by deferral hurts landmark and RCT-ITT most. | Fig 1 (threat panels), Fig 2, Table 3 |
+| **Aim 3. Do negative controls detect the problem?** | When a method's MALE estimate is biased by more than 1 percentage point, how often does each negative control flag it? | H3a: NCO-1 (patient-level) is sensitive to unmeasured severity or frailty. H3b: NCO-2 (physician-level) flags style variation even when style is not confounding, so its specificity is low. H3c: Substantial bias can sit behind a clean negative control. | Fig 3, Table 4 |
+
+*Secondary (supplement): each MALE component separately; total procedures per patient; Sorber "late PVI" definition; HR scale; MALE-free survival; 3-month window; sequential trials per-protocol; preference IV; unmeasured supervised exercise; n = 200,000; full factorial if the one-at-a-time results suggest interactions.*
+
+# 4. Target trial specification
+
+| **Component** | **Target trial** | **Emulation in the simulated claims** |
+|---|---|---|
+| Eligibility | Medicare beneficiaries ≥66 y with a new claudication diagnosis, no PVI, bypass, CLTI or amputation in the prior 12 months | Same, at the diagnosis month (month 0) |
+| Strategies | **A:** PVI within 6 months of diagnosis. **B1 (primary comparator):** no PVI in months 0–5, then usual care. **B2:** no elective PVI at any time (PVI only after a MALE event) | Observed PVI timing compared with each strategy |
+| Assignment | Randomized at diagnosis | Exchangeable conditional on measured baseline and time-varying covariates (and, by variant, on the physician) |
+| Time zero | Diagnosis = eligibility = assignment | Same |
+| Follow-up | Months 0–59. Ends at MALE, death or month 60 | Same; no loss to follow-up (as in Medicare FFS without disenrollment) |
+| Outcome | MALE: first of CLTI, major amputation, open bypass, or any revascularization after the first (index) one (D7). Death is a competing event (D8) | Same, from simulated claims |
+| Causal contrast | Per-protocol (sustained strategies) | CCW with artificial-censoring weights |
+| Measure | 5-year cumulative incidence of MALE; risk difference A − B1 (and A − B2) | Weighted cause-specific hazards → cumulative incidence |
+
+**Clinically relevant estimand (the reference, D11):** RD(A − B1) at month 60, among all new claudicants, with PVI done by each patient's own physician (D21), under full adherence.
+
+# 5. Data-generating mechanism
+
+## 5.1 Structure
+
+Time is in months 0–59 (D16). Each replicate draws **J ≈ 1,500 physicians** and **N = 50,000 patients** (D25). Patients are assigned to physicians with lognormal panel sizes (median ≈ 20, long right tail).
+
+**Physicians (j).**
+- Latent **skill** K_j ~ N(0, σ_K²) and latent **style** S_j ~ N(0, σ_S²), with corr(K, S) = ρ.
+- **Practice setting** OBL_j ∈ {0,1} is observed, with logit P(OBL_j = 1) = o₀ + o_S·S_j (high-style physicians are more often in office-based labs or ASCs).
+
+**Patients (i), at diagnosis.**
+- Measured X_i: age, female, Black race, diabetes, CKD/dialysis, CHF, COPD, prior CVD, coded smoking, calendar year.
+- Latent **frailty** F_i, correlated with the measured comorbidities.
+- Latent **true smoking** (coded smoking is a misclassified version of it).
+- Latent **initial severity** Z_i0 = γ₀ + γ_X·X_i + e_i.
+- Latent **progression rate** r_i = r₀ + r_X·(diabetes, smoking) + u_i.
+
+**Latent severity process.** Z_it = Z_i,t−1 + r_i + ε_it, so severity drifts upward at a patient-specific rate.
+
+**What a PVI does (the treatment mechanism, D20).** At the month of PVI:
+- **Relief.** Severity drops by Δ (symptom relief, which lowers the near-term CLTI hazard).
+- **Restenosis and progression.** After PVI, the progression rate increases by κ·g(K_j), where g decreases with skill. This adds a skill-dependent hazard of repeat revascularization and of later CLTI.
+- Relief early and faster progression later produce the *benefit-then-harm* pattern (D4, factor 7). Δ and κ set the truth scenario (§5.3).
+
+**PVI decision, before the first PVI.** The monthly hazard is
+
+logit h_PVI(t) = α₀ + α_Z·Z_it + α_X·X_i + α_F·F_i + α_S·S_j + α_win·1[t ≥ 6]·(β_def + β_Z,def·Z_it)
+
+- Inside the window (months 0–5) the decision is driven by severity, which is confounding by indication.
+- After month 6, deferred PVI is governed by β_def and β_Z,def. This is factor 4.
+
+**Outcomes (monthly cause-specific hazards; pooled logistic form).**
+
+| **Process** | **Depends on** |
+|---|---|
+| True CLTI | Z_it, X_i, F_i |
+| Coded CLTI | True CLTI + extra coding hazard ∝ S_j·1[factor 3] among patients with moderate Z_it |
+| Major amputation | Mainly after CLTI; small direct hazard from Z_it, diabetes, CKD |
+| Open bypass (elective) | Z_it; higher if a previous PVI failed |
+| Repeat revascularization (only after the index PVI) | Restenosis process (κ, K_j); lower threshold ∝ S_j·1[factor 3] |
+| Death | Age, X_i, F_i, Z_it (weak) |
+| NCO-1: hip-fracture or pneumonia admission | F_i, true smoking, age. **Never** PVI |
+| NCO-2: varicose-vein ablation | S_j, latent venous disease. **Never** PVI |
+| NCE: low-value procedure (vein ablation) in months 0–5 | S_j, F_i. **Never** affects any outcome |
+
+**Measured vs latent (D17).**
+- What the analyst sees: X_i, a claims frailty index, the physician ID, OBL_j, all coded events, and a time-varying **severity proxy**. The proxy is a monthly count of vascular visits and noninvasive arterial tests, P_it ~ Poisson(exp(p₀ + λ·Z_it)), plus its running sum.
+- What is latent: Z_it, r_i, F_i (the frailty index is its proxy), true smoking, K_j, S_j.
+
+**Table A. Variable roles**
+
+| **Block** | **Causal role** | **Clinical analogue** | **Measured?** |
+|---|---|---|---|
+| X | Baseline confounders | Age, sex, race, diabetes, CKD, CHF, COPD, prior CVD, coded smoking, year | Yes |
+| F | Patient confounder (death, PVI candidacy, NCO-1) | Frailty | Through a frailty-index proxy |
+| Z_t | Time-varying confounder (PVI decision, CLTI) | Walking limitation / ABI | Through a visit/testing proxy, P_t (strength λ) |
+| r | Progression rate (prognostic) | Diabetes/smoking-driven disease progression | No |
+| K | Physician effect modifier (and confounder if ρ ≠ 0) | Technical skill, durability of PVI | No (physician ID only) |
+| S | Physician treatment driver; with factor 3 on, also an outcome driver | Practice style / "high user" | No (physician ID, OBL) |
+| NCO-1, NCO-2, NCE | Negative controls (true effect exactly 0) | Hip fracture/pneumonia; vein ablation | Yes |
+| E (supplement only) | Unmeasured co-intervention | Supervised exercise uptake | No |
+
+## 5.2 Worlds (one factor at a time, D24)
+
+Every world is run under each of the three truth scenarios in §5.3.
+
+| **World** | **Skill (1)** | **Style (2)** | **Style → outcome (3)** | **Deferral/crossover (4)** | **Unmeasured severity (5)** |
+|---|---|---|---|---|---|
+| **W0 Base (clean)** | σ_K = 0 | Modest σ_S | Off | Low, unrelated to severity | Off: λ high, proxy ≈ Z |
+| W1 Skill | σ_K > 0, ρ = −0.3 (and ρ = 0) | as W0 | Off | as W0 | Off |
+| W2 Style | as W0 | Strong σ_S, calibrated to the Hicks 2020 distribution | Off | as W0 | Off |
+| W3 Style → outcome | as W0 | Strong | **On** | as W0 | Off |
+| W4 Deferral | as W0 | as W0 | Off | **High, severity-driven after month 6** | Off |
+| W5 Unmeasured severity | as W0 | as W0 | Off | as W0 | **On: λ low (proxy–Z correlation ≈ 0.4), frailty index noisy** |
+| **WR Realistic** | On | Strong | On | High | On |
+
+W1 is run at both ρ = −0.3 and ρ = 0 (D22). WR is calibrated to the Medicare targets in §5.4.
+
+## 5.3 Truth scenarios
+
+| **Scenario** | **Mechanism** | **Target** |
+|---|---|---|
+| **Null (sharp disease null, D34)** | Δ = 0, κ = 0: PVI changes nothing about severity, CLTI, amputation, bypass or death | RD(A − B1) = 0 **exactly** for the disease composite (CLTI or major amputation) and every disease component. MALE is still > 0, because only patients who have had a PVI can have a repeat revascularization: this is the mechanical effect of how the outcome is defined, and it is reported as a finding |
+| **Harm** | Δ small, κ moderate | RD(A − B1) ≈ +5 percentage points |
+| **Benefit then harm** | Δ large, κ moderate | Cumulative incidence curves cross: A lower before about month 18, higher at 60 |
+
+## 5.4 Calibration targets for the usual-care world (WR)
+
+| **Target** | **Value** | **Source** |
+|---|---|---|
+| Early PVI (≤6 months) | 3.2% | Hicks 2020 (31208940); Sorber 2023 (37276171) |
+| Physician early-PVI rate | median 0%, range 0–58%, about 5.6% of physicians ≥2 SD outliers | Hicks 2020 |
+| OBL/ASC share of PVI | about 50% | Dun 2024 (38958132) |
+| Late PVI (>6 months, ~4.4 years) | 22.5% after early PVI vs 3.6% without | Sorber 2023 |
+| CLTI (~4.4 years) | 16.4% vs 7.8% | Sorber 2023 |
+| Mean age | about 76 | Dun 2025 (40782724) |
+| All-cause mortality | Interim target about 35% at 5 years (Swedish population cohort: 10-year mortality in IC 63% vs 27% reference). **Replace with the 5-year mortality of the authors' own Medicare claudication cohort** (Sorber 2023 / VESS 2026 data) | Sigvant 2018 *EJVES* (29478910); cardiovascular-only 5-year mortality 14–16%: Baubeta Fridh 2018 *BJS* (29116654) |
+| Amputation | about 0.4%/year after revascularization for IC | Baubeta Fridh 2017 (28797662) |
+
+Calibration solves for the intercepts and the key coefficients (α₀, α_S, σ_S, β_def, the outcome intercepts) by stochastic approximation on 10⁵ patients, as in `tte-papers/simulations/01_calibrate.R`.
+
+## 5.5 Computing the truths
+
+Truths come from simulating 10⁶ patients, reusing the physicians, under each intervention. The latent processes stay the same; only the PVI decision is overridden.
+
+| **Estimand** | **Intervention used for the truth** |
+|---|---|
+| A | Natural PVI hazard in months 0–4. Everyone still untreated, alive and MALE-free gets PVI at month 5 (the "natural grace" strategy that CCW with end-of-window censoring targets; Wanis 2024) |
+| B1 | PVI hazard = 0 in months 0–5, then natural |
+| B2 | PVI hazard = 0 at all times before MALE |
+| Sequential-trial ITT | For each m = 0–5, take the natural-course population eligible at m. Force PVI at m vs no PVI at m (natural course afterwards in both). Pool over m, weighted by the number eligible |
+| Landmark | No interventional definition, because the comparison is among survivors of the observed regime. Its own estimand is the **oracle large-sample limit**: the landmark IPW estimator on 10⁷ patients, with latent Z and F included in the weight model |
+| Naive | No own estimand. Only total error against the reference is reported |
+| Realistic RCT ITT | Simulate the trial population (hospital-based, top-third skill, low-style physicians; eligible patients), with assignment to A or B under natural adherence (refusal in A; severity-triggered crossover in B) |
+| Ideal pragmatic RCT ITT | Same, with all physicians |
+
+Truths are stored in `simulations/truths/` as small committed CSVs, and each has its own Monte Carlo SE.
+
+# 6. Estimands, in plain words (for the main text)
+
+| **Method** | **The question it answers** |
+|---|---|
+| CCW, A vs B1 (the reference) | "If every newly diagnosed claudicant had PVI within 6 months, compared with if every one waited 6 months and was then treated as usual, how different would 5-year MALE be?" |
+| CCW, A vs B2 | "…compared with if no one ever had an elective PVI?" (abstention, not deferral) |
+| Sequential trials, ITT | "For a patient who is still eligible this month, what does doing the PVI now change, compared with not now (perhaps later)?" |
+| Landmark at 6 months | "Among patients still alive and MALE-free at 6 months, how do those who happened to have a PVI differ from those who didn't?" |
+| Naive "ever early PVI" | No well-defined question (immortal time) |
+| Realistic RCT, ITT | "Among patients that selected, skilled, low-intervention surgeons are willing to randomize, what does *assigning* early PVI change, given that some patients refuse it and some controls cross over?" |
+| Ideal pragmatic RCT, ITT | Same, for all surgeons |
+
+# 7. Methods grid
+
+**Adjustment (D26).** The same machinery is used everywhere. Inverse-probability weights come from pooled logistic models on the measured covariates (X, frailty index, P_t history, OBL, year). Weighted cause-specific pooled logistic hazards for MALE and death are fitted with flexible time terms and combined into standardized cumulative incidence at month 60. Weights are truncated at the 99.5th percentile in a supplementary analysis only.
+
+| **Method** | **Weights** | **Notes** |
+|---|---|---|
+| Naive | Baseline PS for "PVI by month 6" | Follow-up starts at month 0 for both groups, so immortal time is present |
+| Landmark (D13) | PS for "PVI in months 0–5" among patients at risk at month 6, using covariates up to month 6 | Outcome counted through month 60 |
+| Sequential trials (D14) | A PS for treatment at each trial baseline m = 0–5 | Person-trials are stacked. ITT is primary; per-protocol (with IPCW for deviation) is in the supplement |
+| CCW-B1, CCW-B2 | Clones are artificially censored at deviation, with IPCW from the PVI-hazard model | A clones untreated at month 5 are censored at the end of the window |
+| Realistic and ideal RCT | None | ITT by assigned arm; naive per-protocol by treatment received |
+
+**Physician variants (D15), crossed with every observational method.** The physician enters through the **weight models**, because that is where confounding is controlled.
+
+| **Variant** | **Implementation** |
+|---|---|
+| V1 Ignore | No physician terms; model-based SE |
+| V2 Cluster-robust | As V1, with the sandwich SE clustered on the physician |
+| V3 Random intercept | Physician random intercept in the weight models (`glmmTMB`/`lme4`); cluster-robust SE |
+| V4 Within-physician | Physician fixed effects in the weight models. Physicians with no variation in early PVI are dropped, and the change to the estimand that this causes is reported. Mundlak (physician means) is the fallback if FE is too slow. Cluster-robust SE |
+| V5 Preference IV (supplement) | Instrument = the physician's early-PVI rate among the physician's *other* patients; two-stage residual inclusion |
+
+**Negative controls (D32).** NCO-1 and NCO-2 replace MALE as the outcome in every method × variant. The NCE replaces early PVI as the exposure in the CCW and landmark designs. The NCE models adjust for early-PVI status (D36). Without that adjustment, the path NCE ← style → PVI → MALE would make the NCE "detect" a real effect of PVI. The unadjusted NCE analysis is in the supplement.
+
+# 8. Simulated RCTs (D18)
+
+| | **Realistic RCT** | **Ideal pragmatic RCT** |
+|---|---|---|
+| Physicians | Hospital-based (OBL = 0), top third of skill, bottom half of style | All physicians |
+| Patients | Eligible patients of those physicians; a constant 30% consent | All eligible patients |
+| n | 240 (120 per arm) | 50,000 |
+| Arm A | PVI at the next scheduled month in 0–2. Refusal 10% | Same |
+| Arm B | No PVI in months 0–5 unless severity crosses the "lifestyle-limiting" threshold (crossover generated by the DGP); usual care afterwards | Same |
+| Analyses | ITT (unadjusted cumulative incidence) and naive as-treated | Same |
+
+Running replicates of the 240-patient trial also gives its power and CI width, which is the "small trial" result.
+
+# 9. Performance measures and simulation size (D27–D29)
+
+**Notation.** For method m with estimate θ̂_m, its own estimand θ_m, and the clinically relevant estimand θ*:
+
+| **Measure** | **Definition** |
+|---|---|
+| Estimation bias | E[θ̂_m] − θ_m |
+| Question gap | θ_m − θ* |
+| Total error | E[θ̂_m] − θ* (the sum of the two above) |
+| Empirical SE | SD of θ̂_m across replicates |
+| Model SE | Mean of the estimated SEs |
+| Coverage | Coverage of θ_m by the 95% CI |
+| RMSE | Against θ_m and against θ* |
+| **False-alarm rate** (headline) | In the null scenario, the proportion of replicates whose 95% CI excludes 0 **for the disease composite (CLTI or major amputation)**, whose true effect is exactly 0 (D34) |
+
+MCSE is reported for every measure (Morris 2019).
+
+- **Variance (D28).** A sandwich SE clustered on the physician, with the delta method for the cumulative-incidence RD. Validated by a physician-cluster bootstrap (200 resamples) on 100 replicates in W0 and WR.
+- **Replicates (D29).** 100 at draft scale, 500 at final. With 500 replicates, the MCSE of coverage is about 1% and the MCSE of bias is about 0.045 × empirical SE.
+
+# 10. Negative-control diagnostic analysis (Aim 3, D33)
+
+For each world × scenario × method × variant × replicate, record:
+
+- **Flag**: the negative control's 95% CI excludes 0.
+- **Truly biased**: |θ̂_m − θ_m| > 1 percentage point. Also |θ̄_m − θ_m| > 1 pp at the cell level.
+
+**Reported:**
+- the sensitivity and specificity of each negative control, and of "any negative control", for detecting bias;
+- the distribution of |bias| among replicates with a clean negative control ("how much can hide");
+- in the supplement, each negative control's false-positive rate.
+
+# 11. Mock tables
+
+- **Table 1 (introduction). Randomized trials of revascularization for claudication.** Trial, year, n (planned vs actual), comparison, primary endpoint, follow-up, crossover, whether the trial stopped early, whether there was a hard limb outcome. See Appendix A.
+
+- **Table 2. What each method answers (W0, all three scenarios).** Rows: methods. Columns: plain-language estimand, own truth θ_m, question gap, estimation bias, coverage.
+
+- **Table 3. Bias by threat.** Rows: method × physician variant. Columns: W0–W5, WR. Cells: estimation bias (MCSE), shaded by |bias| > 1 pp.
+
+- **Table 4. Negative controls as bias detectors.** Rows: NCO-1, NCO-2, NCE, "any". Columns: sensitivity, specificity, and the median |bias| when the negative control is clean, by world.
+
+# 12. Mock figures
+
+- **Fig 1. Estimand ladder (D30).** One row per method. Each row shows the method's own truth (hollow dot), its mean estimate and 95% interval across replicates (filled dot, coloured by physician variant), and the clinically relevant truth (vertical line). Panels: W0, W1–W5, WR. Two versions: null and harm.
+
+- **Fig 2. False-alarm rate in the null scenario** (disease composite). Heat map of method × variant (rows) against world (columns). The 5% nominal level is marked. An inset shows the MALE estimate in the same null world: "harm" produced purely by how the outcome is defined.
+
+- **Fig 3. What a clean negative control is worth.** For each method, |bias| plotted for flagged and clean replicates, by world.
+
+- **Fig S. The DGP as a DAG**, and cumulative-incidence curves under A, B1 and B2 for each scenario.
+
+# 13. Position relative to the literature
+
+| **Prior work** | **What it showed** | **What this paper adds** |
+|---|---|---|
+| CLEVER, ERASE, IRONIC, SUPER; Fakhry 2018 | Small trials of walking and QoL, with crossover | The trial that could realistically be run, simulated and compared with the question surgeons ask |
+| Sorber 2023; Hicks 2020; Dun 2025/2026 | Medicare associations; physician outliers; CCW emulation | What those designs actually estimated, and which threats remain |
+| Duchesneau 2022 CEBP | Landmark vs CCW against a plasmode truth, in cancer surgery | Adds sequential trials, RCTs and the physician level; estimand-specific truths |
+| Keogh 2023; Hernán 2008/2018; Wanis 2024 | Sequential trials, CCW, grace periods in theory and simulation | The difference between deferral and abstention shown in a surgical question |
+| Li 2020 SMMR; Cafri 2020 Biom J | Unmeasured cluster confounding; within-surgeon matching | Inside TTE designs, with the surgeon acting on treatment, effect and outcome |
+| Myers 2011; Brookhart 2006 | Z-bias; preference IV | Shown for physician style in claims |
+| Paper 1 (tte-papers) | Design matters more than the adjustment method (ICU) | A surgical setting; estimand gaps; negative controls as diagnostics |
+
+# 14. Out of scope
+
+- PS matching and doubly robust estimators (Paper 1 territory).
+- Data-quality problems (Paper 2 territory): coding error beyond factor 3, missing data, disenrollment.
+- Patient–physician sorting, where sicker patients see particular physicians. Patients are assigned to physicians at random (D37); sorting is a possible single sensitivity world.
+- Cost-effectiveness.
+
+# 15. Practicalities
+
+| **Item** | **Plan** |
+|---|---|
+| Software | R (4.6.1). A new package `pvisim/` in this repo with parameters in YAML, a simulator, truths, designs, estimators and metrics (D12). Uses `fastglm`/`speedglm` for pooled logistic, `sandwich`, `glmmTMB`, `data.table`, `arrow`, `furrr`. Quarto reports; renv |
+| Validation | Before any runs: (a) every method is unbiased for its own estimand in W0 at n = 10⁶ (within MCSE); (b) a closed-form check of the CCW truth in a special case with no time-varying confounding; (c) negative-control truths are exactly 0 by construction (tests) |
+| Compute | 7 worlds (+ W1 at ρ = 0) × 3 scenarios × 500 replicates. Each replicate fits about 7 methods × 4 variants × (MALE + 3 negative controls). V3/V4 on person-month data are the bottleneck. Profile at the draft stage and switch to Mundlak if needed |
+| Sequence | 1) parameters YAML → 2) simulator + tests → 3) calibration → 4) truths → 5) estimators, validated at large n → 6) draft runs (100 replicates) → 7) report → 8) final runs → 9) manuscript |
+| Reporting | ADEMP (Morris 2019); TARGET 2025 items for the emulated trial; code and seeds public |
+| Target journal | JAMA Surgery (Original Investigation). The mathematics and the full grid go in the supplement |
+
+# 16. Decisions resolved after the first draft (D34–D37)
+
+1. **Null (D34).** Sharp disease null. The headline false-alarm rate is reported for the disease composite (CLTI or major amputation). MALE in the null world shows "harm" that comes purely from the outcome definition.
+2. **Strength of each threat when switched on (D35).** Use the empirical quantity where one exists (physician-style spread from the Hicks 2020 outlier distribution; deferred PVI from Sorber's late-PVI rates; early-PVI rate). Where none exists (skill, coding intensity), set the threat so that on its own it moves the naive or landmark estimate by about 2–4 percentage points in the harm scenario.
+3. **NCE (D36).** The NCE analysis is adjusted for early-PVI status. The unadjusted version goes in the supplement as a warning.
+4. **Patient–physician sorting (D37).** Not in the main grid. One sensitivity world, if time allows.
+
+Still to do: replace the interim mortality target (§5.4) with the 5-year mortality of the authors' own Medicare cohort.
+
+---
+
+# Appendix A. Draft Table 1: randomized trials of revascularization for intermittent claudication
+
+| **Trial** | **n (planned → randomized)** | **Comparison** | **Primary endpoint** | **Longest follow-up** | **Crossover / adherence** | **Stopped early / under-recruited** | **Hard limb outcome?** |
+|---|---|---|---|---|---|---|---|
+| CLEVER (Murphy 2012, 2015) | 252 → 111 | OMC vs OMC+SE vs OMC+stent (aortoiliac) | Peak walking time, 6 months | 18 months (79 completed) | — | Under-recruited (reduced to 3 arms) | No |
+| ERASE (Fakhry 2015; Klaphake 2022) | 212 | Endovascular + SET vs SET | Maximum walking distance, 12 months | ~5.4 years (post hoc) | SET arm: later revascularization HR 2.5 | No | No (procedures reported) |
+| IRONIC (Nordanstig 2014; Djerf 2020) | 158 | Revascularization + BMT + exercise vs BMT + exercise | SF-36, 1 year | 5 years (benefit lost) | — | Single centre | No |
+| MIMIC (Greenhalgh 2008) | 93 + 34 | PTA added to SE + BMT | Walking distance, 24 months | 24 months | — | Small | No |
+| SUPER (Koelemay 2022) | 240 | SET vs endovascular (iliac) | MWD + VascuQol, 1 year | ~5.5 years | **29% crossover by year 1**; 50% SET compliance | **Stopped early** | No |
+| Mazari (2012, 2017) | 178 | PTA vs SEP vs both (fem-pop) | Walking/QoL, 12 months | 5.2 years | — | — | No (reinterventions reported) |
+| Spronk (2009) | 151 | Endovascular vs SET | Rutherford, 6–12 months | 12 months | — | — | No |
+| OBACT (Nylaende 2007) | 434 screened → 56 | OMT ± PTA | QoL, ABI, walking, 2 years | 2 years | — | **56 of 434 screened** | No |
+
+*Missing cells are to be filled from the full texts (crossover is not reported in several abstracts).*
@@ -0,0 +1,104 @@
+# Decisions log
+
+Each decision is listed with the date and the reason for it. Settled decisions are changed only by adding a new entry below them, never by editing them.
+
+## Round 1: framing (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 1 | **Paper type: a "teaching" paper combined with a "threats" paper.** The backbone is that landmark, CCW and sequential trials answer three different clinical questions. The experimental factors are the threats specific to surgery. We do not claim that "TTE works". | Earlier simulations already show that the methods work. The new contribution is the surgical setting, plus showing what each estimand means in it. |
+| 2 | **Primary outcome: MALE composite. Secondary outcomes: each component on its own.** | Matches the VESS 2026 abstract and practice in vascular surgery. The components show how outcomes that physicians generate (reintervention) behave differently from outcomes that come from disease (CLTI, amputation). The exact definition of MALE is still open (round 2). |
+| 3 | **Include a simulated RCT arm in the same DGP** (ITT and naive per-protocol, with crossover, surgeon selection and small n). | Turns "why RCTs are hard in surgery" into numbers. |
+| 4 | **Experimental factors:** (1) surgeon skill, (2) surgeon preference or style, (3) style acting on the outcome directly, (4) deferral and crossover driven by symptom progression, (5) unmeasured symptom severity. **Built in through the RCT arm:** (6) trial selection of high-skill surgeons. **Fixed feature of the DGP:** (7) treatment effect that changes over time. **Supplement / sensitivity analysis:** (8) unmeasured co-intervention (uptake of supervised exercise). | Keeps the main grid readable. |
+| 5 | **Audience: surgeons. Target journal: JAMA Surgery.** Each method in the main text has a one-sentence clinical definition. The mathematics goes in the supplement. | The user is a surgical health-services researcher. Readers are clinicians. |
+| 6 | **The RCT limitations are covered as narrative in the introduction plus one evidence table** (claudication RCTs: n, endpoint, crossover, recruitment shortfall). There is no separate systematic review. | Scope. |
+
+## Round 2: the target trial (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 7 | **MALE = CLTI, major amputation, open bypass, or any revascularization after the first.** The index PVI does not count as an event in either arm, whenever it happens. Total procedures per patient is a secondary outcome. **Supplement:** the Sorber "late PVI" definition (any PVI after month 6), to show how the outcome definition alone changes the answer. | Avoids an outcome that is asymmetric between arms by construction. Not counting the index procedure is defined the same way in both arms. |
+| 8 | **Death is a competing event.** The primary measure is the 5-year cumulative incidence of MALE, reported as a risk difference. HR, and MALE-free survival (MALE or death), go in the supplement. | Medicare population with a mean age of about 76. Censoring at death would answer an unrealistic hypothetical question. |
+| 9 | **Time zero is the new claudication diagnosis. Strategy A is PVI within 6 months.** A 3-month window is a supplement-only sensitivity analysis. | Matches Sorber 2023, VESS 2026, Dun 2025, and the guidelines. |
+| 10 | **Strategy B, primary: no PVI in the first 6 months, then usual care.** **Strategy B, second (main text): no elective PVI unless a MALE-qualifying event occurs.** | Deferral vs abstention: "conservative" is not a single strategy (the ADAM/CAESAR lesson). |
+| 11 | **Each method gets two errors: estimation bias (estimate vs its own estimand) and question gap (its own estimand vs the clinically relevant estimand). Total error is the sum of the two.** The clinically relevant reference is A vs B(primary) in all new claudicants under full adherence, which is the CCW estimand. | The teaching paper needs to separate "the method fails" from "the method answers another question". |
+| 12 | **R, with a new self-contained package in this repo.** It follows the `tte-papers` conventions (YAML parameters, truths from simulating the interventions in large samples, Quarto, renv) but does not reuse its code. | The DGP differs (months, not hours; physicians as clusters). |
+
+## Round 3: comparators and the simulated world (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 13 | **Landmark:** month 6, among patients alive and MALE-free. Exposure is PVI in months 0–6. Outcome is MALE by month 60 *since diagnosis*. Adjust for baseline covariates plus measured covariates at month 6. | Same horizon as CCW. Mirrors the Sorber 2023 design, so the paper can say what that design estimated. |
+| 14 | **Sequential trials:** monthly trials in months 0–5 among patients still eligible. **ITT is primary**; per-protocol (weighted) goes in the supplement. | Sequential ITT estimates "now vs maybe later" (deferral), in contrast to CCW's "within 6 months vs not". |
+| 15 | **Physician-level variants, crossed with every method:** (1) ignore, (2) cluster-robust SE, (3) physician random intercept, (4) within-physician (FE or Mundlak). **Supplement:** (5) physician-preference IV. | Clustering fixes SEs, not bias. A random intercept does not remove confounding. FE can, but it risks Z-bias. The IV analysis is dubious when the surgeon also decides on reintervention. |
+| 16 | **Monthly discrete time over 60 months. Pooled logistic hazard models.** | Standard for CCW and sequential trials, and simple. |
+| 17 | **Measured:** age, sex, race, diabetes, CKD/dialysis, CHF, COPD, prior CVD, coded smoking, physician ID, practice setting (OBL/ASC vs hospital), calendar year, and a time-varying **proxy for symptom progression** (repeat vascular visits / noninvasive arterial testing). **Latent:** true symptom severity, ABI, true smoking, patient preference, surgeon skill, surgeon style. **Proxy strength is the dial for factor 5.** | What is "measured" means what is in Medicare claims. |
+| 18 | **Two simulated RCTs:** a *realistic* one (n ≈ 240, hospital-based, top-third skill, low-preference surgeons) and an *ideal pragmatic* one (all physicians, large n). Crossover is **generated by the DGP** in both. Each is analysed by ITT and naive per-protocol. | Separates bias from surgeon selection and small n from bias due to crossover. |
+| 19 | **Calibrate usual care to Medicare figures:** early PVI ≈ 3.2%; the physician distribution from Hicks 2020 (median 0%, range 0–58%, about 5.6% outliers); Sorber's crude late PVI and CLTI rates; age-appropriate mortality. **Truth scenarios:** null; harm (e.g. RD +5 points at 5 years); early benefit then late harm. | The null scenario shows pure bias. Calibration keeps the simulated world credible to surgeons. |
+
+## Round 4: DGP core (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 20 | **Mechanistic latent-severity engine.** Each patient has a latent severity that drifts with a patient-specific rate. Severity drives the CLTI hazard and the PVI decision. PVI lowers severity right away but adds a restenosis/progression hazard that depends on skill. | Benefit followed by harm, crossover in the conservative arm, and confounding by indication all come out of one mechanism. Truths are computed by simulating the interventions, so the extra complexity costs nothing there. |
+| 21 | **Estimand: PVI done by the patient's own physician, averaged over the real physician mix.** The gap between the realistic RCT (top-third skill) and the population effect is reported as a finding (external validity), not as bias. | Skill modifies the effect of PVI. |
+| 22 | **Skill and style are separate latent traits. Base correlation −0.3; comparison level 0.** | Sorber's outliers had more late PVI and CLTI. |
+| 23 | **Factor 3 = high users have a lower threshold for repeat revascularization AND code CLTI more intensively.** The two are switched on together. | These are paths outside the treatment decision. Style acting through deferred PVI after month 6 is part of strategy B, not bias. |
+| 24 | **Main analysis: one factor at a time from a clean base world, plus one "realistic, all threats on" world, in each of the 3 truth scenarios. Supplement: full factorial**, only if the one-at-a-time results suggest interactions. | Readability: one panel per threat. |
+| 25 | **Each replicate has 50,000 patients and about 1,500 physicians, with lognormal panel sizes (median about 20). Supplement: 200,000 patients.** Number of replicates is set later, with ADEMP. | About 1,600 early-PVI patients per replicate, enough for CCW and FE positivity. |
+
+## Round 5: analysis structure (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 26 | **All observational methods use the same adjustment.** IP weights come from pooled logistic models on the measured covariates. Weighted cause-specific hazards for MALE and death are then combined into a standardized 5-year cumulative incidence. Landmark uses IPW at month 6. Sequential trials use IPW at each trial baseline. CCW uses artificial-censoring weights. The RCTs are analysed unadjusted (ITT) and as-treated (naive per-protocol). No matching and no doubly robust variants. | Differences between methods then reflect design and estimand only. Paper 1 already covered the choice of adjustment method. |
+| 27 | **Performance measures:** estimation bias (against the method's own estimand), question gap, total error, empirical SE, model SE, coverage of the method's own estimand, RMSE, and MCSE for each. **The headline measure for clinicians is the false-alarm rate under the null** (CI excludes 0). | Separates two problems: whether the method fails, and whether it answers a different question. |
+| 28 | **Variance:** cluster-robust sandwich SE plus the delta method in every replicate. A physician-cluster bootstrap (200 resamples) is run on about 100 replicates in the base and realistic worlds. | A full bootstrap at scale costs too much. Same approach as Paper 1. |
+| 29 | **Replicates:** 100 at draft, 500 at final. MCSE is reported. | MCSE for coverage is about 1% at 500 replicates. |
+| 30 | **Figure 1, the "estimand ladder":** one row per method. Each row shows a hollow dot for the method's own truth, a filled dot with interval for its estimate, and a line for the clinically relevant truth. Panels are worlds; colour is the physician variant. Shown for the null and harm scenarios. **Figure 2:** false-alarm rate by threat and physician variant. | One picture that shows both the question gap and the estimation bias. |
+| 31 | **Next:** the analysis plan (md plus docx), then the table of claudication RCTs, then the build. The docs are committed on the branch `analysis-plan`. | |
+
+## Round 6: negative controls (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 32 | **Three negative controls.** **NCO-1:** hip fracture / pneumonia admission. Shares the latent patient frailty, severity and smoking, so it targets factor 5. **NCO-2:** varicose vein ablation by the same physician. Shares physician procedural style, so it targets factors 2–3. **NCE:** early low-value procedure (e.g. vein ablation) within 6 months, which shares both. All three are generated only from latent traits and never depend on PVI, so their true effect is exactly 0. They go through the same method × physician-variant pipeline as MALE. | Negative controls are the practical bias check available in real claims. CCW and landmark condition on survival, so the controls also need testing under each design. |
+| 33 | **Main text (Figure 3 / Table 3):** how accurately the negative controls detect bias. A control is "flagged" if its CI excludes 0. Bias in the MALE risk difference counts if |bias| > 1 percentage point. Report sensitivity and specificity, plus how much bias can sit behind a clean negative control. **Supplement:** each control's own false-positive rate. | Answers a new, practical question: what does a clean negative control actually tell a surgeon? |
+
+## Round 7: issues raised by drafting the plan (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 34 | **Null scenario = sharp disease null** (Δ = κ = 0). The headline false-alarm rate is reported for the disease composite (CLTI or major amputation), whose true effect is exactly 0. MALE in the null world is reported as "harm" that comes from the outcome definition alone. | Repeat revascularization can only happen after an index PVI, so MALE cannot be exactly null. Showing that is a teaching point. |
+| 35 | **Strength of each threat:** use an empirical anchor where one exists (style spread from Hicks 2020; deferred PVI from Sorber late-PVI rates; the early-PVI rate). Where none exists (skill, coding), set the threat so that on its own it moves the naive or landmark estimate by about 2–4 pp in the harm scenario. | Credible to surgeons while keeping each threat's effect visible. |
+| 36 | **The NCE analysis is adjusted for early-PVI status.** The unadjusted version goes in the supplement. | Otherwise NCE ← style → PVI → MALE flags a real effect of PVI. |
+| 37 | **Patient–physician sorting is not in the main grid.** It becomes one sensitivity world if time allows. | Scope. |
+| — | **Fact (not a decision):** interim calibration target for all-cause mortality is about 35% at 5 years (Sigvant 2018 *EJVES*, PMID 29478910: 10-year mortality in IC 63%). To be replaced by the 5-year mortality in the authors' own Medicare cohort. | No published 5-year all-cause figure for Medicare claudicants was found. |
+
+## Round 8: scope (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 38 | **No data-quality aim.** The paper keeps the three aims of the plan. Claims data problems such as time-zero misdating, physician attribution, laterality and Medicare Advantage switching are not modelled. | Focus. Papers 1 and 2 of `tte-papers` have stopped. |
+| 39 | **Paper 1's comparison of adjustment methods and covariate vibration stay out** (see D26). | This paper is about which question each design answers. |
+
+## Build notes: simulator and calibration (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 40 | **Factor 5 (unmeasured confounding by indication): "off" means the analyst sees true severity Z_t, true frailty and true smoking. "On" means the analyst sees only claims proxies:** a 3-month visit/testing count for severity, a frailty index with reliability 0.4, and coded smoking (sensitivity 0.6). | In the clean worlds every confounder must be measured, so each method can be checked as unbiased for its own estimand. |
+| 41 | **Each world's intercepts are re-solved by calibration (early PVI, deferred PVI, repeat revascularization, CLTI, death) so that every world matches the same Medicare targets.** A world differs from W0 only in the threat it switches on. | A fair comparison between worlds. |
+| 42 | **Deferral:** in W0–W3, W5 and WR, deferred PVI matches Sorber (3.6% late PVI among patients without early PVI). **W4 is calibrated to 15%**, which is trial-like crossover, and in W4 deferral is driven by severity. In WR, deferral is at the Sorber rate but driven by severity. PVI after CLTI was lowered to about 1% per month, so most late PVI is elective. | At the Medicare rate, the deferral threat would be almost invisible. Trials (SUPER 29% in year 1) show much more crossover. |
+| 43 | **Physician style coefficient 1.6** (W2, W3, WR). Among physicians with ≥11 patients this gives a median early-PVI rate of 0.8%, a maximum of 75%, and 4.7% of physicians ≥2 SD outliers who perform 38% of early PVIs. Hicks 2020 reported 0%, 58% and 5.6%. | Close enough to Hicks 2020. Not tuned further. |
+| 44 | **Within a month the order is:** severity update → outcome events → death → elective PVI decision. A MALE event and death in the same month count the MALE event first. | Matches discrete-time pooled logistic analysis. |
+
+## Build notes: truths, estimators, validation (2026-10-27)
+
+| # | Decision | Why |
+|---|---|---|
+| 45 | **Scenario effects, tuned in W0 against the disease composite:** harm (relief 0.3, kappa 0.039) gives RD(A−B1) **+5.1 pp**. Benefit-then-harm (relief 2.0, kappa 0.075) gives **−1.2 pp at 24 months and +3.1 pp at 60 months**. The same mechanism is used in every world. Where skill varies (W1, WR) the true effects are larger (+7 pp), because exp(−skill) is convex: on average a PVI is done worse than one at average skill. | D34 sets the headline on the disease composite. MALE carries an extra +7 to +15 pp that comes purely from the outcome definition. |
+| 46 | **Truths for the negative controls are taken from the truths table, not assumed to be 0.** In the non-null scenarios PVI changes death, so a negative-control outcome's 5-year risk moves slightly through the competing risk (e.g. NCO-1 −0.08 pp). In the null scenario the effect is exactly 0. | Honest scoring. |
+| 47 | **The own estimand of landmark and naive is estimated within each replicate by an oracle run:** the same estimator with exact severity, true frailty and smoking, and the latent skill and style. Estimation bias is the paired difference (estimate − oracle). | A large-n limit was too noisy: the landmark MCSE at n = 10⁵ was 2.8 pp. |
+| 48 | **The PVI model is fitted separately inside the window (month indicators) and in the deferred period (spline in month).** | A pooled model misspecified the window probabilities that the CCW and sequential weights depend on. |
+| 49 | **V2 is the V1 fit with a physician-cluster-robust SE.** Under V4, physicians without an estimable fixed effect in a phase get probability 0, and the V4 population is restricted to physicians with variation in early PVI. | Same point estimate by construction (D15). |
+| 50 | **Validation (W0, harm, 40 replicates of 50,000):** CCW A−B1 is +5.12 vs a truth of +5.11 on the disease composite (z = 0.0); sequential ITT is +5.83 vs +5.71 (z = 0.6); MALE and NCO-1 have |z| ≤ 0.7. The model SE is 0–20% below the empirical SE for CCW, so the bootstrap check in D28 still needs doing. | Plan §15 validation (a). |
@@ -0,0 +1,168 @@
+# Literature review 1: limitations of RCTs in surgery, mapped to TTE on claims data
+
+All citations were checked against PubMed through NCBI E-utilities, 2026-10-27. The numbers come from PubMed abstracts; full texts were not read. Papers with no PubMed abstract are described only in general terms.
+
+**Corrections and caveats**
+- Chapman 2014 was published in *BMJ*, not *Surgery*.
+- The page numbers for McCulloch 2002 are from memory, not from PubMed.
+- BASIL-3 recruited 481 of the planned 861 patients. That the paclitaxel signal caused the shortfall is a guess, not a verified finding.
+
+---
+
+## 1. Classic critiques and the IDEAL framework
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| McCulloch P et al. 2002, *BMJ* "Randomised trials in surgery: problems and possible solutions" | 12065273 | The founding critique. Barriers: equipoise, learning curve, quality control, blinding, surgeon preference. |
+| McLeod RS 1996, *Surgery* 119:483 | 8619200 | Problems specific to RCTs in surgery. |
+| Solomon MJ, McLeod RS 1995, *Surgery* | 7652679 | Many surgical questions could be randomized but are not. |
+| Solomon MJ, McLeod RS 1998, *Med J Aust* | 9803252 | Overview of RCTs in surgery: past, present and future. |
+| Barkun JS et al. 2009, *Lancet* 374:1089 | 19782874 | IDEAL series, paper 1: the stages of surgical innovation. |
+| Ergina PL et al. 2009, *Lancet* 374:1097 | 19782875 | Few of the challenges are unique to surgery, but in surgery they occur together: complexity, surgeon factors, and the range of outcomes. |
+| McCulloch P et al. 2009, *Lancet* 374:1105 | 19782876 | IDEAL recommendations: registries and prospective databases, with RCTs "whenever possible". |
+| McCulloch 2013 / Ergina 2013 / Cook 2013, *BMJ* (IDEAL series) | 23778427, 23778426, 23778425 | Observational studies in the Exploration and Assessment stages; RCTs and long-term registry evaluation. |
+| Hirst A et al. 2019, *Ann Surg* | 29697448 | IDEAL update: adds a Pre-IDEAL stage and covers devices. |
+| Cook JA 2009, *Trials* 10:9 | 19200379 | Checklist for the design, conduct and analysis of surgical RCTs. |
+| Pronk AJM et al. 2023, *Br J Surg* 110:1300 | 37379487 | Of 1,188 surgical RCTs, only **22.1%** published in 2019 were at low risk of bias (14.7% in 2009). |
+| Ahmed Ali U et al. 2018, *J Surg Res* 228:1 | 29907196 | Of 228 negative surgical RCTs, the median sample size was about 76–80, so most were underpowered. |
+
+## 2. Recruitment failure, discontinuation, and equipoise
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| **Chapman SJ et al. 2014, *BMJ* 349:g6870** | 25491195 | Of 395 surgical RCTs, **21% stopped early**. The main reason was poor recruitment (44% of stopped trials). Of completed trials, 34% were unpublished at a median of 4.9 years. |
+| **Rosenthal R et al. 2015, *Ann Surg* 262:68** | 24979608 | Discontinued: **43%** of surgical vs **27%** of medical RCTs. Stopped for slow recruitment: 18% vs 11%. Discontinuation predicted non-publication (OR 4.18). |
+| Kasenda B et al. 2014, *JAMA* 311:1045 | 24618966 | Across all fields, **24.9%** of RCTs were discontinued, most often for poor recruitment. |
+| Cook JA et al. 2008, *Contemp Clin Trials* 29:631 | 18397843 | **Trials of surgery vs medical treatment recruited about 2× slower** than trials of surgery vs surgery (median 68 vs 35 days per recruit). This applies directly to PVI vs exercise. |
+| Donovan J et al. 2002, *BMJ* 325:766 (ProtecT) | 12364308 | Qualitative work on how recruiters presented the options made an "impossible" trial feasible. |
+| Donovan JL et al. 2014, *J Clin Epidemiol* 67:912 | 24811157 | Clinicians' "hunches" undermined equipoise, and **surgeons felt this most intensely**. |
+| Rooshenas L et al. 2016, *PLoS Med* 13:e1002147 | 27755555 | In recorded consultations, the equipoise clinicians stated differed from how they actually presented the options. |
+| Donovan JL et al. 2022, *Trials* 23:258 (QRI-Two) | 35379301 | A rapid recruitment intervention for trials that are failing to recruit. |
+| Hamdy FC et al. 2023, *NEJM* 388:1547 (ProtecT at 15 years) | 36912538 | A monitor-vs-intervene trial can be completed. |
+
+## 3. Surgeon skill, learning curve, expertise-based trials, and clustering
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| **Devereaux PJ et al. 2005, *BMJ* 330:88** | 15637373 | Proposes expertise-based RCTs: each surgeon performs only the procedure they are expert in, to avoid differential-expertise bias. |
+| Cook JA et al. 2015, *Trials* 16:241 | 26025450 | Review of 43 expertise-based RCTs: 92% compliance, but expertise criteria were poorly reported. |
+| Cook JA et al. 2004, *Clin Trials* 1:421 | 16279280 | Statistical methods for learning curves in surgical trials. |
+| Ramsay CR et al. 2000 / 2001 | 11155830, 11319991 | Assessing learning curves in health technology. |
+| **Lee KJ, Thompson SG 2005, *BMJ* 330:142**; *Clin Trials* 2:163 | 15649931, 16279138 | Outcomes cluster by clinician even under individual randomization. Ignoring it makes the estimates look more precise than they are. |
+| Roberts C, Roberts SA 2005, *Clin Trials* 2:152 | 16279137 | Clustering induced by treatment, which can differ between arms (partially nested designs). |
+| Walters SJ 2010, *J Clin Nurs* 19:1102 | 20492055 | Therapist effects: cluster-level, random-effects, and marginal models compared. |
+| **Birkmeyer JD et al. 2013, *NEJM* 369:1434** | 24106936 | Peer ratings of surgical skill from video: the bottom quartile vs the top had **14.5% vs 5.2%** complications and **0.26% vs 0.05%** mortality. |
+| Birkmeyer JD et al. 2003, *NEJM* 349:2117 | 14645640 | Surgeon volume and mortality. |
+
+## 4. Blinding and sham surgery
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| **Wartolowska K et al. 2014, *BMJ* 348:g3253** | 24850821 | Across 53 placebo-controlled surgical RCTs, the placebo arm improved in 74%. Surgery was **no better than placebo in 51%**. |
+| Wartolowska K et al. 2016, *BMJ Open* | 27008687 | Placebo-controlled surgical trials are feasible but recruit slowly. |
+| Boutron I et al. 2017, *Ann Intern Med* 167:40 | 28630973 | CONSORT extension for non-pharmacologic treatments: expertise, blinding, and standardization. |
+
+## 5. Standardization, crossover, and contamination
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| **Garas G et al. 2018, *Ann Surg Oncol* 25:221** | 29110271 | Across 40 RCTs of minimally invasive vs open surgery, crossover averaged **12.6% (0–45%)** and worsened the ITT results. |
+| **Weinstein JN et al. 2006, *JAMA* 296:2441 (SPORT)** | 17119140 | **Only 50%** of the surgery arm had surgery, and **30%** of the nonoperative arm did, so the ITT comparison could not be interpreted. |
+| Hernán MA, Robins JM 2017, *NEJM* 377:1391 | 28976864 | Per-protocol effects in pragmatic trials. |
+
+## 6. Generalizability and external validity
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| **Rothwell PM 2005, *Lancet* 365:82** | 15639683 | "To whom do the results of this trial apply?" |
+| Rothwell PM 2006, *PLoS Clin Trials* | 16871331 | Factors that affect external validity. |
+| Cutteridge J et al. 2025, *BJS Open* | 40105904 | Patients in 307 vascular RCTs were **younger** than those in the UK National Vascular Registry. |
+| Jones HG et al. 2026, *EJVES* | 42526742 | Endovascular lower-limb revascularization is under-researched relative to how often it is done. Its trials are 72% efficacy trials and 39% industry-funded. |
+| Loudon K et al. 2015, *BMJ* (PRECIS-2) | 25956159 | Tool for rating how pragmatic or explanatory a trial is. |
+
+## 7. Trials of timing, "when to operate" (the key analogue)
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| **Fitzgibbons RJ et al. 2006, *JAMA* 295:285** | 16418463 | Hernia watchful waiting: 23% of the waiting arm crossed over to surgery, and 17% of the surgery arm crossed the other way. |
+| **Fitzgibbons RJ et al. 2013, *Ann Surg* 258:508** | 24022443 | Long-term crossover was **68%** (79% in men over 65). "Watchful waiting" works out as "delayed surgery". |
+| O'Dwyer PJ et al. 2006, *Ann Surg* | 16858177 | UK hernia trial: crossover was higher than predicted. |
+| **UK Small Aneurysm Trial 1998, *Lancet*** | 9853436 | HR 0.94. Operative mortality of 5.8% caused an early deficit in survival. |
+| **Lederle FA et al. 2002, *NEJM* (ADAM)** | 12000813 | **61.6% of the surveillance arm were eventually repaired.** Surveillance is a dynamic strategy, not "no treatment". |
+| Cao P et al. 2011, *EJVES* (CAESAR) | 20869890 | Delayed repair under surveillance: 59.7% by 36 months and 84.5% by 54 months. One in six aneurysms became unsuitable for EVAR, a harm of delay that the ITT effect on mortality does not capture. |
+| Ouriel K et al. 2010, *JVS* (PIVOTAL) | 20304589 | 31% of the surveillance arm were repaired by 20 months. HR 1.01. |
+| Ulug P et al. 2020, *Cochrane* | 32609382 | Early repair of small AAAs has no advantage. |
+
+## 8. Pragmatic trials, registry trials, and TTE
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| **Hernán MA, Robins JM 2016, *AJE* 183:758** | 26994063 | The foundational TTE paper. |
+| **Hernán MA et al. 2016, *J Clin Epidemiol* 79:70** | 27237061 | Specifying the target trial avoids self-inflicted biases such as immortal time. |
+| Hernán MA et al. 2022, *JAMA* 328:2446 | 36508210 | JAMA Guide to TTE. |
+| Cain LE et al. 2010, *Int J Biostat* | 21972433 | "When to start treatment?": dynamic regimes. |
+| **Maringe C et al. 2020, *Int J Epidemiol* 49:1719** | 32386426 | CCW tutorial on surgery in older patients with lung cancer. The naive 22% benefit fell to 11%. |
+| Wang SV et al. 2023, *JAMA* 329:1376 (RCT-DUPLICATE) | 37097356 | Of 32 RCTs emulated in claims data, agreement was r = 0.93 when the emulation was close and r = 0.53 when it was not. All were drug trials. |
+| Cashin AG et al. 2025, *JAMA* / *BMJ* (TARGET) | 40899949 / 40903028 | Reporting guideline for TTE. |
+| Scola G et al. 2023, *BMC Med Res Methodol* | 37587484 | Only 16% of TTEs used cloning to handle immortal time. |
+| Katsoulis M et al. 2024, *Epidemiology* 35:730 | 39024012 | TTE of bariatric surgery. |
+| Hoskins W et al. 2026, *JBJS Rev* | 42896011 | Orthopaedic primer. TTE **does not** address unmeasured confounding. |
+| Lauer MS, D'Agostino RB 2013, *NEJM* | 23991657 | The randomized registry trial. |
+| Zwarenstein M et al. 2008, *BMJ* | 19001484 | CONSORT extension for pragmatic trials. |
+| Brookhart MA et al. 2006, *Epidemiology* 17:268 | 16617275 | Physician preference as an instrumental variable. |
+
+**Gap:** we found no surgery-specific systematic review of TTE methods (JAMA Surg / Ann Surg, 2022–2026). Applied examples: Sakowitz 2025, *Ann Surg* (TCAR vs CEA), PMID 40607703; Talathi 2026, *Surgery* (STRASS 2 emulated), PMID 41936769.
+
+## 9. Vascular surgery
+
+| Paper | PMID | Takeaway |
+|---|---|---|
+| Farber A et al. 2022, *NEJM* (BEST-CLI) | 36342173 | Cohort 1: MALE or death 42.6% with surgery vs 57.4% with endovascular (HR 0.68). Whether a vein conduit was available decided which cohort a patient entered. |
+| Bradbury AW et al. 2023, *Lancet* (BASIL-2) | 37116524 | n = 345. Vein bypass first was worse (aHR 1.35), which conflicts with BEST-CLI. |
+| Hunt 2017 *Trials*; Bradbury 2025 *BMJ* (BASIL-3) | 28526046, 39993822 | **481 of 861 planned** randomized. No benefit from drug-coated balloons or drug-eluting stents. |
+| **Katsanos K et al. 2018, *JAHA*** | 30561254 | Paclitaxel: 5-year mortality 14.7% vs 8.1% (RR 1.93). Only 3 RCTs had 5-year data. |
+| Rocha-Singh KJ et al. 2020, *Circulation* | 32370548 | Individual patient data: **23–24% lost to follow-up or withdrew**. After recovering vital status, HR 1.38 → 1.27. |
+| Secemsky EA et al. 2019, *JAMA Cardiol* | 30747949 | Medicare: no excess mortality. |
+| Freisinger E et al. 2020, *Eur Heart J* | 31593987 | German claims data: no excess mortality. |
+| Nordanstig J et al. 2020, *NEJM* (SWEDEPAD, interim) | 33296560 | Registry RCT: HR 1.06. |
+| **Nordanstig J et al. 2025, *Lancet* (SWEDEPAD 2, claudication)** | 40902614 | No benefit in quality of life. 5-year mortality HR 1.47 (1.09–1.98). |
+| Falkenberg M et al. 2025, *Lancet* (SWEDEPAD 1, CLTI) | 40902617 | The CLTI companion trial. |
+| Fakhry F et al. 2015, *JAMA* (ERASE) | 26547465 | Revascularization plus exercise beat exercise alone at 12 months. |
+| Murphy TP et al. 2012, *Circulation* (CLEVER) | 22090168 | Exercise was better for walking time and stenting for quality of life. n = 111. |
+| Nordanstig 2014, *Circulation*; Djerf 2020, *Circ CV Interv* (IRONIC) | 25095886, 31937137 | The quality-of-life benefit seen at 1 year was **lost by 5 years**, at about twice the cost. |
+| **Sorber R et al. 2023, *JVS*** | 37276171 | Medicare. Early PVI 3.2%. Late PVI 22.5% vs 3.6% (aHR 6.89). High-use physicians. |
+| Dun C et al. 2025, *Ann Vasc Surg* 110:424 | 39419322 | Market competition predicts early PVI. |
+| Weaver ML et al. 2023, *JVS* | 37141948 | Market competition and claudication management in the VQI. |
+| Madabhushi V et al. 2021, *JVS* | 33775749 | Revascularizing claudicants: CLTI HR 2.9, amputation HR 4.5. |
+| Conte MS et al. 2025, *JVS* (SVS focused update) | 40316185 | Major evidence gaps for claudication. |
+
+---
+
+## Taxonomy: RCT limitations in surgery → what TTE on claims can do
+
+IV = internal validity, EV = external validity, F = feasibility.
+
+| # | Limitation | Threatens | Can TTE on claims fix it? | Can it go in the simulation DGP? |
+|---|---|---|---|---|
+| 1 | Recruitment failure and discontinuation | F (and EV) | **Yes.** No consent is needed. The price is confounding instead of randomization. | Indirectly: compare a "would-consent" subset with the whole population. |
+| 2 | Patient and surgeon equipoise and preference | F, EV, IV | **Partly.** Preference is what makes treatment vary in practice, and can serve as an IV (high-use physicians). Patient preference is unmeasured. | **Yes.** A latent physician "intervention propensity". |
+| 3 | Surgeon skill and learning curve | IV, EV | **Partly.** Skill is not in claims, but the operator NPI is, so random or fixed effects are possible. Skill and preference cannot be separated if they are correlated. | **Yes.** Latent skill, correlated with preference. |
+| 4 | Clustering and provider heterogeneity | IV (precision), EV | **Mostly.** Mixed models, cluster-robust SEs. | **Yes.** Random intercepts and random slopes by physician. |
+| 5 | Blinding and sham | IV | **No.** No patient-reported outcomes in claims, and outcome ascertainment is not blinded. | Optional: differential ascertainment of CLTI. |
+| 6 | Standardization and co-interventions | IV, EV | **Partly.** Procedure type and Part D fills are visible. Exercise, smoking, and adherence are not. | **Yes.** Unmeasured exercise uptake. |
+| 7 | Crossover and deferral | IV | **Yes, TTE's strongest case**, using dynamic strategies with CCW or sequential trials. It needs the triggers for crossover to be measured. | **Yes, central.** |
+| 8 | Immortal time and misaligned time zero | IV | **Yes, by design.** | **Yes.** The naive design as a comparator. |
+| 9 | Selected surgeons, centers, and patients | EV | **Yes, largely.** Population-based data, plus transportability. | **Yes.** A "trial of high-skill operators" vs the whole population. |
+| 10 | Long horizon, rare outcomes, changing technology | F, IV | **Partly.** Follow-up is long and complete, but device codes are coarse and calendar time confounds. | **Yes.** Time-varying effect (early benefit, late harm). |
+
+**In one line:** TTE fixes feasibility, generalizability, and (with dynamic strategies) crossover. It cannot fix unmeasured confounding by indication or the lack of patient-reported outcomes. Surgeon skill sits in between: it cannot be measured, but its structure can be modelled through the NPI.
+
+## Read first
+1. Sorber 2023 *JVS* (37276171)
+2. ADAM 2002 *NEJM* (12000813) and CAESAR 2011 (20869890)
+3. Fitzgibbons 2006 *JAMA* and 2013 *Ann Surg* (16418463, 24022443)
+4. Hernán & Robins 2016 *AJE*; Hernán 2016 *J Clin Epidemiol* (26994063, 27237061)
+5. Maringe 2020 *Int J Epidemiol* (32386426)
+6. McCulloch 2002 *BMJ*; Ergina 2009 *Lancet* (12065273, 19782875)
+7. Chapman 2014 *BMJ*; Rosenthal 2015 *Ann Surg* (25491195, 24979608)
+8. Katsanos 2018; Nordanstig 2020 and 2025 (30561254, 33296560, 40902614)
@@ -0,0 +1,208 @@
+# Literature review 2: claudication evidence, and what landmark, CCW and sequential trials each estimate
+
+PMIDs were checked through NCBI E-utilities on 2026-10-27. The numbers come from PubMed abstracts. The guideline recommendation wording is from memory and needs checking against the full texts.
+
+---
+
+## Part A. Clinical evidence
+
+### A1. Randomized trials. Every one used walking distance or QoL as the endpoint; none used hard limb outcomes.
+
+| Trial | N and comparison | Result | Limitations |
+|---|---|---|---|
+| **CLEVER**: Murphy 2012 *Circulation* (22090168); 18-month results 2015 *JACC* (25766947); design (18440181); update (19660897) | 111 aortoiliac; OMC vs OMC+SE vs OMC+stent | Peak walking time at 6 months: SE better than stent (p=0.04). At 18 months: no difference (p=0.16). QoL better with stent. | **Under-recruited**: planned 252 in 4 arms, cut to 217 in 3 arms, finished with 111. Only 79 had 18-month data. |
+| **ERASE**: Fakhry 2015 *JAMA* (26547465); long-term Klaphake 2022 *Ann Surg* (33378308) | 212; endovascular + SET vs SET | MWD at 12 months +282 m with the combination | At 5.4 years the walking gap had closed (p=0.62). The SET arm had more later revascularizations (HR 2.50) but **fewer procedures in total (65 vs 149)**. |
+| **IRONIC**: Nordanstig 2014 *Circulation* (25095886); 5-year Djerf 2020 *Circ CV Interv* (31937137) | 158; revascularization + BMT + exercise vs BMT + exercise | SF-36 better with revascularization at 1 year | **No difference at 5 years.** Cost about 2× higher. |
+| **MIMIC**: Greenhalgh 2008 *EJVES* (19022184) | 93 femoropopliteal + 34 aortoiliac; PTA on top of SE + BMT | Walking distance +38% / +78% | Very small. No QoL benefit. |
+| **SUPER**: Koelemay 2022 *EJVES* (35151572) | 240 iliac; SET vs endovascular | MWD the same at 1 year | **Stopped early. 29% of the SET arm crossed over to revascularization within 1 year.** Only 50% complied with SET. |
+| **Mazari**: 2012 *BJS* (22021102); 2017 *BJS* (27763685) | 178; PTA vs SEP vs both | The same at 12 months and 5.2 years | |
+| **Spronk**: 2009 *Radiology* (19188327) | 151; endovascular vs SET | An early advantage for endovascular was gone by 6–12 months | |
+| **OBACT**: Nylaende 2007 *EJVES* (17055756) | 56; OMT vs OMT + PTA | PTA better at 2 years | **434 screened to randomize 56.** |
+| **Cochrane**: Fakhry 2018 (29518253) | 10 RCTs, 1,087 patients | Endovascular vs SET: no clear difference | Small, short trials, and **no hard limb outcomes**. |
+
+**Across these trials:**
+- Revascularization improves walking and QoL early, but the advantage fades by 2–5 years.
+- Crossover and non-adherence to SET are large.
+- No trial was powered for CLTI, amputation or cumulative reintervention.
+
+### A2. Observational evidence that intervention is associated with progression
+
+- **Madabhushi 2021 *JVS* (33775749).** n = 1,051. Revascularization was associated with CLTI (HR 2.9) and amputation (HR 4.5). The exposure is "ever revascularized", so the result is exposed to immortal time bias.
+- **Golledge 2018 *BJS* (29566427).** n = 456. Treatment depended on *which specialist the patient saw*, which is a natural physician-preference contrast. 5-year amputation was 6.2% with early revascularization vs 0.7% with conservative care.
+- **Sorber 2023 *JVS* (37276171).** aHR 6.89 for late PVI. Patients of outlier physicians had more late PVI and more CLTI.
+- **Dun 2025 *Surgery* (40782724).** Markov decision model.
+- **Alabi 2026 *JAMA Surg* (41405976).** VQI–Medicare, 10,012 patients. CLTI within 180 days of revascularization ranged up to 11.8% across subgroups.
+- **Bose 2025 *Ann Surg* (38841837) and *EJVES* (38906366).** Infrapopliteal PVI for claudication had poor outcomes.
+- **Natural history:**
+  - Sigvant 2016 *EJVES* (26777541): about 21% progress to CLI.
+  - Baubeta Fridh 2017 *EJVES* (28797662): amputation after revascularization for IC is 0.4% per year.
+- **Review:** Golledge 2022 *Nat Rev Cardiol* (34997200).
+
+### A3. Variation between physicians
+
+- **Hicks 2020 *JVS* (31208940), "Overuse of early PVI for claudication."** n = 194,974 patients, 3.2% with early PVI. The median physician's rate was 0%, with a range of 0–58%. **320 physicians (5.6%) were ≥2 SD outliers.** Office-based lab (OBL) or ambulatory surgery center (ASC) practice predicted outlier status. **This paper is the main empirical source for the physician component of the DGP.**
+- Dun 2024 *JAHA* (38958132): the OBL/ASC share of PVI went from 12.4% to 55.7% between 2011 and 2022.
+- Dun 2025 *Ann Vasc Surg* (39419322): market competition.
+- Weaver 2023 *JVS* (37141948).
+- Tsou 2024 *JVS* (38908807).
+- Smith 2019 *Ann Vasc Surg* (30684609).
+- Brown 2021 *JVS* (33617980).
+- Iannuzzi 2025 *JVS* (39824263).
+
+### A4. Guidelines
+
+- SVS 2015 (25638515; introduction 25721065)
+- ESVS 2024 (37949800)
+- ESC 2024 (39210722)
+- ACC/AHA 2024 (38743805; *JACC* 38752899). Advises against revascularizing only to prevent progression to CLTI.
+
+### A5. Why there is no RCT of early PVI vs initial conservative therapy on hard outcomes
+
+1. **Event rates are low, so the trial would have to be large and long.** These are rough two-proportion calculations with 80% power, ≥5-year follow-up and no crossover:
+   - Amputation 2% vs 1% at 5 years: about 2,300 per arm.
+   - CLTI 8% vs 5%: about 1,050 per arm.
+2. **Crossover is built into the comparator.** The conservative arm includes deferred PVI, so an ITT analysis dilutes the contrast. Examples: SUPER had 29% crossover in year 1; ERASE's SET arm had a revascularization HR of 2.5; CLEVER finished with 111 of 252 planned; OBACT randomized 56 of 434 screened.
+3. **Equipoise is asymmetric.** The guidelines already recommend against early PVI, and the high-use physicians who do most early PVI would not enroll their patients.
+4. **The hard endpoints are generated by physicians.** Reintervention and the coding of CLTI are decisions made by unblinded clinicians, and patients who get early PVI are watched more closely.
+5. **No sponsor and a lot of heterogeneity.** No device company will fund a "fewer devices" hypothesis. Anatomy and devices vary, and Medicare has covered SET only since 2017.
+
+---
+
+## Part B. Methods: which question each one answers
+
+### B1. Immortal time bias
+
+**Plain version:** if "early PVI" means "had PVI at month 4", that patient had to survive event-free to month 4 to be labelled early PVI. The fix is to line up eligibility, assignment and time zero.
+
+- Suissa 2008 *AJE* (18056625)
+- Lévesque 2010 *BMJ* (20228141)
+- Hernán 2016 *J Clin Epidemiol* (27237061)
+- Hernán 2025 *Epidemiology* (39494894). Reframes the problem: immortal time is a symptom of misclassification or selection.
+
+### B2. Landmark analysis
+
+**Estimand:** among patients alive and CLTI-free at 6 months, those who happened to have PVI in months 0–6 vs those who did not.
+
+- It is conditional and descriptive.
+- The comparator mixes patients who never have PVI with patients who have it at month 7.
+- Events in months 0–6 are dropped.
+- It conditions on a post-treatment variable.
+- **Sorber 2023 is effectively a 6-month landmark design.**
+
+References:
+- Anderson, Cain, Gelber 1983 *JCO* (6668489); revisited 2008 (18711176)
+- Dafni 2011 *Circ CV Qual Outcomes* (21586725)
+- van Houwelingen 2008 (18836831)
+- Putter & van Houwelingen 2017 (29225713)
+
+### B3. Clone–censor–weight (CCW)
+
+**Estimand:** the risk if *everyone* followed strategy A, "PVI within 6 months of diagnosis", vs strategy B, "no PVI in the first 6 months, then usual care (or PVI only for CLTI)". This is a per-protocol effect of sustained strategies starting at diagnosis.
+
+- The 6-month grace period means an event inside the window, before any PVI, counts in both arms.
+- **The answer depends on how strategy B is written.**
+- Natural and stochastic grace-period strategies are different estimands (Wanis 2024).
+
+References:
+- Hernán 2018 *BMJ* (29419381)
+- Cain 2010 *Int J Biostat* (21972433)
+- **Maringe 2020 *IJE* (32386426)**; correspondence 33349843, 33524137
+- Gaber 2024 *Cancer Med* (39642890)
+- **Wanis 2024 *JRSS-A* (39281783)**, on grace periods
+- Hernán & Robins 2017 *NEJM* (28976864)
+- Cole & Hernán 2008 *AJE* (18682488)
+
+### B4. Sequential trial emulation
+
+**Estimand:** among patients still eligible *this month* (no PVI yet, no CLTI, alive), PVI now vs not now.
+
+- "Not now" includes PVI next month, so it is a *deferral* effect, not an *abstention* effect.
+- The pooled population leans toward slowly progressing patients, because they stay eligible for more trials.
+
+References:
+- Hernán 2008 *Epidemiology* (18854702)
+- Danaei 2013 *SMMR* (22016461)
+- Gran 2010 *Stat Med* (20803557)
+- **Keogh 2023 *Stat Med* (37086186)**
+- Limozin/Seaman 2025 *SMMR* (40631395): CIs
+- Robins 2000 *Epidemiology* (10955408): MSMs
+- The TrialEmulation R package is on CRAN and has no PubMed citation.
+
+### B5. Studies comparing methods
+
+- **Duchesneau 2022 *CEBP* (35984990).** Surgery in metastatic breast cancer. Compares landmark with CCW against a plasmode known truth. This is **the closest template.**
+- Mi 2016 *Stat Med* (27350312). Landmark is unbiased under the null but biased when there is an effect.
+- Karim 2016 *AJE* (27455963)
+- Zhou 2005 *AJE* (16192344)
+- Pilleron 2023 *BJC* (36759725). Immortal time bias is larger in older patients.
+- Hoffman 2022 *JAMA Netw Open* (36190729)
+- Moreno-Betancur 2026 *Epidemiology* (41490058). "The Ideal Trial": the estimand you want vs the estimand your data allow.
+
+**Gap:** we found no simulation that compares landmark, CCW and sequential trials side by side as answers to *different* estimands.
+
+### B6. Physician clustering and unmeasured confounding at the surgeon level
+
+There are two separate problems.
+
+**1. Inference.** Clustered SEs or random effects fix only the standard errors.
+- Panageas 2003 *Ann Intern Med* (14568854); 2007 *Stat Med* (17016864)
+- Kahan & Morris 2013 *Stat Med* (23112128)
+
+**2. Confounding by physician.** Physician style affects both early PVI and the outcome. Here the outcome (late PVI, coding of CLTI) is *itself* a physician decision.
+- **Within-physician comparisons** (fixed effects, Mundlak, matching within surgeon) remove time-invariant confounding between physicians:
+  - Neuhaus & Kalbfleisch 1998 (9629647)
+  - Begg & Parides 2003 (12898546)
+  - **Li Y 2020 *SMMR* (31694489)**: closed-form bias results
+  - **Cafri 2020 *Biom J* (32419247)**: matching within surgeon on device data
+  - Li/Zaslavsky/Landrum 2013 (23526267)
+  - Thoemmes & West 2011 (26735885)
+  - Arpino & Cannas 2016 (26833893)
+  - Fuentes 2022 (34128730)
+  - Austin 2017 (27885709): median HR
+- **Caveats:**
+  - A random effect does **not** remove confounding when it is correlated with PVI use.
+  - Frailty HRs are conditional on physician and non-collapsible.
+  - Fixed effects need physicians who sometimes do and sometimes do not intervene early (positivity).
+- **Preference as an instrumental variable:**
+  - Brookhart 2006 (16617275); Brookhart & Schneeweiss 2007 (19655038)
+  - Rassen 2009 (19345561)
+  - Brookhart 2010 (20354968)
+  - Ertefaie 2017 (28239929)
+  - Hutchings 2022 (35588540)
+  - The exclusion restriction is doubtful when the same physician also decides on reintervention.
+- **Bias amplification (Z-bias):** Myers 2011 *AJE* (22025356); Pearl 2011 *AJE* (22034488).
+  - If physician preference is purely instrument-like, conditioning on the physician *amplifies* remaining unmeasured patient confounding.
+  - The simulation can show where the crossover lies.
+
+**Gap:** we found no paper on unmeasured surgeon-level confounding *inside* CCW or sequential trial emulation. The nearest are Cafri 2020 and Li 2020.
+
+### B7. Simulation methods
+
+Morris, White, Crowther 2019 *Stat Med* (30652356), the ADEMP framework. **Score each method against its own true estimand.**
+
+### B8. Reporting
+
+TARGET 2025: *BMJ* 40903028, *JAMA* 40899949; explanation and elaboration in *PLoS Med* 41091811.
+
+---
+
+## Summary table (early PVI ≤6 months)
+
+| Method | Estimand | Requires | Failure mode |
+|---|---|---|---|
+| Naive "ever early PVI" | None that is well defined | — | Immortal time and confounding by indication |
+| Landmark at 6 months | Among patients event-free at 6 months: had PVI in months 0–6 vs not | Confounders measured up to the landmark | Drops early harms. Mixed comparator. Selects on survival after treatment. |
+| CCW with 6-month grace | Everyone follows A vs everyone follows B | Time-varying predictors of PVI timing; positivity; a precise definition of strategy B | Unmeasured symptom progression or physician style. Extreme weights. Sensitive to how B is defined. |
+| Sequential trials (ITT) | PVI now vs not now, among patients eligible this month | Eligibility each month; covariates at each baseline | A deferral effect, not abstention. Leans toward slow progressors. |
+| Physician-clustered models | Clustered SEs: honest SEs. FE: within-physician effect. Frailty: conditional HR | FE needs within-physician variation | RE does not remove confounding. FE can show Z-bias. Conditional vs marginal mismatch. |
+
+## Read first
+
+1. Hernán 2018 *BMJ*, CCW (29419381)
+2. Hernán 2008 *Epidemiology*, sequential trials (18854702)
+3. Keogh 2023 *Stat Med* (37086186)
+4. Maringe 2020 *IJE* (32386426)
+5. Duchesneau 2022 *CEBP* (35984990)
+6. Morris 2019 *Stat Med*, ADEMP (30652356)
+7. Cafri 2020 *Biom J* (32419247)
+8. Myers 2011 *AJE* (22025356)
@@ -0,0 +1,24 @@
+Package: pvisim
+Title: Known-Truth Simulator for Target Trial Emulation of Early Peripheral Vascular Intervention
+Version: 0.1.0
+Authors@R: person("Chen", "Dun", email = "[email protected]", role = c("aut", "cre"))
+Description: Data-generating process (claudicants nested in physicians, latent
+    severity, physician skill and style), truths by simulating interventions,
+    target trial emulation designs (naive, landmark, sequential trials,
+    clone-censor-weight), simulated RCTs, physician-level variants, negative
+    controls and performance metrics. Not for CRAN.
+License: MIT
+Encoding: UTF-8
+Depends: R (>= 4.3)
+Imports:
+    stats,
+    utils,
+    yaml,
+    data.table
+Suggests:
+    testthat (>= 3.0.0),
+    fastglm,
+    sandwich,
+    glmmTMB,
+    pkgload
+Config/testthat/edition: 3
@@ -0,0 +1,3 @@
+exportPattern("^[a-z]")
+import(data.table)
+importFrom(stats, plogis, qlogis, qnorm, qpois, rnorm, runif, rbinom, rlnorm, rgamma, setNames)
@@ -0,0 +1,321 @@
+# Estimators (plan sections 7 and 9; D13-D15, D26, D28).
+#
+# One machinery for every observational design: inverse-probability weights from
+# pooled logistic models of elective PVI, then a weighted Aalen-Johansen estimator of
+# the cumulative incidence of the outcome (death competing), nonparametric in time.
+# Standard errors come from the influence function of the weighted Aalen-Johansen
+# estimator with the weights treated as known (usually conservative for IPW; checked
+# against a physician-cluster bootstrap, D28). Summed by patient for the naive SE, by
+# physician for the cluster-robust SE. Clones and person-trials of one patient are summed
+# together, so cloning and stacking are accounted for.
+
+BASE_COVS <- c("age", "female", "black", "year", "diabetes", "ckd", "chf", "copd", "prior_cvd",
+               "smoking", "frailty_index", "obl")
+VARIANTS <- c("V1_ignore", "V2_cluster", "V3_random", "V4_within")
+
+# ---- inputs -------------------------------------------------------------------------
+
+#' Person-month analysis table: eligibility for elective PVI and the PVI decision.
+#'
+#' A row is month t for a patient who survived month t's events and death and had no
+#' MALE and no PVI before the decision in month t: the risk set of the PVI model.
+pvi_risk_set <- function(rec) {
+  p <- rec$patients
+  extra <- intersect("t_pvi_orig", names(p))
+  pm <- merge(rec$pm, p[, c("id", "phys_id", BASE_COVS, "t_pvi", "pvi_elective", "t_male", "t_death", extra),
+                        with = FALSE], by = "id")
+  pm <- pm[(is.na(t_death) | t_death > t) & (is.na(t_male) | t_male > t) & (is.na(t_pvi) | t_pvi >= t)]
+  pm[, pvi := as.integer(!is.na(t_pvi) & t_pvi == t)]
+  sev0 <- rec$pm[t == 0, .(id, sev0 = sev_obs)]
+  pm <- merge(pm, sev0, by = "id")
+  pm[, `:=`(win = as.integer(t < rec$par$size$window), tt = pmin(t, 12))]
+  if (length(extra)) pm[, pvi_before := as.integer(!is.na(t_pvi_orig) & t_pvi_orig < t)]
+  pm[]
+}
+
+# ---- treatment models ------------------------------------------------------------------
+
+#' Fit P(PVI at t | eligible, history) and return predictions for `newdata`.
+#'
+#' The decision window (months 0..W-1) and the deferred period are fitted as separate
+#' models: practice inside and after the window differs, and a pooled model would
+#' misspecify the window probabilities that the CCW and sequential weights depend on.
+#' variant: V1/V2 no physician terms; V3 physician random intercept (glmer, nAGQ = 0);
+#' V4 physician fixed effects (fixest). Under V4 a physician with no PVI in a phase has
+#' no estimable intercept for it; the maximum-likelihood probability there is 0.
+fit_pvi <- function(data, newdata, variant, extra_covs = NULL) {
+  out <- numeric(nrow(newdata))
+  for (w in c(1L, 0L)) {
+    i_new <- newdata$win == w
+    if (!any(i_new)) next
+    d <- data[win == w]
+    nd <- newdata[i_new]
+    tt_levels <- if (w == 1L) "factor(tt)" else "splines::ns(t, df = 4)"
+    # no exposure at all in a phase (the NCE after the window): probability 0
+    out[i_new] <- if (sum(d$pvi) == 0) 0 else fit_pvi_phase(d, nd, variant, tt_levels, extra_covs)
+  }
+  out
+}
+
+fit_pvi_phase <- function(d, nd, variant, time_terms, extra_covs = NULL) {
+  covs <- c(BASE_COVS, "sev_obs", "sev0", extra_covs)
+  if (variant == "V4_within") covs <- setdiff(covs, "obl")
+  rhs <- paste(c(covs, time_terms), collapse = " + ")
+  if (variant %in% c("V1_ignore", "V2_cluster")) {
+    f <- as.formula(paste("~", rhs))
+    X <- model.matrix(f, d)
+    fit <- fastglm::fastglm(X, d$pvi, family = binomial(), method = 2)
+    b <- coef(fit); b[is.na(b)] <- 0
+    Xn <- model.matrix(f, nd)
+    return(as.vector(plogis(Xn %*% b)))
+  }
+  if (variant == "V3_random") {
+    fit <- lme4::glmer(as.formula(paste("pvi ~", rhs, "+ (1 | phys_id)")), data = d,
+                       family = binomial(), nAGQ = 0, control = lme4::glmerControl(calc.derivs = FALSE))
+    return(as.vector(predict(fit, newdata = nd, type = "response", allow.new.levels = TRUE)))
+  }
+  if (variant == "V4_within") {
+    fit <- fixest::feglm(as.formula(paste("pvi ~", rhs, "| phys_id")), data = d,
+                         family = binomial(), notes = FALSE, warn = FALSE)
+    pr <- as.vector(predict(fit, newdata = nd, type = "response"))
+    pr[is.na(pr)] <- 0
+    return(pr)
+  }
+  stop("unknown variant ", variant)
+}
+
+#' Physicians with variation in early (window) PVI: the V4 analysis population.
+physicians_with_variation <- function(rec) {
+  W <- rec$par$size$window
+  x <- rec$patients[, .(e = sum(pvi_elective %in% TRUE & t_pvi < W), n = .N), by = phys_id]
+  x[e > 0 & e < n, phys_id]
+}
+
+# ---- weighted Aalen-Johansen with influence functions --------------------------------
+
+#' Cumulative incidence at `horizon` from person-month rows.
+#'
+#' @param d data.table with id, phys_id, arm, t, w (weight for the row's risk set),
+#'   ev (1 if the outcome occurs in month t), dth (1 if death without the outcome in t).
+#' @return list(F = named vector of cumulative incidence by arm at the horizon,
+#'   inf = data.table(id, phys_id, arm, inf) of influence-function contributions summed
+#'   by patient, curve = cumulative incidence by arm and month).
+aj_cif <- function(d, horizon = 60L) {
+  d <- d[t < horizon & w > 0]
+  s <- d[, .(R = sum(w), E = sum(w * ev), D = sum(w * dth)), keyby = .(arm, t)]
+  s[, `:=`(he = E / R, hd = D / R)]
+  s[, h := he + hd]
+  s[, S_prev := shift(cumprod(1 - h), fill = 1), by = arm]
+  s[, Fcum := cumsum(S_prev * he), by = arm]
+  s[, FH := Fcum[.N], by = arm]
+  s[, `:=`(a = S_prev - (FH - Fcum) / (1 - h), b = -(FH - Fcum) / (1 - h))]
+  s[!is.finite(a), a := 0]; s[!is.finite(b), b := 0]
+  d <- s[d, on = .(arm, t), .(id, phys_id, arm, x = w * ((ev - he) * a + (dth - hd) * b) / R)]
+  inf <- d[, .(inf = sum(x)), by = .(id, phys_id, arm)]
+  FH <- s[, .(FH = FH[1]), by = arm]
+  list(F = setNames(FH$FH, FH$arm), inf = inf, curve = s[, .(arm, t, cuminc = Fcum)])
+}
+
+#' Risk difference arm1 - arm0 with naive (patient) and cluster-robust (physician) SEs.
+rd_from_aj <- function(aj, arm1, arm0) {
+  x <- dcast(aj$inf, id + phys_id ~ arm, value.var = "inf", fill = 0)
+  diff <- x[[arm1]] - x[[arm0]]
+  se_pat <- sqrt(sum(diff^2))
+  se_phys <- sqrt(sum(tapply(diff, x$phys_id, sum)^2))
+  data.table(est = aj$F[[arm1]] - aj$F[[arm0]], se_patient = se_pat, se_cluster = se_phys,
+             risk1 = aj$F[[arm1]], risk0 = aj$F[[arm0]])
+}
+
+# outcome and death indicators for a person-month table that has columns t_out, t_death
+mark_events <- function(d) {
+  d[, ev := as.integer(!is.na(t_out) & t_out == t)]
+  d[, dth := as.integer(ev == 0L & !is.na(t_death) & t_death == t)]
+  d
+}
+
+# person-months from `start` to the outcome, death or the horizon, one row per month
+expand_followup <- function(p, start, horizon = 60L) {
+  end <- pmin(ifelse(is.na(p$t_out), horizon - 1L, p$t_out),
+              ifelse(is.na(p$t_death), horizon - 1L, p$t_death), horizon - 1L)
+  keep <- end >= start
+  p <- p[keep]; st <- if (length(start) > 1) start[keep] else rep(start, nrow(p)); end <- end[keep]
+  nrep <- end - st + 1L
+  d <- p[rep(seq_len(nrow(p)), nrep)]
+  d[, t := st[rep(seq_along(nrep), nrep)] + sequence(nrep) - 1L]
+  mark_events(d)
+}
+
+# ---- designs --------------------------------------------------------------------------
+
+#' For the within-physician variant, restrict to physicians whose patients include both
+#' early PVI and no early PVI (positivity); the estimand changes accordingly (plan 7).
+restrict_variant <- function(rec, variant) {
+  if (variant != "V4_within") return(rec)
+  keep <- physicians_with_variation(rec)
+  rec$patients <- rec$patients[phys_id %in% keep]
+  rec$pm <- rec$pm[id %in% rec$patients$id]
+  rec
+}
+
+#' Run one design under one physician variant for several outcomes.
+#' Returns one row per outcome (and contrast), with the SE that the variant reports:
+#' patient-level for V1, physician-cluster-robust for V2-V4.
+run_design <- function(rec, design, variant, outcomes = c(male = "t_male", disease = "t_disease",
+                                                           nco1 = "t_nco1", nco2 = "t_nco2")) {
+  rec <- restrict_variant(rec, variant)
+  rs <- if (design %in% c("ccw", "seq")) pvi_risk_set(rec) else NULL
+  phat <- if (design == "ccw") fit_pvi(rs, rs, variant, rec$extra_covs) else NULL
+  cmp <- if (is.null(rec$comparators)) c("B1", "B2") else rec$comparators
+  out <- lapply(names(outcomes), function(o) {
+    r <- switch(design,
+      naive = cbind(contrast = "early_none", est_naive(rec, outcomes[[o]], variant)),
+      landmark = cbind(contrast = "early_none", est_landmark(rec, outcomes[[o]], variant, extra_covs = rec$extra_covs)),
+      ccw = est_ccw(rec, outcomes[[o]], variant, rs = copy(rs), phat = phat, comparators = cmp),
+      seq = cbind(contrast = "now_not", est_seq(rec, outcomes[[o]], variant, rs = copy(rs))))
+    r[, `:=`(outcome = o, design = design, variant = variant, n = nrow(rec$patients))]
+    use_patient <- variant == "V1_ignore"
+    r[, se := if (use_patient) se_patient else se_cluster]
+    r
+  })
+  rbindlist(out, fill = TRUE)
+}
+
+#' Naive "ever early PVI" design: exposure = PVI in months 0..W-1, follow-up from month 0,
+#' baseline propensity score. Immortal time is present by construction.
+est_naive <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60L, extra_covs = NULL) {
+  W <- rec$par$size$window
+  p <- copy(rec$patients)
+  p[, `:=`(t_out = get(outcome), arm = ifelse(!is.na(t_pvi) & t_pvi < W, "early", "none"))]
+  p <- merge(p, rec$pm[t == 0, .(id, sev_obs)], by = "id")
+  p[, pvi := as.integer(arm == "early")]
+  ps <- fit_pvi_ps(p, variant, covs = c(BASE_COVS, "sev_obs", extra_covs))   # sev_obs is the month-0 value here
+  p[, w := ifelse(arm == "early", 1 / ps, 1 / (1 - ps))]
+  d <- expand_followup(p[!is.na(w)], 0L, horizon)
+  rd_from_aj(aj_cif(d, horizon), "early", "none")
+}
+
+# a point-treatment propensity score with the same covariates and physician variants
+fit_pvi_ps <- function(p, variant, covs = c(BASE_COVS, "sev_obs", "sev0")) {
+  if (variant == "V4_within") covs <- setdiff(covs, "obl")
+  rhs <- paste(covs, collapse = " + ")
+  if (variant %in% c("V1_ignore", "V2_cluster")) {
+    X <- model.matrix(as.formula(paste("~", rhs)), p)
+    fit <- fastglm::fastglm(X, p$pvi, family = binomial(), method = 2)
+    return(as.vector(plogis(X %*% coef(fit))))
+  }
+  if (variant == "V3_random") {
+    fit <- lme4::glmer(as.formula(paste("pvi ~", rhs, "+ (1 | phys_id)")), data = p,
+                       family = binomial(), nAGQ = 0, control = lme4::glmerControl(calc.derivs = FALSE))
+    return(as.vector(predict(fit, type = "response")))
+  }
+  fit <- fixest::feglm(as.formula(paste("pvi ~", rhs, "| phys_id")), data = p, family = binomial(),
+                       notes = FALSE, warn = FALSE)
+  as.vector(predict(fit, newdata = p, type = "response"))
+}
+
+#' Landmark at month W (D13): patients alive and MALE-free at the start of month W;
+#' exposure = PVI in months 0..W-1; covariates at diagnosis and at month W-1;
+#' outcome counted in months W..horizon-1.
+est_landmark <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60L, extra_covs = NULL) {
+  W <- rec$par$size$window
+  p <- copy(rec$patients)[(is.na(t_death) | t_death >= W) & (is.na(t_male) | t_male >= W)]
+  p[, `:=`(t_out = get(outcome), arm = ifelse(!is.na(t_pvi) & t_pvi < W, "early", "none"))]
+  p <- p[is.na(t_out) | t_out >= W]
+  p <- merge(p, rec$pm[t == W - 1L, .(id, sev_obs)], by = "id")
+  p <- merge(p, rec$pm[t == 0, .(id, sev0 = sev_obs)], by = "id")
+  p[, pvi := as.integer(arm == "early")]
+  ps <- fit_pvi_ps(p, variant, covs = c(BASE_COVS, "sev_obs", "sev0", extra_covs))
+  p[, w := ifelse(arm == "early", 1 / ps, 1 / (1 - ps))]
+  d <- expand_followup(p[!is.na(w)], W, horizon)
+  rd_from_aj(aj_cif(d, horizon), "early", "none")
+}
+
+#' Clone-censor-weight (D10, D26). Arms: A (PVI by month W-1), B1 (no PVI in months
+#' 0..W-1), B2 (no elective PVI ever). A clone is censored after the PVI decision of the
+#' month in which it deviates; weights are 1 / P(not censored through the previous month).
+est_ccw <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60L,
+                    comparators = c("B1", "B2"), rs = NULL, phat = NULL) {
+  W <- rec$par$size$window
+  if (is.null(rs)) rs <- pvi_risk_set(rec)
+  if (is.null(phat)) phat <- fit_pvi(rs, rs, variant)
+  rs[, ph := phat]
+  rs <- rs[!is.na(ph)]
+  # everyone eligible at diagnosis is cloned, including patients with a MALE event or
+  # death in month 0 (events in the grace period count in both arms)
+  p <- copy(rec$patients)
+  p[, t_out := get(outcome)]
+  base <- expand_followup(p, 0L, horizon)
+  # P(PVI) for each eligible month; months where a patient is not at risk of PVI have
+  # probability 0 of deviating
+  base <- rs[, .(id, t, ph)][base, on = .(id, t)]
+  base[is.na(ph), ph := 0]
+  base[, `:=`(tp = fifelse(is.na(t_pvi), Inf, as.numeric(t_pvi)), el = pvi_elective %in% TRUE)]
+
+  out <- list()
+  # arm A: censored after month W-1 if no PVI by then (and still at risk of PVI then)
+  a <- copy(base)[, arm := "A"]
+  cens_A <- a[t == W - 1L & ph > 0 & tp > W - 1L, id]          # deviated at month W-1
+  a <- a[!(id %in% cens_A & t > W - 1L)]
+  pA <- base[t == W - 1L & ph > 0, .(id, pA = ph)]
+  a <- pA[a, on = "id"]
+  a[, w := ifelse(t > W - 1L & !is.na(pA) & tp == W - 1L, 1 / pA, 1)]
+  out$A <- a[, .(id, phys_id, arm, t, w, ev, dth)]
+  for (cmp in comparators) {
+    b <- copy(base)[, arm := cmp]
+    last <- if (cmp == "B1") W - 1L else horizon
+    b <- b[!(el & tp <= last & t > tp)]                         # censored after month tp
+    b[, q := ifelse(t <= last, 1 - ph, 1)]
+    b[, w := 1 / shift(cumprod(q), fill = 1), by = id]
+    out[[cmp]] <- b[, .(id, phys_id, arm, t, w, ev, dth)]
+  }
+  aj <- aj_cif(rbindlist(out), horizon)
+  rbindlist(lapply(comparators, function(cmp) cbind(contrast = paste0("A_", cmp), rd_from_aj(aj, "A", cmp))))
+}
+
+#' Sequential trials, ITT (D14): a trial at each month m = 0..W-1 among patients eligible
+#' for PVI at m; treatment = PVI at m; propensity score from covariates at m; outcome
+#' counted in months m+1..horizon-1; trial RDs pooled with weights = number eligible.
+est_seq <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60L, rs = NULL) {
+  W <- rec$par$size$window
+  if (is.null(rs)) rs <- pvi_risk_set(rec)
+  trials <- rs[t < W]
+  trials[, tt := t]                       # trial-month indicators in the PS model
+  ps <- fit_pvi(trials, trials, variant)
+  trials[, ps := ps]
+  trials <- trials[!is.na(ps)]
+  trials[, `:=`(m = t, arm = ifelse(pvi == 1L, "pvi", "no"),
+                w = ifelse(pvi == 1L, 1 / ps, 1 / (1 - ps)))]
+  p <- rec$patients[, .(id, t_out = get(outcome))]
+  trials <- p[trials, on = "id"]
+  trials <- trials[is.na(t_out) | t_out > m]
+  res <- list(); infs <- list()
+  for (mm in 0:(W - 1)) {
+    tr <- trials[m == mm]
+    d <- expand_followup(tr[, .(id, phys_id, arm, w, t_out, t_death)], mm + 1L, horizon)
+    aj <- aj_cif(d, horizon)
+    x <- dcast(aj$inf, id + phys_id ~ arm, value.var = "inf", fill = 0)
+    res[[mm + 1L]] <- data.table(m = mm, n = nrow(tr), est = aj$F[["pvi"]] - aj$F[["no"]])
+    infs[[mm + 1L]] <- data.table(id = x$id, phys_id = x$phys_id, m = mm, inf = x$pvi - x$no)
+  }
+  res <- rbindlist(res); infs <- rbindlist(infs)
+  res[, wt := n / sum(n)]
+  infs <- res[, .(m, wt)][infs, on = "m"]
+  pi <- infs[, .(inf = sum(wt * inf)), by = .(id, phys_id)]
+  data.table(est = sum(res$wt * res$est), se_patient = sqrt(sum(pi$inf^2)),
+             se_cluster = sqrt(sum(tapply(pi$inf, pi$phys_id, sum)^2)), risk1 = NA_real_, risk0 = NA_real_)
+}
+
+#' Negative-control exposure (D32, D36): the early low-value procedure (NCE, months
+#' 0..W-1) takes the place of early PVI as the exposure; early-PVI status is adjusted
+#' for (time-varying "PVI before t" in CCW, "PVI in the window" at the landmark).
+#' Only the landmark and CCW (A vs B1) designs are run.
+nce_records <- function(rec, adjust_pvi = TRUE) {
+  W <- rec$par$size$window
+  p <- copy(rec$patients)
+  p[, `:=`(t_pvi_orig = t_pvi, pvi_before = as.integer(!is.na(t_pvi) & t_pvi < W))]
+  p[, `:=`(t_pvi = t_nce, pvi_elective = !is.na(t_nce))]
+  rec$patients <- p
+  rec$extra_covs <- if (adjust_pvi) "pvi_before" else NULL
+  rec$comparators <- "B1"
+  rec
+}
@@ -0,0 +1,29 @@
+# Own estimands of the non-interventional designs (plan 5.5). The landmark and naive
+# contrasts have no intervention that defines them; their own estimand is what the same
+# estimator gives with every confounder measured (oracle): exact severity, true frailty
+# and smoking, and the latent physician traits. A large-n limit is too noisy (the
+# landmark estimator has few treated patients), so the oracle estimator is run on every
+# replicate's own data: estimation bias = mean(estimate - oracle estimate), paired, and
+# the oracle mean estimates the own estimand.
+
+#' The replicate's records with the hidden truth in place of the claims proxies.
+oracle_records <- function(rec) {
+  o <- rec
+  o$pm <- merge(rec$pm[, .(id, t)], rec$hidden_pm, by = c("id", "t"))[, .(id, t, sev_obs = z)]
+  h <- rec$hidden
+  o$patients <- merge(copy(rec$patients)[, c("frailty_index", "smoking") := NULL],
+                      h[, .(id, frailty_index = frailty, smoking = smoking_true, skill, style)], by = "id")
+  o$extra_covs <- c("skill", "style")
+  o
+}
+
+#' Oracle versions of the landmark and naive designs (V1, patient-level SE).
+run_oracle <- function(rec, outcomes = c(male = "t_male", disease = "t_disease", nco1 = "t_nco1", nco2 = "t_nco2")) {
+  o <- oracle_records(rec)
+  rbindlist(lapply(names(outcomes), function(nm) {
+    lm <- est_landmark(o, outcomes[[nm]], "V1_ignore", extra_covs = o$extra_covs)
+    nv <- est_naive(o, outcomes[[nm]], "V1_ignore", extra_covs = o$extra_covs)
+    rbind(cbind(design = "landmark", lm), cbind(design = "naive", nv))[
+      , `:=`(outcome = nm, contrast = "early_none", variant = "oracle", se = se_patient)]
+  }))
+}
@@ -0,0 +1,51 @@
+# Parameters: read the YAML, then apply a world and a truth scenario as overrides.
+
+#' Recursively merge `override` into `base` (lists only; scalars replace).
+deep_merge <- function(base, override) {
+  for (nm in names(override)) {
+    if (is.list(override[[nm]]) && is.list(base[[nm]])) {
+      base[[nm]] <- deep_merge(base[[nm]], override[[nm]])
+    } else {
+      base[[nm]] <- override[[nm]]
+    }
+  }
+  base
+}
+
+#' Path of the default parameter file, inside the installed or loaded package.
+default_params_path <- function() {
+  p <- system.file("params", "dgp_default.yml", package = "pvisim")
+  if (p == "") stop("dgp_default.yml not found; load the package with pkgload::load_all('pvisim')")
+  p
+}
+
+#' Load parameters for one world and one truth scenario.
+#'
+#' @param world name in `worlds:` (W0, W1, W1_rho0, W2, W3, W4, W5, WR).
+#' @param scenario name in `scenarios:` (null, harm, benefit_harm).
+#' @param path YAML file; defaults to the package's dgp_default.yml.
+#' @param overrides optional list merged last (used by calibration and tests).
+#' @param calibrated apply the world's entry in `calibration_file` (intercepts solved by
+#'   simulations/01_calibrate.R under the null scenario), after the world and before the
+#'   scenario. Calibration itself runs with FALSE.
+load_params <- function(world = "W0", scenario = "null", path = default_params_path(),
+                        overrides = NULL, calibrated = TRUE,
+                        calibration_file = "simulations/config/calibrated.yml") {
+  raw <- yaml::read_yaml(path)
+  if (!world %in% names(raw$worlds)) stop("unknown world: ", world)
+  if (!scenario %in% names(raw$scenarios)) stop("unknown scenario: ", scenario)
+  par <- raw[setdiff(names(raw), c("worlds", "scenarios"))]
+  par <- deep_merge(par, raw$worlds[[world]])
+  if (calibrated && file.exists(calibration_file)) {
+    cal <- yaml::read_yaml(calibration_file)[[world]]
+    if (!is.null(cal)) par <- deep_merge(par, cal)
+    par$calibrated <- !is.null(cal)
+  } else {
+    par$calibrated <- FALSE
+  }
+  par <- deep_merge(par, raw$scenarios[[scenario]])
+  if (!is.null(overrides)) par <- deep_merge(par, overrides)
+  par$world <- world
+  par$scenario <- scenario
+  par
+}
@@ -0,0 +1,90 @@
+# Simulated randomized trials (plan section 8, D18).
+#
+# Realistic RCT: hospital-based (OBL = 0) physicians in the top third of skill and the
+# bottom half of style; 30% of their eligible patients consent; n = 240.
+# Ideal pragmatic RCT: all physicians, everyone consents, large n.
+# Assignment A: PVI at a scheduled month 0..2 (uniform); 10% refuse and follow usual care.
+# Assignment B: no protocol PVI in the window; crossover happens through the DGP's own
+# PVI process (usual care), and usual care continues after the window.
+
+RCT_DESIGNS <- list(
+  realistic = list(n = 240L, consent = 0.30, select = TRUE),
+  ideal = list(n = NA, consent = 1.00, select = FALSE)
+)
+
+#' Physicians a realistic trial would recruit.
+trial_physicians <- function(phys) {
+  phys[obl == 0 & skill >= quantile(skill, 2 / 3) & style <= median(style), phys_id]
+}
+
+#' Draw a trial population from `pop` (patients at diagnosis) and its random assignment.
+#' @return list(pop = population restricted to enrolled patients, arm, refuse, sched).
+enrol_trial <- function(pop, design = c("realistic", "ideal"), n = NULL, seed = 1L) {
+  design <- match.arg(design)
+  spec <- RCT_DESIGNS[[design]]
+  set.seed(seed)
+  pat <- pop$pat
+  if (spec$select) pat <- pat[phys_id %in% trial_physicians(pop$phys)]
+  pat <- pat[runif(.N) < spec$consent]
+  n <- if (!is.null(n)) n else if (is.na(spec$n)) nrow(pat) else spec$n
+  if (nrow(pat) < n) stop("not enough consenting patients: ", nrow(pat), " < ", n)
+  pat <- pat[sort(sample.int(.N, n))]
+  pat[, id := seq_len(.N)]
+  arm <- sample(rep(c("A", "B"), length.out = n))
+  list(pop = list(phys = pop$phys, pat = pat), arm = arm,
+       refuse = runif(n) < 0.10, sched = sample(0:2, n, replace = TRUE))
+}
+
+#' The trial protocol as a patient-specific policy.
+trial_policy <- function(trial, assign = NULL) {
+  arm <- if (is.null(assign)) trial$arm else rep(assign, length(trial$arm))
+  force <- arm == "A" & !trial$refuse
+  policy("custom", decide = function(t, natural, idx) {
+    go <- natural
+    f <- force[idx]
+    go[f] <- t >= trial$sched[idx][f]       # PVI at the scheduled month, or the first month eligible after it
+    go
+  })
+}
+
+#' Run one simulated trial and analyse it: ITT (by assignment, unadjusted) and naive
+#' as-treated (PVI in the window vs not, ignoring assignment).
+run_rct <- function(par, design = "realistic", rep = 1L, outcomes = c(male = "t_male", disease = "t_disease"),
+                    source_n = 60000L) {
+  seed <- replicate_seed(par$world, paste0(par$scenario, "|rct|", design), rep)
+  pop <- make_population(par, n = source_n, seed = seed)
+  trial <- enrol_trial(pop, design, n = if (design == "ideal") par$size$n_patients else NULL, seed = seed + 1L)
+  ev <- simulate_course(trial$pop, par, trial_policy(trial), seed = seed + 2L, record = FALSE)$events
+  W <- par$size$window
+  ev[, `:=`(phys_id = trial$pop$pat$phys_id, assigned = trial$arm,
+            treated_win = ifelse(!is.na(t_pvi) & t_pvi < W, "early", "none"))]
+  rbindlist(lapply(names(outcomes), function(o) {
+    p <- ev[, .(id, phys_id, assigned, treated_win, t_out = get(outcomes[[o]]), t_death)]
+    itt <- rd_from_aj(aj_cif(expand_followup(copy(p)[, `:=`(arm = assigned, w = 1)], 0L)), "A", "B")
+    at <- rd_from_aj(aj_cif(expand_followup(copy(p)[, `:=`(arm = treated_win, w = 1)], 0L)), "early", "none")
+    rbind(cbind(analysis = "itt", itt), cbind(analysis = "as_treated", at))[
+      , `:=`(outcome = o, design = paste0("rct_", design), n = nrow(ev),
+             crossover_B = mean(ev$treated_win[ev$assigned == "B"] == "early"),
+             nonadherence_A = mean(ev$treated_win[ev$assigned == "A"] == "none"))]
+  }))
+}
+
+#' ITT truth for a trial design: everyone in a large trial-eligible population assigned
+#' A vs everyone assigned B, same protocol (refusal, scheduling, natural crossover).
+#' Consent is random, so it does not change the trial population and is skipped here.
+rct_truth <- function(par, design = "realistic", n = 1e6, seed = 20261029L, horizon = 60L) {
+  pop <- truth_population(par, n = n, seed = seed)
+  pat <- if (design == "realistic") pop$pat[phys_id %in% trial_physicians(pop$phys)] else pop$pat
+  set.seed(seed + 1L)
+  k <- nrow(pat)
+  trial <- list(pop = list(phys = pop$phys, pat = copy(pat)[, id := seq_len(.N)]),
+                arm = rep("A", k), refuse = runif(k) < 0.10, sched = sample(0:2, k, replace = TRUE))
+  runs <- lapply(c("A", "B"), function(a)
+    simulate_course(trial$pop, par, trial_policy(trial, assign = a), seed = seed + 2L, record = FALSE)$events)
+  rbindlist(lapply(names(OUTCOMES), function(o) {
+    d <- as.numeric(by_h(runs[[1]][[OUTCOMES[[o]]]], horizon)) - as.numeric(by_h(runs[[2]][[OUTCOMES[[o]]]], horizon))
+    data.table(world = par$world, scenario = par$scenario, outcome = o,
+               estimand = paste0("rct_", design, "_itt"), value = mean(d), mcse = sd(d) / sqrt(length(d)),
+               n = length(d))
+  }))
+}
@@ -0,0 +1,70 @@
+# What the analyst sees: simulated Medicare claims (plan D17).
+# Latent quantities (Z_t, true smoking, frailty, skill, style, progression rate) never
+# leave this function except in `hidden`, which analysis code must not touch.
+
+HIDDEN_COLS <- c("skill", "style", "frailty", "fi_noise", "smoking_true", "z0", "rate", "venous", "z")
+
+#' Turn a population and its simulated course into claims-like records.
+#'
+#' @return list(patients = one row per patient: measured baseline and coded event times;
+#'   pm = person-months (id, t, sev_obs) from month 0 to death or month 59;
+#'   hidden = latent traits, for oracle analyses and checks only).
+make_records <- function(pop, course, par) {
+  pat <- pop$pat; ev <- course$events; ms <- par$measurement
+  M <- par$size$months
+  smoking <- if (isTRUE(ms$smoking) || identical(ms$smoking, "true")) pat$smoking_true else pat$smoking_coded
+  rel <- ms$frailty_reliability
+  fi <- sqrt(rel) * pat$frailty + sqrt(1 - rel) * pat$fi_noise
+
+  patients <- data.table(id = pat$id, phys_id = pat$phys_id, obl = pat$obl, age = pat$age,
+                         female = pat$female, black = pat$black, year = pat$year,
+                         diabetes = pat$diabetes, ckd = pat$ckd, chf = pat$chf, copd = pat$copd,
+                         prior_cvd = pat$prior_cvd, smoking = smoking, frailty_index = fi)
+  patients <- cbind(patients, ev[, !"id"])
+  # t_clti_true is latent: only the coded date is in claims
+  patients[, t_clti_true := NULL]
+
+  # observed severity: exact Z_t (clean worlds) or a claims proxy (world W5, WR)
+  last <- pmin(ifelse(is.na(ev$t_death), M - 1L, ev$t_death), M - 1L)
+  nrep <- last + 1L
+  id <- rep(pat$id, nrep)
+  t <- sequence(nrep) - 1L
+  idx <- cbind(id, t + 1L)
+  if (identical(ms$severity, "exact")) {
+    sev <- course$z[idx]
+  } else {
+    w <- ms$proxy_window
+    cs <- t(apply(course$counts, 1, function(x) { x[is.na(x)] <- 0L; cumsum(x) }))
+    lagged <- cbind(matrix(0, nrow(cs), w), cs)[, seq_len(M)]
+    sev <- log1p(cs[idx] - lagged[idx])
+  }
+  pm <- data.table(id = id, t = t, sev_obs = sev)
+
+  hidden <- data.table(id = pat$id, skill = pat$skill, style = pat$style, frailty = pat$frailty,
+                       smoking_true = pat$smoking_true, rate = pat$rate,
+                       t_clti_true = ev$t_clti_true)
+  hidden_pm <- data.table(id = id, t = t, z = course$z[idx])
+  list(patients = patients, pm = pm, hidden = hidden, hidden_pm = hidden_pm, par = par)
+}
+
+#' Stop if an analysis table carries a latent column.
+assert_no_hidden <- function(dt) {
+  bad <- intersect(names(dt), HIDDEN_COLS)
+  if (length(bad)) stop("latent columns in analysis data: ", paste(bad, collapse = ", "))
+  invisible(TRUE)
+}
+
+#' One replicate's records: population, natural course, claims view.
+simulate_records <- function(par, rep = 1L, n = par$size$n_patients) {
+  seed <- replicate_seed(par$world, par$scenario, rep)
+  pop <- make_population(par, n = n, seed = seed)
+  course <- simulate_course(pop, par, policy("natural"), seed = seed + 104729L, record = TRUE)
+  make_records(pop, course, par)
+}
+
+#' Reproducible seed for a world x scenario x replicate.
+replicate_seed <- function(world, scenario, rep) {
+  key <- paste(world, scenario, rep, sep = "|")
+  h <- sum(utf8ToInt(key) * seq_along(utf8ToInt(key))^2)
+  as.integer((h * 2654435761) %% 2147483647)
+}
@@ -0,0 +1,221 @@
+# The simulated world: physicians, patients, and the monthly course of disease and care.
+# Analysis plan v1.0, section 5.1.
+#
+# Common random numbers. simulate_course() draws the same uniforms for every patient
+# in every month, in a fixed order, whatever the patient's state or the policy. So two
+# runs on the same population with the same seed differ only through the policy. This
+# keeps truths under different strategies paired, and makes the sharp null exact at the
+# patient level (tested).
+
+STREAMS <- c("eps", "proxy", "clti", "false_code", "amp", "bypass", "rr", "clti_revasc",
+             "nco1", "nco2", "nce", "pvi", "death")
+
+#' Physicians: latent skill k and style s (standard normal, correlated), setting, panel weight.
+make_physicians <- function(par, seed) {
+  set.seed(seed)
+  J <- par$size$n_physicians
+  rho <- par$physician$rho_skill_style
+  z1 <- rnorm(J); z2 <- rnorm(J)
+  s <- z1
+  k <- rho * z1 + sqrt(1 - rho^2) * z2
+  obl <- rbinom(J, 1, plogis(par$physician$obl_intercept + par$physician$obl_style * s))
+  w <- rlnorm(J, 0, par$size$panel_sdlog)
+  data.table(phys_id = seq_len(J), skill = k, style = s, obl = obl, weight = w / sum(w))
+}
+
+#' Patients at diagnosis: measured baseline, latent traits, and assigned physician.
+make_patients <- function(par, n, phys, seed) {
+  set.seed(seed)
+  pp <- par$patient
+  phys_id <- sample(phys$phys_id, n, replace = TRUE, prob = phys$weight)
+  age <- pp$age_min + rgamma(n, shape = pp$age_shape, scale = pp$age_scale)
+  frailty <- rnorm(n)
+  a10 <- (age - 76) / 10
+  comorb <- lapply(pp$comorbidity, function(cm)
+    rbinom(n, 1, plogis(qlogis(cm$prev) + cm$age_coef * a10 + cm$frailty_coef * frailty)))
+  smoking <- rbinom(n, 1, plogis(qlogis(pp$smoking$prev) + pp$smoking$age_coef * a10))
+  sv <- par$severity
+  z0 <- sv$z0_diabetes * comorb$diabetes + sv$z0_smoking * smoking + sv$z0_ckd * comorb$ckd +
+    rnorm(n, 0, sv$z0_sd)
+  rate <- sv$r0 + sv$r_diabetes * comorb$diabetes + sv$r_smoking * smoking + rnorm(n, 0, sv$r_sd)
+  venous <- rnorm(n)
+  year <- sample(pp$year_min:pp$year_max, n, replace = TRUE)
+  # what claims show for smoking and frailty (the analyst's versions; plan D17)
+  ms <- par$measurement
+  smoking_coded <- ifelse(smoking == 1, rbinom(n, 1, ms$smoking_sens), rbinom(n, 1, 1 - ms$smoking_spec))
+  fi_noise <- rnorm(n)
+  pat <- data.table(id = seq_len(n), phys_id = phys_id, age = age,
+                    female = rbinom(n, 1, pp$p_female), black = rbinom(n, 1, pp$p_black),
+                    year = year)
+  for (nm in names(comorb)) set(pat, j = nm, value = comorb[[nm]])
+  pat[, `:=`(smoking_true = smoking, smoking_coded = smoking_coded, frailty = frailty,
+             fi_noise = fi_noise, z0 = z0, rate = rate, venous = venous)]
+  pat <- merge(pat, phys[, .(phys_id, skill, style, obl)], by = "phys_id", sort = FALSE)
+  setorder(pat, id)
+  pat[]
+}
+
+#' A full population for one replicate (or for the truths).
+make_population <- function(par, n = par$size$n_patients, seed = 1L, phys = NULL) {
+  if (is.null(phys)) phys <- make_physicians(par, seed)
+  pat <- make_patients(par, n, phys, seed + 7919L)
+  list(phys = phys, pat = pat)
+}
+
+# ---- policies -----------------------------------------------------------------------
+# A policy overrides the natural elective-PVI decision among patients eligible in month t
+# (alive, no PVI yet, no MALE yet). Non-elective revascularization after CLTI is never
+# overridden: it happens after a MALE event, which every strategy allows.
+#   natural            observed practice
+#   A                  natural in months 0..(w-2); everyone still eligible gets PVI at w-1
+#   B1                 no PVI in months 0..(w-1); natural afterwards
+#   B2                 no elective PVI ever
+#   seq_force / seq_forbid (m)  natural except PVI forced / forbidden at month m
+
+#   custom             pol$decide(t, natural, idx): idx are the eligible patients' row
+#                      numbers, so the decision can differ by patient (simulated RCTs)
+
+policy <- function(type = "natural", m = NA_integer_, decide = NULL) {
+  stopifnot(type %in% c("natural", "A", "B1", "B2", "seq_force", "seq_forbid", "custom"))
+  list(type = type, m = m, decide = decide)
+}
+
+apply_policy <- function(pol, t, natural, window, idx = NULL) {
+  if (pol$type == "custom") return(pol$decide(t, natural, idx))
+  switch(pol$type,
+    natural = natural,
+    A = if (t < window - 1) natural else if (t == window - 1) rep(TRUE, length(natural)) else natural,
+    B1 = if (t < window) rep(FALSE, length(natural)) else natural,
+    B2 = rep(FALSE, length(natural)),
+    seq_force = if (t == pol$m) rep(TRUE, length(natural)) else natural,
+    seq_forbid = if (t == pol$m) rep(FALSE, length(natural)) else natural)
+}
+
+# ---- the monthly course -------------------------------------------------------------
+
+#' Simulate months 0..(months-1) for every patient under a policy.
+#'
+#' Within a month the order is: severity update; outcome events (CLTI, coding, amputation,
+#' bypass, repeat revascularization, revascularization after CLTI, negative controls);
+#' death; then the elective PVI decision among survivors still eligible. A MALE event
+#' and death in the same month count the MALE event first.
+#'
+#' @return list(pat = patient table with event times, z = N x months matrix of Z_t
+#'   (record = TRUE only), counts = N x months matrix of proxy counts, pvi_month = vector).
+simulate_course <- function(pop, par, pol = policy("natural"), seed = 1L, record = TRUE) {
+  pat <- pop$pat
+  N <- nrow(pat)
+  M <- par$size$months
+  W <- par$size$window
+  tr <- par$treatment; oc <- par$outcomes; nc <- par$negative_controls; ef <- par$effect
+  a10 <- (pat$age - 76) / 10
+
+  # time-fixed parts of each linear predictor
+  lp_win <- tr$window$intercept + tr$window$frailty * pat$frailty + tr$window$age10 * a10 +
+    tr$window$diabetes * pat$diabetes + tr$window$ckd * pat$ckd + tr$window$smoking * pat$smoking_true +
+    tr$window$style * pat$style + tr$window$obl * pat$obl + tr$window$year * (pat$year - par$patient$year_min)
+  lp_def <- tr$deferred$intercept + tr$deferred$style * pat$style + tr$deferred$obl * pat$obl
+  d <- oc$death
+  lp_death0 <- d$intercept + d$chf * pat$chf + d$ckd * pat$ckd + d$copd * pat$copd + d$frailty * pat$frailty
+  cl <- oc$clti
+  lp_clti0 <- cl$intercept + cl$diabetes * pat$diabetes + cl$ckd * pat$ckd + cl$smoking * pat$smoking_true +
+    cl$frailty * pat$frailty
+  lp_false <- oc$clti_false_code$intercept + oc$clti_false_code$style * pat$style
+  am <- oc$amputation
+  lp_amp0 <- am$intercept + am$diabetes * pat$diabetes + am$ckd * pat$ckd
+  rr <- oc$repeat_revasc
+  lp_rr0 <- rr$intercept + rr$skill * pat$skill + rr$style * pat$style
+  lp_nco1 <- nc$nco1$intercept + nc$nco1$frailty * pat$frailty + nc$nco1$smoking * pat$smoking_true +
+    nc$nco1$age10 * a10
+  lp_nco2 <- nc$nco2$intercept + nc$nco2$style * pat$style + nc$nco2$venous * pat$venous
+  lp_nce <- nc$nce$intercept + nc$nce$style * pat$style + nc$nce$frailty * pat$frailty
+  kappa_post <- ef$kappa * exp(-ef$skill_kappa * pat$skill)
+
+  # state
+  Z <- pat$z0
+  alive <- rep(TRUE, N); treated <- rep(FALSE, N); clti_true <- rep(FALSE, N)
+  clti_code <- rep(FALSE, N); male <- rep(FALSE, N)
+  t_pvi <- t_clti_true <- t_clti_code <- t_amp <- t_bypass <- t_rr <- t_male <- t_death <-
+    t_nco1 <- t_nco2 <- t_nce <- rep(NA_integer_, N)
+  pvi_elective <- rep(NA, N)
+  n_proc <- integer(N)
+  zmat <- if (record) matrix(NA_real_, N, M) else NULL
+  cmat <- if (record) matrix(NA_integer_, N, M) else NULL
+  ns <- length(STREAMS)
+
+  set.seed(seed)
+  for (t in 0:(M - 1)) {
+    U <- matrix(runif(N * ns), N, ns, dimnames = list(NULL, STREAMS))
+    a <- alive
+    if (t > 0) Z[a] <- Z[a] + pat$rate[a] + kappa_post[a] * treated[a] +
+      par$severity$eps_sd * qnorm(U[a, "eps"])
+    if (record) {
+      zmat[a, t + 1] <- Z[a]
+      cmat[a, t + 1] <- qpois(U[a, "proxy"], exp(par$measurement$proxy_intercept + par$measurement$proxy_z * Z[a]))
+    }
+
+    # true CLTI, and its coding (true onset always coded; false codes by physician style)
+    new_clti <- a & !clti_true & U[, "clti"] < plogis(lp_clti0 + cl$z * Z)
+    clti_true[new_clti] <- TRUE; t_clti_true[new_clti] <- t
+    new_code <- a & !clti_code & (new_clti | U[, "false_code"] < plogis(lp_false))
+    clti_code[new_code] <- TRUE; t_clti_code[new_code] <- t
+
+    new_amp <- a & is.na(t_amp) & U[, "amp"] < plogis(lp_amp0 + am$clti * clti_true + am$z * Z)
+    t_amp[new_amp] <- t
+    new_byp <- a & is.na(t_bypass) & U[, "bypass"] < plogis(oc$bypass$intercept + oc$bypass$z * Z)
+    t_bypass[new_byp] <- t
+    n_proc[new_byp] <- n_proc[new_byp] + 1L
+
+    # repeat revascularization: only after the index PVI (from the month after it)
+    new_rr <- a & treated & !is.na(t_pvi) & t > t_pvi & U[, "rr"] < plogis(lp_rr0 + rr$z * Z)
+    t_rr[new_rr & is.na(t_rr)] <- t
+    n_proc[new_rr] <- n_proc[new_rr] + 1L
+
+    # revascularization after a CLTI diagnosis (not elective; allowed by every strategy)
+    cr <- a & clti_code & U[, "clti_revasc"] < plogis(oc$clti_revasc$intercept)
+    first <- cr & !treated
+    rep_ <- cr & treated & !new_rr
+    t_rr[rep_ & is.na(t_rr)] <- t
+    n_proc[cr & !new_rr] <- n_proc[cr & !new_rr] + 1L
+    if (any(first)) {
+      treated[first] <- TRUE; t_pvi[first] <- t; pvi_elective[first] <- FALSE
+      Z[first] <- Z[first] - ef$relief
+    }
+
+    # negative controls (never depend on PVI)
+    x <- a & is.na(t_nco1) & U[, "nco1"] < plogis(lp_nco1); t_nco1[x] <- t
+    x <- a & is.na(t_nco2) & U[, "nco2"] < plogis(lp_nco2); t_nco2[x] <- t
+    if (t < W) { x <- a & is.na(t_nce) & U[, "nce"] < plogis(lp_nce); t_nce[x] <- t }
+
+    # MALE: coded CLTI, amputation, bypass, or repeat revascularization
+    new_male <- a & !male & (clti_code | !is.na(t_amp) | !is.na(t_bypass) | !is.na(t_rr))
+    male[new_male] <- TRUE; t_male[new_male] <- t
+
+    # death
+    dies <- a & U[, "death"] < plogis(lp_death0 + d$age10 * a10 + d$age10 * t / 120 + d$z * Z + d$clti * clti_true)
+    alive[dies] <- FALSE; t_death[dies] <- t
+
+    # elective PVI among survivors still eligible
+    elig <- alive & !treated & !male
+    if (any(elig)) {
+      lp <- if (t < W) lp_win[elig] + tr$window$z * Z[elig] else lp_def[elig] + tr$deferred$z * Z[elig]
+      natural <- U[elig, "pvi"] < plogis(lp)
+      go <- apply_policy(pol, t, natural, W, which(elig))
+      idx <- which(elig)[go]
+      if (length(idx)) {
+        treated[idx] <- TRUE; t_pvi[idx] <- t; pvi_elective[idx] <- TRUE
+        n_proc[idx] <- n_proc[idx] + 1L
+        Z[idx] <- Z[idx] - ef$relief
+      }
+    }
+  }
+
+  # disease composite (D34): coded CLTI or major amputation
+  t_disease <- pmin(t_clti_code, t_amp, na.rm = TRUE)
+  out <- data.table(id = pat$id, t_pvi, pvi_elective, t_clti_true, t_clti_code, t_amp, t_bypass,
+                    t_rr, t_male, t_disease, t_death, t_nco1, t_nco2, t_nce, n_proc)
+  list(events = out, z = zmat, counts = cmat, policy = pol)
+}
+
+#' Cumulative incidence by month `h` (event in months 0..h-1, before death or in the same month).
+cum_inc <- function(t_event, h = 60L) mean(!is.na(t_event) & t_event < h)
@@ -0,0 +1,92 @@
+# Truths by simulating the interventions in a large population (plan section 5.5).
+#
+# The truth population keeps the replicate's panel structure (patients per physician)
+# but has proportionally more physicians, so the truth averages over the physician
+# distribution rather than over one draw of 1,500 physicians. Every policy is run on the
+# same population with the same course seed (common random numbers), so risk
+# differences are paired and their Monte Carlo SE is sd(paired difference) / sqrt(n).
+#
+# Landmark, naive and RCT estimands are not interventional in this sense; they are added
+# with the estimators (oracle large-sample limits).
+
+OUTCOMES <- c(male = "t_male", disease = "t_disease", clti = "t_clti_code", amp = "t_amp",
+              bypass = "t_bypass", rr = "t_rr", death = "t_death", nco1 = "t_nco1", nco2 = "t_nco2")
+
+truth_population <- function(par, n, seed) {
+  scale <- n / par$size$n_patients
+  par$size$n_physicians <- max(1L, round(par$size$n_physicians * scale))
+  make_population(par, n = n, seed = seed)
+}
+
+# event indicator by month h (event in months 0..h-1)
+by_h <- function(t, h) !is.na(t) & t < h
+
+#' Truths for one world x scenario.
+#'
+#' @return list(truths = long table: outcome, estimand, value, mcse, n;
+#'   curves = cumulative incidence by month under natural, A, B1, B2 for MALE and disease).
+compute_truths <- function(world, scenario, n = 1e6, seed = 20261027L, horizon = 60L,
+                           window = NULL, verbose = TRUE) {
+  par <- load_params(world, scenario)
+  W <- if (is.null(window)) par$size$window else window
+  par$size$window <- W
+  t0 <- Sys.time()
+  pop <- truth_population(par, n, seed)
+  cseed <- seed + 1L
+  run <- function(p) simulate_course(pop, par, p, seed = cseed, record = FALSE)$events
+
+  ev <- list(natural = run(policy("natural")), A = run(policy("A")),
+             B1 = run(policy("B1")), B2 = run(policy("B2")))
+  rows <- list()
+  add <- function(outcome, estimand, d, nn = length(d)) {
+    rows[[length(rows) + 1L]] <<- data.table(world = world, scenario = scenario, outcome = outcome,
+                                             estimand = estimand, value = mean(d),
+                                             mcse = sd(d) / sqrt(nn), n = nn)
+  }
+  for (o in names(OUTCOMES)) {
+    col <- OUTCOMES[[o]]
+    y <- lapply(ev, function(e) as.numeric(by_h(e[[col]], horizon)))
+    for (p in names(y)) add(o, paste0("risk_", p), y[[p]])
+    add(o, "rd_A_B1", y$A - y$B1)
+    add(o, "rd_A_B2", y$A - y$B2)
+  }
+  # total procedures per patient (secondary outcome, D7)
+  np <- lapply(ev, function(e) e$n_proc)
+  for (p in names(np)) add("n_proc", paste0("mean_", p), np[[p]])
+  add("n_proc", "diff_A_B1", np$A - np$B1)
+
+  # sequential trials, ITT (D14): among patients eligible at month m under the natural
+  # course, PVI forced vs forbidden at m, natural afterwards; risk by month `horizon`
+  # since diagnosis; pooled over m = 0..W-1 weighted by the number eligible.
+  nat <- ev$natural
+  seq_rows <- list()
+  for (m in 0:(W - 1)) {
+    elig <- (is.na(nat$t_death) | nat$t_death > m) & (is.na(nat$t_pvi) | nat$t_pvi >= m) &
+      (is.na(nat$t_male) | nat$t_male > m)
+    f <- run(policy("seq_force", m)); b <- run(policy("seq_forbid", m))
+    for (o in names(OUTCOMES)) {
+      col <- OUTCOMES[[o]]
+      d <- as.numeric(by_h(f[[col]][elig], horizon) & f[[col]][elig] > m) -
+        as.numeric(by_h(b[[col]][elig], horizon) & b[[col]][elig] > m)
+      seq_rows[[length(seq_rows) + 1L]] <- data.table(outcome = o, m = m, n_elig = sum(elig),
+                                                      rd = mean(d), var = var(d))
+    }
+  }
+  seq_dt <- rbindlist(seq_rows)
+  pooled <- seq_dt[, .(value = sum(rd * n_elig) / sum(n_elig),
+                       mcse = sqrt(sum(n_elig^2 * var / n_elig)) / sum(n_elig),
+                       n = sum(n_elig)), by = outcome]
+  rows[[length(rows) + 1L]] <- pooled[, .(world, scenario, outcome, estimand = "seq_itt", value, mcse, n)]
+  for (mm in 0:(W - 1)) {
+    x <- seq_dt[m == mm]
+    rows[[length(rows) + 1L]] <- x[, .(world, scenario, outcome, estimand = paste0("seq_itt_m", mm),
+                                       value = rd, mcse = sqrt(var / n_elig), n = n_elig)]
+  }
+
+  curves <- rbindlist(lapply(names(ev), function(p) rbindlist(lapply(c("male", "disease"), function(o)
+    data.table(world = world, scenario = scenario, policy = p, outcome = o, month = seq_len(horizon),
+               cuminc = vapply(seq_len(horizon), function(h) mean(by_h(ev[[p]][[OUTCOMES[[o]]]], h)), 0))))))
+  if (verbose) message(sprintf("%s/%s: %.1f min", world, scenario,
+                               as.numeric(Sys.time() - t0, units = "mins")))
+  list(truths = rbindlist(rows), curves = curves)
+}
@@ -0,0 +1,126 @@
+# Every number that defines the simulated world. Analysis plan v1.0, section 5.
+# Time unit: month. Hazards are monthly, on the logit scale.
+# Starting values are placeholders, to be tuned by calibration
+# (simulations/01_calibrate.R, targets in plan section 5.4).
+# Physician traits (skill k, style s) and patient frailty F are standard normal;
+# their influence is set by the coefficients that multiply them.
+
+size:
+  n_patients: 50000
+  n_physicians: 1500
+  panel_sdlog: 1.0          # lognormal panel weights; median panel about 20 at n = 50,000
+  months: 60                # follow-up months 0-59
+  window: 6                 # strategy A: PVI in months 0-5
+
+physician:
+  rho_skill_style: -0.3     # corr(k, s). Has no effect unless skill matters (world W1)
+  obl_intercept: 0.0        # logit P(office-based lab / ASC); about 50% of physicians
+  obl_style: 0.8            # high-style physicians more often work in OBL/ASC
+
+patient:
+  age_min: 66
+  age_shape: 2.5            # age = 66 + Gamma(shape, scale); mean about 76
+  age_scale: 4.0
+  p_female: 0.50
+  p_black: 0.10
+  # comorbidities: logit P = qlogis(prev) + age_coef * (age - 76) / 10 + frailty_coef * F
+  comorbidity:
+    diabetes:  {prev: 0.40, age_coef: -0.10, frailty_coef: 0.4}
+    ckd:       {prev: 0.20, age_coef:  0.30, frailty_coef: 0.5}
+    chf:       {prev: 0.20, age_coef:  0.40, frailty_coef: 0.6}
+    copd:      {prev: 0.25, age_coef:  0.10, frailty_coef: 0.4}
+    prior_cvd: {prev: 0.45, age_coef:  0.20, frailty_coef: 0.3}
+  smoking: {prev: 0.25, age_coef: -0.40}
+  year_min: 2015
+  year_max: 2019
+
+severity:                    # latent Z_t: walking limitation / ABI on a standard scale
+  z0_diabetes: 0.3
+  z0_smoking: 0.3
+  z0_ckd: 0.2
+  z0_sd: 1.0
+  r0: 0.02                   # mean drift per month
+  r_diabetes: 0.01
+  r_smoking: 0.01
+  r_sd: 0.01                 # patient-specific progression rate
+  eps_sd: 0.05               # month-to-month noise
+
+treatment:
+  window:                    # elective PVI hazard, months 0-5, among the eligible
+    intercept: -5.4
+    z: 0.9                   # confounding by indication
+    frailty: -0.3            # frail patients are less often offered PVI
+    age10: -0.3
+    diabetes: 0.1
+    ckd: -0.1
+    smoking: 0.1
+    style: 0.5               # physician style; world W2 sets it to the strong level
+    obl: 0.3
+    year: 0.05               # per calendar year after year_min
+  deferred:                  # elective PVI hazard from month 6 on (deferral / crossover)
+    intercept: -7.2
+    z: 0.0                   # world W4 makes deferral severity-driven
+    style: 0.5
+    obl: 0.3
+
+effect:                      # what a PVI does. Set by the truth scenario
+  relief: 0.0                # immediate drop in Z at PVI
+  kappa: 0.0                 # extra Z drift per month after PVI
+  skill_kappa: 0.0           # kappa * exp(-skill_kappa * k): skill slows post-PVI progression
+
+outcomes:
+  death:  {intercept: -5.0, age10: 0.7, chf: 0.4, ckd: 0.4, copd: 0.3, frailty: 0.4, z: 0.1, clti: 0.8}
+  clti:   {intercept: -6.6, z: 0.6, diabetes: 0.4, ckd: 0.3, smoking: 0.3, frailty: 0.2}
+  clti_false_code: {intercept: -9.0, style: 0.0}     # factor 3 sets style
+  amputation: {intercept: -8.5, clti: 3.0, diabetes: 0.4, ckd: 0.4, z: 0.2}
+  bypass: {intercept: -7.8, z: 0.5}
+  repeat_revasc: {intercept: -5.3, skill: 0.0, style: 0.0, z: 0.2}   # only after the index PVI
+  clti_revasc: {intercept: -4.5}                      # PVI after CLTI (not elective); most CLTI patients get bypass, tibial work, amputation or nothing
+
+negative_controls:           # never depend on PVI
+  nco1: {intercept: -6.4, frailty: 0.5, smoking: 0.4, age10: 0.3}   # hip fracture / pneumonia
+  nco2: {intercept: -7.1, style: 0.8, venous: 0.5}                  # varicose-vein ablation
+  nce:  {intercept: -5.7, style: 0.8, frailty: -0.2}                # low-value procedure, months 0-5
+
+measurement:                 # what the analyst sees. World W5 switches to claims proxies
+  severity: exact            # exact: analyst sees Z_t. proxy: only visit/testing counts
+  proxy_intercept: -1.0      # monthly count ~ Poisson(exp(proxy_intercept + proxy_z * Z_t))
+  proxy_z: 0.5
+  proxy_window: 3            # observed severity = log1p(sum of counts over the last 3 months)
+  frailty_reliability: 1.0   # frailty index = sqrt(rel) * F + sqrt(1 - rel) * noise
+  smoking: true              # true: analyst sees true smoking. coded: sensitivity / specificity below
+  smoking_sens: 0.6
+  smoking_spec: 0.98
+
+worlds:                      # one factor at a time (plan section 5.2, D24). Overrides on the default.
+                             # Intercepts are then re-solved per world by calibration, so every world
+                             # matches the same Medicare targets unless `calibration_targets` says otherwise.
+  W0: {}
+  W1:
+    effect: {skill_kappa: 0.5}
+    outcomes: {repeat_revasc: {skill: -0.6}}
+  W1_rho0:
+    physician: {rho_skill_style: 0.0}
+    effect: {skill_kappa: 0.5}
+    outcomes: {repeat_revasc: {skill: -0.6}}
+  W2:
+    treatment: {window: {style: 1.6}, deferred: {style: 1.6}}
+  W3:
+    treatment: {window: {style: 1.6}, deferred: {style: 1.6}}
+    outcomes: {clti_false_code: {style: 0.8}, repeat_revasc: {style: 0.6}}
+  W4:                        # deferral common (as in trial conservative arms) and severity-driven
+    treatment: {deferred: {z: 0.8}}
+    calibration_targets: {late_not: 0.15}
+  W5:
+    measurement: {severity: proxy, frailty_reliability: 0.4, smoking: coded}
+  WR:
+    effect: {skill_kappa: 0.5}
+    treatment: {window: {style: 1.6}, deferred: {style: 1.6, z: 0.8}}   # deferral at the Sorber rate, but severity-driven
+    outcomes: {repeat_revasc: {skill: -0.6, style: 0.6}, clti_false_code: {style: 0.8}}
+    measurement: {severity: proxy, frailty_reliability: 0.4, smoking: coded}
+
+scenarios:                   # truth scenarios (plan section 5.3, D34)
+  "null": {effect: {relief: 0.0, kappa: 0.0}}           # sharp disease null
+  harm: {effect: {relief: 0.3, kappa: 0.039}}         # W0: disease-composite RD(A-B1) +5.1 pp at 60 months
+  benefit_harm: {effect: {relief: 2.0, kappa: 0.075}} # W0: RD -1.2 pp at 24 months, +3.1 pp at 60
+  # tuned in W0 with simulations/explore_scenarios.R (2e5 patients); the same mechanism is used in every world
@@ -0,0 +1,55 @@
+test_that("every world and scenario loads", {
+  raw <- yaml::read_yaml(default_params_path())
+  for (w in names(raw$worlds)) for (s in names(raw$scenarios)) {
+    par <- load_params(w, s)
+    expect_equal(par$world, w)
+    expect_true(is.numeric(par$effect$relief))
+  }
+  expect_equal(load_params("W2")$treatment$window$style, 1.6)
+  expect_equal(load_params("W2")$treatment$window$z, 0.9)   # untouched by the override
+})
+
+small_world <- function(world = "WR", scenario = "null", n = 4000) {
+  par <- load_params(world, scenario, overrides = list(size = list(n_physicians = 150)))
+  list(par = par, pop = make_population(par, n = n, seed = 11L))
+}
+
+test_that("under the sharp null, disease, death and negative controls are identical under every policy", {
+  sw <- small_world("WR", "null")
+  runs <- lapply(list(policy("natural"), policy("A"), policy("B1"), policy("B2"),
+                      policy("seq_force", 3L)),
+                 function(p) simulate_course(sw$pop, sw$par, p, seed = 5L, record = FALSE)$events)
+  for (r in runs[-1]) {
+    for (col in c("t_disease", "t_clti_code", "t_amp", "t_death", "t_nco1", "t_nco2", "t_nce")) {
+      expect_identical(r[[col]], runs[[1]][[col]], info = col)
+    }
+  }
+  # MALE is not null: repeat revascularization needs an index PVI (D34)
+  expect_gt(cum_inc(runs[[2]]$t_male), cum_inc(runs[[3]]$t_male))
+})
+
+test_that("policies are followed", {
+  sw <- small_world("W0", "harm")
+  W <- sw$par$size$window
+  a <- simulate_course(sw$pop, sw$par, policy("A"), seed = 5L, record = FALSE)$events
+  # under A, everyone who survives month W-1 without a MALE event has had a PVI by then
+  # (death and MALE in a month come before that month's PVI decision)
+  free <- (is.na(a$t_death) | a$t_death > W - 1) & (is.na(a$t_male) | a$t_male > W - 1 | (!is.na(a$t_pvi) & a$t_pvi < a$t_male))
+  expect_true(all(!is.na(a$t_pvi[free]) & a$t_pvi[free] <= W - 1))
+  b1 <- simulate_course(sw$pop, sw$par, policy("B1"), seed = 5L, record = FALSE)$events
+  expect_false(any(b1$pvi_elective %in% TRUE & b1$t_pvi < W))
+  b2 <- simulate_course(sw$pop, sw$par, policy("B2"), seed = 5L, record = FALSE)$events
+  expect_false(any(b2$pvi_elective %in% TRUE))
+})
+
+test_that("records carry no latent columns and person-months end at death", {
+  par <- load_params("W5", "harm", overrides = list(size = list(n_physicians = 100)))
+  rec <- simulate_records(par, rep = 1L, n = 2000)
+  expect_true(assert_no_hidden(rec$patients))
+  expect_true(assert_no_hidden(rec$pm))
+  last <- rec$pm[, .(tmax = max(t)), by = id]
+  dead <- rec$patients[!is.na(t_death), .(id, t_death)]
+  chk <- merge(last, dead, by = "id")
+  expect_true(all(chk$tmax == chk$t_death))
+  expect_true(all(is.finite(rec$pm$sev_obs)))
+})
@@ -0,0 +1,80 @@
+# Calibrate each world's natural course to the Medicare targets (plan section 5.4).
+#
+# For every world, under the null scenario, solve one intercept per target, cycling
+# until all targets hold. Every world gets the same targets, so worlds differ only in
+# the threat they switch on, not in base rates.
+#
+# Usage: Rscript simulations/01_calibrate.R [n] [world ...]
+# Writes simulations/config/calibrated.yml (merged over existing entries).
+
+pkgload::load_all("pvisim", quiet = TRUE)
+suppressPackageStartupMessages(library(data.table))
+
+args <- commandArgs(TRUE)
+n <- if (length(args) && grepl("^[0-9e.]+$", args[1])) as.numeric(args[1]) else 1e5
+worlds <- setdiff(args, as.character(args[1]))
+if (!length(worlds)) worlds <- names(yaml::read_yaml(default_params_path())$worlds)
+H <- 53L   # Sorber 2023 median follow-up, about 4.4 years
+
+# Natural-course statistics that the targets are defined on.
+stats_of <- function(par) {
+  pop <- make_population(par, n = n, seed = 20261027L)
+  ev <- simulate_course(pop, par, policy("natural"), seed = 20261028L, record = FALSE)$events
+  early <- ev$pvi_elective %in% TRUE & ev$t_pvi < par$size$window
+  late_early <- !is.na(ev$t_rr) & ev$t_rr < H
+  late_not <- !is.na(ev$t_pvi) & ev$t_pvi >= par$size$window & ev$t_pvi < H
+  list(early = mean(early),
+       late_not = mean(late_not[!early]),
+       late_early = mean(late_early[early]),
+       clti_not = mean((!is.na(ev$t_clti_code) & ev$t_clti_code < H)[!early]),
+       death = cum_inc(ev$t_death))
+}
+
+# target name -> parameter path and value (plan 5.4)
+targets <- list(
+  early      = list(path = c("treatment", "window", "intercept"), value = 0.032),
+  late_not   = list(path = c("treatment", "deferred", "intercept"), value = 0.036),
+  late_early = list(path = c("outcomes", "repeat_revasc", "intercept"), value = 0.225),
+  clti_not   = list(path = c("outcomes", "clti", "intercept"), value = 0.078),
+  death      = list(path = c("outcomes", "death", "intercept"), value = 0.35)
+)
+
+set_path <- function(lst, path, value) {
+  if (length(path) == 1) { lst[[path]] <- value; return(lst) }
+  lst[[path[1]]] <- set_path(if (is.null(lst[[path[1]]])) list() else lst[[path[1]]], path[-1], value)
+  lst
+}
+get_path <- function(lst, path) Reduce(function(x, p) x[[p]], path, lst)
+
+# Each target depends mainly on its own intercept, and for rare monthly events
+# logit(rate) moves about one-for-one with the intercept. So update all intercepts
+# together: b <- b - damp * (logit(stat) - logit(target)), until every target is within
+# `tol` on the logit scale. Common random numbers make each evaluation deterministic.
+calibrate_world <- function(world, max_iter = 25, tol = 0.02, damp = 0.9) {
+  base <- load_params(world, "null", calibrated = FALSE)
+  for (tg in names(base$calibration_targets)) targets[[tg]]$value <- base$calibration_targets[[tg]]
+  ov <- list()
+  for (tg in names(targets)) ov <- set_path(ov, targets[[tg]]$path, get_path(base, targets[[tg]]$path))
+  for (it in seq_len(max_iter)) {
+    s <- stats_of(load_params(world, "null", overrides = ov, calibrated = FALSE))
+    gap <- sapply(names(targets), function(tg) qlogis(s[[tg]]) - qlogis(targets[[tg]]$value))
+    cat(sprintf("%s iter %2d: %s\n", world, it,
+                paste(sprintf("%s %.3f", names(s), unlist(s)), collapse = "  ")))
+    if (all(abs(gap) < tol)) break
+    for (tg in names(targets)) {
+      b <- get_path(ov, targets[[tg]]$path)
+      ov <- set_path(ov, targets[[tg]]$path, round(b - damp * gap[[tg]], 4))
+    }
+  }
+  if (any(abs(gap) >= tol)) warning(world, ": not converged; largest logit gap ", round(max(abs(gap)), 3))
+  ov
+}
+
+out_path <- "simulations/config/calibrated.yml"
+cal <- if (file.exists(out_path)) yaml::read_yaml(out_path) else list()
+for (w in worlds) {
+  t0 <- Sys.time()
+  cal[[w]] <- calibrate_world(w)
+  cat(sprintf("%s done in %.1f min\n", w, as.numeric(Sys.time() - t0, units = "mins")))
+  yaml::write_yaml(cal, out_path)
+}
@@ -0,0 +1,20 @@
+# Interventional truths for every world x scenario (plan section 5.5).
+# Usage: Rscript simulations/02_truths.R [n] [cores]
+# Writes simulations/truths/truths.csv and simulations/truths/risk_curves.csv.
+
+pkgload::load_all("pvisim", quiet = TRUE)
+suppressPackageStartupMessages(library(data.table))
+args <- commandArgs(TRUE)
+n <- if (length(args) >= 1) as.numeric(args[1]) else 1e6
+cores <- if (length(args) >= 2) as.integer(args[2]) else 6L
+
+raw <- yaml::read_yaml(default_params_path())
+jobs <- CJ(world = names(raw$worlds), scenario = names(raw$scenarios), sorted = FALSE)
+res <- parallel::mclapply(seq_len(nrow(jobs)), function(i)
+  compute_truths(jobs$world[i], jobs$scenario[i], n = n), mc.cores = cores)
+bad <- vapply(res, inherits, TRUE, "try-error")
+if (any(bad)) stop("failed: ", paste(jobs$world[bad], jobs$scenario[bad], collapse = "; "), "\n", res[bad][[1]])
+dir.create("simulations/truths", showWarnings = FALSE)
+fwrite(rbindlist(lapply(res, `[[`, "truths")), "simulations/truths/truths.csv")
+fwrite(rbindlist(lapply(res, `[[`, "curves")), "simulations/truths/risk_curves.csv")
+cat("wrote simulations/truths/truths.csv\n")
@@ -0,0 +1,20 @@
+# RCT ITT truths for every world x scenario (the landmark and naive own estimands are
+# estimated within replicates by run_oracle()). Writes simulations/truths/truths_rct.csv.
+# Usage: Rscript simulations/03_limits.R [n] [cores]
+pkgload::load_all("pvisim", quiet = TRUE)
+suppressPackageStartupMessages(library(data.table))
+a <- commandArgs(TRUE)
+n <- if (length(a) >= 1) as.numeric(a[1]) else 5e5
+cores <- if (length(a) >= 2) as.integer(a[2]) else 6L
+raw <- yaml::read_yaml(default_params_path())
+jobs <- CJ(world = names(raw$worlds), scenario = names(raw$scenarios), sorted = FALSE)
+res <- parallel::mclapply(seq_len(nrow(jobs)), function(i) {
+  w <- jobs$world[i]; s <- jobs$scenario[i]
+  par <- load_params(w, s)
+  rbind(rct_truth(par, "realistic", n = 1e6)[outcome %in% c("male", "disease", "nco1", "nco2")],
+        rct_truth(par, "ideal", n = 1e6)[outcome %in% c("male", "disease", "nco1", "nco2")])
+}, mc.cores = cores)
+bad <- vapply(res, inherits, TRUE, "try-error")
+if (any(bad)) stop(res[bad][[1]])
+fwrite(rbindlist(res), "simulations/truths/truths_rct.csv")
+cat("wrote simulations/truths/truths_rct.csv\n")
@@ -0,0 +1,104 @@
+W0:
+  treatment:
+    window:
+      intercept: -6.3552
+    deferred:
+      intercept: -7.7172
+  outcomes:
+    repeat_revasc:
+      intercept: -5.6375
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6504
+W4:
+  treatment:
+    window:
+      intercept: -6.3535
+    deferred:
+      intercept: -6.7744
+  outcomes:
+    repeat_revasc:
+      intercept: -5.6429
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6542
+W1:
+  treatment:
+    window:
+      intercept: -6.3552
+    deferred:
+      intercept: -7.7172
+  outcomes:
+    repeat_revasc:
+      intercept: -5.7951
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6504
+W1_rho0:
+  treatment:
+    window:
+      intercept: -6.3552
+    deferred:
+      intercept: -7.7172
+  outcomes:
+    repeat_revasc:
+      intercept: -5.7187
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6504
+W2:
+  treatment:
+    window:
+      intercept: -7.3097
+    deferred:
+      intercept: -8.5615
+  outcomes:
+    repeat_revasc:
+      intercept: -5.5855
+    clti:
+      intercept: -7.4913
+    death:
+      intercept: -5.6505
+W3:
+  treatment:
+    window:
+      intercept: -7.3134
+    deferred:
+      intercept: -8.5848
+  outcomes:
+    repeat_revasc:
+      intercept: -6.4032
+    clti:
+      intercept: -7.5232
+    death:
+      intercept: -5.649
+W5:
+  treatment:
+    window:
+      intercept: -6.3552
+    deferred:
+      intercept: -7.7172
+  outcomes:
+    repeat_revasc:
+      intercept: -5.6375
+    clti:
+      intercept: -7.4903
+    death:
+      intercept: -5.6504
+WR:
+  treatment:
+    window:
+      intercept: -7.3131
+    deferred:
+      intercept: -9.598
+  outcomes:
+    repeat_revasc:
+      intercept: -6.8121
+    clti:
+      intercept: -7.5233
+    death:
+      intercept: -5.6489
@@ -0,0 +1,12 @@
+# Experiment grids (plan sections 5.2, 7, 9). Run with simulations/run_experiment.R <grid> [--scale draft|final]
+main:
+  worlds: [W0, W1, W1_rho0, W2, W3, W4, W5, WR]
+  scenarios: ["null", harm, benefit_harm]
+  designs: [naive, landmark, ccw, seq]
+  variants: [V1_ignore, V3_random, V4_within]     # V1 also yields V2 (cluster-robust SE)
+  nce: true                                        # negative-control exposure (landmark, CCW), D32/D36
+  rct: [realistic, ideal]
+scales:
+  smoke: {n_patients: 20000, reps: 1}
+  draft: {n_patients: 50000, reps: 100}
+  final: {n_patients: 50000, reps: 500}                     # to be revisited after the draft (sample size, D25/D29)
@@ -0,0 +1,17 @@
+# Explore effect sizes: RD(A - B1) by month for MALE and the disease composite in W0.
+pkgload::load_all("pvisim", quiet = TRUE)
+suppressMessages(library(data.table))
+args <- commandArgs(TRUE)
+grid <- if (length(args)) eval(parse(text = args[1])) else
+  list(c(0, 0), c(0.3, 0.03), c(1.2, 0.04))
+par0 <- load_params("W0", "null")
+pop <- make_population(par0, n = 2e5, seed = 9L)
+ci <- function(t, h) mean(!is.na(t) & t < h)
+for (g in grid) {
+  par <- load_params("W0", "null", overrides = list(effect = list(relief = g[1], kappa = g[2])))
+  ev <- lapply(c("A", "B1"), function(p) simulate_course(pop, par, policy(p), seed = 10L, record = FALSE)$events)
+  rd <- function(col, h) 100 * (ci(ev[[1]][[col]], h) - ci(ev[[2]][[col]], h))
+  cat(sprintf("relief %.2f kappa %.3f | disease RD @12 %+.1f @24 %+.1f @60 %+.1f | MALE RD @60 %+.1f (A %.1f%% B1 %.1f%%) | death RD %+.2f\n",
+              g[1], g[2], rd("t_disease", 12), rd("t_disease", 24), rd("t_disease", 60), rd("t_male", 60),
+              100 * ci(ev[[1]]$t_male, 60), 100 * ci(ev[[2]]$t_male, 60), rd("t_death", 60)))
+}
@@ -0,0 +1,26 @@
+# Quick look: the natural course of each world against the calibration targets (plan 5.4).
+pkgload::load_all("pvisim", quiet = TRUE)
+library(data.table)
+args <- commandArgs(TRUE); worlds <- if (length(args)) args else c("W0", "WR")
+for (w in worlds) {
+  par <- load_params(w, "harm")
+  t0 <- Sys.time()
+  pop <- make_population(par, seed = 1L)
+  ev <- simulate_course(pop, par, policy("natural"), seed = 2L, record = FALSE)$events
+  ev <- merge(ev, pop$pat[, .(id, phys_id)], by = "id")
+  early <- !is.na(ev$t_pvi) & ev$t_pvi < 6 & ev$pvi_elective %in% TRUE
+  h <- 53  # Sorber median follow-up about 4.4 years
+  late <- !is.na(ev$t_pvi) & ev$t_pvi >= 6 | !is.na(ev$t_rr) & ev$t_rr >= 6
+  late <- late & pmin(ev$t_pvi * NA^(ev$t_pvi < 6), ifelse(!is.na(ev$t_rr) & ev$t_rr >= 6, ev$t_rr, NA), na.rm = TRUE) < h
+  clti <- !is.na(ev$t_clti_code) & ev$t_clti_code < h
+  phys <- ev[, .(n = .N, rate = mean(!is.na(t_pvi) & t_pvi < 6 & pvi_elective %in% TRUE)), by = phys_id]
+  sdr <- sd(phys$rate); out <- mean(phys$rate >= mean(phys$rate) + 2 * sdr)
+  cat(sprintf("\n%s (%.0fs)\n  early PVI %.1f%% (target 3.2)\n  physician rate: median %.1f%%, max %.0f%%, >=2SD outliers %.1f%% (0, 58, 5.6)\n",
+              w, as.numeric(Sys.time() - t0, units = "secs"), 100 * mean(early),
+              100 * median(phys$rate), 100 * max(phys$rate), 100 * out))
+  cat(sprintf("  late PVI by 4.4y: early %.1f%% vs not %.1f%% (22.5 vs 3.6)\n", 100 * mean(late[early]), 100 * mean(late[!early])))
+  cat(sprintf("  coded CLTI by 4.4y: early %.1f%% vs not %.1f%% (16.4 vs 7.8)\n", 100 * mean(clti[early]), 100 * mean(clti[!early])))
+  cat(sprintf("  death by 5y %.1f%% (about 35)   MALE by 5y %.1f%%   amputation %.2f%%   disease composite %.1f%%\n",
+              100 * cum_inc(ev$t_death), 100 * cum_inc(ev$t_male), 100 * cum_inc(ev$t_amp), 100 * cum_inc(ev$t_disease)))
+  cat(sprintf("  NCO1 %.1f%%  NCO2 %.1f%%  NCE %.1f%%\n", 100 * cum_inc(ev$t_nco1), 100 * cum_inc(ev$t_nco2), 100 * cum_inc(ev$t_nce)))
+}
@@ -0,0 +1,60 @@
+# Run an experiment grid; one CSV per world x scenario x replicate, so runs are resumable.
+# Usage: Rscript simulations/run_experiment.R main [--scale draft] [--cores 10] [--worlds W0,W2]
+pkgload::load_all("pvisim", quiet = TRUE)
+suppressPackageStartupMessages(library(data.table))
+a <- commandArgs(TRUE)
+opt <- function(flag, default) { i <- match(flag, a); if (is.na(i)) default else a[i + 1] }
+grid_name <- a[1]
+scale <- opt("--scale", "draft")
+cores <- as.integer(opt("--cores", "10"))
+cfg <- yaml::read_yaml("simulations/config/experiments.yml")
+g <- cfg[[grid_name]]; sc <- cfg$scales[[scale]]
+worlds <- if (is.null(opt("--worlds", NULL))) g$worlds else strsplit(opt("--worlds", NULL), ",")[[1]]
+outdir <- file.path("results", grid_name, scale)
+jobs <- CJ(world = worlds, scenario = unlist(g$scenarios), rep = seq_len(sc$reps), sorted = FALSE)
+jobs[, file := file.path(outdir, world, scenario, sprintf("rep%03d.csv", rep))]
+todo <- jobs[!file.exists(file)]
+cat(sprintf("%s/%s: %d jobs, %d to do, %d cores\n", grid_name, scale, nrow(jobs), nrow(todo), cores))
+
+one_job <- function(i) {
+  j <- todo[i]
+  par <- load_params(j$world, j$scenario)
+  rec <- simulate_records(par, rep = j$rep, n = sc$n_patients)
+  out <- list()
+  for (v in g$variants) {
+    for (d in g$designs) {
+      r <- run_design(rec, d, v)
+      if (v == "V1_ignore") {   # same fit, two SEs
+        r2 <- copy(r)[, `:=`(variant = "V2_cluster", se = se_cluster)]
+        r <- rbind(r, r2)
+      }
+      out[[length(out) + 1]] <- r[, analysis := "main"]
+    }
+    if (isTRUE(g$nce)) for (d in c("landmark", "ccw")) {
+      r <- run_design(nce_records(rec), d, v, outcomes = c(male = "t_male", disease = "t_disease"))
+      if (v == "V1_ignore") r <- rbind(r, copy(r)[, `:=`(variant = "V2_cluster", se = se_cluster)])
+      out[[length(out) + 1]] <- r[, analysis := "nce"]
+    }
+  }
+  for (d in c("landmark", "ccw")) {     # unadjusted NCE (supplement, D36), V1/V2 only
+    r <- run_design(nce_records(rec, adjust_pvi = FALSE), d, "V1_ignore", outcomes = c(male = "t_male", disease = "t_disease"))
+    out[[length(out) + 1]] <- r[, analysis := "nce_unadjusted"]
+  }
+  out[[length(out) + 1]] <- run_oracle(rec)[, analysis := "main"]   # own estimand of landmark / naive
+  for (dz in g$rct) {
+    r <- run_rct(par, dz, rep = j$rep)
+    out[[length(out) + 1]] <- r[, `:=`(contrast = analysis, analysis = "rct", variant = "none", se = se_patient)]
+  }
+  res <- rbindlist(out, fill = TRUE)[, `:=`(world = j$world, scenario = j$scenario, rep = j$rep)]
+  dir.create(dirname(j$file), recursive = TRUE, showWarnings = FALSE)
+  fwrite(res, j$file)
+  NULL
+}
+t0 <- Sys.time()
+res <- parallel::mclapply(seq_len(nrow(todo)), function(i) tryCatch(one_job(i), error = function(e)
+  sprintf("%s/%s/%d: %s", todo$world[i], todo$scenario[i], todo$rep[i], conditionMessage(e))),
+  mc.cores = cores, mc.preschedule = FALSE)
+err <- unlist(res[!vapply(res, is.null, TRUE)])
+cat(sprintf("done in %.1f h; %d errors\n", as.numeric(Sys.time() - t0, units = "hours"), length(err)))
+dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
+if (length(err)) { writeLines(err, file.path(outdir, sprintf("errors_%s.txt", format(Sys.time(), "%Y%m%d_%H%M")))); cat(head(err, 3), sep = "\n") }
@@ -0,0 +1,34 @@
+# Validation (plan section 15): in a clean world, each design averaged over replicates
+# should match its own truth within Monte Carlo error.
+# Usage: Rscript simulations/validate_designs.R [world] [scenario] [reps] [cores] [n]
+
+pkgload::load_all("pvisim", quiet = TRUE)
+suppressPackageStartupMessages(library(data.table))
+a <- commandArgs(TRUE)
+world <- if (length(a) >= 1) a[1] else "W0"
+scenario <- if (length(a) >= 2) a[2] else "harm"
+reps <- if (length(a) >= 3) as.integer(a[3]) else 40L
+cores <- if (length(a) >= 4) as.integer(a[4]) else 10L
+n <- if (length(a) >= 5) as.numeric(a[5]) else 5e4
+
+par <- load_params(world, scenario)
+res <- rbindlist(parallel::mclapply(seq_len(reps), function(r) {
+  rec <- simulate_records(par, rep = 1000L + r, n = n)
+  rbindlist(lapply(c("naive", "landmark", "ccw", "seq"), function(d)
+    cbind(rep = r, run_design(rec, d, "V1_ignore", outcomes = c(male = "t_male", disease = "t_disease", nco1 = "t_nco1")))),
+    fill = TRUE)
+}, mc.cores = cores), fill = TRUE)
+
+w_ <- world; s_ <- scenario
+tr <- fread("simulations/truths/truths.csv")[world == w_ & scenario == s_]
+map <- data.table(design = c("ccw", "ccw", "seq"), contrast = c("A_B1", "A_B2", "now_not"),
+                  estimand = c("rd_A_B1", "rd_A_B2", "seq_itt"))
+s <- res[, .(mean_est = mean(est), emp_se = sd(est), mean_se = mean(se), reps = .N),
+         by = .(design, contrast, outcome)]
+s <- merge(s, map, by = c("design", "contrast"), all.x = TRUE)
+s <- merge(s, tr[, .(outcome, estimand, truth = value)], by = c("outcome", "estimand"), all.x = TRUE)
+s[, `:=`(bias = mean_est - truth, mcse_bias = emp_se / sqrt(reps))]
+s[, z := bias / mcse_bias]
+print(s[order(design, outcome), .(design, contrast, outcome, truth = round(truth, 4), mean_est = round(mean_est, 4),
+          bias = round(bias, 4), z = round(z, 1), emp_se = round(emp_se, 4), mean_se = round(mean_se, 4))])
+fwrite(res, sprintf("results/validate_%s_%s.csv", world, scenario))
Author	SHA1	Message	Date
chendun1995andClaude Opus 5.5	e1457ccc1e	Decisions 45-50: scenario tuning, NC truths, oracle own estimand, separate PVI phases, validation Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 09:20:22 -04:00
chendun1995andClaude Opus 5.5	383ef89c2a	Validation, NCE analysis, per-replicate oracle for landmark/naive, RCT truths script, resumable experiment runner - W0 harm validation (40 reps): CCW and sequential ITT unbiased for their truths (\|z\| <= 0.7) - PVI models fitted separately inside and after the window - experiments.yml grid (8 worlds x 3 scenarios; designs x V1-V4; NCE; RCTs); smoke run clean Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 09:19:59 -04:00
chendun1995andClaude Opus 5.5	fc1461d5b5	Truths, estimators (naive, landmark, CCW, sequential; physician variants V1-V4) and simulated RCTs - truths.R: interventional truths by common random numbers (A, B1, B2, sequential ITT) - estimators.R: weighted Aalen-Johansen with influence-function SEs (patient / physician cluster) - rct.R: realistic and ideal trials, ITT and as-treated; RCT ITT truths - scenario effects tuned in W0: harm +5.1 pp, benefit-then-harm -1.2 pp at 24 mo, +3.1 at 60 Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 09:01:13 -04:00
chendun1995andClaude Opus 5.5	5ae8874693	Decisions 40-44: build notes for simulator and calibration Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 08:15:30 -04:00
chendun1995andClaude Opus 5.5	5c4407f876	pvisim: parameters, simulator, claims records, tests; calibration of all worlds to Medicare targets - latent-severity DGP with physicians (skill, style), policies A/B1/B2/sequential, common random numbers - sharp null exact at the patient level (tested) - 01_calibrate.R: joint logit-scale update; all worlds match early PVI 3.2%, late PVI 22.5/3.6%, CLTI 7.8%, death 35% - W4 deferral target 15% (trial-like crossover); post-CLTI PVI lowered so late PVI is mostly elective Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 08:15:18 -04:00
chendun1995andClaude Opus 5.5	eb886d5505	Decisions 38-39: keep the three aims Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 07:54:19 -04:00
chendun1995andClaude Opus 5.5	a59363a4b9	Decisions 34-37 (sharp disease null, threat strengths, NCE adjustment, sorting) applied to the analysis plan Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 07:01:05 -04:00
chendun1995andClaude Opus 5.5	d2d1e3a1d4	Analysis plan v1.0 (md + docx) with draft Table 1 of claudication RCTs Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 06:58:58 -04:00
chendun1995andClaude Opus 5.5	645713f7ec	Decisions 32-33: negative controls Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 06:52:22 -04:00
chendun1995andClaude Opus 5.5	c3c22fdcdf	Literature reviews and decisions log (rounds 1-5) for the early-PVI TTE simulation Co-Authored-By: Claude Opus 5.5 <[email protected]>	2026-10-27 06:51:17 -04:00