# The review's registered analysis: does ketamine lower depressive symptoms when it is given only while
# patients are unconscious under general anesthesia, and by how much less than when they are awake?
#
#   Rscript code/analyze.R [--blind] [data dir] [results dir]
#
# Reads data/trials.json (each trial's design, coded from its Methods, and its integrity facts),
# data/extraction.json (the numbers each trial reported), data/parameters.json, and data/registry.json.
# Writes results/R*.json (declared results), results/effects.csv (every trial's effect in every window),
# and results/analyses.csv (every model fitted). --blind swaps each trial's arms at random (seed in
# parameters.json), so the code can be run and checked before the outcome data 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)
names(extraction) <- vapply(extraction, function(e) e$trial, "")
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)) NA_real_ else as.numeric(x)
set.seed(P$seed)

# A trial that randomized ketamine both while patients were unconscious (given after induction of general
# anesthesia and ended by the end of surgery) and while they were conscious, beside an inert arm, gives two
# comparisons: each kind of ketamine arm against the control arm, whose participants are split between the
# two (Cochrane Handbook 23.3.4), the first masked and the second awake. Its two kinds of ketamine arm are
# also compared with each other (`within`).
arm_unconscious <- function(a, t) {
  general <- (t$design$anesthesia %||% "") %in% c("general", "general_and_regional")
  general && (a$start %||% "") %in% c("at_induction", "after_induction") && (a$end %||% "") %in% c("before_end_of_surgery", "end_of_surgery") && !isTRUE(a$pump)
}
within <- list()
original <- trials
expanded <- list()
for (t in trials) {
  ket <- Filter(function(a) identical(a$role, "ketamine"), t$design$arms %||% list())
  ctl <- Filter(function(a) identical(a$role, "control"), t$design$arms %||% list())
  kinds <- vapply(ket, function(a) arm_unconscious(a, t), TRUE)
  if (isTRUE(t$design$eligible) && isTRUE(t$design$timing_comparison) && length(ctl) > 0 && any(kinds) && any(!kinds)) {
    others <- Filter(function(a) !(a$role %in% c("ketamine", "control")), t$design$arms %||% list())
    for (side in c("m", "a")) {
      u <- t
      u$trial <- paste0(t$trial, side)
      u$parent <- t$trial
      u$control_share <- 0.5
      u$design$arms <- c(if (side == "m") ket[kinds] else ket[!kinds], ctl, others)
      u$design$conscious_dosing <- if (side == "m") "none" else "some"
      u$design$late_dosing <- if (side == "m") any(vapply(ket[kinds], function(a) identical(a$end, "end_of_surgery"), TRUE)) else NULL
      expanded[[length(expanded) + 1]] <- u
    }
    within[[t$trial]] <- list(unconscious = vapply(ket[kinds], function(a) a$name, ""), conscious = vapply(ket[!kinds], function(a) a$name, ""))
  } else {
    expanded[[length(expanded) + 1]] <- t
  }
}
trials <- expanded

# ---------------------------------------------------------------------------------------------------
# Arms: which arms are ketamine and which the inert control. --blind swaps them in a random half.

arm_roles <- function(t) {
  roles <- list(ketamine = character(), control = character())
  for (a in t$design$arms) {
    if (identical(a$role, "ketamine")) roles$ketamine <- c(roles$ketamine, a$name)
    if (identical(a$role, "control")) roles$control <- c(roles$control, a$name)
  }
  roles
}
swap <- setNames(runif(length(trials)) < 0.5, vapply(trials, function(t) t$trial, ""))

# ---------------------------------------------------------------------------------------------------
# Converting one arm's reported values to a mean, SD, and n.

arm_n <- function(v, arms, name) {
  num(v$n) %|na|% num(arms[[name]]$analysed) %|na|% num(arms[[name]]$randomized)
}
`%|na|%` <- function(a, b) if (length(a) == 0 || is.na(a)) b else a

