# The review's registered analysis: does mass azithromycin lower mortality among young children when all
# children 1 to 59 months are treated, when only infants 1 to 11 months are, and which lowers infant mortality more?
#
#   Rscript code/analyze.R [--blind] [data dir] [results dir]
#
# Reads data/trials.json (each trial's decided design: its comparisons and their strategies, masking, and
# integrity facts), data/extraction.json (every mortality estimate the trials' sources print, as decided, with
# each entry's age band in months and period), data/parameters.json, and data/registry.json. Writes
# results/R*.json (declared results), results/estimates.csv (every estimate with its log ratio), and
# results/analyses.csv (every model fitted). --blind flips the sign of each comparison's log ratios at random
# (seed in parameters.json), so the code can be run and checked before the numbers mean anything.

suppressPackageStartupMessages({
  library(metafor)
  library(jsonlite)
})

args <- commandArgs(trailingOnly = TRUE)
blind <- "--blind" %in% args
args <- setdiff(args, "--blind")
data_dir <- if (length(args) >= 1) args[1] else "data"
out_dir <- if (length(args) >= 2) args[2] else "results"
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

P <- fromJSON(file.path(data_dir, "parameters.json"), simplifyVector = FALSE)
trials <- fromJSON(file.path(data_dir, "trials.json"), simplifyVector = FALSE)
extraction <- fromJSON(file.path(data_dir, "extraction.json"), simplifyVector = FALSE)
registry <- fromJSON(file.path(data_dir, "registry.json"), simplifyVector = FALSE)

`%||%` <- function(a, b) if (is.null(a) || length(a) == 0) b else a
num <- function(x) if (is.null(x) || length(x) == 0) NA_real_ else suppressWarnings(as.numeric(x))
`%|na|%` <- function(a, b) if (length(a) == 0 || is.na(a[1])) b else a
fill_na <- function(x, b) { x[is.na(x)] <- b; x }
set.seed(P$seed)

# ---------------------------------------------------------------------------------------------------
# Comparisons: one row per randomized comparison of an eligible trial.

comparisons <- do.call(rbind, lapply(Filter(function(t) isTRUE(t$eligible), trials), function(t) {
  do.call(rbind, lapply(t$comparisons, function(c) data.frame(
    trial = t$trial, name = t$name, comparison = c$id, strategy = c$strategy, control = c$control %||% NA_character_,
    frequency = c$frequency %||% NA_character_, label = c$label, subgroup = c$subgroup %||% NA_character_,
    treated_lo = num(c$treated_ages[[1]]), treated_hi = num(c$treated_ages[[2]]),
    shares_control_with = c$shares_control_with %||% NA_character_,
    masked = identical(t$masking, "double") && identical(c$control, "placebo"),
    integrity_ok = !isFALSE(t$integrity$ok), stringsAsFactors = FALSE)))
}))
sign_of <- setNames(ifelse(blind & runif(nrow(comparisons)) < 0.5, -1, 1), comparisons$comparison)

# A trial randomized within countries also has a comparison for all its countries together (`fallback_for`
# lists the country comparisons it stands for): it is used for an age band or period only where none of its
# countries has an estimate of its own, and it never enters beside them.
fallback <- list()
for (t in trials) for (c in t$comparisons) if (!is.null(c$fallback_for)) fallback[[c$id]] <- unlist(c$fallback_for)
add_fallback <- function(cmps, has) {
  out <- cmps
  for (parent in names(fallback)) {
    members <- intersect(fallback[[parent]], cmps)
    if (length(members) && !any(vapply(members, has, TRUE)) && has(parent)) out <- c(setdiff(out, members), parent)
  }
  out
}

# ---------------------------------------------------------------------------------------------------
# Estimates: one row per extracted entry, with its log ratio (azithromycin against control) and variance.

z_of <- function(level) qnorm(1 - (1 - (num(level) %|na|% 0.95)) / 2)

