Multiple Comparisons

The code that produced every figure in this report, exactly as it ran. Fingerprint d28e7aa14b6bcfe5. The same code on the same data gives the same report.

`standard_multiple_comparisons_v2` <- function(pf) {
  `%||%` <- function(a, b) if (!is.null(a)) a else b
  #' Readable figures (LAT-3180, LAT-3181): whole numbers from a thousand up, one decimal from a hundred, two from
  #' one, three significant figures below one. A cell carries what the value needs, not what R prints.
  tidy <- function(x) {
    x <- as.numeric(x)
    ifelse(is.na(x), NA_real_,
      ifelse(abs(x) >= 1000, round(x, 0),
        ifelse(abs(x) >= 100, round(x, 1),
          ifelse(abs(x) >= 1, round(x, 2), signif(x, 3)))))
  }
  #' P-VALUES BELOW 0.0001 LEAVE AS A FLOOR (LAT-3181). Results serialize at four decimal digits, which rounds 1.4e-05
  #' to 0; 1e-12 survives and the mapper shows it as "<0.0001". Table cells only; the answer keeps the computed p.
  p_cell <- function(p) { p <- as.numeric(p); ifelse(is.na(p), NA_real_, ifelse(p < 1e-4, 1e-12, signif(p, 3))) }
  inputs <- pf$taskList$inputs
  params <- inputs$module_parameters %||% list()
  # THE QUESTION this tool answers: the customer's objective, verbatim, when given.
  question <- (inputs$userContext %||% list())$objective %||%
    "I ran many tests: which results are still significant after correcting for multiple testing?"

  #' ## Column mapping
  #' The customer maps one `p_value` column and, optionally, a `test_label` column naming each test. Each row is one
  #' hypothesis test. Semantic names are used inside; the customer's own headers live in `col_map`.
  col_map <- inputs$column_mapping %||% list()
  df <- renderObject.taskFunction.init(inputs, col_map)   # df has SEMANTIC names
  human <- function(sem) {
    v <- col_map[[sem]]
    if (is.null(v) || !nzchar(as.character(v))) sem else as.character(v)
  }
  p_name <- human("p_value")
  has_label <- "test_label" %in% names(df)
  label_name <- if (has_label) human("test_label") else "test"
  if (!("p_value" %in% names(df)))
    stop(sprintf("column_mapping must map p_value to your column of raw p-values ('%s' was not found).", p_name))
  n_in <- nrow(df)

  #' ## The columns this tool did not look at
  #' Read from the RAW rows still on `inputs`, because `renderObject.taskFunction.init` has already narrowed `df` to the
  #' mapped columns. When the raw shape cannot be read the method says so instead of implying there were none.
  raw_names <- local({
    ds <- inputs$dataset %||% inputs$df
    if (is.data.frame(ds)) return(names(ds))
    if (is.list(ds) && length(ds) > 0) {
      rows <- ds[seq_len(min(length(ds), 50))]
      nm <- unique(unlist(lapply(rows, function(r) if (is.list(r)) names(r) else NULL)))
      if (length(nm)) return(nm)
      if (!is.null(names(ds)) && all(nzchar(names(ds)))) return(names(ds))
    }
    character(0)
  })
  mapped_actual <- unique(as.character(unlist(col_map)))
  ignored_cols <- setdiff(raw_names, unique(c(mapped_actual, make.names(mapped_actual))))

  #' ## Parameters
  #' `alpha`: the significance level every method is judged at, one of 0.01, 0.05 (the default) or 0.10.
  alpha <- suppressWarnings(as.numeric(params$alpha %||% 0.05))
  if (is.na(alpha) || !(alpha %in% c(0.01, 0.05, 0.1)))
    stop("module_parameters$alpha must be 0.01, 0.05 or 0.10")
  #' `family_size` (LAT-4597): how many tests were run in all, when the file holds only some of them (36 run, the 10
  #' interesting ones uploaded). Every correction then counts the whole family, which is what stops cherry-picking.
  #' Optional; empty means the file is the whole family. Checked against the tests read once they are counted below.
  family_raw <- params$family_size
  #' A list or vector of several numbers is not "not given": it is refused below, never silently dropped.
  family_given <- !is.null(family_raw) && (length(family_raw) != 1 || (!is.na(family_raw) && nzchar(trimws(as.character(family_raw)))))
  FAMILY_MAX <- 1000000L   # p.adjust takes n as an integer; a larger total is a typo, not a family

  #' ## Reading the p-values
  #' Numeric as given; text is read with a leading comparison stripped, so "<0.0001" (how p-values arrive from stats
  #' software and published tables) is analysed at its bound, which is conservative for significance. The 95% rule
  #' applies: a column where fewer than 95% of values read as numbers is refused. Missing, unreadable and out-of-range
  #' values are excluded and counted by reason.
  v <- df$p_value
  n_censored <- 0L
  if (is.numeric(v)) {
    conv <- as.numeric(v); unreadable <- rep(FALSE, length(v)); blank <- is.na(v)
  } else {
    raw_chr <- trimws(as.character(v)); raw_chr[is.na(raw_chr)] <- ""
    blank <- raw_chr == ""
    stripped <- sub("^(<=|>=|\u2264|\u2265|<|>)\\s*", "", raw_chr)
    conv <- suppressWarnings(as.numeric(stripped))
    n_nonblank <- sum(!blank)
    if (n_nonblank == 0 || sum(!is.na(conv[!blank])) < 0.95 * n_nonblank)
      stop(sprintf("The mapped p-value column ('%s') is not numeric: map a column of raw p-values between 0 and 1 (censored forms like '<0.0001' are accepted).", p_name))
    n_censored <- sum(raw_chr != stripped & !is.na(conv))
    unreadable <- !blank & is.na(conv)
  }
  out_of_range <- !is.na(conv) & (conv < 0 | conv > 1)
  keep <- !is.na(conv) & !out_of_range
  n_blank <- sum(blank); n_unreadable <- sum(unreadable); n_out <- sum(out_of_range)
  p <- as.numeric(conv[keep])
  m <- length(p)
  if (m < 5) stop(sprintf("Only %d valid p-values remained in '%s' (%d missing, %d unreadable, %d outside 0 to 1); a multiple-testing correction needs at least 5 tests.",
                          m, p_name, n_blank, n_unreadable, n_out))
  M <- m
  if (family_given) {
    M <- if (length(family_raw) == 1) suppressWarnings(as.numeric(unlist(family_raw))) else NA_real_
    if (is.na(M) || !is.finite(M) || M != round(M))
      stop("The total tests you ran must be a whole number.")
    if (M > FAMILY_MAX)
      stop("The total tests you ran must be 1,000,000 or fewer.")   # pipeline review LAT-4597: 1e10 and Inf died in as.integer
    if (M < m)
      stop("The total cannot be smaller than the number of tests in your file.")   # marketing, LAT-4597
    M <- as.integer(M)
  }
  partial <- M > m   # only some of the family was entered

  #' ## Test names
  #' The label column when mapped, else "Test 1", "Test 2" in file order. A blank label takes its row number. Repeated
  #' labels are kept apart with a counter and counted, because the same test entered twice enlarges the family.
  lab <- if (has_label) trimws(as.character(df$test_label)) else paste0("Test ", seq_len(n_in))
  lab[is.na(lab) | lab == ""] <- paste0("Test ", which(is.na(lab) | lab == ""))
  lab <- lab[keep]
  n_dup_labels <- sum(duplicated(lab))
  if (n_dup_labels > 0) lab <- make.unique(lab, sep = " #")

  #' ## The corrections
  #' Bonferroni and Holm control the family-wise error rate (the chance of even one false positive); Holm is uniformly
  #' more powerful than Bonferroni. Benjamini-Hochberg controls the false discovery rate (the expected share of false
  #' positives among the results called significant). All three from base R's p.adjust, over the whole family M
  #' (its n argument): equal to m unless family_size says more tests were run than entered.
  bonf <- p.adjust(p, "bonferroni", n = M); holm <- p.adjust(p, "holm", n = M); bh <- p.adjust(p, "BH", n = M)
  sig_raw <- p < alpha; sig_bonf <- bonf < alpha; sig_holm <- holm < alpha; sig_bh <- bh < alpha
  n_sig_raw <- sum(sig_raw); n_sig_bonf <- sum(sig_bonf); n_sig_holm <- sum(sig_holm); n_sig_bh <- sum(sig_bh)
  expected_false <- alpha * M
  #' Storey's estimate of the share of tests that are truly null (lambda = 0.5): p-values of null tests are uniform, so
  #' twice the share above one half estimates it. Capped at 100%; rough below about 50 tests, which the method says.
  #' Not estimated when only part of the family was entered: the tests left out are likely the null ones.
  null_share <- if (partial) NA_real_ else min(1, sum(p > 0.5) / (0.5 * m))

  survives <- ifelse(sig_bonf, "every correction", ifelse(sig_holm, "Holm and BH", ifelse(sig_bh, "BH only",
              ifelse(sig_raw, "uncorrected only", "none"))))
  ord <- order(p, seq_along(p))
  res <- data.frame(test = lab[ord], raw_p = p_cell(p[ord]), bonferroni_p = p_cell(bonf[ord]),
                    holm_p = p_cell(holm[ord]), bh_p = p_cell(bh[ord]), survives = survives[ord], stringsAsFactors = FALSE)
  rownames(res) <- NULL
  TABLE_SHOWN <- 60L
  res_shown <- head(res, TABLE_SHOWN)

  #' ## The frames for the places
  counts_df <- data.frame(standard = c("Uncorrected", "Benjamini-Hochberg", "Holm", "Bonferroni"),
                          significant = c(n_sig_raw, n_sig_bh, n_sig_holm, n_sig_bonf), stringsAsFactors = FALSE)

  #' Evidence by rank: each test's -log10 p against its rank, with the line each step-wise method compares it to. A
  #' p-value axis crushes every small p against zero; on the evidence scale a larger value is stronger evidence and the
  #' lines are the bars to clear. Holm's bar at rank k is alpha / (m - k + 1); Benjamini-Hochberg's is k * alpha / m.
  #' The first 100 ranks are drawn, which is where every decision in a family of any size is made.
  R <- min(m, 100L); k <- seq_len(R); ps <- sort(p)[k]
  ev <- function(x) tidy(-log10(pmax(x, 1e-300)))
  rank_df <- rbind(
    data.frame(rank = k, evidence = ev(ps), series = "Tests, smallest p first", stringsAsFactors = FALSE),
    data.frame(rank = k, evidence = ev(k * alpha / M), series = "Benjamini-Hochberg bar", stringsAsFactors = FALSE),
    data.frame(rank = k, evidence = ev(alpha / (M - k + 1)), series = "Holm bar", stringsAsFactors = FALSE))

  #' The p-value histogram against what an all-null family would give (a flat m / 10 per bin): a spike in the lowest bin
  #' is the signature of real effects; a flat shape says most tests are null.
  breaks <- seq(0, 1, by = 0.1)
  bin <- pmin(pmax(findInterval(p, breaks, rightmost.closed = TRUE), 1), 10)
  counts <- as.integer(table(factor(bin, levels = 1:10)))
  bin_lab <- sprintf("%.1f-%.1f", breaks[-11], breaks[-1])
  dist_df <- rbind(
    data.frame(p_range = bin_lab, count = counts, series = "Observed", stringsAsFactors = FALSE),
    data.frame(p_range = bin_lab, count = rep(tidy(m / 10), 10), series = "Expected by chance", stringsAsFactors = FALSE))

  #' The tests the methods disagree about: significant uncorrected but not under every correction. These are the
  #' decisions the choice of method actually changes. Conditional: written as a dropped place with its reason when
  #' there are none (LAT-3102).
  disputed <- res[res$survives %in% c("Holm and BH", "BH only", "uncorrected only"), c("test", "raw_p", "bh_p", "holm_p", "survives"), drop = FALSE]
  names(disputed)[names(disputed) == "survives"] <- "reading"
  disputed$reading <- c("Holm and BH" = "survives Holm and BH, not Bonferroni", "BH only" = "survives BH, not Holm",
                        "uncorrected only" = "significant only before correction")[disputed$reading]
  rownames(disputed) <- NULL

  #' ## Assumption checks (LAT-3138): one verdict per condition, each with the statistic a reader can re-derive.
  verdict_p <- function(pv, soft = 0.05, hard = 0.001) if (is.na(pv)) "unknown" else if (pv >= soft) "holds" else if (pv >= hard) "strained" else "violated"
  fmt_p <- function(pv) if (is.na(pv)) "" else if (pv < 1e-4) "<0.0001" else as.character(signif(pv, 3))
  n_excluded <- n_blank + n_unreadable + n_out
  excl_share <- n_excluded / max(1, n_in)
  upper <- p[p > 0.5]
  unif_p <- if (!partial && length(upper) >= 8) tryCatch(suppressWarnings(stats::ks.test((upper - 0.5) / 0.5, "punif")$p.value), error = function(e) NA_real_) else NA_real_
  rounded_share <- (sum(duplicated(p)) + n_censored) / m
  checks_df <- data.frame(
    check = c("Enough tests", "Valid p-values", "Each test entered once", "Null end looks uniform", "Few rounded or censored values"),
    statistic = c(paste0(m, " tests"),
                  paste0(n_excluded, " of ", n_in, " rows excluded"),
                  paste0(n_dup_labels, " repeated test names"),
                  if (partial) "not tested: only part of the family was entered"
                  else if (is.na(unif_p)) paste0(length(upper), " p-values above 0.5, too few to test") else paste0("uniformity of the ", length(upper), " p-values above 0.5"),
                  paste0(round(100 * rounded_share, 1), "% repeated or given as a bound")),
    p_value = c("", "", "", fmt_p(unif_p), ""),
    verdict = c(if (m >= 10) "holds" else "strained",
                if (excl_share < 0.05) "holds" else if (excl_share < 0.2) "strained" else "violated",
                if (n_dup_labels == 0) "holds" else if (n_dup_labels / m < 0.1) "strained" else "violated",
                verdict_p(unif_p),
                if (rounded_share < 0.1) "holds" else if (rounded_share < 0.3) "strained" else "violated"),
    note = c("with few tests the corrections change little and the null-share estimate is rough",
             "an excluded test is left out of the family, which makes every correction slightly less strict",
             "the same test entered twice enlarges the family and makes the corrections stricter than they should be",
             "p-values of null tests are uniform; a pile-up near one suggests rounding or tests that are not continuous",
             "rounded or censored p-values tie, and a bound such as <0.0001 is analysed at the bound"),
    stringsAsFactors = FALSE)

  #' ## Method text, assumptions and the answer
  excluded_rows <- c(if (n_blank > 0) sprintf("%d row%s with no p-value", n_blank, if (n_blank > 1) "s" else ""),
                     if (n_unreadable > 0) sprintf("%d row%s whose p-value could not be read as a number", n_unreadable, if (n_unreadable > 1) "s" else ""),
                     if (n_out > 0) sprintf("%d row%s with a p-value outside 0 to 1", n_out, if (n_out > 1) "s" else ""))
  method <- paste0(
    "Bonferroni, Holm and Benjamini-Hochberg corrections (base R p.adjust) of ", m, " p-values from '", p_name, "'",
    if (has_label) paste0(", each test named by '", label_name, "'") else ", tests named by row order",
    ": ", m, " tests used of ", n_in, " rows",
    if (partial) sprintf("; corrected for a family of %d tests (family_size), of which %d were entered", M, m) else "",
    if (length(excluded_rows)) paste0(" (excluded: ", paste(excluded_rows, collapse = ", "), ")") else "",
    if (n_censored > 0) sprintf("; %d censored value%s such as <0.0001 analysed at the bound, which is conservative", n_censored, if (n_censored > 1) "s" else "") else "",
    if (n_dup_labels > 0) sprintf("; %d repeated test name%s kept as separate tests", n_dup_labels, if (n_dup_labels > 1) "s" else "") else "",
    if (length(ignored_cols)) paste0("; not used: ", paste(utils::head(ignored_cols, 12), collapse = ", "), " (a correction needs only the p-values and their names)") else "",
    "; every method judged at alpha = ", alpha,
    if (partial) "; the share of truly null tests not estimated, because only part of the family was entered"
    else paste0("; the share of truly null tests estimated by Storey's method at lambda 0.5", if (m < 50) " (rough with fewer than 50 tests)" else ""),
    "; the evidence chart draws the first ", R, " ranks",
    if (nrow(res) > TABLE_SHOWN) sprintf("; the table lists the %d smallest p-values of %d", TABLE_SHOWN, nrow(res)) else "", ".")
  assumptions <- list(
    "Bonferroni and Holm control the chance of any false positive whatever the dependence between tests; Benjamini-Hochberg controls the share of false positives among discoveries for independent or positively related tests.",
    if (partial) sprintf("The family is the %d tests that were run (family_size), not only the %d entered; the corrections count all %d. Benjamini-Hochberg over a partial file is at least as strict as over the full family.", M, m, M)
    else "The family is the set of p-values given: tests that were run and not included are not corrected for.",
    "A corrected p-value says whether a result survives the correction; it is not the probability that the effect is real.",
    "Censored p-values such as <0.0001 are analysed at their bound, which can only make a result harder to keep.",
    "The estimated share of null tests assumes null p-values are uniform and is rough below about 50 tests.")
  strongest <- res$test[1]
  answer <- list(n_tests = m, alpha = alpha, n_significant_uncorrected = n_sig_raw, n_significant_bonferroni = n_sig_bonf,
                 n_significant_holm = n_sig_holm, n_significant_bh = n_sig_bh, expected_false_positives = tidy(expected_false),
                 estimated_null_share_pct = round(100 * null_share, 1), strongest_test = strongest,
                 smallest_p = if (min(p) < 1e-4) "<0.0001" else signif(min(p), 3), n_censored = n_censored, n = m)   # four-decimal serializer
  if (partial) answer$family_size <- M   # only when given and larger, so a default run's output is unchanged

  results <- list()
  #' The verdict and the headline are NOT places of a library tool (LAT-3130): the last mile writes them.
  results$summary_metrics <- place_metric(list(
    n_tests = m, n_sig_raw = n_sig_raw, n_sig_bh = n_sig_bh, n_sig_holm = n_sig_holm,
    expected_false_raw = tidy(expected_false), est_null_share_pct = round(100 * null_share, 1)),
    lead = "n_tests", place = "summary_metrics")
  if (partial) results$summary_metrics$values$family_size <- M
  results$survival_counts <- place_comparison(counts_df, category = "standard", value = "significant", place = "survival_counts")
  results$evidence_by_rank <- place_trend(rank_df, x = "rank", y = "evidence", series = "series", place = "evidence_by_rank")
  results$p_value_distribution <- place_comparison(dist_df, category = "p_range", value = "count", series = "series",
    place = "p_value_distribution")
  if (nrow(disputed) > 0) {
    results$disputed_tests <- place_table(disputed, place = "disputed_tests")
  } else {
    results$disputed_tests <- place_dropped(if (n_sig_raw == 0)
      "no test was significant even before correction, so there is no result for the methods to disagree about"
      else "every test that was significant before correction survives every correction, so the methods agree on all of them",
      place = "disputed_tests")
  }
  results$results_table <- place_table(res_shown, place = "results_table")
  results$assumption_checks <- place_table(checks_df, place = "assumption_checks")
  results$comparison_method <- list(kind = "metric", values = list(
    method = method, n_in = n_in, n_used = m, excluded = as.list(excluded_rows), assumptions = assumptions,
    # the method card prints x_column -> y_column: the columns read, and what the correction produced
    x_column = if (has_label) paste(label_name, p_name, sep = ", ") else p_name,
    y_column = sprintf("%d tests corrected", m)),
    value_order = list("n_used", "n_in"))

  objects <- list()   # filled by the object layer, not here
  list(answer = answer, method = method, n = m, results = results, objects = objects,
       json_output = list(answer = answer, method = method, n = m))
}