# 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"
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)]

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"),
                      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)
# 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", "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")
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 MALE RD
# against the method's own estimand; a control "flags" when its 95% CI excludes 0.
m <- d[analysis == "main" & 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 == "male", 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 == "male" & 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)),
                              p90_absbias_clean = quantile(abs(cell_bias), 0.9)), 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")