from_interval <- function(est) {
  if (is.null(est) || is.null(est$value)) return(NULL)
  v <- num(est$value); lo <- num(est$lower); hi <- num(est$upper)
  if (identical(est$measure, "percent reduction")) {
    v <- 1 - v / 100; l2 <- 1 - hi / 100; h2 <- 1 - lo / 100; lo <- l2; hi <- h2
  }
  if (any(is.na(c(v, lo, hi))) || lo <= 0 || hi <= 0 || v <= 0) return(NULL)
  c(yi = log(v), sei = (log(hi) - log(lo)) / (2 * z_of(est$level)))
}

from_counts <- function(t, c) {
  dt <- num(t$deaths); dc <- num(c$deaths)
  pt <- num(t$person_years); pc <- num(c$person_years)
  nt <- num(t$children); nc <- num(c$children)
  if (!any(is.na(c(dt, dc, pt, pc))) && dt > 0 && dc > 0) {
    return(c(yi = log((dt / pt) / (dc / pc)), var_poisson = 1 / dt + 1 / dc, dt = dt, dc = dc))
  }
  if (!any(is.na(c(dt, dc, nt, nc))) && dt > 0 && dc > 0) {
    return(c(yi = log((dt / nt) / (dc / nc)), var_poisson = 1 / dt - 1 / nt + 1 / dc - 1 / nc, dt = dt, dc = dc))
  }
  NULL
}

rows <- list()
for (tr in extraction) {
  for (i in seq_along(tr$entries)) {
    e <- tr$entries[[i]]
    ci <- from_interval(e$estimate)
    ct <- from_counts(e$treated %||% list(), e$control %||% list())
    rate_c <- num(e$control$rate_per_1000_person_years) %|na|% (1000 * num(e$control$deaths) / num(e$control$person_years))
    rows[[length(rows) + 1]] <- data.frame(
      trial = tr$trial, entry = paste0(tr$trial, "#", i), comparison = e$comparison_id %||% NA_character_,
      ages_lo = num(e$ages_lo), ages_hi = num(e$ages_hi), period = e$period %||% "whole",
      period_mid = num(e$period_mid_years), subgroup = e$subgroup %||% NA_character_,
      primary = isTRUE(e$primary), clustered = !isFALSE(e$adjusted_for_clustering), role = e$source_role %||% "secondary",
      yi_ci = if (is.null(ci)) NA_real_ else ci[["yi"]], sei_ci = if (is.null(ci)) NA_real_ else ci[["sei"]],
      yi_ct = if (is.null(ct)) NA_real_ else ct[["yi"]], var_poisson = if (is.null(ct)) NA_real_ else ct[["var_poisson"]],
      dt = if (is.null(ct)) NA_real_ else ct[["dt"]], dc = if (is.null(ct)) NA_real_ else ct[["dc"]],
      control_rate = rate_c, stringsAsFactors = FALSE)
  }
}
est <- do.call(rbind, rows)

# A trial's design effect: the variance of its primary estimate as printed with its interval, over the Poisson
# variance of the same estimate from its counts; the median over trials where a trial gives no such pair.
de_trial <- sapply(split(est, est$trial), function(d) {
  ok <- d$primary & !is.na(d$sei_ci) & !is.na(d$var_poisson) & d$clustered
  if (any(ok)) median(d$sei_ci[ok]^2 / d$var_poisson[ok]) else NA_real_
})
de_median <- median(de_trial, na.rm = TRUE) %|na|% 1
design_effect <- function(trial) if (trial %in% names(de_trial) && !is.na(de_trial[[trial]])) de_trial[[trial]] else de_median

est$yi <- ifelse(!is.na(est$yi_ci), est$yi_ci, est$yi_ct)
est$vi <- ifelse(!is.na(est$sei_ci), est$sei_ci^2, est$var_poisson * vapply(est$trial, design_effect, 1))
est$how <- ifelse(!is.na(est$sei_ci), "interval", ifelse(!is.na(est$yi_ct), "counts", "none"))
est <- est[est$how != "none", ]
est$yi <- est$yi * fill_na(unname(sign_of[est$comparison]), 1)
write.csv(est, file.path(out_dir, "estimates.csv"), row.names = FALSE)