to_mean_sd <- function(form, v, n) {
  # Returns c(mean, sd, n, converted) or NULL. `converted` is 1 when the mean or SD was estimated
  # from medians, quartiles, or ranges (Luo et al. 2018 for the mean, Wan et al. 2014 for the SD).
  # A change from baseline (change_...) is converted the same way.
  form <- sub("^change_", "", form)
  if (form == "mean_sd") return(c(num(v$mean), num(v$sd), n, 0))
  if (form == "mean_se") return(c(num(v$mean), num(v$se) * sqrt(n), n, 0))
  if (form == "mean_range") {
    a <- num(v$min); b <- num(v$max)
    return(c(num(v$mean), (b - a) / (2 * qnorm((n - 0.375) / (n + 0.25))), n, 1))
  }
  if (form == "mean_ci") {
    level <- num(v$level) %|na|% 0.95
    tcrit <- qt(1 - (1 - level) / 2, df = max(n - 1, 1))
    return(c(num(v$mean), sqrt(n) * (num(v$upper) - num(v$lower)) / (2 * tcrit), n, 0))
  }
  if (form == "median_iqr") {
    m <- num(v$median); q1 <- num(v$q1); q3 <- num(v$q3)
    w <- 0.7 + 0.39 / n
    mean <- w * (q1 + q3) / 2 + (1 - w) * m
    sd <- (q3 - q1) / (2 * qnorm((0.75 * n - 0.125) / (n + 0.25)))
    return(c(mean, sd, n, 1))
  }
  if (form == "median_range") {
    m <- num(v$median); a <- num(v$min); b <- num(v$max)
    w <- 4 / (4 + n^0.75)
    mean <- w * (a + b) / 2 + (1 - w) * m
    sd <- (b - a) / (2 * qnorm((n - 0.375) / (n + 0.25)))
    return(c(mean, sd, n, 1))
  }
  NULL
}

# Combining several ketamine arms into one (Cochrane Handbook, section 23.3.4).
combine_arms <- function(rows) {
  if (length(rows) == 1) return(rows[[1]])
  acc <- rows[[1]]
  for (r in rows[-1]) {
    n1 <- acc[3]; n2 <- r[3]
    m <- (n1 * acc[1] + n2 * r[1]) / (n1 + n2)
    s <- sqrt(((n1 - 1) * acc[2]^2 + (n2 - 1) * r[2]^2 + n1 * n2 / (n1 + n2) * (acc[1] - r[1])^2) / (n1 + n2 - 1))
    acc <- c(m, s, n1 + n2, max(acc[4], r[4]))
  }
  acc
}

# ---------------------------------------------------------------------------------------------------
# One effect, ketamine against control, from one extracted entry: Hedges' g and its variance, with a
# note of how it was obtained. Negative g means less depression (or pain) with ketamine.

