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
parent 5ae8874693
commit fc1461d5b5
8 changed files with 566 additions and 7 deletions
@@ -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)
+}