From fc1461d5b58bbcd9625561c1895419a308af9aa0 Mon Sep 17 00:00:00 2001 From: Chen Date: Tue, 27 Oct 2026 09:01:13 -0400 Subject: [PATCH] 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 --- pvisim/R/estimators.R | 302 +++++++++++++++++++++++++++++ pvisim/R/rct.R | 90 +++++++++ pvisim/R/simulate.R | 14 +- pvisim/R/truths.R | 92 +++++++++ pvisim/inst/params/dgp_default.yml | 5 +- simulations/02_truths.R | 20 ++ simulations/explore_scenarios.R | 17 ++ simulations/validate_designs.R | 33 ++++ 8 files changed, 566 insertions(+), 7 deletions(-) create mode 100644 pvisim/R/estimators.R create mode 100644 pvisim/R/rct.R create mode 100644 pvisim/R/truths.R create mode 100644 simulations/02_truths.R create mode 100644 simulations/explore_scenarios.R create mode 100644 simulations/validate_designs.R diff --git a/pvisim/R/estimators.R b/pvisim/R/estimators.R new file mode 100644 index 0000000..a7c829b --- /dev/null +++ b/pvisim/R/estimators.R @@ -0,0 +1,302 @@ +# 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 + pm <- merge(rec$pm, p[, c("id", "phys_id", BASE_COVS, "t_pvi", "pvi_elective", "t_male", "t_death"), + 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))] + 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) { + 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)" + out[i_new] <- fit_pvi_phase(d, nd, variant, tt_levels) + } + out +} + +fit_pvi_phase <- function(d, nd, variant, time_terms) { + covs <- c(BASE_COVS, "sev_obs", "sev0") + 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) else NULL + 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)), + ccw = est_ccw(rec, outcomes[[o]], variant, rs = copy(rs), phat = phat), + 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) { + 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")) # 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) { + 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) + 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_) +} diff --git a/pvisim/R/rct.R b/pvisim/R/rct.R new file mode 100644 index 0000000..d1c3cae --- /dev/null +++ b/pvisim/R/rct.R @@ -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)) + })) +} diff --git a/pvisim/R/simulate.R b/pvisim/R/simulate.R index c2d6f62..947b4b3 100644 --- a/pvisim/R/simulate.R +++ b/pvisim/R/simulate.R @@ -72,12 +72,16 @@ make_population <- function(par, n = par$size$n_patients, seed = 1L, phys = NULL # B2 no elective PVI ever # seq_force / seq_forbid (m) natural except PVI forced / forbidden at month m -policy <- function(type = "natural", m = NA_integer_) { - stopifnot(type %in% c("natural", "A", "B1", "B2", "seq_force", "seq_forbid")) - list(type = type, m = 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) { +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, @@ -196,7 +200,7 @@ simulate_course <- function(pop, par, pol = policy("natural"), seed = 1L, record 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) + 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 diff --git a/pvisim/R/truths.R b/pvisim/R/truths.R new file mode 100644 index 0000000..b51eca1 --- /dev/null +++ b/pvisim/R/truths.R @@ -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) +} diff --git a/pvisim/inst/params/dgp_default.yml b/pvisim/inst/params/dgp_default.yml index b6a74e3..c94f50e 100644 --- a/pvisim/inst/params/dgp_default.yml +++ b/pvisim/inst/params/dgp_default.yml @@ -121,5 +121,6 @@ worlds: # one factor at a time (plan section 5.2, D24). Ove 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.03}} # tuned to RD about +5 pp (calibration) - benefit_harm: {effect: {relief: 1.2, kappa: 0.04}} # curves cross; tuned in calibration + 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 diff --git a/simulations/02_truths.R b/simulations/02_truths.R new file mode 100644 index 0000000..e27cbb6 --- /dev/null +++ b/simulations/02_truths.R @@ -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") diff --git a/simulations/explore_scenarios.R b/simulations/explore_scenarios.R new file mode 100644 index 0000000..159aa63 --- /dev/null +++ b/simulations/explore_scenarios.R @@ -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))) +} diff --git a/simulations/validate_designs.R b/simulations/validate_designs.R new file mode 100644 index 0000000..82b880b --- /dev/null +++ b/simulations/validate_designs.R @@ -0,0 +1,33 @@ +# 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) + +tr <- fread("simulations/truths/truths.csv")[world == ..world & scenario == ..scenario] +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))