effect_of <- function(entry, t, e) {
  share <- t$control_share %||% 1
  roles <- arm_roles(t)
  if (blind && isTRUE(swap[t$trial])) roles <- list(ketamine = roles$control, control = roles$ketamine)
  form <- entry$form
  vals <- entry$values
  arms <- e$arms
  # Counts of improvement (responders, remitters) run the other way from counts of participants with
  # depression: their sign is turned so negative g still means fewer depressive symptoms.
  flip <- if (identical(entry$counts, "improvement")) -1 else 1
  if (form == "percent") {
    # A percentage of an arm with depression, as a count of that arm's n.
    vals <- lapply(vals, function(v) { n <- num(v$n); v$events <- num(v$percent) / 100 * n; v })
    for (a in names(vals)) if (is.na(num(vals[[a]]$n))) {
      n <- arm_n(list(), arms, a); vals[[a]]$n <- n; vals[[a]]$events <- num(vals[[a]]$percent) / 100 * n
    }
    form <- "count"
  }
  if (form == "count") {
    ev <- function(names) sum(vapply(names, function(a) num(vals[[a]]$events) %|na|% NA_real_, 0))
    nn <- function(names) sum(vapply(names, function(a) num(vals[[a]]$n) %|na|% arm_n(list(), arms, a), 0))
    a <- ev(roles$ketamine); n1 <- nn(roles$ketamine); c <- ev(roles$control) * share; n2 <- nn(roles$control) * share
    if (anyNA(c(a, n1, c, n2)) || n1 <= 0 || n2 <= 0) return(NULL)
    es <- escalc(measure = "OR2DL", ai = a, n1i = n1, ci = c, n2i = n2, add = 0.5, to = "only0")
    return(list(yi = flip * as.numeric(es$yi), vi = as.numeric(es$vi), n = n1 + n2, how = "count", converted = 1))
  }
  if (form %in% c("difference", "change_difference")) {
    # The estimate and its interval, given on the entry (or in its values), with the arm subtracted from
    # the other as its source says (contrast_sign); a difference whose contrast isn't known isn't used.
    get <- function(k) num(entry[[k]] %||% vals[[k]])
    level <- get("level") %|na|% 0.95
    z <- qnorm(1 - (1 - level) / 2)
    n1 <- sum(vapply(roles$ketamine, function(a) arm_n(vals[[a]] %||% list(), arms, a), 0))
    n2 <- sum(vapply(roles$control, function(a) arm_n(vals[[a]] %||% list(), arms, a), 0)) * share
    est <- get("estimate"); lo <- get("lower"); hi <- get("upper")
    scale <- entry$scale %||% vals$scale
    sign <- switch(entry$contrast_sign %||% "", ketamine_minus_control = 1, ketamine_versus_control = 1, control_minus_ketamine = -1, NA)
    if (anyNA(c(est, lo, hi, n1, n2, sign))) return(NULL)
    how <- if (form == "change_difference") "change" else "difference"
    if (identical(scale, "mean_difference")) {
      se <- (hi - lo) / (2 * z)
      sdp <- se / sqrt(1 / n1 + 1 / n2)
      es <- escalc(measure = "SMD", m1i = sign * est, sd1i = sdp, n1i = n1, m2i = 0, sd2i = sdp, n2i = n2)
      return(list(yi = as.numeric(es$yi), vi = as.numeric(es$vi), n = n1 + n2, how = how, converted = 1))
    }
    if (identical(scale, "odds_ratio")) {
      lor <- sign * log(est); se <- (log(hi) - log(lo)) / (2 * z)
      return(list(yi = flip * lor * sqrt(3) / pi, vi = se^2 * 3 / pi^2, n = n1 + n2, how = how, converted = 1))
    }
    return(NULL)
  }
  side <- function(names) {
    # The arms of one side that report this measure at this time, combined; an arm with no value here
    # (a trial reporting one of its ketamine arms on a day it doesn't report the other) is left out.
    rows <- list()
    for (a in names) {
      v <- vals[[a]]
      if (is.null(v)) next
      r <- to_mean_sd(form, v, arm_n(v, arms, a))
      if (is.null(r) || anyNA(r[1:3]) || r[2] <= 0 || r[3] <= 1) return(NULL)
      rows[[length(rows) + 1]] <- r
    }
    if (length(rows) == 0) return(NULL)
    combine_arms(rows)
  }
  k <- side(roles$ketamine); c <- side(roles$control)
  if (is.null(k) || is.null(c)) return(NULL)
  c[3] <- c[3] * share
  es <- escalc(measure = "SMD", m1i = k[1], sd1i = k[2], n1i = k[3], m2i = c[1], sd2i = c[2], n2i = c[3])
  how <- if (startsWith(form, "change_")) "change" else if (k[4] + c[4] > 0) "median" else "mean"
  list(yi = as.numeric(es$yi), vi = as.numeric(es$vi), n = k[3] + c[3], how = how, converted = as.numeric(how %in% c("median", "change")))
}

# ---------------------------------------------------------------------------------------------------
# Choosing each trial's outcome in a window: the measure the trial named primary, else the first of
# the hierarchy it reports there; the earliest day in the window; the best form for that measure and day.

hierarchy <- unlist(P$depression_hierarchy)
form_order <- unlist(P$form_order)
# The code the design coding gave a measure, by its name: depression codes (MADRS, HAMD, ...) or pain
# codes (rest, movement, opioid, other). A name the design coding doesn't list is OTHER.
code_of <- function(t, kind, m) {
  for (x in t$design[[paste0(kind, "_measures")]] %||% list()) if (identical(x$measure, m)) return(x$code %||% "OTHER")
  if (kind == "depression") "OTHER" else "other"
}

pick <- function(entries, t, window, kind) {
  lo <- window[[1]]; hi <- window[[2]]
  inwin <- Filter(function(x) !is.null(x$day) && x$day >= lo && x$day <= hi, entries)
  if (length(inwin) == 0) return(NULL)
  measures <- unique(vapply(inwin, function(x) x$measure, ""))
  order <- character()
  if (kind == "depression") {
    prim <- Filter(function(m) isTRUE(m$primary), t$design$depression_measures %||% list())
    if (length(prim)) order <- c(order, vapply(prim, function(m) m$measure, ""))
    ranked <- measures[order(match(vapply(measures, function(m) code_of(t, "depression", m), ""), hierarchy, nomatch = length(hierarchy) + 1))]
    order <- unique(c(intersect(order, measures), ranked))
  } else {
    codes <- vapply(measures, function(m) code_of(t, "pain", m), "")
    keep <- codes %in% unlist(P$pain_codes)
    order <- measures[keep][order(match(codes[keep], unlist(P$pain_codes)))]
  }
  for (m in order) {
    cand <- Filter(function(x) x$measure == m, inwin)
    days <- sort(unique(vapply(cand, function(x) as.numeric(x$day), 0)))
    for (d in days) {
      same <- Filter(function(x) as.numeric(x$day) == d, cand)
      # By form, and among counts, participants with depression before responders or remitters.
      same <- same[order(match(vapply(same, function(x) x$form, ""), form_order, nomatch = 99),
                         vapply(same, function(x) identical(x$counts, "improvement"), TRUE))]
      for (x in same) {
        ef <- effect_of(x, t, extraction[[t$parent %||% t$trial]])
        if (!is.null(ef)) return(c(ef, list(measure = m, day = d, form = x$form)))
      }
    }
  }
  NULL
}

