# Writes the per-association results the paper cites, from results/associations.json (code/main.R)
# and the sampling frame's labels (data/frame.csv): results/associations.csv, one row per
# association, results/R2.json, the same values by association for the paper's placeholders, and
# results/R3.json, the size of the sampling frame. It computes nothing new; it rounds what main.R
# computed to the precision the paper reports.

results <- jsonlite::fromJSON("results/associations.json", simplifyVector = FALSE)
frame <- read.csv("data/frame.csv", colClasses = "character", encoding = "UTF-8")

# Four significant digits for estimates and bounds (an odds ratio per milligram such as 1.002
# keeps its information), two for p-values, and two decimals for ratios and power.
est <- function(x) if (is.null(x) || is.na(x)) NA else signif(x, 4)
pv <- function(x) if (is.null(x) || is.na(x)) NA else signif(x, 2)
dec <- function(x) if (is.null(x) || is.na(x)) NA else round(x, 2)
bound <- function(fit, i) if (is.null(fit$ci) || length(fit$ci) < i) NA else est(fit$ci[[i]])
headline <- function(r) Filter(function(d) isTRUE(d$affects_headline), r$departures)
# A 2021-2023 refit of a sensitivity analysis, where the association has one.
refit <- function(fit) if (is.null(fit)) NULL else list(estimate = est(fit$estimate), low = bound(fit, 1), high = bound(fit, 2), n = fit$n, p = pv(fit$p))

entry <- function(r) {
  label <- frame[frame$row == as.character(r$row), ]
  departures <- headline(r)
  list(
    row = r$row,
    predictor = label$predictor, condition = label$condition,
    measure = r$measure,
    published = list(estimate = est(r$published$estimate), low = est(r$published$low), high = est(r$published$high),
                     z = dec(r$published_z), interval_ratio = dec(r$interval_ratio_original)),
    original = list(estimate = est(r$original$estimate), low = bound(r$original, 1), high = bound(r$original, 2), n = r$original$n),
    harmonized = list(estimate = est(r$harmonized$estimate), low = bound(r$harmonized, 1), high = bound(r$harmonized, 2), n = r$harmonized$n),
    replication = list(estimate = est(r$replication$estimate), low = bound(r$replication, 1), high = bound(r$replication, 2),
                       n = r$replication$n, p = pv(r$replication$p), q = pv(r$q), power = dec(r$power)),
    replicated = isTRUE(r$replicated), informative = isTRUE(r$informative),
    difference = r$difference, difference_q = pv(r$published_difference_q), difference_q_t = pv(r$published_difference_q_t),
    power_bonferroni = dec(r$power_bonferroni), power_half = dec(r$power_half),
    # The sensitivity analyses' refits in 2021-2023 (code/lib/pipeline.R), added in answer to the
    # second version's reviews: the survey design for an analysis re-run unweighted, weighted or
    # with its unweighted estimate's design-based variance, and the phlebotomy weight.
    weighted = isTRUE(r$replication$weighted),
    replication_weighted = refit(r$replication_weighted), replicated_weighted = isTRUE(r$replicated_weighted),
    replication_design = refit(r$replication_design), replicated_design = isTRUE(r$replicated_design),
    replication_phlebotomy = refit(r$replication_phlebotomy), difference_q_phlebotomy = pv(r$published_difference_q_phlebotomy),
    difference_q_t_phlebotomy = pv(r$published_difference_q_t_phlebotomy),
    difference_q_weighted = pv(r$published_difference_q_weighted), difference_q_t_weighted = pv(r$published_difference_q_t_weighted),
    # The 2021-2023 estimate against our harmonized one on the paper's own cycles (C8's tests).
    heterogeneity_p = pv(r$heterogeneity_p), heterogeneity_q = pv(r$heterogeneity_q), heterogeneity_q_t = pv(r$heterogeneity_q_t),
    # The analysis as the paper describes it, on its own cycles (pipeline.R's described_version).
    described = c(refit(r$described), list(applies = length(r$described$applies), not_applied = length(r$described$not_applied),
                                            same_as_paper_version = isTRUE(r$described$same_as_paper_version))),
    described_in_published_ci = isTRUE(r$described_in_published_ci), described_ratio = dec(r$described_ratio),
    ratio = list(overall = dec(r$ratio), original = dec(r$ratio_original), harmonized = dec(r$ratio_harmonized),
                 between_cycles = dec(r$ratio_own)),
    departures_headline = length(departures),
    departure_kinds = paste(sort(unique(vapply(departures, function(d) d$kind, ""))), collapse = ", "),
    departures_not_followed = sum(vapply(departures, function(d) isFALSE(d$followed), TRUE)),
    departure_evidence = paste(vapply(departures, function(d) d$evidence, ""), collapse = ", "),
    # The variants tried on the paper's own cycles, in the order its file lists them.
    variants = lapply(r$variants, function(v) list(label = v$label, estimate = est(v$estimate), low = bound(v, 1), high = bound(v, 2), n = v$n))
  )
}

entries <- lapply(results, entry)
names(entries) <- vapply(results, function(r) r$id, "")
writeLines(jsonlite::toJSON(entries, auto_unbox = TRUE, digits = NA, pretty = TRUE, null = "null", na = "null"), "results/R2.json")

flat <- do.call(rbind, lapply(names(entries), function(id) {
  e <- entries[[id]]
  data.frame(
    id = id, row = e$row, doi = results[[id]]$doi, predictor = e$predictor, condition = e$condition, measure = e$measure,
    contrast = results[[id]]$contrast,
    published = e$published$estimate, published_low = e$published$low, published_high = e$published$high,
    published_z = e$published$z, interval_ratio = e$published$interval_ratio,
    original = e$original$estimate, original_low = e$original$low, original_high = e$original$high, original_n = e$original$n,
    harmonized = e$harmonized$estimate, harmonized_low = e$harmonized$low, harmonized_high = e$harmonized$high,
    replication = e$replication$estimate, replication_low = e$replication$low, replication_high = e$replication$high,
    replication_n = e$replication$n, replication_p = e$replication$p, replication_q = e$replication$q, power = e$replication$power,
    replicated = e$replicated, informative = e$informative, difference = e$difference, difference_q = e$difference_q,
    difference_q_t = e$difference_q_t, power_bonferroni = e$power_bonferroni,
    ratio = e$ratio$overall, ratio_original = e$ratio$original, ratio_harmonized = e$ratio$harmonized, ratio_between_cycles = e$ratio$between_cycles,
    departures_headline = e$departures_headline, departure_kinds = e$departure_kinds, departures_not_followed = e$departures_not_followed,
    departure_evidence = e$departure_evidence,
    replication_weights = paste(unlist(results[[id]]$replication$weights), collapse = " "),
    dropped = paste(unlist(results[[id]]$replication$dropped), collapse = " "),
    error = if (is.null(results[[id]]$replication$error)) "" else results[[id]]$replication$error,
    weighted = e$weighted, power_half = e$power_half,
    replicated_weighted = e$replicated_weighted, replicated_design = e$replicated_design,
    described = e$described$estimate, described_low = e$described$low, described_high = e$described$high,
    described_departures = e$described$applies, described_in_published_ci = e$described_in_published_ci,
    stringsAsFactors = FALSE
  )
}))
write.csv(flat, "results/associations.csv", row.names = FALSE, fileEncoding = "UTF-8")

writeLines(jsonlite::toJSON(list(frame_papers = nrow(frame)), auto_unbox = TRUE, pretty = TRUE), "results/R3.json")