D53: landmark on baseline covariates (month-5 version as landmark_m5), NCO-1 depends on PVI-free severity, physician-cluster bootstrap script

Co-Authored-By: Claude Opus 5.5 <[email protected]>

docs/decisions.md

@@ -109,3 +109,4 @@ Each decision is listed with the date and the reason for it. Settled decisions a
 | # | Decision | Why |
 |---|---|---|
 | 52 | **Draft run done:** 8 worlds × 3 scenarios × 100 replicates of 50,000 patients; 9.4 h on 11 cores; 0 errors. 16 of 288,000 estimates are NA (sequential trials, V4, 4 replicates) and are dropped. Summary: `simulations/04_summarise.R` writes `results/summary/`. Report: `reports/draft-report.Rmd`. Two open questions: (Q1) the landmark propensity model adjusts for month-5 severity, a post-PVI mediator, which drives the landmark question gap (W1 own truth about 20 pp vs a reference of 7.2 pp); with baseline covariates only, the 20-replicate oracle mean is 8.1. (Q2) NCO-1 has no severity term in the DGP, contrary to D32, so it cannot detect W5/WR bias. The bootstrap check (D28) waits on Q1/Q2. | Surfaced by the draft results. |
+| 53 | **Answers to D52 (user, 2026-10-28).** (Q1) **The primary landmark propensity score uses baseline covariates only** (diagnosis covariates and month-0 severity). The version that also adjusts for month-5 severity becomes `landmark_m5` (supplement). It conditions on a mediator, because a PVI in months 0–5 has already changed that severity. Both versions get their own oracle. (Q2) **NCO-1 now depends on severity (coefficient 0.3)**, using the patient's **PVI-free severity path** (the same baseline, rate and monthly noise, with no relief and no kappa). That way NCO-1 shares the confounder but still cannot be affected by PVI (only through death as a competing risk, D46). Intercept −6.75 keeps 5-year NCO-1 incidence at about 9.3%. This is the best case for a patient-level negative control, and the paper says so. **Bootstrap (D28, reduced):** W0 and WR × null and harm, replicates 1–20 of the main run, 50 physician-cluster resamples each, V1 point estimates, about 8 s per resample. Compare the bootstrap SD with the model SE (patient and cluster) and the empirical SE. **The rerun goes to `results/main/n50k_v2`;** the first draft stays in `n50k`. Truths are recomputed (NCO-1 changed). | Landmark analyses in practice adjust for baseline covariates. D32 intended NCO-1 to share severity. |

pvisim/R/estimators.R

@@ -197,6 +197,9 @@ run_designs <- function(rec, designs, variant, outcomes = c(male = "t_male", dis
         cbind(contrast = "early_none", est_naive(rec, outcomes[[o]], variant), outcome = o))),
       landmark = rbindlist(lapply(names(outcomes), function(o)
         cbind(contrast = "early_none", est_landmark(rec, outcomes[[o]], variant, extra_covs = rec$extra_covs), outcome = o))),
+      landmark_m5 = rbindlist(lapply(names(outcomes), function(o)
+        cbind(contrast = "early_none", est_landmark(rec, outcomes[[o]], variant, extra_covs = rec$extra_covs,
+                                                    at_landmark = TRUE), outcome = o))),
       ccw = est_ccw(rec, outcomes, variant, rs = copy(rs), phat = phat, comparators = cmp),
       seq = cbind(contrast = "now_not", est_seq(rec, outcomes, variant, rs = copy(rs), phat = phat)))
     r[, `:=`(design = design, variant = variant, n = nrow(rec$patients))]
@@ -239,9 +242,11 @@ fit_pvi_ps <- function(p, variant, covs = c(BASE_COVS, "sev_obs", "sev0")) {
 }
 
 #' 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;
+#' exposure = PVI in months 0..W-1; covariates at diagnosis; outcome counted in months
-#' outcome counted in months W..horizon-1.
+#' W..horizon-1. With at_landmark = TRUE the propensity score also uses severity at month
-est_landmark <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60L, extra_covs = NULL) {
+#' W-1, which a PVI has already changed (supplement, D53).
+est_landmark <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60L, extra_covs = NULL,
+                         at_landmark = FALSE) {
   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"))]
@@ -249,7 +254,7 @@ est_landmark <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon
   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))
+  ps <- fit_pvi_ps(p, variant, covs = c(BASE_COVS, "sev0", if (at_landmark) "sev_obs", 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")

pvisim/R/limits.R

@@ -17,13 +17,14 @@ oracle_records <- function(rec) {
   o
 }
 
-#' Oracle versions of the landmark and naive designs (V1, patient-level SE).
+#' Oracle versions of the landmark (both versions) 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)
+    l5 <- est_landmark(o, outcomes[[nm]], "V1_ignore", extra_covs = o$extra_covs, at_landmark = TRUE)
     nv <- est_naive(o, outcomes[[nm]], "V1_ignore", extra_covs = o$extra_covs)
-    rbind(cbind(design = "landmark", lm), cbind(design = "naive", nv))[
+    rbind(cbind(design = "landmark", lm), cbind(design = "landmark_m5", l5), cbind(design = "naive", nv))[
       , `:=`(outcome = nm, contrast = "early_none", variant = "oracle", se = se_patient)]
   }))
 }