# ---------------------------------------------------------------------------------------------------
# Integrity: facts about each trial that bear on whether its data can be trusted.

# GRIM (Brown and Heathers 2017): a mean of n whole-number scores, reported to d decimals, must be some
# whole number divided by n. Checked only where the n behind the mean is known exactly (given with it,
# or an arm with no loss between randomization and analysis) and where n < 10^d, so the test can fail.
grim_failures <- function(t, e) {
  fails <- 0L; checked <- 0L
  for (x in e$depression %||% list()) {
    if (!(x$form %in% c("mean_sd", "mean_se", "mean_ci"))) next
    if (!(code_of(t, "depression", x$measure) %in% unlist(P$integer_scales))) next
    for (a in names(x$values)) {
      v <- x$values[[a]]
      raw <- if (is.null(v$mean)) NULL else as.character(v$mean)
      n <- num(v$n)
      if (is.na(n)) {
        r <- num(e$arms[[a]]$randomized); an <- num(e$arms[[a]]$analysed)
        if (!is.na(r) && !is.na(an) && r == an) n <- r
      }
      if (is.na(n) || is.null(raw) || !grepl("^-?[0-9]+\\.[0-9]+$", raw)) next
      dec <- nchar(sub("^[^.]*\\.", "", raw))
      if (n >= 10^dec) next
      mean <- as.numeric(raw)
      s <- round(mean * n)
      ok <- any(abs(c(s - 1, s, s + 1) / n - mean) <= 0.5 * 10^(-dec) + 1e-9)
      checked <- checked + 1L
      if (!ok) fails <- fails + 1L
    }
  }
  c(checked = checked, fails = fails)
}

month_of <- function(s) {
  if (is.null(s) || is.na(s) || !nzchar(s)) return(NA_real_)
  p <- as.numeric(strsplit(substr(s, 1, 7), "-")[[1]])
  if (length(p) < 2 || anyNA(p)) return(NA_real_)
  p[1] * 12 + p[2]
}

integrity <- function(t, e) {
  reg <- t$design$registration
  regid <- reg$id %||% ""
  regdate <- reg$date %||% (t$registry_dates$first_submitted %||% NA)
  start <- t$design$enrollment$start %||% e$enrollment$start %||% NA
  retro <- if (!is.na(month_of(regdate)) && !is.na(month_of(start))) month_of(regdate) > month_of(start) else NA
  if (is.na(retro) && nzchar(regid) && !is.na(t$registration_year %||% NA) && !is.na(month_of(start))) {
    retro <- if (t$registration_year > floor((month_of(start) - 1) / 12)) TRUE else NA
  }
  g <- grim_failures(t, e)
  list(
    notice = length(t$notices %||% list()) > 0 && any(unlist(t$notices) %in% c("ExpressionOfConcernIn", "RetractionIn")),
    design_conflict = isTRUE(t$design_conflict),
    registration_shared = length(t$registration_shared_with %||% list()) > 0,
    registered = nzchar(regid),
    retrospective = isTRUE(retro),
    prospective_verified = identical(retro, FALSE),
    grim_checked = g[["checked"]],
    grim_fails = g[["fails"]]
  )
}

# ---------------------------------------------------------------------------------------------------
# Every trial's effects.

