chendun1995andClaude Opus 5.5
fc1461d5b5
- 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]>
91 lines
4.8 KiB
R
91 lines
4.8 KiB
R
# 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))
|
|
}))
|
|
}
|