Categorical Tests

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

`standard_categorical_tests_v2` <- function(pf) {
  `%||%` <- function(a, b) if (!is.null(a)) a else b
  #' Readable figures (LAT-3181): whole numbers from a thousand up, one decimal from a hundred, two from one,
  #' three significant figures below one.
  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_FLOOR <- 2e-16
  p_text <- function(p) if (is.na(p)) "could not be computed" else if (p < P_FLOOR) "< 2e-16" else paste0("= ", format(signif(p, 3)))
  #' P-VALUES BELOW 0.0001 LEAVE AS A FLOOR (LAT-3181): the results serialize at four decimal digits, so 1.4e-05 would
  #' arrive as 0. 1e-12 survives and the mapper shows any p under 0.0001 as "<0.0001". Table cells only.
  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()
  question <- (inputs$userContext %||% list())$objective %||%
    "Are the two categories related, which combinations are over-represented, and how far apart are the rates?"

  col_map <- inputs$column_mapping %||% list()
  raw <- inputs$dataset
  raw_names <- if (is.data.frame(raw)) names(raw) else if (is.list(raw) && length(raw)) names(raw[[1]]) else character(0)
  unmapped <- setdiff(raw_names, as.character(unlist(col_map)))
  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) }
  group_name <- human("group"); outcome_name <- human("outcome")

  max_levels <- suppressWarnings(as.integer(params$max_levels %||% 10L))
  if (is.na(max_levels) || max_levels < 2 || max_levels > 12) stop("module_parameters$max_levels must be an integer between 2 and 12")
  positive_level <- trimws(as.character(params$positive_level %||% ""))

  n_in <- nrow(df)
  for (need in c("group", "outcome")) {
    if (!need %in% names(df)) stop(sprintf("column_mapping must map a '%s' column; this analysis compares two categorical columns", need))
  }
  g <- trimws(as.character(df$group)); o <- trimws(as.character(df$outcome))
  keep <- !is.na(g) & nzchar(g) & !is.na(o) & nzchar(o)
  #' ## A file that is ALREADY a summary carries one line per combination and a count of the rows it stands for.
  w <- rep(1, length(g))
  weighted <- FALSE
  if ("count" %in% names(df)) {
    cv <- suppressWarnings(as.numeric(df$count))
    if (sum(!is.na(cv)) >= 0.95 * sum(keep)) {
      if (any(cv[keep] < 0, na.rm = TRUE)) stop(sprintf("The count column '%s' has negative values; it must be a whole count of rows.", human("count")))
      w <- ifelse(is.na(cv), 0, cv); weighted <- TRUE
      keep <- keep & !is.na(cv) & cv > 0
    }
  }
  n_bad <- sum(!keep)
  g <- g[keep]; o <- o[keep]; w <- w[keep]
  if (!length(g)) stop(sprintf("No rows with a readable value in both '%s' and '%s' remained after cleaning.", group_name, outcome_name))

  #' ## Levels: the largest are kept and the rest folded into Other, on BOTH columns, and the folding is named.
  fold <- function(x, wt, what) {
    tab <- sort(tapply(wt, x, sum), decreasing = TRUE)
    if (length(tab) <= max_levels) return(list(x = x, folded = character(0)))
    folded <- setdiff(names(tab), names(tab)[seq_len(max_levels)])
    x[x %in% folded] <- "Other"
    list(x = x, folded = sprintf("%d %s level(s) folded into Other", length(folded), what))
  }
  fg <- fold(g, w, group_name); g <- fg$x
  fo <- fold(o, w, outcome_name); o <- fo$x
  n_used <- sum(w)
  tb <- tapply(w, list(g, o), sum); tb[is.na(tb)] <- 0
  tb <- as.table(as.matrix(tb))
  k_g <- nrow(tb); k_o <- ncol(tb)
  if (k_g < 2) stop(sprintf("'%s' has only one level after cleaning; there is nothing to compare.", group_name))
  if (k_o < 2) stop(sprintf("'%s' has only one level after cleaning; every row has the same outcome.", outcome_name))
  if (n_used < 20) stop(sprintf("Only %s rows remained after cleaning; at least 20 are needed to test a table.", format(n_used)))

  #' ## The tests. Chi-square is the primary reading; Fisher leads when the table is 2x2 or an expected count is small,
  #' because the chi-square p-value is unreliable there.
  cs <- suppressWarnings(stats::chisq.test(tb, correct = FALSE))
  cs_corr <- if (k_g == 2 && k_o == 2) suppressWarnings(stats::chisq.test(tb, correct = TRUE)) else NULL
  expected <- cs$expected
  min_expected <- min(expected)
  share_small <- mean(expected < 5)
  fisher <- NULL
  if ((k_g == 2 && k_o == 2) || min_expected < 5) {
    fisher <- tryCatch(stats::fisher.test(tb, simulate.p.value = (k_g * k_o > 4), B = 20000), error = function(e) NULL)
  }
  #' Cramer's V, and Bergsma's bias-corrected form. V is NOT a correlation: it has no sign, and what counts as large
  #' depends on the table's shape, which is why the method text carries the table's dimensions beside it.
  chi <- unname(cs$statistic); dfree <- unname(cs$parameter)
  v <- sqrt((chi / n_used) / min(k_g - 1, k_o - 1))
  phi2 <- chi / n_used
  phi2c <- max(0, phi2 - (k_g - 1) * (k_o - 1) / (n_used - 1))
  kc <- k_g - (k_g - 1)^2 / (n_used - 1); oc <- k_o - (k_o - 1)^2 / (n_used - 1)
  v_corr <- sqrt(phi2c / max(1e-12, min(kc - 1, oc - 1)))
  leading <- if (!is.null(fisher) && (min_expected < 5 || (k_g == 2 && k_o == 2))) "Fisher's exact test" else "chi-square test of independence"
  lead_p <- if (leading == "Fisher's exact test") fisher$p.value else unname(cs$p.value)

  #' ## The rate cards track ONE outcome level: the one asked for, else the most common.
  o_levels <- colnames(tb)
  #' WHICH LEVEL THE RATE CARDS TRACK. The most common level is the wrong default on a two-level outcome: a churn
  #' column is mostly "No", so the page reads as the share who did NOT cancel and the question was about cancelling
  #' (the first live package run on telco did exactly that). For a two-level outcome the level that READS as the event
  #' leads; otherwise the most common. `positive_level` overrides either.
  EVENT_LEVELS <- c("yes", "y", "true", "t", "1", "churn", "churned", "cancelled", "canceled", "failed", "fail",
                    "converted", "positive", "pass", "lost", "left", "dead", "relapsed", "default", "defaulted")
  pos <- if (nzchar(positive_level) && positive_level %in% o_levels) positive_level else {
    ev <- o_levels[tolower(trimws(o_levels)) %in% EVENT_LEVELS]
    if (k_o == 2 && length(ev) == 1) ev else o_levels[which.max(colSums(tb))]
  }
  binary <- k_o == 2
  totals <- rowSums(tb)
  hits <- tb[, pos]
  rate_ci <- function(x, n) {   # Wilson, which behaves at 0 and 1 where the textbook interval does not
    if (n == 0) return(c(NA_real_, NA_real_))
    z <- 1.959964; p <- x / n; d <- 1 + z^2 / n
    centre <- (p + z^2 / (2 * n)) / d; half <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2)) / d
    c(max(0, centre - half), min(1, centre + half))
  }
  ci <- t(vapply(seq_len(k_g), function(i) rate_ci(hits[i], totals[i]), numeric(2)))
  comp_df <- data.frame(group = rownames(tb), rate_pct = tidy(100 * hits / totals),
                        low = tidy(100 * ci[, 1]), high = tidy(100 * ci[, 2]), stringsAsFactors = FALSE)
  comp_df <- comp_df[order(-comp_df$rate_pct), , drop = FALSE]

  #' ## The counts, as the reader checks everything else against.
  #' LONG, not a cross-tab: a wide table carries one column per outcome level, so its columns cannot be declared in
  #' the spec before the data arrives, and a place whose columns the spec cannot name lets the estimate door guess them.
  cont_df <- do.call(rbind, lapply(rownames(tb), function(gi) data.frame(
    group = gi, outcome = o_levels, count = as.integer(tb[gi, ]),
    share_of_group_pct = tidy(100 * as.numeric(tb[gi, ]) / max(1, totals[[gi]])), stringsAsFactors = FALSE)))
  rownames(cont_df) <- NULL

  #' ## Where independence breaks: the standardized (adjusted) residual per cell, which is on a standard-normal scale,
  #' so beyond about 2 is the usual mark of a cell worth naming. POSITIVE means the combination happens MORE than
  #' chance. The frame is long, one row per cell, and the group is x and the outcome is y in that order.
  resid <- cs$stdres
  res_df <- do.call(rbind, lapply(rownames(tb), function(gi) data.frame(
    group = gi, outcome = o_levels, residual = tidy(as.numeric(resid[gi, ])), stringsAsFactors = FALSE)))
  rownames(res_df) <- NULL
  #' ORDERED ON ABSOLUTE SIZE, and the label says so, because a bar to the left is as far from chance as the same bar
  #' to the right (LAT-3231: an ordering claim must say what it orders on).
  top_df <- res_df[order(-abs(res_df$residual)), , drop = FALSE]
  top_df <- utils::head(top_df, min(10, nrow(top_df)))
  top_df <- data.frame(cell = paste0(top_df$group, " and ", top_df$outcome), residual = top_df$residual, stringsAsFactors = FALSE)

  #' ## How far apart the groups are: every pair, as a difference in rates with a 95% interval and a two-proportion
  #' test, Holm-adjusted across the pairs. Only when the outcome has exactly two levels, where a rate means one thing.
  pair_df <- NULL
  if (binary && k_g >= 2) {
    combos <- utils::combn(seq_len(k_g), 2)
    rows <- lapply(seq_len(ncol(combos)), function(j) {
      a <- combos[1, j]; b <- combos[2, j]
      x <- c(hits[a], hits[b]); nn <- c(totals[a], totals[b])
      if (any(nn == 0)) return(NULL)
      pt <- suppressWarnings(stats::prop.test(x, nn, correct = TRUE))
      data.frame(comparison = paste0(rownames(tb)[a], " vs ", rownames(tb)[b]),
                 difference = tidy(100 * (x[1] / nn[1] - x[2] / nn[2])),
                 low = tidy(100 * pt$conf.int[1]), high = tidy(100 * pt$conf.int[2]),
                 p_raw = unname(pt$p.value), stringsAsFactors = FALSE)
    })
    pair_df <- do.call(rbind, Filter(Negate(is.null), rows))
    if (!is.null(pair_df) && nrow(pair_df)) {
      pair_df$p_value <- p_cell(stats::p.adjust(pair_df$p_raw, method = "holm"))
      pair_df$p_raw <- NULL
      pair_df <- pair_df[order(-abs(pair_df$difference)), , drop = FALSE]
      rownames(pair_df) <- NULL
    }
  }

  #' ## Every test in one table. The rows are in different units, which is why this is a table and not a chart.
  test_rows <- list()
  add <- function(test, statistic, p, estimate, note)
    test_rows[[length(test_rows) + 1]] <<- data.frame(test = test, statistic = if (is.na(statistic)) NA_real_ else tidy(statistic),
      p_value = if (is.na(p)) NA_real_ else p_cell(p), estimate = estimate, note = note, stringsAsFactors = FALSE)
  add(sprintf("Chi-square test of independence (%d by %d)", k_g, k_o), chi, unname(cs$p.value), sprintf("%d degrees of freedom", dfree),
      "whether the two columns are related at all")
  if (!is.null(cs_corr)) add("Chi-square with continuity correction", unname(cs_corr$statistic), unname(cs_corr$p.value), "1 degree of freedom",
      "the conservative form for a two by two table")
  if (!is.null(fisher)) add("Fisher's exact test", NA_real_, fisher$p.value,
      if (!is.null(fisher$estimate)) sprintf("odds ratio %s", format(tidy(unname(fisher$estimate)))) else "exact",
      "exact, and the reading to trust when an expected count is below 5")
  add("Cramer's V", v, NA_real_, sprintf("%s on a %d by %d table", format(tidy(v)), k_g, k_o),
      "strength of the relation, 0 to 1, with no direction; it is not a correlation")
  add("Cramer's V, bias corrected", v_corr, NA_real_, format(tidy(v_corr)), "the form to quote for a small table or a small sample")
  if (binary && k_g == 2) {
    pt2 <- suppressWarnings(stats::prop.test(c(hits[1], hits[2]), c(totals[1], totals[2]), correct = TRUE))
    add("Two-proportion test", unname(pt2$statistic), unname(pt2$p.value),
        sprintf("%s points apart", format(tidy(100 * (hits[1] / totals[1] - hits[2] / totals[2])))),
        sprintf("the gap in the share of %s between the two groups", pos))
  }
  tests_df <- do.call(rbind, test_rows)

  #' ## Assumption checks (LAT-3138)
  empty_cells <- sum(tb == 0)
  checks_df <- data.frame(
    check = c("Expected counts of 5 or more", "No empty combinations", "Levels not capped away", "Enough rows", "Rows complete"),
    statistic = c(sprintf("smallest expected count %s; %s%% of cells below 5", format(tidy(min_expected)), format(tidy(100 * share_small))),
                  sprintf("%d of %d combinations have no rows", empty_cells, k_g * k_o),
                  if (length(c(fg$folded, fo$folded))) paste(c(fg$folded, fo$folded), collapse = "; ") else sprintf("all %d by %d levels kept", k_g, k_o),
                  sprintf("%s rows across %d combinations", format(n_used), k_g * k_o),
                  sprintf("%d of %d rows excluded", n_bad, n_in)),
    p_value = NA_real_,
    verdict = c(if (min_expected >= 5) "holds" else if (!is.null(fisher)) "strained" else "violated",
                if (empty_cells == 0) "holds" else if (empty_cells <= 0.1 * k_g * k_o) "strained" else "violated",
                if (!length(c(fg$folded, fo$folded))) "holds" else "strained",
                if (n_used >= 20 * k_g * k_o) "holds" else if (n_used >= 5 * k_g * k_o) "strained" else "violated",
                if (n_in == 0 || n_bad / n_in <= 0.05) "holds" else if (n_bad / n_in <= 0.15) "strained" else "violated"),
    note = c("below 5 the chi-square p-value is unreliable, and Fisher's exact test is the reading to trust",
             "a combination with no rows makes the residual for that cell meaningless",
             "folding levels into Other hides differences between the folded levels",
             "a sparse table gives every cell a wide margin of error",
             "rows missing either column are left out"),
    stringsAsFactors = FALSE)

  #' ## Method, assumptions, answer
  excluded <- character(0)
  if (n_bad > 0) excluded <- c(excluded, sprintf("%d row%s missing one of the two columns", n_bad, if (n_bad > 1) "s" else ""))
  if (length(c(fg$folded, fo$folded))) excluded <- c(excluded, paste(c(fg$folded, fo$folded), collapse = "; "))
  if (length(unmapped)) excluded <- c(excluded, sprintf("columns not used by this analysis: %s", paste(unmapped, collapse = ", ")))
  method <- paste0(
    "A ", k_g, " by ", k_o, " table of '", group_name, "' against '", outcome_name, "' over ", format(n_used), " rows",
    if (weighted) sprintf(" (a pre-aggregated file: '%s' carries the rows each line stands for)", human("count")) else "",
    ". Chi-square test of independence without the continuity correction",
    if (!is.null(cs_corr)) ", the corrected form beside it" else "",
    if (!is.null(fisher)) paste0("; Fisher's exact test", if (k_g * k_o > 4) " (p simulated over 20,000 tables)" else "",
                                 ", which leads because ", if (min_expected < 5) sprintf("the smallest expected count is %s", format(tidy(min_expected))) else "the table is two by two") else "",
    "; strength as Cramer's V (", format(tidy(v)), " on this ", k_g, " by ", k_o, " table, which has no direction and is not a correlation)",
    "; the share of '", pos, "' per group with Wilson 95% intervals (the level the question is about",
    if (nzchar(positive_level) && positive_level %in% o_levels) ", as asked for" else if (k_o == 2) ", read from the level names" else ", the most common level", ")",
    if (!is.null(pair_df) && nrow(pair_df)) "; every pair of groups compared with a two-proportion test, Holm-adjusted across the pairs" else "",
    "; standardized residuals per cell, where a positive value means the combination happens more often than independence predicts",
    if (length(excluded)) paste0("; excluded: ", paste(excluded, collapse = "; ")) else "", ".")
  assumptions <- list(
    "Every row is assumed to be one independent case; a file with repeated measurements of the same subject breaks this and the p-values are then too small.",
    "Cramer's V has no direction and its size depends on the table's shape, so it is not comparable across tables of different sizes.",
    "A standardized residual says how far a combination sits from independence, not how common it is.",
    if (min_expected < 5) "An expected count below 5 makes the chi-square p-value unreliable; Fisher's exact test is the reading to trust here." else "Expected counts are large enough for the chi-square reading.",
    "The table shows association, not cause: something else may drive both columns.")
  strongest <- top_df[1, ]   # ordered on ABSOLUTE residual, so this is the cell furthest from chance in either direction
  answer <- paste0(
    "'", group_name, "' and '", outcome_name, "' ",
    if (lead_p < 0.05) "are related" else "cannot be shown to be related",
    " (", leading, ", p ", p_text(lead_p), "; Cramer's V ", format(tidy(v)), " on a ", k_g, " by ", k_o, " table). ",
    "The share of '", pos, "' runs from ", format(comp_df$rate_pct[nrow(comp_df)]), "% in ", comp_df$group[nrow(comp_df)],
    " to ", format(comp_df$rate_pct[1]), "% in ", comp_df$group[1],
    if (abs(strongest$residual) < 2)
      paste0(", and no combination departs materially from chance: the largest standardized residual is ",
             format(strongest$residual), " (", strongest$cell, "), within the range ordinary sampling produces.")
    else paste0(", and the combination furthest from chance is ", strongest$cell, " (standardized residual ",
                format(strongest$residual), ")."))

  #' ## The places. The verdict and the headline are NOT places of a library tool (LAT-3130).
  results <- list()
  summary_vals <- list(largest_gap = tidy(max(comp_df$rate_pct) - min(comp_df$rate_pct)), cramers_v = tidy(v),
                       p_value = p_cell(lead_p), leading_test = leading, n = n_used,
                       group = sprintf("%d levels", k_g), outcome = sprintf("%d levels", k_o))
  results$summary_metrics <- place_metric(summary_vals, lead = "largest_gap", place = "summary_metrics")
  results$composition <- place_comparison(comp_df, category = "group", value = "rate_pct", low = "low", high = "high", place = "composition")
  results$contingency_table <- place_table(cont_df, place = "contingency_table")
  results$residual_matrix <- place_matrix(res_df, x = "group", y = "outcome", z = "residual", place = "residual_matrix")
  results$top_deviations <- place_comparison(top_df, category = "cell", value = "residual", place = "top_deviations")
  #' A CONDITIONAL PLACE IS WRITTEN EITHER WAY (LAT-3102).
  if (!is.null(pair_df) && nrow(pair_df)) {
    results$pairwise_proportions <- place_interval(pair_df, term = "comparison", value = "difference", low = "low", high = "high",
                                                   place = "pairwise_proportions")
  } else {
    results$pairwise_proportions <- place_dropped(sprintf(
      "'%s' has %d levels, and a difference in rates only means one thing when it has exactly two, so no pairwise gap is shown; the table, the residuals and the tests above cover every level",
      outcome_name, k_o), place = "pairwise_proportions")
  }
  results$test_results <- place_table(tests_df, place = "test_results")
  results$assumption_checks <- place_table(checks_df, place = "assumption_checks")
  results$categorical_method <- list(kind = "metric", values = list(
    method = method, n_in = n_in, n_used = n_used, excluded = as.list(unname(excluded)), assumptions = assumptions,
    x_column = group_name, y_column = outcome_name), value_order = list("n_used", "n_in"))

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