Multiple Comparisons Correction — What Survives?

suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))

compute_shared <- function(df, params, col_map = list()) {
  # === SHARED EXPORTS ===
  #   initial_rows/final_rows/rows_removed  $ row accounting (rows_removed = invalid p dropped)
  #   p_name / label_name    $ humanized user column names
  #   n_tests                $ integer — valid tests entering the correction
  #   n_na / n_out_of_range  $ integer — dropped: missing vs outside [0,1]
  #   alpha                  $ 0.05
  #   results_all            $ data.frame(test, raw_p, holm_p, bh_p, sig_raw, sig_holm, sig_bh) sorted by raw_p
  #   adjusted_results_df    $ results_all capped at 50 rows
  #   results_truncated      $ logical — TRUE when results_all > 50 rows
  #   n_sig_raw/n_sig_holm/n_sig_bh $ integer — significant counts per method
  #   significance_summary_df $ data.frame(method, significant_count)
  #   p_value_distribution_df $ data.frame(p_range, count) — 10 bins of width 0.1
  #   dist_diagnosis          $ character — computed narrative for the histogram shape
  #   metrics / json_output
  # === /SHARED EXPORTS ===

initial_rows <- nrow(df)
  p_name <- humanize_semantic("p_value", col_map)
  label_name <- humanize_semantic("test_label", col_map)
  if (!("p_value" %in% names(df))) {
    stop(sprintf("column_mapping must map p_value to your p-value column(&#x27;%s' not found).", p_name))
  }
  if (!("test_label" %in% names(df))) {
    stop(sprintf("column_mapping must map test_label to your test-name column(&#x27;%s' not found).", label_name))
  }

lab <- as.character(df$test_label)
  blank <- is.na(lab) | trimws(lab) == ""
  if (any(blank)) lab[blank] <- paste0("Test ", which(blank))

v <- df$p_value
  if (!is.numeric(v)) {
    conv <- suppressWarnings(as.numeric(as.character(v)))
    n_orig <- sum(!is.na(v) & as.character(v) != "")
    if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) {
      v <- conv
    } else {
      stop(sprintf("The mapped p-value column(&#x27;%s') is not numeric — map a column of raw p-values between 0 and 1.", p_name))
    }
  }
  n_na <- sum(is.na(v))
  n_out_of_range <- sum(!is.na(v) & (v < 0 | v > 1))
  keep <- !is.na(v) & v >= 0 & v <= 1
  p <- as.numeric(v[keep])
  labels <- lab[keep]
  final_rows <- length(p)
  rows_removed <- initial_rows - final_rows
  if (final_rows < 5) {
    stop(sprintf(
      "Only %d valid p-values remained in &#x27;%s' after dropping %d invalid (missing or outside 0-1) — need at least 5 tests to correct.",
      final_rows, p_name, rows_removed))
  }
  n_tests <- final_rows
  alpha <- 0.05

holm <- p.adjust(p, method = "holm")
  bh   <- p.adjust(p, method = "BH")
  sig_raw  <- p    < alpha
  sig_holm <- holm < alpha
  sig_bh   <- bh   < alpha
  n_sig_raw  <- sum(sig_raw)
  n_sig_holm <- sum(sig_holm)
  n_sig_bh   <- sum(sig_bh)

  results_all <- data.frame(
    test     = labels,
    raw_p    = signif(p, 4),
    holm_p   = signif(holm, 4),
    bh_p     = signif(bh, 4),
    sig_raw  = ifelse(sig_raw,  "Yes", "No"),
    sig_holm = ifelse(sig_holm, "Yes", "No"),
    sig_bh   = ifelse(sig_bh,   "Yes", "No"),
    stringsAsFactors = FALSE
  )
  results_all <- results_all[order(results_all$raw_p), , drop = FALSE]
  rownames(results_all) <- NULL
  results_truncated <- nrow(results_all) > 50
  adjusted_results_df <- head(results_all, 50)

  significance_summary_df <- data.frame(
    method = c("Uncorrected", "Holm-Bonferroni", "Benjamini-Hochberg"),
    significant_count = c(n_sig_raw, n_sig_holm, n_sig_bh),
    stringsAsFactors = FALSE
  )

