GLM and Count Regression

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

`standard_glm_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.
  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 the mapper shows as "<0.0001" (LAT-3181).
  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 %||%
    "What drives this outcome, by how much, and which model fits the way the outcome is measured?"

  #' ## Column mapping
  #' One `outcome`, one or more `driver_N` columns (a series, any number), and an optional `exposure` (how long or how
  #' much each row was observed, for counts). Semantic names inside; the customer's own headers in `col_map`.
  col_map <- inputs$column_mapping %||% list()
  df <- renderObject.taskFunction.init(inputs, col_map)
  human <- function(sem) {
    v <- col_map[[sem]]
    if (is.null(v) || !nzchar(as.character(v))) sem else as.character(v)
  }
  n_in <- nrow(df)
  if (!"outcome" %in% names(df)) stop("column_mapping must map an 'outcome' column (the value to explain)")
  driver_cols <- grep("^driver_[0-9]+$", names(df), value = TRUE)
  driver_cols <- driver_cols[order(as.integer(sub("^driver_", "", driver_cols)))]
  if (length(driver_cols) == 0) stop("column_mapping must map at least one driver column (driver_1)")
  outcome_name <- human("outcome")

  family_param <- tolower(trimws(as.character(params$family %||% "auto")))
  valid_families <- c("auto", "gaussian", "binomial", "poisson", "negative_binomial", "gamma")
  if (!family_param %in% valid_families)
    stop(sprintf("module_parameters$family must be one of %s; got '%s'.", paste(valid_families, collapse = ", "), family_param))
  positive_param <- tolower(trimws(as.character(params$positive %||% "")))

  #' ## The outcome and its family
  #' Two distinct values are a yes/no outcome (logistic). Otherwise the outcome must be numeric for 95% of its values:
  #' non-negative whole numbers are counts (Poisson, or negative binomial when the counts vary more than Poisson allows);
  #' strictly positive values with a long right tail (skewness above 1) take a gamma model with a log link; anything else
  #' is gaussian. A named family is honoured when the outcome can carry it and refused with the reason when it cannot.
  raw_chr <- trimws(as.character(df$outcome))
  nonblank <- !is.na(df$outcome) & raw_chr != "" & tolower(raw_chr) != "na"
  if (sum(nonblank) < 30) stop(sprintf("Only %d rows have a value in %s; at least 30 are required.", sum(nonblank), outcome_name))
  lv <- unique(tolower(raw_chr[nonblank]))
  if (length(lv) == 1) stop(sprintf("%s has one value, so there is nothing to explain.", outcome_name))
  positive <- NA_character_
  skew <- function(x) { m <- mean(x); s <- stats::sd(x); if (!is.finite(s) || s == 0) 0 else mean(((x - m) / s)^3) }
  if (family_param == "binomial" || (family_param == "auto" && length(lv) == 2)) {
    if (length(lv) != 2) stop(sprintf("A logistic model needs an outcome with exactly two values; %s has %d.", outcome_name, length(lv)))
    yes_words <- c("1", "true", "yes", "y", "t")
    positive <- if (nzchar(positive_param) && positive_param %in% lv) positive_param
                else if (sum(lv %in% yes_words) == 1) lv[lv %in% yes_words]
                else names(sort(table(tolower(raw_chr[nonblank]))))[1]
    y <- ifelse(nonblank, as.integer(tolower(raw_chr) == positive), NA_integer_)
    family_key <- "binomial"
    positive_display <- raw_chr[nonblank][match(positive, tolower(raw_chr[nonblank]))]
  } else {
    conv <- suppressWarnings(as.numeric(gsub("[$,£€ ]", "", raw_chr)))
    if (sum(nonblank & !is.na(conv)) < 0.95 * sum(nonblank))
      stop(sprintf("%s is not numeric for 95%% of its values and has %d distinct values; this tool needs a number, a count or a yes/no outcome.", outcome_name, length(lv)))
    y <- ifelse(nonblank, conv, NA_real_)
    yv <- y[!is.na(y)]
    if (isTRUE(stats::var(yv) == 0)) stop(sprintf("%s is constant, so there is nothing to explain.", outcome_name))
    countlike <- all(yv >= 0) && all(abs(yv - round(yv)) < 1e-9)
    if (family_param == "auto") {
      family_key <- if (countlike) "count" else if (all(yv > 0) && skew(yv) > 1) "gamma" else "gaussian"
    } else if (family_param %in% c("poisson", "negative_binomial")) {
      if (!countlike) stop(sprintf("A %s model needs non-negative whole numbers; %s has fractions or negatives.", sub("_", " ", family_param), outcome_name))
      family_key <- "count"
    } else if (family_param == "gamma") {
      if (!all(yv > 0)) stop(sprintf("A gamma model needs values above zero; %s has zeros or negatives.", outcome_name))
      family_key <- "gamma"
    } else family_key <- "gaussian"
  }
  df$outcome <- y
  n_blank_outcome <- sum(is.na(df$outcome))
  df <- df[!is.na(df$outcome), , drop = FALSE]

  #' ## Exposure
  #' Counts are modelled per unit of exposure through a log offset. Rows with a blank, zero or negative exposure are
  #' left out and counted; exposure mapped on a model that is not a count is not used, and the method says so.
  has_exposure <- "exposure" %in% names(df)
  n_bad_exposure <- 0L
  exposure_ignored <- FALSE
  if (has_exposure) {
    ex <- suppressWarnings(as.numeric(as.character(df$exposure)))
    if (family_key == "count") {
      bad <- is.na(ex) | ex <= 0
      n_bad_exposure <- sum(bad)
      df <- df[!bad, , drop = FALSE]; ex <- ex[!bad]
      df$.log_exposure <- log(ex)
    } else { exposure_ignored <- TRUE; has_exposure <- FALSE }
  }

  #' ## Drivers
  #' Numeric by the 95% rule, blanks filled with the median and every fill counted (LAT-3210); text becomes a factor with
  #' blanks as "Missing" and levels beyond twelve lumped into "Other"; identifier-like and constant drivers are excluded.
  dropped <- character(0); why <- character(0); filled <- list()
  for (dc in driver_cols) {
    v <- df[[dc]]
    if (!is.numeric(v)) {
      cv <- suppressWarnings(as.numeric(as.character(v)))
      n_orig <- sum(!is.na(v) & as.character(v) != "")
      if (n_orig > 0 && sum(!is.na(cv)) >= 0.95 * n_orig) df[[dc]] <- cv
    }
    v <- df[[dc]]
    if (is.numeric(v)) {
      med <- stats::median(v, na.rm = TRUE)
      if (is.na(med)) { dropped <- c(dropped, dc); why <- c(why, "empty"); next }
      n_na <- sum(is.na(v)); v[is.na(v)] <- med; df[[dc]] <- v
      if (isTRUE(stats::var(v) == 0)) { dropped <- c(dropped, dc); why <- c(why, "constant"); next }
      if (length(unique(v)) == length(v) && all(v == round(v)) && length(v) > 50 && isTRUE(all(diff(sort(v)) == 1))) {
        dropped <- c(dropped, dc); why <- c(why, "identifier-like (a running index)"); next }
      if (n_na > 0) filled[[dc]] <- list(n = n_na, value = med)
    } else {
      v <- as.character(v); v[is.na(v) | trimws(v) == ""] <- "Missing"
      if (length(unique(v)) > nrow(df) / 2) { dropped <- c(dropped, dc); why <- c(why, "identifier-like (almost every row its own value)"); next }
      tab <- sort(table(v), decreasing = TRUE)
      if (length(tab) > 12) v[!(v %in% names(tab)[1:12])] <- "Other"
      if (length(unique(v)) <= 1) { dropped <- c(dropped, dc); why <- c(why, "constant"); next }
      df[[dc]] <- factor(v, levels = names(sort(table(v), decreasing = TRUE)))
    }
  }
  model_drivers <- setdiff(driver_cols, dropped)
  if (length(model_drivers) == 0) stop("No usable driver columns remained after cleaning (all constant, empty, or identifier-like).")
  fmt_num <- function(x) format(signif(x, 6), big.mark = ",", scientific = FALSE, trim = TRUE)
  filled <- filled[names(filled) %in% model_drivers]
  filled_items <- vapply(names(filled), function(dc) paste0(human(dc), " (", filled[[dc]]$n, " cell",
    if (filled[[dc]]$n > 1) "s" else "", ", median ", fmt_num(filled[[dc]]$value), ")"), character(1), USE.NAMES = FALSE)
  filled_txt <- paste(filled_items, collapse = ", ")
  n_used <- nrow(df)
  n_terms <- sum(vapply(model_drivers, function(dc) if (is.factor(df[[dc]])) nlevels(df[[dc]]) - 1 else 1, numeric(1)))
  if (n_used < n_terms + 10) stop(sprintf("Only %d usable rows for %d model terms; at least %d are needed.", n_used, n_terms, n_terms + 10))
  if (family_key == "binomial" && min(table(df$outcome)) < 5)
    stop(sprintf("%s needs at least 5 rows in each of its two values; the rarer has %d.", outcome_name, min(table(df$outcome))))

  #' ## Fit
  model_df <- df[, c(model_drivers, "outcome", if (has_exposure) ".log_exposure"), drop = FALSE]
  fml <- stats::as.formula(paste("outcome ~", paste(model_drivers, collapse = " + "), if (has_exposure) "+ offset(.log_exposure)" else ""))
  dispersion <- NA_real_; disp_p <- NA_real_; theta <- NA_real_; compare_df <- NULL; zero_obs <- NA_real_; zero_exp <- NA_real_
  if (family_key == "binomial") {
    model <- stats::glm(fml, family = stats::binomial(), data = model_df)
    family_label <- "Logistic (binomial)"; effect_label <- "Odds ratio"; exp_scale <- TRUE
  } else if (family_key == "gaussian") {
    model <- stats::glm(fml, family = stats::gaussian(), data = model_df)
    family_label <- "Linear (gaussian)"; effect_label <- "Coefficient"; exp_scale <- FALSE
  } else if (family_key == "gamma") {
    model <- tryCatch(stats::glm(fml, family = stats::Gamma(link = "log"), data = model_df, control = stats::glm.control(maxit = 100)),
                      error = function(e) stop(sprintf("The gamma model did not converge on %s: %s", outcome_name, conditionMessage(e))))
    family_label <- "Gamma (log link)"; effect_label <- "Mean ratio"; exp_scale <- TRUE
    gauss <- stats::glm(fml, family = stats::gaussian(), data = model_df)
    compare_df <- data.frame(model = c("Gamma (log link)", "Linear (gaussian)"),
      aic = tidy(c(stats::AIC(model), stats::AIC(gauss))), dispersion = tidy(c(summary(model)$dispersion, NA)),
      chosen = c("yes", "no"),
      note = c(if (family_param == "gamma") "the family requested" else "strictly positive values with a long right tail",
               "assumes a symmetric spread that does not grow with the mean"), stringsAsFactors = FALSE)
  } else {
    fit_p <- stats::glm(fml, family = stats::poisson(), data = model_df)
    pearson <- sum(stats::residuals(fit_p, type = "pearson")^2)
    dispersion <- pearson / stats::df.residual(fit_p)
    disp_p <- stats::pchisq(pearson, stats::df.residual(fit_p), lower.tail = FALSE)
    fit_nb <- tryCatch(suppressWarnings(MASS::glm.nb(fml, data = model_df)), error = function(e) NULL)
    use_nb <- if (family_param == "poisson") FALSE
              else if (family_param == "negative_binomial") TRUE
              else isTRUE(dispersion > 1.5 && disp_p < 0.05 && !is.null(fit_nb))
    if (use_nb && is.null(fit_nb)) stop(sprintf("The negative binomial model did not converge on %s.", outcome_name))
    model <- if (use_nb) fit_nb else fit_p
    theta <- if (!is.null(fit_nb)) fit_nb$theta else NA_real_
    family_key <- if (use_nb) "negative_binomial" else "poisson"
    family_label <- if (use_nb) "Negative binomial" else "Poisson"
    effect_label <- "Rate ratio"; exp_scale <- TRUE
    compare_df <- data.frame(model = c("Poisson", "Negative binomial"),
      aic = tidy(c(stats::AIC(fit_p), if (!is.null(fit_nb)) stats::AIC(fit_nb) else NA)),
      dispersion = tidy(c(dispersion, if (!is.null(fit_nb)) sum(stats::residuals(fit_nb, type = "pearson")^2) / stats::df.residual(fit_nb) else NA)),
      chosen = c(if (use_nb) "no" else "yes", if (use_nb) "yes" else "no"),
      note = c("assumes the counts vary as much as their mean",
               if (is.null(fit_nb)) "did not converge" else paste0("allows extra spread (theta ", tidy(theta), ")")), stringsAsFactors = FALSE)
    mu <- stats::fitted(model)
    zero_obs <- sum(model_df$outcome == 0)
    zero_exp <- if (use_nb) sum(stats::dnbinom(0, size = theta, mu = mu)) else sum(stats::dpois(0, mu))
  }

  #' ## Terms, effects and importance
  #' A term's driver is matched exactly, longest driver name first (LAT-3210); aliased terms are named and left out.
  cf_all <- stats::coef(model)
  aliased_raw <- names(cf_all)[is.na(cf_all)]
  owners_order <- model_drivers[order(-nchar(model_drivers))]
  term_owner <- function(t) {
    for (dc in owners_order) {
      v <- model_df[[dc]]
      if (if (is.factor(v)) t %in% paste0(dc, levels(v)) else identical(t, dc)) return(dc)
    }
    NA_character_
  }
  human_term <- function(t) {
    if (t == "(Intercept)") return("(Intercept)")
    dc <- term_owner(t); if (is.na(dc)) return(t)
    lvl <- substring(t, nchar(dc) + 1)
    if (nzchar(lvl)) paste0(human(dc), " = ", lvl) else human(dc)
  }
  aliased <- vapply(aliased_raw, human_term, character(1), USE.NAMES = FALSE)
  sm <- summary(model)$coefficients
  terms <- rownames(sm)
  est <- sm[, 1]; se <- sm[, 2]; zv <- sm[, 3]; pv <- sm[, 4]
  lo <- est - stats::qnorm(0.975) * se; hi <- est + stats::qnorm(0.975) * se
  scale_fn <- if (exp_scale) exp else identity
  effects_df <- data.frame(term = vapply(terms, human_term, character(1), USE.NAMES = FALSE),
    effect = tidy(scale_fn(est)), low = tidy(scale_fn(lo)), high = tidy(scale_fn(hi)),
    estimate = tidy(est), std_error = tidy(se), z_value = tidy(zv), p_value = p_cell(pv), stringsAsFactors = FALSE)
  rownames(effects_df) <- NULL
  interval_df <- effects_df[terms != "(Intercept)", c("term", "effect", "low", "high"), drop = FALSE]
  owners <- vapply(terms, term_owner, character(1), USE.NAMES = FALSE)
  imp <- do.call(rbind, lapply(model_drivers, function(dc) {
    idx <- which(owners == dc & is.finite(zv)); if (!length(idx)) return(NULL)
    data.frame(driver = human(dc), abs_z = max(abs(zv[idx])), p = min(pv[idx]), stringsAsFactors = FALSE)
  }))
  imp <- imp[order(-imp$abs_z), , drop = FALSE]
  imp_df <- data.frame(driver = imp$driver, importance = round(100 * imp$abs_z / max(imp$abs_z), 1), stringsAsFactors = FALSE)
  n_sig <- sum(imp$p < 0.05)

  #' ## Fit, predictions and the outcome's shape
  dev_expl <- if (is.finite(model$null.deviance) && model$null.deviance > 0) 1 - model$deviance / model$null.deviance else NA_real_
  pred <- as.numeric(stats::fitted(model))
  set.seed(42)
  if (family_key == "binomial") {
    bins <- cut(rank(pred, ties.method = "first"), breaks = 10, labels = FALSE)
    pva_df <- data.frame(predicted = tidy(tapply(pred, bins, mean)), actual = tidy(tapply(model_df$outcome, bins, mean)), stringsAsFactors = FALSE)
    dist_df <- data.frame(value = c(positive_display, paste("not", positive_display)),
                          rows = as.integer(c(sum(model_df$outcome == 1), sum(model_df$outcome == 0))), stringsAsFactors = FALSE)
  } else {
    idx <- if (n_used > 1500) sort(sample(n_used, 1500)) else seq_len(n_used)
    pva_df <- data.frame(predicted = tidy(pred[idx]), actual = tidy(model_df$outcome[idx]), stringsAsFactors = FALSE)
    if (family_key %in% c("poisson", "negative_binomial")) {
      top <- 10
      lab <- ifelse(model_df$outcome >= top, paste0(top, " or more"), as.character(model_df$outcome))
      lvls <- c(as.character(0:(top - 1)), paste0(top, " or more"))
      dist_df <- data.frame(value = lvls, rows = as.integer(table(factor(lab, levels = lvls))), stringsAsFactors = FALSE)
    } else {
      br <- pretty(model_df$outcome, n = 10)
      cutv <- cut(model_df$outcome, breaks = br, include.lowest = TRUE)
      lab <- paste(format(br[-length(br)], big.mark = ",", trim = TRUE), "to", format(br[-1], big.mark = ",", trim = TRUE))
      dist_df <- data.frame(value = lab, rows = as.integer(table(cutv)), stringsAsFactors = FALSE)
    }
  }

  #' ## Assumption checks (LAT-3138)
  vif_by <- tryCatch({
    X <- stats::model.matrix(model)
    X <- X[, setdiff(colnames(X), c("(Intercept)", aliased_raw)), drop = FALSE]
    if (ncol(X) < 2) NULL else {
      v <- vapply(seq_len(ncol(X)), function(j) {
        r2 <- summary(stats::lm.fit(cbind(1, X[, -j, drop = FALSE]), X[, j]) |> (\(f) list(r.squared = 1 - sum(f$residuals^2) / sum((X[, j] - mean(X[, j]))^2)))())$r.squared
        if (!is.finite(r2) || r2 >= 1) Inf else 1 / (1 - r2) }, numeric(1))
      own <- vapply(colnames(X), function(c) { o <- term_owner(c); if (is.na(o)) c else o }, character(1))
      tapply(v, own, max)
    }
  }, error = function(e) NULL)
  max_vif <- if (length(vif_by)) max(vif_by) else NA_real_
  max_vif_driver <- if (length(vif_by)) human(names(vif_by)[which.max(vif_by)]) else ""
  cooks <- tryCatch(stats::cooks.distance(model), error = function(e) rep(NA_real_, n_used))
  share_infl <- mean(cooks > 4 / n_used, na.rm = TRUE)
  rows_per_term <- n_used / max(1, n_terms)
  fill_share <- if (length(filled)) max(vapply(filled, function(f) f$n, numeric(1))) / n_used else 0
  chk <- list()
  add_chk <- function(check, statistic, p, verdict, note) chk[[length(chk) + 1]] <<- data.frame(check = check, statistic = statistic, p_value = p, verdict = verdict, note = note, stringsAsFactors = FALSE)
  add_chk("Enough rows per term", paste0(round(rows_per_term, 1), " rows per term"), "",
          if (rows_per_term >= 10) "holds" else if (rows_per_term >= 5) "strained" else "violated",
          "fewer than ten rows per term makes the effects unstable")
  if (family_key %in% c("poisson", "negative_binomial")) {
    add_chk("Spread of the counts", paste0("Pearson dispersion ", tidy(dispersion), " under Poisson"), format(signif(p_cell(disp_p), 3)),
            if (family_key == "negative_binomial" || dispersion <= 1.5) "holds" else "violated",
            if (family_key == "negative_binomial") "the counts vary more than Poisson allows; the negative binomial model absorbs the extra spread"
            else "counts that vary more than their mean understate every interval under a Poisson model")
    zr <- if (zero_exp > 0) zero_obs / zero_exp else NA_real_
    add_chk("Zeros as expected", paste0(zero_obs, " zeros observed against ", tidy(zero_exp), " expected"), "",
            if (!is.finite(zr) || zr <= 1.25) "holds" else if (zr <= 1.6) "strained" else "violated",
            "far more zeros than the model expects points to a separate group that never has an event")
  }
  if (family_key == "binomial") {
    epv <- min(table(model_df$outcome)) / max(1, n_terms)
    sep <- any(abs(est[terms != "(Intercept)"]) > 10 | se[terms != "(Intercept)"] > 50, na.rm = TRUE)
    add_chk("Enough events per term", paste0(round(epv, 1), " rarer-class rows per term"), "",
            if (epv >= 10) "holds" else if (epv >= 5) "strained" else "violated", "few events per term makes the odds ratios unstable")
    add_chk("No separation", if (sep) "a term's estimate runs off to infinity" else "every term estimated", "",
            if (sep) "violated" else "holds", "a driver that splits the two values perfectly gives an odds ratio with no finite value")
  }
  if (family_key == "gaussian") {
    rs <- stats::residuals(model)
    add_chk("Residuals roughly symmetric", paste0("skewness ", tidy(skew(rs))), "",
            if (abs(skew(rs)) < 1) "holds" else if (abs(skew(rs)) < 2) "strained" else "violated",
            "a long tail in the residuals makes the intervals too narrow on one side")
  }
  add_chk("No strong collinearity",
          paste0(if (is.na(max_vif)) "VIF not assessed" else paste0("max VIF = ", tidy(max_vif), " (", max_vif_driver, ")"),
                 if (length(aliased)) paste0("; ", length(aliased), " term", if (length(aliased) > 1) "s" else "", " exactly collinear and dropped") else ""), "",
          if (!is.na(max_vif) && max_vif > 10) "violated" else if ((!is.na(max_vif) && max_vif > 5) || length(aliased)) "strained" else "holds",
          paste0("a VIF above 5 inflates that driver's interval", if (length(aliased)) paste0("; dropped because other terms determine them: ", paste(aliased, collapse = ", ")) else ""))
  add_chk("No dominant rows", paste0(round(100 * share_infl, 1), "% of rows over Cook's 4/n"), "",
          if (share_infl <= 0.05) "holds" else if (share_infl <= 0.1) "strained" else "violated",
          "a few rows steering the fit move the effects on their own")
  add_chk("Few filled cells", if (!length(filled)) "no blank numeric cells" else paste0(sum(vapply(filled, function(f) f$n, numeric(1))), " blank cells set to the median"), "",
          if (fill_share < 0.01) "holds" else if (fill_share < 0.05) "strained" else "violated",
          if (length(filled)) paste0(filled_txt, "; a filled row sits at the middle of that driver") else "a filled row sits at the middle of that driver")
  checks_df <- do.call(rbind, chk)

  #' ## Method and answer
  excluded <- c(if (n_blank_outcome > 0) sprintf("%d row%s with no %s", n_blank_outcome, if (n_blank_outcome > 1) "s" else "", outcome_name),
                if (n_bad_exposure > 0) sprintf("%d row%s with a blank, zero or negative %s", n_bad_exposure, if (n_bad_exposure > 1) "s" else "", human("exposure")),
                if (length(dropped)) paste0(vapply(dropped, human, character(1)), " (", why, ")"))
  family_why <- switch(family_key,
    binomial = paste0(outcome_name, " has two values, so a logistic model gives the odds that it is '", positive_display, "'"),
    gaussian = paste0(outcome_name, " is a number that is not a count and not strictly positive with a long tail, so a linear model applies"),
    gamma = paste0(outcome_name, " is strictly positive with a long right tail, so a gamma model with a log link gives mean ratios"),
    poisson = paste0(outcome_name, " is a count whose spread matches its mean (dispersion ", tidy(dispersion), "), so a Poisson model applies"),
    negative_binomial = paste0(outcome_name, " is a count that varies more than Poisson allows (dispersion ", tidy(dispersion), "), so a negative binomial model applies"))
  if (family_param != "auto") family_why <- paste0(family_why, " (the family requested)")
  method <- paste0(family_label, " regression of ", outcome_name, " on ", length(model_drivers), " driver", if (length(model_drivers) > 1) "s" else "",
    " over ", n_used, " rows: ", family_why,
    if (has_exposure) paste0("; counts are modelled per unit of ", human("exposure"), " through a log offset") else "",
    if (exposure_ignored) paste0("; ", human("exposure"), " was mapped but only counts use an exposure, so it was not used") else "",
    "; effects are ", tolower(effect_label), "s with 95% Wald intervals; a driver's importance is its largest |z| scaled to 100",
    if (length(filled)) paste0("; blank numeric cells filled with the column median: ", filled_txt) else "",
    if (length(excluded)) paste0("; excluded: ", paste(excluded, collapse = "; ")) else "", ".")
  assumptions <- Filter(Negate(is.null), list(
    "Rows are independent; repeated rows per person or place make the intervals too narrow.",
    if (exp_scale) "Each numeric driver changes the outcome by a constant ratio per unit; a curved effect is understated." else "Each numeric driver changes the outcome by a constant amount per unit; a curved effect is understated.",
    "These are associations in the data, not effects of changing a driver.",
    if (has_exposure) paste0("Counts are proportional to ", human("exposure"), ": twice the exposure, twice the expected count.") else NULL,
    if (length(filled)) paste0("Blank numeric cells were set to the column median: ", filled_txt, ".") else NULL))
  answer <- list(family = family_key, n = n_used, deviance_explained = if (is.na(dev_expl)) NULL else round(dev_expl, 3),
                 top_driver = if (nrow(imp_df)) imp_df$driver[1] else NULL, significant_drivers = n_sig)

  results <- list()
  summary_vals <- list(rows = n_used, deviance_explained_pct = if (is.na(dev_expl)) NA_real_ else round(100 * dev_expl, 1),
                       aic = tidy(stats::AIC(model)), significant_drivers = as.integer(n_sig))
  if (family_key %in% c("poisson", "negative_binomial")) summary_vals$dispersion <- tidy(dispersion)
  results$summary_metrics <- place_metric(summary_vals, lead = "deviance_explained_pct", place = "summary_metrics")
  results$effect_interval <- place_interval(interval_df, term = "term", value = "effect", low = "low", high = "high", place = "effect_interval")
  results$driver_importance <- place_comparison(imp_df, category = "driver", value = "importance", place = "driver_importance")
  results$predicted_vs_actual <- place_relationship(pva_df, x = "predicted", y = "actual", place = "predicted_vs_actual")
  results$outcome_distribution <- place_comparison(dist_df, category = "value", value = "rows", place = "outcome_distribution")
  if (!is.null(compare_df)) {
    results$family_comparison <- place_table(compare_df, place = "family_comparison")
  } else {
    results$family_comparison <- place_dropped(if (family_key == "binomial") "a yes/no outcome has one model family, the logistic model, so there is nothing to compare"
                                               else "the outcome is not a count and not strictly positive, so the linear model is the only family that fits it", place = "family_comparison")
  }
  results$effects_table <- place_table(effects_df, place = "effects_table")
  results$assumption_checks <- place_table(checks_df, place = "assumption_checks")
  results$glm_method <- list(kind = "metric", values = list(
    method = method, n_in = n_in, n_used = n_used, excluded = as.list(excluded),
    imputed = if (length(filled)) as.list(filled_items) else list(), assumptions = assumptions,
    x_column = paste(vapply(model_drivers, human, character(1)), collapse = ", "), y_column = outcome_name),
    value_order = list("n_used", "n_in"))

  objects <- list()
  list(answer = answer, method = method, n = n_used, results = results, objects = objects,
       json_output = list(answer = answer, method = method, n = n_used))
}