93 lines
4.6 KiB
R
93 lines
4.6 KiB
R
# 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(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) {
|
|
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("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,
|
|
as.numeric(Sys.time() - t0, units = "mins")))
|
|
list(truths = rbindlist(rows), curves = curves)
|
|
}
|