# Aggregate results over the pre-registered associations. Proportions carry exact
# (Clopper-Pearson) 95% intervals and medians the distribution-free interval from order
# statistics. Values are rounded to the precision their uncertainty supports, so a re-run on
# another platform declares the same numbers.

share <- function(k, n) {
  if (n == 0) return(list(k = k, n = n, share = NA, ci = c(NA, NA)))
  list(k = k, n = n, share = round(k / n, 3), ci = round(as.numeric(stats::binom.test(k, n)$conf.int), 3))
}

median_with_ci <- function(x) {
  x <- x[is.finite(x)]
  if (length(x) == 0) return(list(n = 0, median = NA, ci = c(NA, NA), quartiles = c(NA, NA)))
  list(n = length(x), median = round(stats::median(x), 3), ci = round(median_ci(x), 3),
       quartiles = round(as.numeric(stats::quantile(x, c(0.25, 0.75), type = 7)), 3))
}

# The departures bearing on a headline estimate, keyed as in data/departure_coding.csv: the
# association's id and the departure's place in its file's list.
headline_departures <- function(results) {
  unlist(unname(lapply(results, function(r) {
    if (length(r$departures) == 0) return(list())
    keys <- paste0(r$id, ".", seq_along(r$departures))
    setNames(r$departures, keys)[vapply(r$departures, function(d) isTRUE(d$affects_headline), TRUE)]
  })), recursive = FALSE)
}

# Two independent codings of each headline departure's evidence (data/departure_coding.csv),
# made after registration: how often they agree, and Cohen's kappa. The association files hold
# the label each departure carries, settled where the two differ.
coding_agreement <- function(departures, path = "data/departure_coding.csv") {
  coding <- read.csv(path, colClasses = "character")
  if (!setequal(coding$key, names(departures)) || anyDuplicated(coding$key)) stop(path, " must code each headline departure once")
  if (!all(c(coding$coder_1, coding$coder_2) %in% DEPARTURE_EVIDENCE)) stop(path, ": unknown evidence label")
  m <- nrow(coding)
  observed <- mean(coding$coder_1 == coding$coder_2)
  chance <- sum(vapply(DEPARTURE_EVIDENCE, function(e) mean(coding$coder_1 == e) * mean(coding$coder_2 == e), 0))
  final <- vapply(departures[coding$key], function(d) d$evidence, "")
  list(agree = share(sum(coding$coder_1 == coding$coder_2), m), kappa = round((observed - chance) / (1 - chance), 2),
       settled = sum(coding$coder_1 != coding$coder_2),
       final_differs_from = list(coder_1 = sum(final != coding$coder_1), coder_2 = sum(final != coding$coder_2)))
}

# Added in answer to the second version's reviews: 20 headline departures drawn at random (the seed
# below, from the sorted keys), each checked against its paper's full text by an agent that recorded
# none of them (data/departure_validation.csv): confirmed, partly (the record misstates the paper but
# the departure stands), not confirmed, or not checkable. Its checker_verdict column keeps the
# checking agent's own verdict where the lead agent counted one otherwise under these definitions.
VALIDATION_VERDICTS <- c("confirmed", "partly", "not_confirmed", "could_not_check")

departure_validation <- function(departures, path = "data/departure_validation.csv", n = 20) {
  set.seed(20261019, kind = "Mersenne-Twister", normal.kind = "Inversion", sample.kind = "Rejection")
  drawn <- sample(sort(names(departures)), n)
  checked <- read.csv(path, colClasses = "character")
  if (!setequal(checked$key, drawn) || anyDuplicated(checked$key)) stop(path, " must check each of the ", n, " departures drawn once")
  if (!all(checked$verdict %in% VALIDATION_VERDICTS)) stop(path, ": unknown verdict")
  c(list(drawn = n), setNames(lapply(VALIDATION_VERDICTS, function(v) share(sum(checked$verdict == v), n)), VALIDATION_VERDICTS),
    # Records misstating anything about their paper, whether or not the departure stands.
    list(misstating = share(sum(checked$verdict %in% c("partly", "not_confirmed")), n)))
}