# ---------------------------------------------------------------------------------------------------
# Choosing one estimate per comparison, age band, and period: whole follow-up, no subgroup beyond the
# comparison's own (a country), the primary analysis first, then estimates adjusted for clustering, then the
# main report, then the registry, then other reports.

role_rank <- c(main = 1, extension = 2, registry = 3, secondary = 4)
pick <- function(d) {
  if (nrow(d) == 0) return(NULL)
  d <- d[order(!d$primary, !d$clustered, fill_na(unname(role_rank[d$role]), 5), d$vi), ]
  d[1, ]
}
same_subgroup <- function(d, cmp) {
  sg <- comparisons$subgroup[comparisons$comparison == cmp]
  if (is.na(sg)) is.na(d$subgroup) else (!is.na(d$subgroup) & tolower(d$subgroup) == tolower(sg))
}
band <- function(cmp, lo, hi, period = "whole") {
  d <- est[est$comparison %in% cmp & est$period == period, ]
  d <- d[same_subgroup(d, cmp), ]
  # An estimate stands for the band [lo, hi] when its own band starts within a month below lo and ends within a
  # month above hi (so 1-60 months stands for 1-59, but 1-59 never stands for 12-59).
  full <- d[!is.na(d$ages_lo) & d$ages_lo <= lo & d$ages_lo >= lo - 1 & d$ages_hi >= hi & d$ages_hi <= hi + 1, ]
  got <- pick(full)
  if (!is.null(got)) return(got)
  # Narrower bands that together cover [lo, hi] without overlapping, combined by inverse variance: starting
  # at lo, each next band starts where the last ended, taking the widest such band at each step.
  parts <- d[!is.na(d$ages_lo) & d$ages_lo >= lo - 1 & d$ages_hi <= hi + 1, ]
  if (nrow(parts) < 2) return(NULL)
  chain <- list(); at <- lo
  repeat {
    nxt <- parts[parts$ages_lo >= at - 1 & parts$ages_lo <= at + 0.5, ]
    if (length(chain)) nxt <- nxt[nxt$ages_lo > chain[[length(chain)]]$ages_hi - 0.5, ]
    if (nrow(nxt) == 0) return(NULL)
    nxt <- nxt[order(-nxt$ages_hi, !nxt$primary, nxt$vi), ][1, ]
    chain[[length(chain) + 1]] <- nxt
    if (nxt$ages_hi >= hi - 1) break
    at <- nxt$ages_hi + 1
  }
  if (length(chain) < 2) return(NULL)
  parts <- do.call(rbind, chain)
  w <- 1 / parts$vi
  out <- parts[1, ]
  out$yi <- sum(w * parts$yi) / sum(w)
  out$vi <- 1 / sum(w)
  out$dt <- sum(parts$dt); out$dc <- sum(parts$dc)
  out$how <- "combined bands"
  out$ages_lo <- min(parts$ages_lo); out$ages_hi <- max(parts$ages_hi)
  out
}

# A comparison's own primary estimate, in whatever age band the trial measured (for trials that aren't
# pooled by band: individual dosing of infants, community-wide treatment, schedule comparisons).
primary_entry <- function(cmp) {
  d <- est[est$comparison %in% cmp & est$period == "whole", ]
  d <- d[same_subgroup(d, cmp), ]
  b <- pick(d)
  if (!is.null(b)) b$comparison <- cmp
  b
}

collect <- function(cmps, lo, hi) {
  cmps <- add_fallback(cmps, function(cmp) !is.null(band(cmp, lo, hi)))
  out <- lapply(cmps, function(cmp) {
    b <- band(cmp, lo, hi)
    if (is.null(b)) return(NULL)
    b$comparison <- cmp
    b
  })
  do.call(rbind, Filter(Negate(is.null), out))
}

# ---------------------------------------------------------------------------------------------------
# Models.