breaks <- seq(0, 1, by = 0.1)
  bin_idx <- pmin(pmax(findInterval(p, breaks, rightmost.closed = TRUE), 1), 10)
  bin_labels <- paste0(format(breaks[-11], nsmall = 1), "-", format(breaks[-1], nsmall = 1))
  counts <- as.integer(table(factor(bin_idx, levels = 1:10)))
  p_value_distribution_df <- data.frame(
    p_range = bin_labels,
    count = counts,
    stringsAsFactors = FALSE
  )

uniform_expected <- n_tests / 10
  first_bin <- counts[1]
  dist_diagnosis <- if (first_bin >= 2 * uniform_expected && first_bin >= 3) {
    sprintf(paste0(
      "The lowest bin(0.0-0.1) holds %d of %d p-values — about %.1fx the ~%.0f a ",
      "uniform(all-null) distribution would put there. That spike near zero is the ",
      "classic signature of genuinely non-null effects among your tests."),
      first_bin, n_tests, first_bin / max(uniform_expected, 1e-9), uniform_expected)
  } else if (first_bin > uniform_expected) {
    sprintf(paste0(
      "The lowest bin(0.0-0.1) holds %d of %d p-values, modestly above the ~%.0f a ",
      "uniform(all-null) distribution would produce — weak evidence of a few true ",
      "effects mixed into mostly null tests."),
      first_bin, n_tests, uniform_expected)
  } else {
    sprintf(paste0(
      "The histogram is close to uniform(lowest bin: %d of %d p-values vs ~%.0f ",
      "expected under all-null) — consistent with mostly null effects, so treat any ",
      "uncorrected significant results with suspicion."),
      first_bin, n_tests, uniform_expected)
  }

  metrics <- list(
    `Tests Evaluated`         = n_tests,
    `Invalid P-Values Dropped` = as.integer(rows_removed),
    `Significant(Uncorrected)` = as.integer(n_sig_raw),
    `Significant(Holm)`      = as.integer(n_sig_holm),
    `Significant(BH)`        = as.integer(n_sig_bh),
    `Smallest P-Value`        = signif(min(p), 3)
  )

  json_output <- list(
    answer = paste0(
      "Of ", n_tests, " tests, ", n_sig_raw, " were significant uncorrected(p < 0.05); ",
      "Holm-Bonferroni keeps ", n_sig_holm, " (strict family-wise control) and ",
      "Benjamini-Hochberg keeps ", n_sig_bh, " (false-discovery-rate control). ",
      dist_diagnosis
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "adjusted_results",
        "significance_summary", "p_value_distribution"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows = initial_rows, final_rows = final_rows,
    rows_removed = rows_removed,
    p_name = p_name, label_name = label_name,
    n_tests = n_tests, n_na = n_na, n_out_of_range = n_out_of_range,
    alpha = alpha,
    results_all = results_all,
    adjusted_results_df = adjusted_results_df,
    results_truncated = results_truncated,
    n_sig_raw = n_sig_raw, n_sig_holm = n_sig_holm, n_sig_bh = n_sig_bh,
    significance_summary_df = significance_summary_df,
    p_value_distribution_df = p_value_distribution_df,
    dist_diagnosis = dist_diagnosis,
    metrics = metrics, json_output = json_output
  )
}

# Card: tldr (tldr)
card_tldr <- function(shared, df, params) {
  survival_word <- if (shared$n_sig_raw == 0) {
    "No result was significant even before correction."
  } else if (shared$n_sig_holm == shared$n_sig_raw) {
    "Every uncorrected finding survived even the strictest correction — these results are robust to multiple testing."
  } else if (shared$n_sig_holm == 0 && shared$n_sig_bh == 0) {
    "No result survived either correction — the uncorrected &#x27;significant' findings are consistent with multiple-testing luck."
  } else {
    sprintf("%d uncorrected finding(s) did not survive Holm — likely multiple-testing artifacts.",
            shared$n_sig_raw - shared$n_sig_holm)
  }
  text <- paste0(
    "Of ", shared$n_tests, " tests, ", shared$n_sig_raw,
    " were significant uncorrected(p < ", shared$alpha, "); Holm-Bonferroni keeps ",
    shared$n_sig_holm, " (strict family-wise control) and Benjamini-Hochberg keeps ",
    shared$n_sig_bh, " (false-discovery-rate control). ", survival_word
  )
  list(
    title = "Executive Summary",
    description = paste0("Multiple-testing correction across ", shared$n_tests, " tests"),
    metrics = shared$metrics,
    text = text
  )
}