rows <- list()
for (t in trials) {
  if (!isTRUE(t$design$eligible)) next
  e <- extraction[[t$parent %||% t$trial]]
  if (is.null(e)) next
  integ <- integrity(t, e)
  base <- list(
    trial = t$trial, parent = t$parent %||% t$trial, source = t$source, setting = t$design$setting,
    group = switch(t$design$conscious_dosing %||% "unclear", none = "masked", some = "awake", "unclear"),
    late = isTRUE(t$design$late_dosing), selected = identical(t$design$depression_selection, "selected"),
    pump_only = isTRUE(t$pump_only), drug = t$drug, dose_ke = num(t$dose_ketamine_equivalent_mg_per_kg),
    notice = integ$notice, design_conflict = integ$design_conflict, registration_shared = integ$registration_shared,
    registered = integ$registered, retrospective = integ$retrospective,
    prospective_verified = integ$prospective_verified, grim_checked = integ$grim_checked, grim_fails = integ$grim_fails
  )
  for (w in names(P$windows)) {
    d <- pick(e$depression %||% list(), t, P$windows[[w]], "depression")
    if (!is.null(d)) rows[[length(rows) + 1]] <- c(base, list(outcome = "depression", window = w, yi = d$yi, vi = d$vi, n = d$n, how = d$how, converted = d$converted, measure = d$measure, day = d$day, form = d$form))
  }
  pn <- pick(e$pain %||% list(), t, P$pain_window, "pain")
  if (!is.null(pn)) rows[[length(rows) + 1]] <- c(base, list(outcome = "pain", window = "day1", yi = pn$yi, vi = pn$vi, n = pn$n, how = pn$how, converted = pn$converted, measure = pn$measure, day = pn$day, form = pn$form))
}
eff <- do.call(rbind, lapply(rows, function(r) as.data.frame(r, stringsAsFactors = FALSE)))
# A trial is trustworthy unless a notice, reports that disagree on whether allocation was random, a
# registration number another trial also gives, a missing or retrospective registration, two or more GRIM
# failures, or an early depression effect beyond the plausible bound counts against it.
early <- eff[eff$outcome == "depression" & eff$window == "early", c("trial", "yi")]
implausible <- early$trial[abs(early$yi) > P$implausible_g]
eff$implausible <- eff$trial %in% implausible
eff$trustworthy <- !eff$notice & !eff$design_conflict & !eff$registration_shared & eff$registered & !eff$retrospective &
  eff$grim_fails < P$grim_fail_threshold & !eff$implausible
write.csv(eff, file.path(out_dir, "effects.csv"), row.names = FALSE)

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

fits <- list()
fit_re <- function(d, label) {
  if (nrow(d) < 2) return(NULL)
  m <- rma(yi, vi, data = d, method = "REML", test = "knha")
  pr <- predict(m)
  r <- list(label = label, k = m$k, n = sum(d$n), estimate = as.numeric(m$b), lower = m$ci.lb, upper = m$ci.ub, p = m$pval,
            tau2 = m$tau2, i2 = m$I2, pi_lower = pr$pi.lb, pi_upper = pr$pi.ub)
  fits[[label]] <<- r
  invisible(r)
}
fit_mr <- function(d, label, mods) {
  if (nrow(d) < 4 || length(unique(d$unmasked)) < 2) return(NULL)
  m <- rma(yi, vi, mods = mods, data = d, method = "REML", test = "knha")
  j <- which(rownames(m$b) == "unmasked")
  r <- list(label = label, k = m$k, n = sum(d$n), estimate = as.numeric(m$b[j]), lower = m$ci.lb[j], upper = m$ci.ub[j], p = m$pval[j],
            tau2 = m$tau2, r2 = m$R2)
  fits[[label]] <<- r
  invisible(r)
}

dep <- function(window = "early") eff[eff$outcome == "depression" & eff$window == window, ]
surg <- function(d) {
  d <- d[d$setting == "surgery" & d$group %in% c("masked", "awake"), ]
  d$unmasked <- as.numeric(d$group == "awake")
  d$selected <- as.numeric(d$selected)
  d
}

# Primary 1: the masked effect, every masked trial, days 1 to 3.
p1 <- fit_re(dep()[dep()$group == "masked", ], "P1_masked_early")
# Primary 2: how much larger the effect is when patients are awake, among trials in surgery, adjusted
# for whether trials enrolled only depressed patients.
p2 <- fit_mr(surg(dep()), "P2_awake_minus_masked_surgery_early", ~ unmasked + selected)
holm <- p.adjust(c(p1$p %||% NA, p2$p %||% NA), method = "holm")

# The two groups separately, in surgery.
fit_re(surg(dep())[surg(dep())$unmasked == 0, ], "masked_surgery_early")
fit_re(surg(dep())[surg(dep())$unmasked == 1, ], "awake_surgery_early")
fit_re(dep()[dep()$group == "awake", ], "awake_all_early")

# Positive control: pain at about 24 hours.
pain <- eff[eff$outcome == "pain", ]
fit_re(pain[pain$group == "masked", ], "pain_masked")
fit_re(pain[pain$group == "awake", ], "pain_awake")
fit_mr(surg(pain), "pain_awake_minus_masked_surgery", ~ unmasked + selected)