analyses <- list()
record <- function(name, fit, coef = 1) {
  analyses[[length(analyses) + 1]] <<- data.frame(analysis = name, k = fit$k, estimate = exp(as.numeric(fit$beta[coef])),
    lower = exp(fit$ci.lb[coef]), upper = exp(fit$ci.ub[coef]), p = fit$pval[coef], tau2 = fit$tau2 %||% NA, stringsAsFactors = FALSE)
  fit
}
re <- function(d, name, mods = NULL) {
  if (is.null(mods)) return(record(name, rma(d$yi, d$vi, method = "REML", test = "knha")))
  record(name, rma(d$yi, d$vi, mods = mods, method = "REML", test = "knha"), coef = 2)
}
fe <- function(d, name) record(name, rma(d$yi, d$vi, method = "FE"))

bayes <- function(d, scale = P$bayes_tau_scale) {
  # Random-effects model with a half-normal prior of the given scale on tau and a flat prior on the mean,
  # integrated over a grid of tau: the mean's posterior is a mixture of normals.
  taus <- seq(0, 4 * scale, length.out = 4001)
  lp <- vapply(taus, function(tau) {
    w <- 1 / (d$vi + tau^2)
    m <- sum(w * d$yi) / sum(w)
    0.5 * sum(log(w)) - 0.5 * log(sum(w)) - 0.5 * sum(w * (d$yi - m)^2)
  }, 1) - taus^2 / (2 * scale^2)
  post <- exp(lp - max(lp)); post <- post / sum(post)
  means <- vapply(taus, function(tau) { w <- 1 / (d$vi + tau^2); sum(w * d$yi) / sum(w) }, 1)
  sds <- vapply(taus, function(tau) sqrt(1 / sum(1 / (d$vi + tau^2))), 1)
  cdf <- function(x) sum(post * pnorm(x, means, sds))
  q <- function(p) uniroot(function(x) cdf(x) - p, c(-5, 5), tol = 1e-10)$root
  list(estimate = exp(q(0.5)), lower = exp(q(0.025)), upper = exp(q(0.975)), prob_benefit = cdf(0))
}

loo <- function(d) {
  if (nrow(d) < 3) return(list(max_p = NA_real_, estimates = list()))
  l <- leave1out(rma(d$yi, d$vi, method = "REML", test = "knha"))
  list(max_p = max(l$pval), min_estimate = exp(min(l$estimate)), max_estimate = exp(max(l$estimate)),
       left_out = setNames(as.list(round(exp(l$estimate), 3)), d$comparison))
}

summary_of <- function(fit, d, digits = 3) {
  pi <- predict(fit)
  list(k = fit$k, trials = length(unique(d$trial)), comparisons = d$comparison,
       rate_ratio = round(exp(as.numeric(fit$beta[1])), digits), lower = round(exp(fit$ci.lb[1]), digits), upper = round(exp(fit$ci.ub[1]), digits),
       reduction_pct = round(100 * (1 - exp(as.numeric(fit$beta[1]))), 1),
       reduction_lower_pct = round(100 * (1 - exp(fit$ci.ub[1])), 1), reduction_upper_pct = round(100 * (1 - exp(fit$ci.lb[1])), 1),
       p = signif(fit$pval[1], 3), tau2 = signif(fit$tau2, 3), i2_pct = round(fit$I2, 1),
       pi_lower = round(exp(pi$pi.lb), digits), pi_upper = round(exp(pi$pi.ub), digits),
       deaths = sum(d$dt + d$dc, na.rm = TRUE))
}

p_text <- function(p) if (is.na(p)) NA_character_ else if (p < 0.001) "p < 0.001" else sprintf("p = %.3f", p)

primary_cmp <- function(strategy, freq = NULL) {
  ok <- comparisons$strategy == strategy & comparisons$control %in% "placebo" & !(comparisons$comparison %in% names(fallback))
  if (!is.null(freq)) ok <- ok & comparisons$frequency %in% freq
  comparisons$comparison[ok]
}

