# Exploratory and descriptive results added after registration (listed in deviations.json), beside the
# registered analysis in analyze.R. None of these is a planned test.
#
#   Rscript code/summarize.R [data dir] [results dir]
#
# 1. Two comparisons whose sources print each arm's mortality rate with its 95% confidence interval but no ratio
#    and no counts, so the registered code can't use them: MORDOR II (re-randomized, Niger) and AVENIR's child arm
#    against its infant arm on infant mortality. Here each gets the ratio of its two rates, with the standard error
#    of its log from the two intervals (each arm's log rate taken as normal and the arms as independent). For
#    AVENIR's contrast, also the inverse-variance combination of the two ratios its secondary report prints for
#    infants with and without an older sibling.
# 2. P1 with MORDOR II added that way.
# 3. The two individually randomized infant trials combined by a common-effect model, since Knapp-Hartung
#    intervals from two trials are too wide to inform.
# 4. AVENIR's child arm against placebo at 12 to 59 months from its printed rates, and the age analysis with it.
# 5. PRET's estimate for the band its report names as its mortality outcome (6 months to under 5 years).
# Writes results/S.json.

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

args <- commandArgs(trailingOnly = TRUE)
data_dir <- if (length(args) >= 1) args[1] else "data"
out_dir <- if (length(args) >= 2) args[2] else "results"
extraction <- fromJSON(file.path(data_dir, "extraction.json"), simplifyVector = FALSE)
R1 <- fromJSON(file.path(out_dir, "R1.json"), simplifyVector = FALSE)
R3 <- fromJSON(file.path(out_dir, "R3.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))
z <- qnorm(0.975)

entries <- unlist(lapply(extraction, function(t) lapply(t$entries, function(e) c(e, list(trial = t$trial)))), recursive = FALSE)
find <- function(cmp, lo, hi, subgroup = NULL, roles = c("main", "extension", "secondary", "registry")) {
  ok <- Filter(function(e) identical(e$comparison_id, cmp) && identical(num(e$ages_lo), lo) && identical(num(e$ages_hi), hi) &&
                 identical(e$period, "whole") && (if (is.null(subgroup)) is.null(e$subgroup) else identical(e$subgroup, subgroup)), entries)
  ok[order(match(vapply(ok, function(e) e$source_role %||% "secondary", ""), roles))]
}

from_rates <- function(e) {
  rt <- num(e$treated$rate_per_1000_person_years); rc <- num(e$control$rate_per_1000_person_years)
  st <- (log(num(e$treated$rate_upper)) - log(num(e$treated$rate_lower))) / (2 * z)
  sc <- (log(num(e$control$rate_upper)) - log(num(e$control$rate_lower))) / (2 * z)
  y <- log(rt / rc); se <- sqrt(st^2 + sc^2)
  list(yi = y, sei = se, rate_ratio = round(exp(y), 3), lower = round(exp(y - z * se), 3), upper = round(exp(y + z * se), 3),
       p = signif(2 * pnorm(-abs(y / se)), 3), rates = c(rt, rc), where = e$where)
}
with_rates <- function(es) Filter(function(e) !is.na(num(e$treated$rate_lower)) && !is.na(num(e$control$rate_lower)), es)

out <- list()

# 1. MORDOR II and AVENIR's arm-to-arm contrast.
m2 <- with_rates(find("T02-1", 1, 60))
out$mordor2 <- if (length(m2)) from_rates(m2[[1]]) else NULL
aw <- with_rates(find("T03-within", 1, 11))
out$avenir_within_rates <- if (length(aw)) from_rates(aw[[1]]) else NULL
sib <- Filter(function(e) identical(e$comparison_id, "T03-within") && !is.null(e$subgroup) && grepl("sibling", e$subgroup) && !is.null(e$estimate), entries)
if (length(sib) == 2) {
  y <- vapply(sib, function(e) log(num(e$estimate$value)), 1)
  s <- vapply(sib, function(e) (log(num(e$estimate$upper)) - log(num(e$estimate$lower))) / (2 * z), 1)
  w <- 1 / s^2; m <- sum(w * y) / sum(w); se <- sqrt(1 / sum(w))
  out$avenir_within_siblings <- list(rate_ratio = round(exp(m), 3), lower = round(exp(m - z * se), 3), upper = round(exp(m + z * se), 3),
                                     p = signif(2 * pnorm(-abs(m / se)), 3), subgroups = vapply(sib, function(e) e$subgroup, ""))
}

# 2. P1 with MORDOR II.
p1 <- R1$each
y <- vapply(p1, function(x) log(x$rate_ratio), 1)
v <- vapply(p1, function(x) ((log(x$upper) - log(x$lower)) / (2 * z))^2, 1)
if (!is.null(out$mordor2)) {
  f <- rma(c(y, out$mordor2$yi), c(v, out$mordor2$sei^2), method = "REML", test = "knha")
  out$p1_with_mordor2 <- list(k = f$k, rate_ratio = round(exp(as.numeric(f$beta)), 3), lower = round(exp(f$ci.lb), 3), upper = round(exp(f$ci.ub), 3),
                              p = signif(f$pval, 3), p_text = if (f$pval < 0.001) "p < 0.001" else sprintf("p = %.3f", f$pval))
}

# 3. Individual dosing of infants, common effect: each trial's whole follow-up, and NAITRE's 6-month primary outcome instead.
ind <- function(es) {
  y <- vapply(es, function(e) log(num(e$estimate$value)), 1)
  s <- vapply(es, function(e) (log(num(e$estimate$upper)) - log(num(e$estimate$lower))) / (2 * z), 1)
  f <- rma(y, s^2, method = "FE")
  list(k = f$k, rate_ratio = round(exp(as.numeric(f$beta)), 3), lower = round(exp(f$ci.lb), 3), upper = round(exp(f$ci.ub), 3), p = signif(f$pval, 3))
}
pick_est <- function(cmp, period) Filter(function(e) identical(e$comparison_id, cmp) && identical(e$period, period) && is.null(e$subgroup) && !is.null(e$estimate), entries)
t5 <- pick_est("T05-1", "whole"); t6 <- pick_est("T06-1", "whole"); t6p <- Filter(function(e) isTRUE(e$primary), Filter(function(e) identical(e$comparison_id, "T06-1") && is.null(e$subgroup) && !is.null(e$estimate), entries))
if (length(t5) && length(t6)) out$individual_infant_common <- ind(list(t5[[1]], t6[[1]]))
if (length(t5) && length(t6p)) out$individual_infant_common_6_months <- ind(list(t5[[1]], t6p[[1]]))

# 4. AVENIR's child arm against placebo, 12 to 59 months, from its rates; and the age analysis with it.
a1259 <- with_rates(find("T03-child", 12, 59))
out$avenir_child_12_59 <- if (length(a1259)) from_rates(a1259[[1]]) else NULL
r6 <- fromJSON(file.path(out_dir, "R6.json"), simplifyVector = FALSE)
inf <- R3$each[["T03_child"]]
if (!is.null(out$avenir_child_12_59) && !is.null(inf)) {
  # the infant estimate as printed (23% lower), not the split-control variance P3 used
  ai <- Filter(function(e) !is.null(e$estimate), find("T03-child", 1, 11))[[1]]
  yi <- log(1 - num(ai$estimate$value) / 100); si <- (log(1 - num(ai$estimate$lower) / 100) - log(1 - num(ai$estimate$upper) / 100)) / (2 * z)
  out$avenir_infant_vs_older <- list(ratio = round(exp(yi - out$avenir_child_12_59$yi), 3),
                                     lower = round(exp(yi - out$avenir_child_12_59$yi - z * sqrt(si^2 + out$avenir_child_12_59$sei^2)), 3),
                                     upper = round(exp(yi - out$avenir_child_12_59$yi + z * sqrt(si^2 + out$avenir_child_12_59$sei^2)), 3))
}

# 5. PRET, its report's mortality outcome band.
pret <- Filter(function(e) !is.null(e$estimate), find("T09-1", 6, 59))
out$pret_6m_5y <- if (length(pret)) list(rate_ratio = num(pret[[1]]$estimate$value), lower = num(pret[[1]]$estimate$lower), upper = num(pret[[1]]$estimate$upper)) else NULL

write(toJSON(out, auto_unbox = TRUE, digits = NA, pretty = TRUE, null = "null", na = "null"), file.path(out_dir, "S.json"))
cat("summaries written\n")