101 lines
6.0 KiB
R
101 lines
6.0 KiB
R
# Performance against the known truths (plan sections 9 and 10).
|
|
# Usage: Rscript simulations/04_summarise.R [results_dir]
|
|
# Writes results/summary/{reps,perf,nc_reps,nc_summary}.csv.
|
|
|
|
suppressPackageStartupMessages(library(data.table))
|
|
a <- commandArgs(TRUE)
|
|
indir <- if (length(a) >= 1) a[1] else "results/main/n50k_v2"
|
|
outdir <- "results/summary"
|
|
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
|
|
|
|
files <- list.files(indir, pattern = "^rep[0-9]+\\.csv$", recursive = TRUE, full.names = TRUE)
|
|
d <- rbindlist(lapply(files, fread), fill = TRUE)
|
|
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)]
|
|
|
|
# ---- own estimand and reference for every estimate -----------------------------------
|
|
# reference (clinically relevant) = RD(A - B1); NCE has a true effect of exactly 0 (D32)
|
|
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", 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)
|
|
d[design %in% c("landmark", "landmark_m5", "naive") & variant != "oracle", own := oracle]
|
|
d[analysis %like% "^nce", `:=`(own = 0, ref = 0)]
|
|
d[, oracle := NULL]
|
|
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", "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],
|
|
bias = mean(err), mcse_bias = sd(err) / sqrt(R),
|
|
gap = own_m - ref[1], total = mean(est) - ref[1],
|
|
emp_se = sd(est), mod_se = sqrt(mean(se^2)),
|
|
se_ratio = sqrt(mean(se^2)) / sd(est),
|
|
cover_own = mean(lo <= own_m & hi >= own_m),
|
|
cover_ref = mean(lo <= ref[1] & hi >= ref[1]),
|
|
reject0 = mean(lo > 0 | hi < 0),
|
|
rmse_own = sqrt(mean((est - own_m)^2)), rmse_ref = sqrt(mean((est - ref[1])^2)))
|
|
}, by = keys]
|
|
perf[, `:=`(mcse_cover = sqrt(cover_own * (1 - cover_own) / reps),
|
|
mcse_reject = sqrt(reject0 * (1 - reject0) / 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 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 == "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 == "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)]
|
|
nc[, `:=`(biased_rep = abs(err) > 0.01, biased_cell = abs(cell_bias) > 0.01)]
|
|
fwrite(nc, file.path(outdir, "nc_reps.csv"))
|
|
|
|
acc <- function(x, truth) {
|
|
ok <- !is.na(x)
|
|
list(sens = mean(x[ok & truth]), spec = mean(!x[ok & !truth]),
|
|
n_biased = sum(ok & truth), n_clean = sum(ok & !truth))
|
|
}
|
|
ncs <- rbindlist(lapply(c("nco1", "nco2", "nce", "any_nc"), function(k) rbind(
|
|
nc[, c(control = k, level = "replicate", acc(get(k), biased_rep)), by = world],
|
|
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), 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")
|