# ---------------------------------------------------------------------------------------------------
# P1: all children 1 to 59 months.

p1_cmps <- primary_cmp("all_children_mda", "biannual")
d1 <- collect(p1_cmps, 1, 59)
f1 <- re(d1, "P1")
f1_fe <- fe(d1, "P1 common effect")
b1 <- bayes(d1)
l1 <- loo(d1)
d1m <- d1[comparisons$masked[match(d1$comparison, comparisons$comparison)], ]
f1_m <- if (nrow(d1m) >= 2) re(d1m, "P1 double-masked placebo") else NULL
no_rerand <- d1[!(d1$comparison %in% P$p1_rerandomized), ]
f1_nr <- if (nrow(no_rerand) >= 2 && nrow(no_rerand) < nrow(d1)) re(no_rerand, "P1 without re-randomized") else NULL

# P2: infants only, twice yearly.

p2_cmps <- primary_cmp("infant_only_mda", "biannual")
d2 <- collect(p2_cmps, 1, 11)
f2 <- re(d2, "P2")
f2_fe <- fe(d2, "P2 common effect")
b2 <- bayes(d2)
l2 <- loo(d2)

# P3: infant mortality under each strategy. Comparisons of one trial that share a control arm have their
# variance recomputed with half the control's deaths, times the trial's design effect.

d3a <- collect(p1_cmps, 1, 11)
d3 <- rbind(transform(d3a, infant_only = 0), transform(d2, infant_only = 1))
split_shared <- function(d) {
  partner <- comparisons$shares_control_with[match(d$comparison, comparisons$comparison)]
  for (i in which(!is.na(partner) & partner %in% d$comparison)) {
    d$vi[i] <- if (!is.na(d$dt[i]) && !is.na(d$dc[i])) design_effect(d$trial[i]) * (1 / d$dt[i] + 2 / d$dc[i]) else d$vi[i] * 1.5
  }
  d
}
d3 <- split_shared(d3)
f3 <- re(d3, "P3", mods = ~ d3$infant_only)
f3_fe <- record("P3 common effect", rma(d3$yi, d3$vi, mods = ~ d3$infant_only, method = "FE"), coef = 2)

holm <- p.adjust(c(f1$pval[1], f2$pval[1], f3$pval[2]), method = "holm")

# AVENIR's within-trial contrast: the child arm against the infant arm, infant mortality.
within <- intersect(P$within_trial_comparison %||% character(), comparisons$comparison)
w3 <- if (length(within)) band(within[1], 1, 11) else NULL

# ---------------------------------------------------------------------------------------------------
# Sensitivity analyses of P2.

p2_quarterly <- c(setdiff(p2_cmps, P$p2_twice_yearly_in_quarterly_trial), P$p2_quarterly)
d2q <- collect(p2_quarterly, 1, 11)
f2q <- if (nrow(d2q) >= 2) re(d2q, "P2 with quarterly arm") else NULL
both <- c(p2_cmps, P$p2_quarterly)
d2b <- split_shared(collect(both, 1, 11))
f2b <- if (nrow(d2b) >= 2) re(d2b, "P2 with both LAKANA arms") else NULL
d2nb <- d2[!(d2$comparison %in% P$p2_birth_history), ]
f2nb <- if (nrow(d2nb) >= 2 && nrow(d2nb) < nrow(d2)) re(d2nb, "P2 without birth-history trial") else NULL
ind_cmps <- comparisons$comparison[comparisons$strategy == "individual_infant" & comparisons$control %in% "placebo"]
d_ind <- do.call(rbind, lapply(ind_cmps, primary_entry))
d2i <- if (!is.null(d_ind)) rbind(d2, d_ind[, names(d2)]) else d2
f2i <- if (!is.null(d_ind) && nrow(d2i) > nrow(d2)) re(d2i, "P2 with individual infant dosing") else NULL

# ---------------------------------------------------------------------------------------------------
# Secondary analyses.