# The departures the coding check found bearing on the headline estimate: how many associations
# have one, by kind, and how many have one the re-implementation doesn't follow (it corrects it,
# or the paper's numbers rule out the stated computation without revealing its own); then, added
# after registration, by evidence: the departures of each kind of evidence, and the associations
# with one the paper itself shows, with one only recomputation identifies, and with only
# unresolved ones.
departure_counts <- function(results) {
  n <- length(results)
  count <- function(f) sum(vapply(results, function(r) isTRUE(f(r)), TRUE))
  affecting <- function(r, keep = function(d) TRUE) {
    any(vapply(r$departures, function(d) isTRUE(d$affects_headline) && keep(d), TRUE))
  }
  shown <- function(e) function(d) identical(d$evidence, e)
  departures <- headline_departures(results)
  m <- length(departures)
  evidence <- vapply(departures, function(d) d$evidence, "")
  list(
    affecting_headline = share(count(function(r) affecting(r)), n),
    by_kind = setNames(lapply(DEPARTURE_KINDS, function(k) share(count(function(r) affecting(r, function(d) d$kind == k)), n)), DEPARTURE_KINDS),
    not_followed = share(count(function(r) affecting(r, function(d) isFALSE(d$followed))), n),
    headline_departures = m,
    by_evidence = setNames(lapply(DEPARTURE_EVIDENCE, function(e) share(sum(evidence == e), m)), DEPARTURE_EVIDENCE),
    shown_by_paper = share(count(function(r) affecting(r, shown("paper"))), n),
    shown_by_paper_computation = share(count(function(r) affecting(r, function(d) identical(d$evidence, "paper") && d$kind != "reporting")), n),
    identified_by_data = share(count(function(r) !affecting(r, shown("paper")) && affecting(r, shown("data"))), n),
    only_unresolved = share(count(function(r) affecting(r) && !affecting(r, function(d) d$evidence %in% c("paper", "data"))), n),
    coding = coding_agreement(departures),
    validation = departure_validation(departures)
  )
}

# Added in answer to the second version's reviews: each analysis as its paper describes it
# (described_version in pipeline.R). The associations it changes; the departures it puts back and
# those whose text's version can't be built; how many described estimates fall outside the published
# interval (the registered reproduction criterion failed), over all associations and over those it
# changes; how many of those it changes lack a significant estimate with the published sign, as
# described and as followed; and the ratio of the described estimate to the published one.
described_counts <- function(results) {
  count <- function(rs, f) sum(vapply(rs, function(r) isTRUE(f(r)), TRUE))
  n <- length(results)
  changed <- Filter(function(r) !isTRUE(r$described$same_as_paper_version), results)
  m <- length(changed)
  list(
    changed = m,
    applied = sum(vapply(results, function(r) length(r$described$applies), 0L)),
    not_applied = sum(vapply(results, function(r) length(r$described$not_applied), 0L)),
    not_fitted = count(results, function(r) !r$described$converged),
    outside_published_ci = share(count(results, function(r) !r$described_in_published_ci), n),
    outside_published_ci_changed = share(count(changed, function(r) !r$described_in_published_ci), m),
    not_significant_changed = share(count(changed, function(r) !r$described_significant), m),
    not_significant_changed_followed = share(count(changed, function(r) !r$paper_version_significant), m),
    ratio = median_with_ci(vapply(results, function(r) r$described_ratio, 0)),
    ratio_changed = median_with_ci(vapply(changed, function(r) r$described_ratio, 0)),
    to_followed_changed = median_with_ci(vapply(changed, function(r) if (isTRUE(r$described$converged)) r$described$b / r$original$b else NA, 0)),
    # Those whose described estimate is within 1% of the followed one on the analysis scale.
    within_one_percent = count(changed, function(r) isTRUE(abs(r$described$b / r$original$b - 1) < 0.01))
  )
}