# Sensitivity analyses of both primaries.
sens <- list(
  trustworthy = function(d) d[d$trustworthy, ],
  prospective = function(d) d[d$prospective_verified & !d$notice & !d$design_conflict & !d$registration_shared & d$grim_fails < P$grim_fail_threshold & !d$implausible, ],
  strict_masking = function(d) d[!(d$group == "masked" & d$late), ],
  no_pump_only = function(d) d[!d$pump_only, ],
  full_text = function(d) d[d$source %in% c("fulltext", "preprint", "library"), ],
  open_sources = function(d) d[d$source != "library", ],
  no_conversion = function(d) d[d$converted == 0, ],
  selected_only = function(d) d[d$selected, ],
  unselected_only = function(d) d[!d$selected, ],
  esketamine_only = function(d) d[d$drug == "esketamine", ]
)
for (s in names(sens)) {
  d <- sens[[s]](dep())
  fit_re(d[d$group == "masked", ], paste0("P1_", s))
  sd <- surg(d)
  if (s %in% c("selected_only", "unselected_only")) fit_mr(sd, paste0("P2_", s), ~ unmasked) else fit_mr(sd, paste0("P2_", s), ~ unmasked + selected)
}
for (w in setdiff(names(P$windows), "early")) {
  fit_re(dep(w)[dep(w)$group == "masked", ], paste0("P1_", w))
  fit_mr(surg(dep(w)), paste0("P2_", w), ~ unmasked + selected)
}

# Leave one trial out, for both primaries: the largest p.
loo <- function(d, f) {
  ps <- c()
  for (tr in unique(d$trial)) ps <- c(ps, f(d[d$trial != tr, ])$p %||% NA)
  max(ps, na.rm = TRUE)
}
p1_loo <- loo(dep()[dep()$group == "masked", ], function(d) {
  m <- rma(yi, vi, data = d, method = "REML", test = "knha"); list(p = m$pval)
})
p2_loo <- loo(surg(dep()), function(d) {
  if (length(unique(d$unmasked)) < 2) return(list(p = NA))
  m <- rma(yi, vi, mods = ~ unmasked + selected, data = d, method = "REML", test = "knha"); list(p = m$pval[2])
})

# Small-study effects among masked trials: Egger's test and the PET-PEESE estimate.
masked <- dep()[dep()$group == "masked", ]
egger <- if (nrow(masked) >= 10) regtest(rma(yi, vi, data = masked, method = "REML"), model = "rma", predictor = "sei") else NULL
pet <- if (nrow(masked) >= 10) rma(yi, vi, mods = ~ sqrt(vi), data = masked, method = "REML", test = "knha") else NULL
peese <- if (nrow(masked) >= 10) rma(yi, vi, mods = ~ vi, data = masked, method = "REML", test = "knha") else NULL
petpeese <- if (!is.null(pet)) {
  use <- if (pet$pval[1] < 0.05) peese else pet
  list(estimate = as.numeric(use$b[1]), lower = use$ci.lb[1], upper = use$ci.ub[1], model = if (pet$pval[1] < 0.05) "PEESE" else "PET")
} else NULL

# Dose and response among masked trials.
dose <- masked[!is.na(masked$dose_ke) & masked$dose_ke > 0, ]
dose_fit <- if (nrow(dose) >= 10) rma(yi, vi, mods = ~ log(dose_ke), data = dose, method = "REML", test = "knha") else NULL

# ---------------------------------------------------------------------------------------------------
# Manipulation check: psychotomimetic adverse events, masked against awake.