# Mortality where it is given: log rate ratio on the log of the control arm's rate, with infant_only.
dm <- rbind(transform(d1, infant_only = 0), transform(d2, infant_only = 1))
dm <- dm[!is.na(dm$control_rate) & dm$control_rate > 0, ]
fm <- if (nrow(dm) >= 4) rma(dm$yi, dm$vi, mods = ~ log(dm$control_rate) + dm$infant_only, method = "REML", test = "knha") else NULL

# Durability: period estimates against years since the first distribution.
period_rows <- function(cmp) {
  d <- est[est$comparison == cmp & est$period != "whole" & !is.na(est$period_mid), ]
  d <- d[same_subgroup(d, cmp), ]
  lo <- comparisons$treated_lo[comparisons$comparison == cmp]; hi <- comparisons$treated_hi[comparisons$comparison == cmp]
  d <- d[!is.na(d$ages_lo) & d$ages_lo <= lo & d$ages_lo >= lo - 1 & d$ages_hi >= hi & d$ages_hi <= hi + 1, ]
  if (nrow(d) == 0) return(NULL)
  d <- d[order(d$period_mid, !d$primary, fill_na(unname(role_rank[d$role]), 5)), ]
  d <- d[!duplicated(d$period_mid), ]
  if (nrow(d) == 0) return(NULL)
  d$comparison <- cmp
  d
}
dd <- do.call(rbind, lapply(add_fallback(c(p1_cmps, p2_cmps), function(cmp) !is.null(period_rows(cmp))), period_rows))
fd <- if (!is.null(dd) && nrow(dd) >= 4 && length(unique(dd$comparison)) >= 2) rma.mv(yi, vi, mods = ~ period_mid, random = ~ 1 | comparison, data = dd, method = "REML", test = "t") else NULL

# Age within all-children trials: infant against 12 to 59 months.
da <- do.call(rbind, lapply(add_fallback(p1_cmps, function(cmp) !is.null(band(cmp, 1, 11)) && !is.null(band(cmp, 12, 59))), function(cmp) {
  a <- band(cmp, 1, 11); b <- band(cmp, 12, 59)
  if (is.null(a) || is.null(b)) return(NULL)
  data.frame(comparison = cmp, trial = a$trial, yi = a$yi - b$yi, vi = a$vi + b$vi)
}))
fa <- if (!is.null(da) && nrow(da) >= 2) re(da, "Infants against 12-59 months") else NULL

# Individual dosing of infants.
fi <- if (!is.null(d_ind) && nrow(d_ind) >= 2) re(d_ind, "Individual infant dosing") else NULL

# Frequency: LAKANA's quarterly against twice-yearly arm; PRET as printed.
freq_q <- if (length(P$frequency_comparison)) band(P$frequency_comparison, 1, 11) else NULL
pret <- if (length(P$pret_comparison)) primary_entry(P$pret_comparison) else NULL

# Community-wide and individual-child comparisons, each as printed.
as_printed <- function(strategy) {
  cm <- comparisons[comparisons$strategy == strategy, ]
  lapply(seq_len(nrow(cm)), function(i) {
    b <- primary_entry(cm$comparison[i])
    if (is.null(b)) return(NULL)
    list(comparison = cm$comparison[i], trial = cm$trial[i], ages = c(b$ages_lo, b$ages_hi),
         rate_ratio = round(exp(b$yi), 3), lower = round(exp(b$yi - 1.959964 * sqrt(b$vi)), 3), upper = round(exp(b$yi + 1.959964 * sqrt(b$vi)), 3))
  })
}

# The file drawer.
regs <- registry$registrations
due <- Filter(function(r) isTRUE(r$due), regs)
unreported <- Filter(function(r) !isTRUE(r$reported_mortality %||% r$reported), due)

# ---------------------------------------------------------------------------------------------------
# Declared results.

