Adjusted Group Comparison (ANCOVA)

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

`standard_ancova_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_text <- function(p) if (is.na(p)) "" else if (p < 1e-4) "<0.0001" else format(signif(p, 3))
  inputs <- pf$taskList$inputs
  params <- inputs$module_parameters %||% list()
  question <- (inputs$userContext %||% list())$objective %||%
    "Do these groups differ on the outcome once the covariate is adjusted for, and which groups differ from which?"

  #' ## Column mapping
  #' `outcome` (numeric), `group` (any) and `covariate` (numeric). Semantic names 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)
  }
  for (sem in c("outcome", "group", "covariate"))
    if (!(sem %in% names(df))) stop(sprintf("column_mapping must map '%s' (%s was not found).", sem, human(sem)))
  outcome_h <- human("outcome"); group_h <- human("group"); cov_h <- human("covariate")
  n_in <- nrow(df)

  #' ## The columns this tool did not look at, read from the raw rows (init narrows `df` to the mapped columns)
  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
  max_groups <- suppressWarnings(as.integer(params$max_groups %||% 8L))
  if (is.na(max_groups) || max_groups < 2 || max_groups > 12) stop("module_parameters$max_groups must be an integer from 2 to 12")

  #' ## Data preparation
  #' Outcome and covariate coerced to numeric; a row missing any of the three is excluded and counted. Group labels
  #' trimmed; the largest `max_groups` kept and the rest folded into Other (counted); a group with fewer than 3 rows
  #' is excluded and named. At least 2 groups and 10 rows are required.
  y <- suppressWarnings(as.numeric(as.character(df$outcome)))
  x <- suppressWarnings(as.numeric(as.character(df$covariate)))
  g <- trimws(as.character(df$group)); g[is.na(g) | g == "" | tolower(g) %in% c("na", "n/a", "null", "nan")] <- NA
  keep <- !is.na(y) & !is.na(x) & !is.na(g)
  n_missing <- sum(!keep)
  y <- y[keep]; x <- x[keep]; g <- g[keep]
  if (length(y) < 10) stop(sprintf("Only %d rows carry an outcome, a covariate and a group; the adjusted comparison needs at least 10.", length(y)))
  tab <- sort(table(g), decreasing = TRUE)
  n_folded <- 0L
  if (length(tab) > max_groups) {
    keep_levels <- names(tab)[seq_len(max_groups - 1)]
    n_folded <- length(tab) - length(keep_levels)
    g[!(g %in% keep_levels)] <- "Other"
  }
  small <- names(which(table(g) < 3))
  if (length(small)) { drop_rows <- g %in% small; y <- y[!drop_rows]; x <- x[!drop_rows]; g <- g[!drop_rows] }
  levs <- names(sort(table(g), decreasing = TRUE)); k <- length(levs)
  if (k < 2) stop("Fewer than two groups with at least three rows each; there is nothing to compare.")
  gf <- factor(g, levels = levs); n <- length(y)
  if (stats::sd(x) == 0) stop(sprintf("The covariate (%s) is constant, so there is nothing to adjust for; use the plain group comparison.", cov_h))

  #' ## The model: outcome ~ covariate + group, and the interaction model for the equal-slopes check
  dat <- data.frame(y = y, x = x, gf = gf)
  fit <- stats::lm(y ~ x + gf, data = dat)
  if (is.na(stats::coef(fit)["x"])) stop(sprintf("The covariate (%s) is fixed within each group, so its slope cannot be separated from the groups; use the plain group comparison.", cov_h))
  df_res <- stats::df.residual(fit)
  if (df_res < 3) stop("Too few rows for the number of groups; fewer than three residual degrees of freedom remain.")
  a_tab <- stats::anova(fit)
  cov_F <- a_tab["x", "F value"]; cov_p <- a_tab["x", "Pr(>F)"]
  group_F <- a_tab["gf", "F value"]; group_p <- a_tab["gf", "Pr(>F)"]
  ss_g <- a_tab["gf", "Sum Sq"]; ss_r <- a_tab["Residuals", "Sum Sq"]
  partial_eta2 <- ss_g / (ss_g + ss_r)
  slope <- unname(stats::coef(fit)["x"])
  slope_ci <- suppressWarnings(as.numeric(stats::confint(fit, "x", level = 0.95)))
  slope_p <- summary(fit)$coefficients["x", 4]
  fit_int <- stats::lm(y ~ x * gf, data = dat)
  cmp <- tryCatch(stats::anova(fit, fit_int), error = function(e) NULL)
  slopes_F <- if (!is.null(cmp)) cmp$F[2] else NA_real_
  slopes_p <- if (!is.null(cmp)) cmp$`Pr(>F)`[2] else NA_real_

  #' ## Adjusted means at the covariate mean, raw means, pairwise adjusted differences (Holm)
  xbar <- mean(x); tcrit <- stats::qt(0.975, df_res)
  pr <- stats::predict(fit, newdata = data.frame(x = rep(xbar, k), gf = factor(levs, levels = levs)), se.fit = TRUE)
  adj <- as.numeric(pr$fit); adj_se <- as.numeric(pr$se.fit)
  raw <- as.numeric(tapply(y, gf, mean)); grp_n <- as.integer(table(gf))
  adj_df <- data.frame(group = levs, adjusted_mean = tidy(adj), ci_low = tidy(adj - tcrit * adj_se), ci_high = tidy(adj + tcrit * adj_se), stringsAsFactors = FALSE)
  rva_df <- data.frame(group = rep(levs, 2), mean = tidy(c(raw, adj)), kind = rep(c("Raw", "Adjusted"), each = k), stringsAsFactors = FALSE)
  b <- stats::coef(fit); V <- stats::vcov(fit)
  cvec <- function(l) { cv <- stats::setNames(rep(0, length(b)), names(b)); dn <- paste0("gf", l); if (dn %in% names(b)) cv[dn] <- 1; cv }
  pair_df <- NULL; n_sig <- 0L
  if (k >= 3) {
    rows <- list()
    for (i in seq_len(k - 1)) for (j in (i + 1):k) {
      cw <- cvec(levs[i]) - cvec(levs[j]); est <- sum(cw * b); se <- sqrt(as.numeric(t(cw) %*% V %*% cw))
      tv <- if (se > 0) est / se else NA_real_; pv <- if (!is.na(tv)) 2 * stats::pt(-abs(tv), df_res) else NA_real_
      rows[[length(rows) + 1]] <- data.frame(comparison = paste0(levs[i], " - ", levs[j]), adj_difference = est, ci_low = est - tcrit * se, ci_high = est + tcrit * se, p_raw = pv, stringsAsFactors = FALSE)
    }
    pair_df <- do.call(rbind, rows)
    pair_df$adj_p <- stats::p.adjust(pair_df$p_raw, method = "holm"); pair_df$p_raw_adj <- pair_df$adj_p
    pair_df <- pair_df[order(-abs(pair_df$adj_difference)), c("comparison", "adj_difference", "ci_low", "ci_high", "adj_p", "p_raw_adj")]
    pair_df$adj_difference <- tidy(pair_df$adj_difference); pair_df$ci_low <- tidy(pair_df$ci_low); pair_df$ci_high <- tidy(pair_df$ci_high); pair_df$adj_p <- vapply(pair_df$adj_p, p_text, character(1))
    n_sig <- sum(!is.na(pair_df$p_raw_adj) & pair_df$p_raw_adj < 0.05)
    pair_df$p_raw_adj <- NULL
    rownames(pair_df) <- NULL
  }
  table_df <- data.frame(
    term = c(cov_h, group_h, "Residuals", paste0("Slope of ", cov_h)),
    df = c(a_tab["x", "Df"], a_tab["gf", "Df"], a_tab["Residuals", "Df"], NA),
    sum_sq = tidy(c(a_tab["x", "Sum Sq"], ss_g, ss_r, NA)),
    statistic = tidy(c(cov_F, group_F, NA, slope / summary(fit)$coefficients["x", 2])),
    p_value = c(p_text(cov_p), p_text(group_p), "", p_text(slope_p)),
    estimate = tidy(c(NA, NA, NA, slope)), ci_low = tidy(c(NA, NA, NA, slope_ci[1])), ci_high = tidy(c(NA, NA, NA, slope_ci[2])),
    stringsAsFactors = FALSE)
  scatter <- data.frame(x = x, y = y, group = g, stringsAsFactors = FALSE)
  names(scatter) <- c(cov_h, outcome_h, group_h)
  if (nrow(scatter) > 1000) { set.seed(20260918); scatter <- scatter[sort(sample.int(nrow(scatter), 1000)), , drop = FALSE] }

  #' ## Assumption checks (LAT-3138)
  cov_bal_p <- tryCatch(stats::anova(stats::lm(x ~ gf, data = dat))$`Pr(>F)`[1], error = function(e) NA_real_)
  res <- stats::residuals(fit)
  lin_p <- tryCatch({ f2 <- stats::lm(y ~ x + I(x^2) + gf, data = dat); stats::anova(fit, f2)$`Pr(>F)`[2] }, error = function(e) NA_real_)
  sw_p <- if (n >= 3 && n <= 5000) tryCatch(stats::shapiro.test(res)$p.value, error = function(e) NA_real_) else NA_real_
  lev_p <- tryCatch({ ad <- abs(res - stats::ave(res, gf, FUN = stats::median)); stats::anova(stats::lm(ad ~ gf))$`Pr(>F)`[1] }, error = function(e) NA_real_)
  min_n <- min(grp_n)
  v3 <- function(p, hold = 0.05, strain = 0.01) if (is.na(p)) "strained" else if (p >= hold) "holds" else if (p >= strain) "strained" else "violated"
  checks_df <- data.frame(
    check = c("Equal slopes across groups", sprintf("%s balanced across groups", cov_h), "Linear covariate relation", "Residuals are normal", "Equal residual spread across groups", "Enough rows per group"),
    statistic = c(if (is.na(slopes_F)) "interaction F not estimable" else sprintf("interaction F = %s", format(tidy(slopes_F))),
                  if (is.na(cov_bal_p)) "not estimable" else sprintf("one-way F on %s across %s", cov_h, group_h),
                  if (is.na(lin_p)) "not estimable" else "added squared term of the covariate",
                  if (is.na(sw_p)) sprintf("Shapiro-Wilk not run (n = %d)", n) else "Shapiro-Wilk on the residuals",
                  if (is.na(lev_p)) "not estimable" else "Levene (median) on the residuals",
                  sprintf("smallest group n = %d", min_n)),
    p_value = c(p_text(slopes_p), p_text(cov_bal_p), p_text(lin_p), p_text(sw_p), p_text(lev_p), ""),
    verdict = c(v3(slopes_p), if (is.na(cov_bal_p)) "strained" else if (cov_bal_p >= 0.05) "holds" else "strained", v3(lin_p),
                if (is.na(sw_p)) "strained" else if (sw_p >= 0.05) "holds" else if (n >= 30) "strained" else "violated",
                v3(lev_p), if (min_n >= 20) "holds" else if (min_n >= 10) "strained" else "violated"),
    note = c("when the slopes differ, the group difference depends on the covariate level and one adjusted mean per group misleads",
             "a covariate that differs between groups means the adjustment carries the comparison; the adjusted means then extrapolate",
             "a curved relation makes the linear adjustment incomplete",
             "with 30 or more rows the F test tolerates non-normal residuals; the intervals are approximate",
             "unequal spread makes the pooled residual variance wrong for some groups",
             "few rows in a group widen its interval and weaken the equal-slopes check"),
    stringsAsFactors = FALSE)

  #' ## Method, assumptions, answer
  excluded <- c(if (n_missing > 0) sprintf("%d row%s missing the outcome, the covariate or the group", n_missing, if (n_missing == 1) "" else "s"),
                if (n_folded > 0) sprintf("%d smaller group%s folded into Other", n_folded, if (n_folded == 1) "" else "s"),
                if (length(small)) sprintf("%d group%s with fewer than 3 rows excluded (%s)", length(small), if (length(small) == 1) "" else "s", paste(small, collapse = ", ")))
  method <- paste0(
    "Analysis of covariance of ", outcome_h, " on ", group_h, " (", k, " groups: ", paste(levs, collapse = ", "), ") adjusting for ", cov_h,
    ", fitted as outcome ~ covariate + group by least squares on ", n, " rows of ", n_in, "; group means adjusted to the covariate's mean (", format(tidy(xbar)),
    ") with t intervals; the group effect tested after the covariate (sequential F); partial eta squared from the group and residual sums of squares",
    if (k >= 3) sprintf("; %d pairwise adjusted differences with t intervals and Holm-adjusted p-values, %d differing at p < 0.05", nrow(pair_df), n_sig) else "; with two groups the overall test is the pairwise comparison",
    "; equal slopes checked by the interaction F of the covariate by group model against the additive one",
    if (length(excluded)) paste0("; excluded: ", paste(excluded, collapse = "; ")) else "",
    if (length(ignored_cols)) paste0("; not used: ", paste(utils::head(ignored_cols, 12), collapse = ", "), " (not mapped)") else "", ".")
  assumptions <- list(
    "The covariate relates to the outcome with the same slope in every group; when it does not, the adjusted means mislead.",
    "The covariate was measured before, or independently of, whatever made the groups different; adjusting for a consequence of the group removes the effect.",
    "The groups are compared at the covariate's overall mean; where groups occupy different covariate ranges that point is an extrapolation for some of them.",
    "Residuals are roughly normal with equal spread across groups; with many rows the F test tolerates departures, the intervals less so.",
    "A difference after adjustment is not a causal effect of the group unless the groups were assigned at random.")
  answer <- list(group_f = tidy(group_F), group_p = p_text(group_p), partial_eta2 = round(partial_eta2, 3), slope = tidy(slope),
                 slope_low = tidy(slope_ci[1]), slope_high = tidy(slope_ci[2]), slopes_p = p_text(slopes_p), groups = k, n = n,
                 highest_adjusted = levs[which.max(adj)], lowest_adjusted = levs[which.min(adj)],
                 adjusted_means = as.list(stats::setNames(tidy(adj), levs)), n_pairs_differing = n_sig)

  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(group_f = tidy(group_F), group_p = p_text(group_p), partial_eta2 = round(partial_eta2, 3),
    slope = tidy(slope), slopes_p = p_text(slopes_p), groups = k, n = n), lead = "group_f", place = "summary_metrics")
  results$adjusted_means <- place_comparison(adj_df, category = "group", value = "adjusted_mean", low = "ci_low", high = "ci_high", place = "adjusted_means")
  results$raw_vs_adjusted <- place_comparison(rva_df, category = "group", value = "mean", series = "kind", place = "raw_vs_adjusted")
  results$covariate_scatter <- place_relationship(scatter, x = cov_h, y = outcome_h, series = group_h, draws = "dataset", place = "covariate_scatter")
  if (!is.null(pair_df)) {
    results$pairwise_interval <- place_interval(pair_df, term = "comparison", value = "adj_difference", low = "ci_low", high = "ci_high", place = "pairwise_interval")
  } else {
    results$pairwise_interval <- place_dropped("with two groups the overall adjusted test is the pairwise comparison; see the adjusted means and the ANCOVA table", place = "pairwise_interval")
  }
  results$ancova_table <- place_table(table_df, place = "ancova_table")
  results$assumption_checks <- place_table(checks_df, place = "assumption_checks")
  results$comparison_method <- place_method(method = method, n_in = n_in, n_used = n, assumptions = assumptions, excluded = as.list(excluded), place = "comparison_method")

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