pvisim/R/records.R

@@ -68,3 +68,24 @@ replicate_seed <- function(world, scenario, rep) {
   h <- sum(utf8ToInt(key) * seq_along(utf8ToInt(key))^2)
   as.integer((h * 2654435761) %% 2147483647)
 }
+
+#' Physician-cluster bootstrap resample of a replicate's records (D28): physicians are
+#' drawn with replacement and each draw brings all of its patients, with new patient and
+#' physician ids so repeated physicians count as separate clusters.
+bootstrap_records <- function(rec, seed) {
+  set.seed(seed)
+  phys <- unique(rec$patients$phys_id)
+  draw <- data.table(phys_id = sample(phys, length(phys), replace = TRUE), new_phys = seq_along(phys))
+  map <- rec$patients[, .(id, phys_id)][draw, on = "phys_id", allow.cartesian = TRUE]
+  map[, new_id := seq_len(.N)]
+  remap <- function(x) {
+    x <- x[map[, .(id, new_id)], on = "id", allow.cartesian = TRUE]
+    x[, id := new_id][, new_id := NULL][]
+  }
+  out <- rec
+  out$patients <- remap(rec$patients)[map[, .(new_id, new_phys)], on = c(id = "new_id"), phys_id := i.new_phys]
+  out$pm <- setkey(remap(rec$pm), id, t)
+  out$hidden <- remap(rec$hidden)
+  out$hidden_pm <- remap(rec$hidden_pm)
+  out
+}

pvisim/R/simulate.R

@@ -127,12 +127,14 @@ simulate_course <- function(pop, par, pol = policy("natural"), seed = 1L, record
   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
+  nco1_z <- if (is.null(nc$nco1$z)) 0 else nc$nco1$z
   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
+  Zn <- pat$z0     # severity on the natural course without any PVI (relief or kappa): drives NCO-1 (D53)
   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 <-
@@ -147,8 +149,11 @@ simulate_course <- function(pop, par, pol = policy("natural"), seed = 1L, record
   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] +
+    if (t > 0) {
-      par$severity$eps_sd * qnorm(U[a, "eps"])
+      eps <- par$severity$eps_sd * qnorm(U[a, "eps"])
+      Z[a] <- Z[a] + pat$rate[a] + kappa_post[a] * treated[a] + eps
+      Zn[a] <- Zn[a] + pat$rate[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]))
@@ -183,7 +188,7 @@ simulate_course <- function(pop, par, pol = policy("natural"), seed = 1L, record
     }
 
     # 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_nco1) & U[, "nco1"] < plogis(lp_nco1 + nco1_z * Zn); 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 }
 

pvisim/inst/params/dgp_default.yml

@@ -78,7 +78,7 @@ outcomes:
   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
+  nco1: {intercept: -6.75, frailty: 0.5, smoking: 0.4, age10: 0.3, z: 0.3}   # hip fracture / pneumonia; z = PVI-free severity (D53)
   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
 

reports/draft-report.Rmd