# Each comparison's estimate as an analysis used it (added after registration: output only).
each <- function(d) setNames(lapply(seq_len(nrow(d)), function(i) list(
  trial = d$trial[i], rate_ratio = round(exp(d$yi[i]), 3), lower = round(exp(d$yi[i] - 1.959964 * sqrt(d$vi[i])), 3),
  upper = round(exp(d$yi[i] + 1.959964 * sqrt(d$vi[i])), 3), how = d$how[i])), gsub("-", "_", d$comparison))

out <- function(name, x) write(toJSON(x, auto_unbox = TRUE, digits = NA, pretty = TRUE, null = "null", na = "null"), file.path(out_dir, paste0(name, ".json")))

r1 <- summary_of(f1, d1)
r1$p_text <- p_text(r1$p)
r1$p_holm <- signif(holm[1], 3); r1$p_holm_text <- p_text(holm[1])
r1$common_effect <- list(rate_ratio = round(exp(f1_fe$beta[1]), 3), lower = round(exp(f1_fe$ci.lb), 3), upper = round(exp(f1_fe$ci.ub), 3), p = signif(f1_fe$pval, 3))
r1$bayes <- lapply(b1, function(x) round(x, 3))
r1$leave_one_out <- l1
r1$double_masked <- if (!is.null(f1_m)) list(k = f1_m$k, rate_ratio = round(exp(f1_m$beta[1]), 3), lower = round(exp(f1_m$ci.lb), 3), upper = round(exp(f1_m$ci.ub), 3), p = signif(f1_m$pval, 3)) else NULL
r1$without_rerandomized <- if (!is.null(f1_nr)) list(k = f1_nr$k, rate_ratio = round(exp(f1_nr$beta[1]), 3), lower = round(exp(f1_nr$ci.lb), 3), upper = round(exp(f1_nr$ci.ub), 3), p = signif(f1_nr$pval, 3)) else NULL
r1$each <- each(d1)
r1$lowers_mortality <- holm[1] < 0.05
r1$robust <- holm[1] < 0.05 && f1_fe$pval < 0.05 && isTRUE(l1$max_p < 0.05)
out("R1", r1)

r2 <- summary_of(f2, d2)
r2$p_text <- p_text(r2$p)
r2$p_holm <- signif(holm[2], 3); r2$p_holm_text <- p_text(holm[2])
r2$common_effect <- list(rate_ratio = round(exp(f2_fe$beta[1]), 3), lower = round(exp(f2_fe$ci.lb), 3), upper = round(exp(f2_fe$ci.ub), 3), p = signif(f2_fe$pval, 3))
r2$bayes <- lapply(b2, function(x) round(x, 3))
r2$leave_one_out <- l2
r2$largest_reduction_allowed_pct <- round(100 * (1 - exp(f2$ci.lb[1])), 1)
r2$excludes_p1_reduction <- exp(f2$ci.lb[1]) > exp(as.numeric(f1$beta[1]))
sens <- function(f) if (is.null(f)) NULL else list(k = f$k, rate_ratio = round(exp(f$beta[1]), 3), lower = round(exp(f$ci.lb[1]), 3), upper = round(exp(f$ci.ub[1]), 3), p = signif(f$pval[1], 3))
r2$with_quarterly_arm <- sens(f2q)
r2$with_both_lakana_arms <- sens(f2b)
r2$without_birth_history <- sens(f2nb)
r2$with_individual_infant <- sens(f2i)
r2$each <- each(d2)
r2$lowers_infant_mortality <- holm[2] < 0.05
out("R2", r2)

r3 <- list(k = f3$k, comparisons = d3$comparison, infant_only = d3$infant_only,
           ratio = round(exp(f3$beta[2]), 3), lower = round(exp(f3$ci.lb[2]), 3), upper = round(exp(f3$ci.ub[2]), 3),
           p = signif(f3$pval[2], 3), p_text = p_text(f3$pval[2]), p_holm = signif(holm[3], 3), p_holm_text = p_text(holm[3]),
           all_children_infant_rate_ratio = round(exp(f3$beta[1]), 3), all_children_infant_lower = round(exp(f3$ci.lb[1]), 3), all_children_infant_upper = round(exp(f3$ci.ub[1]), 3),
           common_effect = list(ratio = round(exp(f3_fe$beta[2]), 3), lower = round(exp(f3_fe$ci.lb[2]), 3), upper = round(exp(f3_fe$ci.ub[2]), 3), p = signif(f3_fe$pval[2], 3)))