ae_rows <- list()
types <- unlist(P$psychotomimetic_types)
for (t in trials) {
  if (!isTRUE(t$design$eligible)) next
  e <- extraction[[t$parent %||% t$trial]]
  if (is.null(e) || length(e$adverse_events %||% list()) == 0) next
  roles <- arm_roles(t)
  if (blind && isTRUE(swap[t$trial])) roles <- list(ketamine = roles$control, control = roles$ketamine)
  # Participants with each event in each arm: the largest count over the windows reported (counts in
  # different windows can be the same participants, so they aren't added), at most the arm's size.
  tot <- setNames(rep(0, length(types)), types)
  cnt <- list()
  for (side in c("ketamine", "control")) for (a in roles[[side]]) for (ty in types) {
    best <- NA
    for (w in e$adverse_events) {
      v <- num(w$values[[a]][[ty]])
      if (!is.na(v)) best <- max(best, v, na.rm = TRUE)
    }
    if (!is.na(best)) {
      best <- min(best, arm_n(list(), e$arms, a), na.rm = TRUE)
      key <- paste(side, ty)
      cnt[[key]] <- (cnt[[key]] %||% 0) + best
      tot[[ty]] <- tot[[ty]] + best
    }
  }
  if (all(tot == 0)) next
  dom <- names(which.max(tot))
  nn <- function(names) sum(vapply(names, function(a) arm_n(list(), e$arms, a), 0))
  share <- t$control_share %||% 1
  ae_rows[[length(ae_rows) + 1]] <- data.frame(
    trial = t$trial, group = switch(t$design$conscious_dosing %||% "unclear", none = "masked", some = "awake", "unclear"),
    type = dom, ai = cnt[[paste("ketamine", dom)]] %||% 0, n1i = nn(roles$ketamine),
    ci = (cnt[[paste("control", dom)]] %||% 0) * share, n2i = nn(roles$control) * share, stringsAsFactors = FALSE)
}
ae <- if (length(ae_rows)) do.call(rbind, ae_rows) else data.frame()
mh <- function(d) {
  if (nrow(d) < 2) return(NULL)
  m <- rma.mh(ai = ai, n1i = n1i, ci = ci, n2i = n2i, data = d, measure = "OR")
  list(k = m$k, or = exp(as.numeric(m$b)), lower = exp(m$ci.lb), upper = exp(m$ci.ub), p = m$pval, logor = as.numeric(m$b), se = m$se)
}
ae_masked <- if (nrow(ae)) mh(ae[ae$group == "masked", ]) else NULL
ae_awake <- if (nrow(ae)) mh(ae[ae$group == "awake", ]) else NULL
ae_diff <- if (!is.null(ae_masked) && !is.null(ae_awake)) {
  z <- (ae_awake$logor - ae_masked$logor) / sqrt(ae_awake$se^2 + ae_masked$se^2)
  list(ratio = exp(ae_awake$logor - ae_masked$logor), p = 2 * pnorm(-abs(z)))
} else NULL

# ---------------------------------------------------------------------------------------------------
# Registrations on ClinicalTrials.gov that never reported: the file drawer.

reg <- registry$registrations %||% list()
eligible_regs <- Filter(function(r) isTRUE(r$eligible), reg)
due <- Filter(function(r) isTRUE(r$due), eligible_regs)
unreported <- Filter(function(r) !isTRUE(r$published) && !isTRUE(r$results_posted), due)

# ---------------------------------------------------------------------------------------------------
# Within a trial that randomized ketamine both ways: conscious against unconscious dosing, days 1 to 3
# (negative g: fewer depressive symptoms when the drug was given to patients who were awake).

within_results <- list()
for (tid in names(within)) {
  t <- Find(function(x) identical(x$trial, tid), original)
  e <- extraction[[tid]]
  if (is.null(t) || is.null(e)) next
  u <- t
  u$design$arms <- c(lapply(Filter(function(a) a$name %in% within[[tid]]$conscious, t$design$arms), function(a) { a$role <- "ketamine"; a }),
                     lapply(Filter(function(a) a$name %in% within[[tid]]$unconscious, t$design$arms), function(a) { a$role <- "control"; a }))
  w <- pick(e$depression %||% list(), u, P$windows$early, "depression")
  within_results[[tid]] <- if (is.null(w)) NULL else list(g = round(w$yi, 2), lower = round(w$yi - 1.96 * sqrt(w$vi), 2),
    upper = round(w$yi + 1.96 * sqrt(w$vi), 2), p = signif(2 * pnorm(-abs(w$yi / sqrt(w$vi))), 2), n = w$n, measure = w$measure, day = w$day)
}

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

r3 <- function(x) if (is.null(x) || length(x) == 0 || is.na(x)) NULL else round(x, 2)
r1 <- function(x) if (is.null(x) || is.na(x)) NULL else round(x, 1)
pval <- function(x) if (is.null(x) || is.na(x)) NULL else signif(x, 2)
fit_summary <- function(f) if (is.null(f)) NULL else list(k = f$k, n = f$n, g = r3(f$estimate), lower = r3(f$lower), upper = r3(f$upper), p = pval(f$p))
# Each p value also as the paper prints it ("p = 0.011", or "p < 0.001" below 0.001), beside it as <key>_text.
p_text <- function(p) if (p < 0.001) "p < 0.001" else paste("p =", format(signif(p, 2), scientific = FALSE))
with_p_text <- function(x) {
  if (!is.list(x)) return(x)
  if (is.null(names(x))) return(lapply(x, with_p_text))
  out <- list()
  for (k in names(x)) {
    out[[k]] <- with_p_text(x[[k]])
    v <- x[[k]]
    if ((k == "p" || startsWith(k, "p_") || endsWith(k, "_p")) && is.numeric(v) && length(v) == 1 && !is.na(v)) out[[paste0(k, "_text")]] <- p_text(v)
  }
  out
}
write_result <- function(name, x) write(toJSON(with_p_text(x), auto_unbox = TRUE, null = "null", digits = NA, pretty = TRUE), file.path(out_dir, paste0(name, ".json")))

