diff --git a/docs/decisions.md b/docs/decisions.md index 40c8a1f..4ca90f5 100644 --- a/docs/decisions.md +++ b/docs/decisions.md @@ -116,3 +116,4 @@ Each decision is listed with the date and the reason for it. Settled decisions a | # | Decision | Why | |---|---|---| | 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. | diff --git a/pvisim/R/rct.R b/pvisim/R/rct.R index d1c3cae..048d83b 100644 --- a/pvisim/R/rct.R +++ b/pvisim/R/rct.R @@ -1,16 +1,27 @@ # 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. +# bottom half of style; on average 30% of their eligible patients consent, less often when +# older or frailer (OR 0.7 per 10 years and per SD of frailty; D55); 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) + realistic = list(n = 240L, consent = 0.30, select = TRUE, consent_or = c(age10 = 0.7, frailty = 0.7)), + ideal = list(n = NA, consent = 1.00, select = FALSE, consent_or = NULL) ) +RCT_HORIZONS <- c(12L, 24L, 60L) # trials usually report at 1-2 years (D55) + +#' Probability that each eligible patient consents: mean `spec$consent`, lower for older +#' and frailer patients (Wasmann 2019, King 2005: refusers differ moderately by covariates). +consent_prob <- function(pat, spec) { + if (is.null(spec$consent_or)) return(rep(spec$consent, nrow(pat))) + lp <- log(spec$consent_or[["age10"]]) * (pat$age - 76) / 10 + log(spec$consent_or[["frailty"]]) * pat$frailty + a <- uniroot(function(a) mean(plogis(a + lp)) - spec$consent, c(-10, 10))$root + plogis(a + lp) +} #' Physicians a realistic trial would recruit. trial_physicians <- function(phys) { @@ -25,7 +36,7 @@ enrol_trial <- function(pop, design = c("realistic", "ideal"), n = NULL, seed = set.seed(seed) pat <- pop$pat if (spec$select) pat <- pat[phys_id %in% trial_physicians(pop$phys)] - pat <- pat[runif(.N) < spec$consent] + pat <- pat[runif(.N) < consent_prob(pat, spec)] 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))] @@ -58,33 +69,38 @@ run_rct <- function(par, design = "realistic", rep = 1L, outcomes = c(male = "t_ 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) { + rbindlist(lapply(names(outcomes), function(o) rbindlist(lapply(RCT_HORIZONS, function(h) { 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") + itt <- rd_from_aj(aj_cif(expand_followup(copy(p)[, `:=`(arm = assigned, w = 1)], 0L, h), h), "A", "B") + at <- rd_from_aj(aj_cif(expand_followup(copy(p)[, `:=`(arm = treated_win, w = 1)], 0L, h), h), "early", "none") rbind(cbind(analysis = "itt", itt), cbind(analysis = "as_treated", at))[ - , `:=`(outcome = o, design = paste0("rct_", design), n = nrow(ev), + , `:=`(outcome = o, horizon = h, 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) { +#' Patients are weighted by their probability of consenting, so the truth refers to the +#' population that enrols. Horizons: 60 months (estimand rct__itt) and the +#' shorter trial horizons (rct__itt_h12, _h24). +rct_truth <- function(par, design = "realistic", n = 1e6, seed = 20261029L, horizons = RCT_HORIZONS) { pop <- truth_population(par, n = n, seed = seed) pat <- if (design == "realistic") pop$pat[phys_id %in% trial_physicians(pop$phys)] else pop$pat + pc <- consent_prob(pat, RCT_DESIGNS[[design]]) 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)) + w <- pc / sum(pc) + rbindlist(lapply(names(OUTCOMES), function(o) rbindlist(lapply(horizons, function(h) { + d <- as.numeric(by_h(runs[[1]][[OUTCOMES[[o]]]], h)) - as.numeric(by_h(runs[[2]][[OUTCOMES[[o]]]], h)) + m <- sum(w * d) 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)) - })) + estimand = paste0("rct_", design, "_itt", if (h != 60L) paste0("_h", h)), value = m, + mcse = sqrt(sum(w^2 * (d - m)^2)), n = length(d)) + })))) } diff --git a/pvisim/R/records.R b/pvisim/R/records.R index e67dba3..630d8ea 100644 --- a/pvisim/R/records.R +++ b/pvisim/R/records.R @@ -32,6 +32,8 @@ make_records <- function(pop, course, par) { idx <- cbind(id, t + 1L) if (identical(ms$severity, "exact")) { sev <- course$z[idx] + } else if (identical(ms$severity, "noisy")) { + sev <- course$z[idx] + rnorm(length(id), 0, ms$noise_sd) # EHR/registry: measured with error } else { w <- ms$proxy_window cs <- t(apply(course$counts, 1, function(x) { x[is.na(x)] <- 0L; cumsum(x) })) diff --git a/pvisim/inst/params/dgp_default.yml b/pvisim/inst/params/dgp_default.yml index 4968d1e..20a2a8b 100644 --- a/pvisim/inst/params/dgp_default.yml +++ b/pvisim/inst/params/dgp_default.yml @@ -83,7 +83,8 @@ negative_controls: # never depend on PVI 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 + severity: exact # exact: analyst sees Z_t. proxy: only visit/testing counts. noisy: Z_t + error (EHR/registry) + noise_sd: 0.75 # noisy: SD of the monthly measurement error; Z_t has SD about 1 at diagnosis, so cor(observed, true) is about 0.8 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) @@ -113,6 +114,8 @@ worlds: # one factor at a time (plan section 5.2, D24). Ove calibration_targets: {late_not: 0.15} W5: measurement: {severity: proxy, frailty_reliability: 0.4, smoking: coded} + WE: # EHR/registry-like: severity recorded (ABI, symptom class) but with error (D55) + measurement: {severity: noisy, frailty_reliability: 0.7, smoking: true} 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 diff --git a/reports/draft-report.Rmd b/reports/draft-report.Rmd index 95bc348..cc52b70 100644 --- a/reports/draft-report.Rmd +++ b/reports/draft-report.Rmd @@ -17,9 +17,9 @@ S <- file.path(root, "results", "summary") perf <- fread(file.path(S, "perf.csv")); nc <- fread(file.path(S, "nc_reps.csv")) ncs <- fread(file.path(S, "nc_summary.csv")); reps <- fread(file.path(S, "reps.csv")) -WORLDS <- c("W0", "W1", "W1_rho0", "W2", "W3", "W4", "W5", "WR") +WORLDS <- c("W0", "W1", "W1_rho0", "W2", "W3", "W4", "W5", "WE", "WR") WLAB <- c(W0 = "W0 Clean", W1 = "W1 Skill", W1_rho0 = "W1 Skill (rho 0)", W2 = "W2 Style", - W3 = "W3 Style -> outcome", W4 = "W4 Deferral", W5 = "W5 Unmeasured severity", WR = "WR Realistic") + W3 = "W3 Style -> outcome", W4 = "W4 Deferral", W5 = "W5 Unmeasured severity", WE = "WE EHR/registry-like", WR = "WR Realistic") VAR <- c(V2_cluster = "Ignore physician (cluster SE)", V3_random = "Random intercept", V4_within = "Within-physician (FE)", none = "Trial (no adjustment)") VCOL <- c("Ignore physician (cluster SE)" = "#2a78d6", "Random intercept" = "#eb6834", @@ -35,7 +35,9 @@ pp <- function(x, k = 1) formatC(100 * x, format = "f", digits = k) theme_set(theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank(), panel.grid.major.y = element_blank(), strip.text = element_text(face = "bold", hjust = 0), legend.position = "bottom")) -perf <- perf[analysis %in% c("main", "rct")] # NCE analyses are used only in Table 4 +perf_h <- perf[analysis == "rct"] # trials at 12, 24 and 60 months (D55) +perf <- perf[analysis %in% c("main", "rct") & horizon == 60L] # NCE analyses are used only in Table 4 +reps <- reps[horizon == 60L] perf[, key := paste(design, contrast, sep = "|")] perf <- merge(perf, METH, by = "key") perf[, method := factor(method, levels = rev(METH$method))] diff --git a/simulations/04_summarise.R b/simulations/04_summarise.R index 5fe116f..35ec84f 100644 --- a/simulations/04_summarise.R +++ b/simulations/04_summarise.R @@ -14,6 +14,9 @@ cat(sprintf("%d files, %d rows\n", length(files), nrow(d))) # rare failures (sequential V4: fixed effects not estimable in a trial month) are dropped and counted cat(sprintf("dropped %d estimates with NA\n", d[is.na(est), .N])) d <- d[!is.na(est)] +# trials are analysed at 12, 24 and 60 months (D55); everything else at 60 +if (!"horizon" %in% names(d)) d[, horizon := NA_integer_] +d[is.na(horizon), horizon := 60L] tr <- rbind(fread("simulations/truths/truths.csv"), fread("simulations/truths/truths_rct.csv")) tr <- tr[, .(world, scenario, outcome, estimand, truth = value, truth_mcse = mcse)] @@ -23,12 +26,19 @@ tr <- tr[, .(world, scenario, outcome, estimand, truth = value, truth_mcse = mcs d[, estimand := fcase(design == "ccw" & contrast == "A_B1", "rd_A_B1", design == "ccw" & contrast == "A_B2", "rd_A_B2", design == "seq", "seq_itt", - design %like% "^rct_" & contrast == "itt", paste0(design, "_itt"), + design %like% "^rct_" & contrast == "itt", + paste0(design, "_itt", fifelse(horizon == 60L, "", paste0("_h", horizon))), default = NA_character_)] d <- merge(d, tr[, .(world, scenario, outcome, estimand, own = truth)], by = c("world", "scenario", "outcome", "estimand"), all.x = TRUE) d <- merge(d, tr[estimand == "rd_A_B1", .(world, scenario, outcome, ref = truth)], by = c("world", "scenario", "outcome"), all.x = TRUE) +# at the shorter trial horizons the reference is RD(A - B1) at that month, from the risk curves +cv <- dcast(fread("simulations/truths/risk_curves.csv")[policy %in% c("A", "B1")], + world + scenario + outcome + month ~ policy, value.var = "cuminc") +d <- merge(d, cv[, .(world, scenario, outcome, horizon = month, ref_h = A - B1)], + by = c("world", "scenario", "outcome", "horizon"), all.x = TRUE) +d[horizon != 60L, ref := ref_h][, ref_h := NULL] # 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) @@ -39,7 +49,7 @@ d[, `:=`(lo = est - 1.96 * se, hi = est + 1.96 * se)] fwrite(d, file.path(outdir, "reps.csv")) # ---- performance by cell (plan 9) ----------------------------------------------------- -keys <- c("world", "scenario", "analysis", "design", "contrast", "variant", "outcome") +keys <- c("world", "scenario", "analysis", "design", "contrast", "variant", "outcome", "horizon") perf <- d[variant != "oracle", { err <- est - own; R <- .N; own_m <- mean(own) .(reps = R, mean_est = mean(est), own = own_m, ref = ref[1], @@ -83,8 +93,8 @@ ncs <- rbindlist(lapply(c("nco1", "nco2", "nce", "any_nc"), function(k) rbind( nc[, c(control = k, level = "cell", acc(get(k), biased_cell)), by = world], nc[, c(world = "all", control = k, level = "replicate", acc(get(k), biased_rep))], nc[, c(world = "all", control = k, level = "cell", acc(get(k), biased_cell))]))) -hide <- nc[any_nc == FALSE, .(control = "any_nc", med_absbias_clean = median(abs(cell_bias)), - p90_absbias_clean = quantile(abs(cell_bias), 0.9)), by = world] +hide <- nc[any_nc == FALSE, .(control = "any_nc", med_absbias_clean = median(abs(cell_bias), na.rm = TRUE), + p90_absbias_clean = quantile(abs(cell_bias), 0.9, na.rm = TRUE)), by = world] ncs <- merge(ncs, hide, by = c("world", "control"), all.x = TRUE) fwrite(ncs, file.path(outdir, "nc_summary.csv")) cat("wrote", outdir, "\n") diff --git a/simulations/config/calibrated.yml b/simulations/config/calibrated.yml index 17bad7e..b58b33a 100644 --- a/simulations/config/calibrated.yml +++ b/simulations/config/calibrated.yml @@ -102,3 +102,16 @@ WR: intercept: -7.5233 death: intercept: -5.6489 +WE: + treatment: + window: + intercept: -6.3552 + deferred: + intercept: -7.7172 + outcomes: + repeat_revasc: + intercept: -5.6375 + clti: + intercept: -7.4903 + death: + intercept: -5.6504 diff --git a/simulations/config/experiments.yml b/simulations/config/experiments.yml index 774324d..5a48f24 100644 --- a/simulations/config/experiments.yml +++ b/simulations/config/experiments.yml @@ -1,6 +1,6 @@ # Experiment grids (plan sections 5.2, 7, 9). Run with simulations/run_experiment.R [--scale draft|final] main: - worlds: [W0, W1, W1_rho0, W2, W3, W4, W5, WR] + worlds: [W0, W1, W1_rho0, W2, W3, W4, W5, WE, WR] # WE: EHR/registry-like measurement (D55) scenarios: ["null", harm, benefit_harm] 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)