summarize_results <- function(results, replication) {
  count <- function(rs, f) sum(vapply(rs, function(r) isTRUE(f(r)), TRUE))
  n <- length(results)
  out <- list(reproduction = list(
    associations = n,
    reproduced = share(count(results, function(r) r$reproduced), n),
    same_sign_p05 = share(count(results, function(r) r$same_sign_original && r$original$p < 0.05), n),
    ratio_original = median_with_ci(vapply(results, function(r) r$ratio_original, 0)),
    departures = departure_counts(results),
    described = described_counts(results)
  ))
  if (replication) {
    informative <- Filter(function(r) isTRUE(r$informative), results)
    informative_bonferroni <- Filter(function(r) isTRUE(r$informative_bonferroni), results)
    reproduced <- Filter(function(r) isTRUE(r$reproduced), results)
    unweighted <- Filter(function(r) isFALSE(r$replication$weighted), results)
    informative_half <- Filter(function(r) isTRUE(r$power_half >= 0.8), results)
    out$replication <- list(
      associations = n,
      replicated = share(count(results, function(r) r$replicated), n),
      informative = length(informative),
      replicated_informative = share(count(informative, function(r) r$replicated), length(informative)),
      same_sign_p05 = share(count(results, function(r) r$same_sign && r$replication$p < 0.05), n),
      in_published_ci = share(count(results, function(r) r$in_published_ci), n),
      reversed = share(count(results, function(r) !r$same_sign && r$q < 0.05), n),
      ratio = median_with_ci(vapply(results, function(r) r$ratio, 0)),
      ratio_informative = median_with_ci(vapply(informative, function(r) r$ratio, 0)),
      ratio_own = median_with_ci(vapply(results, function(r) r$ratio_own, 0)),
      heterogeneous = share(count(results, function(r) r$heterogeneity_q < 0.05), n),
      differs_from_published = share(count(results, function(r) r$differs_from_published), n),
      differs_by_direction = lapply(c(smaller = "smaller", larger = "larger", opposite_sign = "opposite sign"),
                                    function(d) share(count(results, function(r) identical(r$difference, d)), n)),
      replicated_reproduced = share(count(reproduced, function(r) r$replicated), length(reproduced)),
      replicated_paper_weight = share(count(results, function(r) r$replicated_paper_weight), n),
      # Sensitivities added after registration, in answer to review.
      informative_bonferroni = length(informative_bonferroni),
      replicated_informative_bonferroni = share(count(informative_bonferroni, function(r) r$replicated), length(informative_bonferroni)),
      differs_from_published_t = share(count(results, function(r) r$differs_from_published_t), n),
      heterogeneous_t = share(count(results, function(r) r$heterogeneity_q_t < 0.05), n),
      # Added in answer to the second version's reviews: the analyses re-run unweighted, as their
      # papers' numbers show they were computed, and the replication with the survey design for
      # them, weighted (registered) or with their unweighted estimates' design-based variance;
      # tests informative with each, and at half the published effect.
      unweighted = share(length(unweighted), n),
      replicated_unweighted = share(count(unweighted, function(r) r$replicated), length(unweighted)),
      design = design_counts(results, unweighted),
      informative_half = length(informative_half),
      replicated_informative_half = share(count(informative_half, function(r) r$replicated), length(informative_half)),
      expected = expected_replications(results),
      expected_design = expected_replications(results, "design"),
      expected_weighted = expected_replications(results, "weighted"),
      # The effect ratio with the registered weighted refits in place of the unweighted fits, and
      # over the analyses re-run unweighted and the others apart.
      ratio_weighted = median_with_ci(vapply(results, function(r) sensitivity_fit(r, "weighted")$b / r$published_b[1], 0)),
      ratio_unweighted = median_with_ci(vapply(unweighted, function(r) r$ratio, 0)),
      ratio_others = median_with_ci(vapply(Filter(function(r) !isFALSE(r$replication$weighted), results), function(r) r$ratio, 0))
    )
  }
  out
}

