diff --git a/docs/decisions.md b/docs/decisions.md index 4ca90f5..5d9b32f 100644 --- a/docs/decisions.md +++ b/docs/decisions.md @@ -117,3 +117,4 @@ Each decision is listed with the date and the reason for it. Settled decisions a |---|---|---| | 54 | **The paper's question: "Can target trial emulation overcome the problems that make randomized trials hard in surgery?"** Early PVI for claudication is the worked example. Medicare only supplies realistic simulation inputs, and the paper makes no claims about Medicare findings or the authors' earlier papers. The structure is a scorecard: each surgical-RCT barrier (small/stopped early, short follow-up/soft endpoints, crossover, selected surgeons, practice style, unmeasured severity) against whether TTE solves it, checked against the known truth. Aims 1–3 stay. Aim 1 = "TTE is not one method", Aim 2 = the scorecard rows, Aim 3 = whether the one unfixable problem can be detected. The MALE/repeat-revascularization point is secondary, general outcome-definition guidance. The Medicare "dial" idea is dropped. **Novelty (quick PubMed check, 2026-10-28):** there are TTE guides for surgeons (e.g. JBJS Rev 2026) and immortal-time design comparisons (Duchesneau 2022 *CEBP*), but no search hit for a known-truth simulation comparing TTE with a feasible surgical RCT under surgeon-level threats. | User: the aim is to test whether TTE can solve RCT challenges in surgery, without arguing against the authors' earlier results. | | 55 | **Three additions from literature review 3 (user, 2026-10-28), for the overnight run.** (1) **New world WE, "EHR/registry-like":** the same DGP and calibration as W0, but the analyst sees severity with measurement error (monthly Z_t + N(0, 0.75²), correlation with true severity about 0.8), a frailty index with reliability 0.7, and true smoking. It sits between W0 (exact) and W5 (claims proxy) and answers the claims-vs-EHR row of Table 1. (2) **Non-random consent in the realistic trial:** mean 30%, OR 0.7 per 10 years of age and per SD of frailty (Wasmann 2019, King 2005). The trial's ITT truth is weighted by consent probability. (3) **Trial horizons 12, 24 and 60 months:** ITT and as-treated at each. Truths `rct__itt_h12/_h24`. The reference at the short horizons is RD(A − B1) at that month, from the risk curves. **Learning curve:** not simulated; limitation (user). | Completes Table 1 rows 3, 5 and 17. | +| 56 | **Primary outcome changed to limb MALE: the first of coded CLTI, major amputation or open bypass** (user, 2026-10-28). D7's MALE (which also counts any revascularization after the first) moves to the supplement as "MALE including repeat revascularization". Repeat revascularization is reported on its own as a physician-driven outcome. Limb MALE also replaces the disease composite (CLTI or amputation) as the headline (D34); the disease composite stays as a secondary outcome. Negative-control bias (Aim 3) is now measured on limb MALE. Clean-world truths: 0 under the null (exactly); +6.2 under harm; −0.7 at 12 months, −1.7 at 24 and +3.6 at 60 under benefit-then-harm. Acute limb ischemia is not simulated (rare in claudication; limitation). | Standard MALE definitions (Conte 2009 *JVS* 19897335; Nordanstig 2024 *EJVES* 39121905) are built for follow-up **after** a revascularization: reintervention exists only for treated patients. D7 argued that excluding the index PVI in both arms made the outcome symmetric. It does not: under A everyone is at risk of a repeat for about 55 months, under B1 about 3.6% are. In the clean world under the null, repeat revascularization is 18.6% vs 0.7%, giving +14 pp of MALE from the definition alone. | diff --git a/pvisim/R/estimators.R b/pvisim/R/estimators.R index 8b5601d..63e442a 100644 --- a/pvisim/R/estimators.R +++ b/pvisim/R/estimators.R @@ -174,14 +174,14 @@ restrict_variant <- function(rec, variant) { #' 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", +run_design <- function(rec, design, variant, outcomes = c(limb = "t_limb", male = "t_male", disease = "t_disease", nco1 = "t_nco1", nco2 = "t_nco2")) { run_designs(rec, design, variant, outcomes) } #' Several designs under one variant, sharing the PVI fit between CCW and sequential #' trials (the sequential propensity score is the window part of the same model). -run_designs <- function(rec, designs, variant, outcomes = c(male = "t_male", disease = "t_disease", +run_designs <- function(rec, designs, variant, outcomes = c(limb = "t_limb", male = "t_male", disease = "t_disease", nco1 = "t_nco1", nco2 = "t_nco2")) { rec <- restrict_variant(rec, variant) need <- intersect(designs, c("ccw", "seq")) diff --git a/pvisim/R/limits.R b/pvisim/R/limits.R index 4e0dc8c..cffbfa9 100644 --- a/pvisim/R/limits.R +++ b/pvisim/R/limits.R @@ -18,7 +18,7 @@ oracle_records <- function(rec) { } #' 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")) { +run_oracle <- function(rec, outcomes = c(limb = "t_limb", 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) diff --git a/pvisim/R/rct.R b/pvisim/R/rct.R index 048d83b..869247f 100644 --- a/pvisim/R/rct.R +++ b/pvisim/R/rct.R @@ -60,7 +60,7 @@ trial_policy <- function(trial, assign = NULL) { #' 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"), +run_rct <- function(par, design = "realistic", rep = 1L, outcomes = c(limb = "t_limb", 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) diff --git a/pvisim/R/simulate.R b/pvisim/R/simulate.R index 2bb059c..7809c98 100644 --- a/pvisim/R/simulate.R +++ b/pvisim/R/simulate.R @@ -217,8 +217,11 @@ simulate_course <- function(pop, par, pol = policy("natural"), seed = 1L, record # disease composite (D34): coded CLTI or major amputation t_disease <- pmin(t_clti_code, t_amp, na.rm = TRUE) + # limb MALE (primary, D56): coded CLTI, major amputation or open bypass. Every component + # can occur with or without a prior PVI; repeat revascularization is reported separately + t_limb <- pmin(t_clti_code, t_amp, t_bypass, 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) + t_rr, t_male, t_disease, t_limb, t_death, t_nco1, t_nco2, t_nce, n_proc) list(events = out, z = zmat, counts = cmat, policy = pol) } diff --git a/pvisim/R/truths.R b/pvisim/R/truths.R index b51eca1..341b002 100644 --- a/pvisim/R/truths.R +++ b/pvisim/R/truths.R @@ -9,7 +9,7 @@ # 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", +OUTCOMES <- c(limb = "t_limb", 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) { @@ -83,7 +83,7 @@ compute_truths <- function(world, scenario, n = 1e6, seed = 20261027L, horizon = value = rd, mcse = sqrt(var / n_elig), n = n_elig)] } - curves <- rbindlist(lapply(names(ev), function(p) rbindlist(lapply(c("male", "disease"), function(o) + curves <- rbindlist(lapply(names(ev), function(p) rbindlist(lapply(c("limb", "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, diff --git a/pvisim/tests/testthat/test-simulate.R b/pvisim/tests/testthat/test-simulate.R index 2732ca6..87e1fa3 100644 --- a/pvisim/tests/testthat/test-simulate.R +++ b/pvisim/tests/testthat/test-simulate.R @@ -20,7 +20,7 @@ test_that("under the sharp null, disease, death and negative controls are identi 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")) { + for (col in c("t_limb", "t_disease", "t_clti_code", "t_amp", "t_bypass", "t_death", "t_nco1", "t_nco2", "t_nce")) { expect_identical(r[[col]], runs[[1]][[col]], info = col) } } diff --git a/reports/draft-report.Rmd b/reports/draft-report.Rmd index cc52b70..40796ab 100644 --- a/reports/draft-report.Rmd +++ b/reports/draft-report.Rmd @@ -52,7 +52,7 @@ perf[, vlab := factor(VAR[variant], levels = VAR)] variants, 3 negative controls, an oracle run, and 2 simulated RCTs. 2,400 replicate files, 0 errors. 16 of 288,000 estimates are missing (sequential trials with physician fixed effects, 4 replicates), and these are dropped. -* **Outcome shown.** The headline outcome is the **disease composite (CLTI or major amputation)**, +* **Outcome shown.** The headline outcome is **limb MALE (CLTI, major amputation or open bypass; D56)**, whose true effect under the null is exactly 0 (D34). MALE is shown where it adds something. * **Two errors, kept separate (D27).** *Estimation bias* = mean estimate − the method's own truth. *Question gap* = the method's own truth − the reference truth (CCW A vs B1). *Total error* is the sum. @@ -64,7 +64,7 @@ perf[, vlab := factor(VAR[variant], levels = VAR)] # 1. Aim 1: what each method answers (clean world W0) ```{r tab2} -t2 <- perf[world == "W0" & outcome == "disease" & variant %in% c("V2_cluster", "none") & +t2 <- perf[world == "W0" & outcome == "limb" & variant %in% c("V2_cluster", "none") & !grepl("as_treated", key)] t2 <- t2[, .(Scenario = scenario, Method = as.character(method), `Reference truth` = pp(ref), `Own truth` = pp(own), `Question gap` = pp(gap), `Mean estimate` = pp(mean_est), @@ -89,7 +89,7 @@ knitr::kable(t2[order(Scenario, Method)], align = "llrrrrrrr", the 5-point harm in only 25% of replicates. ```{r fig1, fig.width = 11, fig.height = 9} -ladder <- function(sc, oc = "disease") { +ladder <- function(sc, oc = "limb") { x <- perf[scenario == sc & outcome == oc & variant != "V1_ignore" & !grepl("as_treated", key)] r <- reps[scenario == sc & outcome == oc & analysis %in% c("main", "rct") & variant != "V1_ignore" & variant != "oracle" & !grepl("as_treated", contrast)] @@ -109,24 +109,24 @@ ladder <- function(sc, oc = "disease") { "Hollow dot: the method's own truth. Solid line: the reference truth (CCW A vs B1).")) + coord_cartesian(xlim = c(-15, 35)) } -ladder("harm") + ggtitle("Figure 1a. Estimand ladder, harm scenario (disease composite)") +ladder("harm") + ggtitle("Figure 1a. Estimand ladder, harm scenario (limb MALE)") ``` ```{r fig1b, fig.width = 11, fig.height = 9} -ladder("null") + ggtitle("Figure 1b. Estimand ladder, null scenario (disease composite)") +ladder("null") + ggtitle("Figure 1b. Estimand ladder, null scenario (limb MALE)") ``` # 2. Aim 2: surgery-specific threats ```{r tab3} -x <- perf[scenario == "harm" & outcome == "disease" & analysis == "main" & +x <- perf[scenario == "harm" & outcome == "limb" & analysis == "main" & variant %in% c("V2_cluster", "V3_random", "V4_within") & key != "ccw|A_B2"] x[, cell := sprintf("%s (%s)", pp(bias), pp(mcse_bias))] t3 <- dcast(x, method + vlab ~ wlab, value.var = "cell") setnames(t3, c("method", "vlab"), c("Method", "Physician variant")) knitr::kable(t3[order(-as.integer(Method))], caption = paste( "Table 3. Estimation bias (MCSE) against each method's own truth, percentage points.", - "Harm scenario, disease composite.")) + "Harm scenario, limb MALE.")) ``` ```{r tab3b} @@ -135,7 +135,7 @@ t3b <- dcast(x, method + vlab ~ wlab, value.var = "cell2") setnames(t3b, c("method", "vlab"), c("Method", "Physician variant")) knitr::kable(t3b[order(-as.integer(Method))], caption = paste( "Table 3b. Total error against the reference truth (CCW A vs B1), percentage points.", - "Harm scenario, disease composite.")) + "Harm scenario, limb MALE.")) ``` **What it shows.** @@ -157,28 +157,28 @@ knitr::kable(t3b[order(-as.integer(Method))], caption = paste( # 3. False alarms under the null (headline for clinicians) ```{r fig2, fig.width = 10, fig.height = 5.2} -fa <- perf[scenario == "null" & outcome == "disease" & analysis == "main" & key != "ccw|A_B2" & variant != "V1_ignore"] +fa <- perf[scenario == "null" & outcome == "limb" & analysis == "main" & key != "ccw|A_B2" & variant != "V1_ignore"] fa[, row := factor(paste0(method, " · ", sub(" \\(.*", "", vlab)))] fa[, row := factor(row, levels = rev(sort(levels(row))))] p1 <- ggplot(fa, aes(wlab, row, fill = 100 * reject0)) + geom_tile(colour = "white", linewidth = 1) + geom_text(aes(label = sprintf("%.0f", 100 * reject0), colour = reject0 > 0.5), size = 3) + scale_colour_manual(values = c(`TRUE` = "white", `FALSE` = "#0b0b0b"), guide = "none") + scale_fill_gradient(low = "#eef4fc", high = "#0d3d78", limits = c(0, 100), name = "% of replicates\nCI excludes 0") + - labs(x = NULL, y = NULL, title = "Figure 2. False-alarm rate, null scenario (disease composite; nominal 5%)") + + labs(x = NULL, y = NULL, title = "Figure 2. False-alarm rate, null scenario (limb MALE; nominal 5%)") + theme(axis.text.x = element_text(angle = 30, hjust = 1), panel.grid = element_blank(), legend.position = "right") p1 ``` ```{r fig2b, fig.width = 8, fig.height = 3} m0 <- perf[scenario == "null" & world == "W0" & variant %in% c("V2_cluster", "none") & - outcome %in% c("male", "disease") & !grepl("as_treated", key) & key != "ccw|A_B2"] -m0[, oc := factor(ifelse(outcome == "male", "MALE", "Disease composite"), levels = c("Disease composite", "MALE"))] + outcome %in% c("male", "limb") & !grepl("as_treated", key) & key != "ccw|A_B2"] +m0[, oc := factor(ifelse(outcome == "male", "MALE incl. repeat revascularization", "Limb MALE"), levels = c("Limb MALE", "MALE incl. repeat revascularization"))] ggplot(m0, aes(100 * mean_est, method, colour = oc, group = oc)) + geom_vline(xintercept = 0, colour = "#c3c2b7") + geom_errorbarh(aes(xmin = 100 * (mean_est - 1.96 * emp_se), xmax = 100 * (mean_est + 1.96 * emp_se)), height = 0, linewidth = 0.6, position = position_dodge(width = 0.5)) + - geom_point(size = 2, position = position_dodge(width = 0.5)) + scale_colour_manual(values = c("Disease composite" = "#2a78d6", MALE = "#eb6834"), name = NULL) + + geom_point(size = 2, position = position_dodge(width = 0.5)) + scale_colour_manual(values = c("Limb MALE" = "#2a78d6", "MALE incl. repeat revascularization" = "#eb6834"), name = NULL) + labs(x = "Risk difference (pp), W0 null", y = NULL, - title = "Figure 2 inset. Sharp null: disease 0, yet MALE about +14 pp") + title = "Figure 2 inset. Sharp null: limb MALE 0, yet MALE incl. repeat revascularization about +14 pp") ``` **What it shows.** @@ -242,10 +242,10 @@ ggplot(cb, aes(wlab, 100 * absbias, colour = flag)) + # 5. Simulated RCTs ```{r tabrct} -r <- perf[analysis == "rct" & outcome == "disease" & scenario != "null"] +r <- perf[analysis == "rct" & outcome == "limb" & scenario != "null"] r <- r[, .(World = wlab, Scenario = scenario, Trial = method, `Reference` = pp(ref), `Own ITT truth` = pp(own), `Mean estimate` = pp(mean_est), `Empirical SE` = pp(emp_se), `Power (CI excludes 0)` = sprintf("%.0f%%", 100 * reject0))] -knitr::kable(r[order(World, Scenario, Trial)], caption = "Table 5. Simulated trials, disease composite (pp).") +knitr::kable(r[order(World, Scenario, Trial)], caption = "Table 5. Simulated trials, limb MALE (pp).") ``` * **Surgeon selection changes the answer.** In WR, a realistic trial recruits hospital-based, @@ -291,7 +291,7 @@ if (length(bf)) { x <- bt[, .(reps = .N, boot = sqrt(mean(boot_se^2)), patient = sqrt(mean(se_patient^2)), cluster = sqrt(mean(se_cluster^2))), by = .(world, scenario, design, contrast, outcome)] x <- merge(x, emp, by = c("world", "scenario", "design", "contrast", "outcome")) - x <- x[outcome == "disease" & contrast != "A_B2"] + x <- x[outcome == "limb" & contrast != "A_B2"] knitr::kable(x[order(world, scenario, design), .(World = world, Scenario = scenario, Design = design, `Empirical SE (100 reps)` = pp(emp_se, 2), `Bootstrap SE` = pp(boot, 2), `Model SE, cluster` = pp(cluster, 2), `Model SE, patient` = pp(patient, 2), `Bootstrap / empirical` = round(boot / emp_se, 2), @@ -303,7 +303,7 @@ if (length(bf)) { # Supplement: full performance table ```{r supp} -x <- perf[analysis == "main" & outcome %in% c("disease", "male")] +x <- perf[analysis == "main" & outcome %in% c("limb", "male", "disease")] knitr::kable(x[order(outcome, scenario, world, method, variant), .(outcome, scenario, world, method = as.character(method), variant, est = pp(mean_est, 2), own = pp(own, 2), ref = pp(ref, 2), bias = pp(bias, 2), mcse = pp(mcse_bias, 2), gap = pp(gap, 2), emp_se = pp(emp_se, 2), diff --git a/simulations/03_limits.R b/simulations/03_limits.R index 0337bbd..478b607 100644 --- a/simulations/03_limits.R +++ b/simulations/03_limits.R @@ -11,8 +11,8 @@ jobs <- CJ(world = names(raw$worlds), scenario = names(raw$scenarios), sorted = 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")]) + rbind(rct_truth(par, "realistic", n = 1e6)[outcome %in% c("limb", "male", "disease", "rr", "nco1", "nco2")], + rct_truth(par, "ideal", n = 1e6)[outcome %in% c("limb", "male", "disease", "rr", "nco1", "nco2")]) }, mc.cores = cores) bad <- vapply(res, inherits, TRUE, "try-error") if (any(bad)) stop(res[bad][[1]]) diff --git a/simulations/04_summarise.R b/simulations/04_summarise.R index 35ec84f..82bc5df 100644 --- a/simulations/04_summarise.R +++ b/simulations/04_summarise.R @@ -67,16 +67,16 @@ perf[, `:=`(mcse_cover = sqrt(cover_own * (1 - cover_own) / reps), 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 +# One unit = world x scenario x design x variant x replicate. Bias is in the limb-MALE RD (D56) # against the method's own estimand; a control "flags" when its 95% CI excludes 0. 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 -nc <- m[outcome == "male", c(..u, "est", "own")][, err := est - own] +nc <- m[outcome == "limb", c(..u, "est", "own")][, err := est - own] nc <- merge(nc, m[outcome == "nco1", c(..u)][, nco1 := flag(m[outcome == "nco1"])], by = u, all.x = TRUE) nc <- merge(nc, m[outcome == "nco2", c(..u)][, nco2 := flag(m[outcome == "nco2"])], by = u, all.x = TRUE) -ne <- d[analysis == "nce" & outcome == "male" & contrast %in% c("early_none", "A_B1")] +ne <- d[analysis == "nce" & outcome == "limb" & contrast %in% c("early_none", "A_B1")] nc <- merge(nc, ne[, c(..u)][, nce := flag(ne)], by = u, all.x = TRUE) nc[, any_nc := nco1 | nco2 | fifelse(is.na(nce), FALSE, nce)] nc[, cell_bias := mean(err), by = .(world, scenario, design, variant)] diff --git a/simulations/05_bootstrap.R b/simulations/05_bootstrap.R index 115d1f8..d9d225d 100644 --- a/simulations/05_bootstrap.R +++ b/simulations/05_bootstrap.R @@ -13,7 +13,7 @@ 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") +outs <- c(limb = "t_limb", 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))] diff --git a/simulations/run_experiment.R b/simulations/run_experiment.R index abb4c8d..aab4dd9 100644 --- a/simulations/run_experiment.R +++ b/simulations/run_experiment.R @@ -23,7 +23,7 @@ one_job <- function(i) { out <- list() with_v2 <- function(r) if (r$variant[1] == "V1_ignore") # same fit, two SEs rbind(r, copy(r)[, `:=`(variant = "V2_cluster", se = se_cluster)]) else r - nce_out <- c(male = "t_male", disease = "t_disease") + nce_out <- c(limb = "t_limb", male = "t_male", disease = "t_disease") for (v in g$variants) { out[[length(out) + 1]] <- with_v2(run_designs(rec, unlist(g$designs), v))[, analysis := "main"] if (isTRUE(g$nce))