# Robust: significant after Holm's correction, among trustworthy trials, with strict masking, and with
# every single trial left out.
robust <- function(prefix, holm_p, loo_p) {
  needed <- c("trustworthy", "strict_masking")
  ok <- !is.na(holm_p) && holm_p < 0.05 && !is.na(loo_p) && loo_p < 0.05
  for (s in needed) ok <- ok && !is.null(fits[[paste0(prefix, s)]]) && fits[[paste0(prefix, s)]]$p < 0.05
  ok
}

write_result("R1", c(fit_summary(p1), list(p_holm = pval(holm[1]), tau2 = r3(p1$tau2), i2_pct = r1(p1$i2),
  pi_lower = r3(p1$pi_lower), pi_upper = r3(p1$pi_upper), loo_max_p = pval(p1_loo), robust = robust("P1_", holm[1], p1_loo))))
write_result("R2", c(fit_summary(p2), list(p_holm = pval(holm[2]), loo_max_p = pval(p2_loo), robust = robust("P2_", holm[2], p2_loo),
  masked = fit_summary(fits$masked_surgery_early), awake = fit_summary(fits$awake_surgery_early))))
write_result("R3", list(masked = fit_summary(fits$pain_masked), awake = fit_summary(fits$pain_awake),
  difference = fit_summary(fits$pain_awake_minus_masked_surgery)))
write_result("R4", list(masked = if (is.null(ae_masked)) NULL else list(k = ae_masked$k, or = r3(ae_masked$or), lower = r3(ae_masked$lower), upper = r3(ae_masked$upper), p = pval(ae_masked$p)),
  awake = if (is.null(ae_awake)) NULL else list(k = ae_awake$k, or = r3(ae_awake$or), lower = r3(ae_awake$lower), upper = r3(ae_awake$upper), p = pval(ae_awake$p)),
  ratio_awake_to_masked = if (is.null(ae_diff)) NULL else r3(ae_diff$ratio), p_difference = if (is.null(ae_diff)) NULL else pval(ae_diff$p)))
write_result("R5", list(
  P1 = lapply(setNames(names(sens), names(sens)), function(s) fit_summary(fits[[paste0("P1_", s)]])),
  P2 = lapply(setNames(names(sens), names(sens)), function(s) fit_summary(fits[[paste0("P2_", s)]]))))
write_result("R6", list(
  P1 = lapply(setNames(setdiff(names(P$windows), "early"), setdiff(names(P$windows), "early")), function(w) fit_summary(fits[[paste0("P1_", w)]])),
  P2 = lapply(setNames(setdiff(names(P$windows), "early"), setdiff(names(P$windows), "early")), function(w) fit_summary(fits[[paste0("P2_", w)]]))))
write_result("R7", list(
  egger_p = if (is.null(egger)) NULL else pval(egger$pval),
  pet_peese = if (is.null(petpeese)) NULL else list(model = petpeese$model, g = r3(petpeese$estimate), lower = r3(petpeese$lower), upper = r3(petpeese$upper)),
  dose_slope_per_log_mg_kg = if (is.null(dose_fit)) NULL else list(estimate = r3(dose_fit$b[2]), lower = r3(dose_fit$ci.lb[2]), upper = r3(dose_fit$ci.ub[2]), p = pval(dose_fit$pval[2])),
  registrations_eligible = length(eligible_regs), registrations_due = length(due), registrations_unreported = length(unreported)))
write_result("R8", list(within_trial = within_results))
write_result("S", list(
  trials_eligible = sum(vapply(original, function(t) isTRUE(t$design$eligible), TRUE)),
  trials_with_early_effect = nrow(dep()), masked_early = sum(dep()$group == "masked"), awake_early = sum(dep()$group == "awake"),
  participants_early = sum(dep()$n), trustworthy_early = sum(dep()$trustworthy), blind = blind))

an <- do.call(rbind, lapply(fits, function(f) data.frame(label = f$label, k = f$k, n = f$n, estimate = f$estimate, lower = f$lower, upper = f$upper, p = f$p, stringsAsFactors = FALSE)))
write.csv(an, file.path(out_dir, "analyses.csv"), row.names = FALSE)
cat("done:", nrow(eff), "effects,", length(fits), "models", if (blind) "(blind)" else "", "\n")