# Added in answer to the second version's reviews: the replication tests simulated as if every
# association were true, its effect `scale` times the published one, each test drawn from its own
# 2021-2023 standard error and degrees of freedom and decided as registered (same sign, Benjamini-
# Hochberg over all of them), the tests independent: the count of replications in each of `draws`.
simulate_replications <- function(results, scale, sensitivity = NULL, draws = 10000) {
  set.seed(20261019, kind = "Mersenne-Twister", normal.kind = "Inversion", sample.kind = "Rejection")
  fits <- lapply(results, function(r) if (is.null(sensitivity)) r$replication else sensitivity_fit(r, sensitivity))
  b <- vapply(results, function(r) r$published_b[1], 0)
  se <- vapply(fits, function(fit) fit$se, 0)
  df <- vapply(fits, function(fit) fit$df, 0)
  finite <- is.finite(df)
  chisq_df <- ifelse(finite, df, 1)
  vapply(seq_len(draws), function(i) {
    t <- (stats::rnorm(length(b)) + scale * b / se) / ifelse(finite, sqrt(stats::rchisq(length(b), chisq_df) / chisq_df), 1)
    q <- stats::p.adjust(2 * stats::pt(-abs(t), df), method = "BH")
    sum(sign(t) == sign(b) & q < 0.05)
  }, 0)
}

# How many replications that predicts at the published effect and at half of it (mean and the central
# range holding 95% of draws), the chance of as few as were observed, and the effect size, as a
# fraction of the published one, at which as many as were observed are expected, with the fractions
# (in steps of 0.05) whose central range holds the observed count: for the registered tests, or, with
# `sensitivity`, for a sensitivity analysis's fits and its count of replications.
expected_replications <- function(results, sensitivity = NULL) {
  replicated <- if (is.null(sensitivity)) "replicated" else paste0("replicated_", sensitivity)
  observed <- sum(vapply(results, function(r) isTRUE(r[[replicated]]), TRUE))
  at <- function(scale) {
    counts <- simulate_replications(results, scale, sensitivity)
    list(mean = round(mean(counts), 1), range = as.numeric(stats::quantile(counts, c(0.025, 0.975), type = 1)),
         at_most_observed = round(mean(counts <= observed), 3), at_least_observed = round(mean(counts >= observed), 3))
  }
  grid <- seq(0.3, 1.5, by = 0.05)
  sims <- lapply(grid, function(scale) simulate_replications(results, scale, sensitivity))
  means <- vapply(sims, mean, 0)
  consistent <- grid[vapply(sims, function(counts) mean(counts <= observed) >= 0.025 && mean(counts >= observed) >= 0.025, TRUE)]
  list(observed = observed, published = at(1), half = at(0.5),
       implied_scale = round(stats::approx(means, grid, xout = observed)$y, 2), implied_scale_range = range(consistent))
}

# The replication count with the survey design for the analyses re-run unweighted (`unweighted`),
# weighted or with their unweighted estimates' design-based variance: replicated over all and among
# the unweighted ones, and the tests informative with it and replicated among them.
design_counts <- function(results, unweighted) {
  count <- function(rs, f) sum(vapply(rs, function(r) isTRUE(f(r)), TRUE))
  n <- length(results)
  lapply(c(weighted = "weighted", design = "design"), function(s) {
    replicated <- function(r) r[[paste0("replicated_", s)]]
    informative <- Filter(function(r) isTRUE(r[[paste0("power_", s)]] >= 0.8), results)
    # The refit's standard error over the registered unweighted fit's, for the analyses re-run unweighted.
    se_ratio <- vapply(unweighted, function(r) sensitivity_fit(r, s)$se / r$replication$se, 0)
    list(
      replicated = share(count(results, replicated), n),
      replicated_unweighted = share(count(unweighted, replicated), length(unweighted)),
      informative = length(informative),
      informative_also_registered = count(informative, function(r) r$informative),
      replicated_informative = share(count(informative, replicated), length(informative)),
      se_ratio = list(median = round(stats::median(se_ratio), 2), range = round(range(se_ratio), 2), smaller = sum(se_ratio < 1))
    )
  })
}