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

w_ <- world; s_ <- scenario
tr <- fread("simulations/truths/truths.csv")[world == w_ & scenario == s_]
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))