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

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

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)
  ))
  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)
    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)
    )
  }
  out
}