r3$within_trial <- if (!is.null(w3)) list(comparison = within[1], rate_ratio = round(exp(w3$yi), 3), lower = round(exp(w3$yi - 1.959964 * sqrt(w3$vi)), 3), upper = round(exp(w3$yi + 1.959964 * sqrt(w3$vi)), 3), p = signif(2 * pnorm(-abs(w3$yi / sqrt(w3$vi))), 3), how = w3$how) else NULL
r3$each <- each(d3)
r3$more_with_all_children <- holm[3] < 0.05 && exp(f3$beta[2]) > 1
r3$robust <- isTRUE(r3$more_with_all_children) && !is.null(w3) && exp(w3$yi) < 1 && 2 * pnorm(-abs(w3$yi / sqrt(w3$vi))) < 0.05
out("R3", r3)

out("R4", if (is.null(fm)) list(k = nrow(dm), fitted = FALSE) else list(k = fm$k, fitted = TRUE,
    slope = round(as.numeric(fm$beta[2]), 3), slope_lower = round(fm$ci.lb[2], 3), slope_upper = round(fm$ci.ub[2], 3), p = signif(fm$pval[2], 3),
    control_rates = round(dm$control_rate, 1), comparisons = dm$comparison))
out("R5", if (is.null(fd)) list(k = if (is.null(dd)) 0 else nrow(dd), fitted = FALSE) else list(k = fd$k, fitted = TRUE,
    slope_per_year = round(as.numeric(fd$beta[2]), 3), slope_lower = round(fd$ci.lb[2], 3), slope_upper = round(fd$ci.ub[2], 3), p = signif(fd$pval[2], 3),
    ratio_per_year = round(exp(as.numeric(fd$beta[2])), 3), comparisons = unique(dd$comparison)))
out("R6", if (is.null(fa)) list(k = if (is.null(da)) 0 else nrow(da), fitted = FALSE) else c(list(fitted = TRUE), summary_of(fa, transform(da, dt = NA, dc = NA))[setdiff(names(summary_of(fa, transform(da, dt = NA, dc = NA))), "deaths")]))
out("R7", if (is.null(fi)) list(k = if (is.null(d_ind)) 0 else nrow(d_ind), fitted = FALSE) else c(list(fitted = TRUE), summary_of(fi, d_ind)))
out("R8", list(
  pret_ages = if (!is.null(pret)) c(pret$ages_lo, pret$ages_hi) else NULL,
  quarterly_vs_twice_yearly = if (!is.null(freq_q)) list(rate_ratio = round(exp(freq_q$yi), 3), lower = round(exp(freq_q$yi - 1.959964 * sqrt(freq_q$vi)), 3), upper = round(exp(freq_q$yi + 1.959964 * sqrt(freq_q$vi)), 3), how = freq_q$how) else NULL,
  pret = if (!is.null(pret)) list(rate_ratio = round(exp(pret$yi), 3), lower = round(exp(pret$yi - 1.959964 * sqrt(pret$vi)), 3), upper = round(exp(pret$yi + 1.959964 * sqrt(pret$vi)), 3), how = pret$how) else NULL,
  community_wide = as_printed("community_wide_mda"),
  individual_child = as_printed("individual_child")))
out("R9", list(eligible = length(regs), due = length(due), reported = length(due) - length(unreported), unreported = length(unreported),
               unreported_ncts = vapply(unreported, function(r) r$nct, "")))
write.csv(do.call(rbind, analyses), file.path(out_dir, "analyses.csv"), row.names = FALSE)
cat("P1", r1$rate_ratio, r1$lower, r1$upper, "| P2", r2$rate_ratio, r2$lower, r2$upper, "| P3", r3$ratio, r3$lower, r3$upper, "\n")