@@ -25,9 +25,9 @@ VAR <- c(V2_cluster = "Ignore physician (cluster SE)", V3_random = "Random inter
 VCOL <- c("Ignore physician (cluster SE)" = "#2a78d6", "Random intercept" = "#eb6834",
           "Within-physician (FE)" = "#1baf7a", "Trial (no adjustment)" = "#52514e")
 METH <- data.table(
-  k = c("ccw|A_B1", "ccw|A_B2", "seq|now_not", "landmark|early_none", "naive|early_none",
+  k = c("ccw|A_B1", "ccw|A_B2", "seq|now_not", "landmark|early_none", "landmark_m5|early_none", "naive|early_none",
           "rct_realistic|itt", "rct_ideal|itt", "rct_realistic|as_treated", "rct_ideal|as_treated"),
-  method = c("CCW, A vs B1 (reference)", "CCW, A vs B2", "Sequential trials, ITT", "Landmark at 6 months",
+  method = c("CCW, A vs B1 (reference)", "CCW, A vs B2", "Sequential trials, ITT", "Landmark at 6 months", "Landmark, + month-5 severity (suppl.)",
              "Naive 'ever early PVI'", "Realistic RCT (n = 240), ITT", "Ideal pragmatic RCT, ITT",
              "Realistic RCT, as-treated", "Ideal RCT, as-treated"))
 setnames(METH, "k", "key")

simulations/04_summarise.R

@@ -4,7 +4,7 @@
 
 suppressPackageStartupMessages(library(data.table))
 a <- commandArgs(TRUE)
-indir <- if (length(a) >= 1) a[1] else "results/main/n50k"
+indir <- if (length(a) >= 1) a[1] else "results/main/n50k_v2"
 outdir <- "results/summary"
 dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
 
@@ -32,7 +32,7 @@ d <- merge(d, tr[estimand == "rd_A_B1", .(world, scenario, outcome, ref = truth)
 # landmark and naive: own estimand = the paired oracle estimate in the same replicate (D47)
 orc <- d[variant == "oracle", .(world, scenario, rep, design, outcome, oracle = est)]
 d <- merge(d, orc, by = c("world", "scenario", "rep", "design", "outcome"), all.x = TRUE)
-d[design %in% c("landmark", "naive") & variant != "oracle", own := oracle]
+d[design %in% c("landmark", "landmark_m5", "naive") & variant != "oracle", own := oracle]
 d[analysis %like% "^nce", `:=`(own = 0, ref = 0)]
 d[, oracle := NULL]
 d[, `:=`(lo = est - 1.96 * se, hi = est + 1.96 * se)]
@@ -59,7 +59,7 @@ fwrite(perf, file.path(outdir, "perf.csv"))
 # ---- negative controls as bias detectors (plan 10) ------------------------------------
 # One unit = world x scenario x design x variant x replicate. Bias is in the MALE RD
 # against the method's own estimand; a control "flags" when its 95% CI excludes 0.
-m <- d[analysis == "main" & variant %in% c("V1_ignore", "V2_cluster", "V3_random", "V4_within") &
+m <- d[analysis == "main" & design != "landmark_m5" & variant %in% c("V1_ignore", "V2_cluster", "V3_random", "V4_within") &
          contrast %in% c("early_none", "A_B1", "now_not")]
 u <- c("world", "scenario", "design", "variant", "rep")
 flag <- function(x) x$lo > 0 | x$hi < 0

simulations/05_bootstrap.R

@@ -0,0 +1,35 @@
+# Physician-cluster bootstrap check of the model-based SEs (D28, reduced in D53).
+# Usage: Rscript simulations/05_bootstrap.R [reps] [B] [cores]
+# Writes results/bootstrap/<world>_<scenario>_rep<k>.csv: the estimate and the
+# bootstrap SD for each design, contrast and outcome (V1 point estimates).
+
+pkgload::load_all("pvisim", quiet = TRUE)
+suppressPackageStartupMessages(library(data.table))
+a <- commandArgs(TRUE)
+reps <- if (length(a) >= 1) as.integer(a[1]) else 20L
+B <- if (length(a) >= 2) as.integer(a[2]) else 50L
+cores <- if (length(a) >= 3) as.integer(a[3]) else 11L
+cfg <- yaml::read_yaml("simulations/config/experiments.yml")
+n <- cfg$scales$draft$n_patients
+outdir <- "results/bootstrap"; dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
+designs <- c("naive", "landmark", "landmark_m5", "ccw", "seq")
+outs <- c(male = "t_male", disease = "t_disease")
+
+jobs <- CJ(world = c("W0", "WR"), scenario = c("null", "harm"), rep = seq_len(reps), sorted = FALSE)
+jobs[, file := file.path(outdir, sprintf("%s_%s_rep%03d.csv", world, scenario, rep))]
+todo <- jobs[!file.exists(file)]
+cat(sprintf("bootstrap: %d jobs, %d to do, B = %d, %d cores\n", nrow(jobs), nrow(todo), B, cores))
+res <- parallel::mclapply(seq_len(nrow(todo)), function(i) tryCatch({
+  j <- todo[i]
+  par <- load_params(j$world, j$scenario)
+  rec <- simulate_records(par, rep = j$rep, n = n)          # the same data as the main run's replicate
+  bs <- rbindlist(lapply(seq_len(B), function(b)
+    cbind(b = b, run_designs(bootstrap_records(rec, seed = 1e6 * j$rep + b), designs, "V1_ignore", outs))), fill = TRUE)
+  out <- bs[, .(boot_se = sd(est, na.rm = TRUE), boot_mean = mean(est, na.rm = TRUE), B = sum(!is.na(est))),
+            by = .(design, contrast, outcome)]
+  fwrite(out[, `:=`(world = j$world, scenario = j$scenario, rep = j$rep)], j$file)
+  NULL
+}, 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; %d errors\n", length(err))); if (length(err)) cat(head(err, 3), sep = "\n")

simulations/config/experiments.yml

@@ -2,11 +2,11 @@
 main:
   worlds: [W0, W1, W1_rho0, W2, W3, W4, W5, WR]
   scenarios: ["null", harm, benefit_harm]
-  designs: [naive, landmark, ccw, seq]
+  designs: [naive, landmark, landmark_m5, ccw, seq]   # landmark_m5: also adjusts for month-5 severity (supplement, D53)
   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, dir: n50k}
+  draft: {n_patients: 50000, reps: 100, dir: n50k_v2}   # v2: D53 (landmark baseline covariates, NCO-1 severity)
-  final: {n_patients: 50000, reps: 200, dir: n50k}   # extends the draft's replicates 1-100 to 200 (D51)
+  final: {n_patients: 50000, reps: 200, dir: n50k_v2}   # extends the draft's replicates 1-100 to 200 (D51)