Classification

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

`standard_classification_v2` <- function(pf) {
  `%||%` <- function(a, b) if (!is.null(a)) a else b
  #' Readable figures (LAT-3181, as LAT-3180): 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)))))
  }
  inputs <- pf$taskList$inputs
  #' P-VALUES BELOW 0.0001 LEAVE AS A FLOOR (LAT-3181). The results serialize at four decimal digits, which rounds a p of
  #' 1.4e-05 to 0 and 6e-05 up to 0.0001; a table printed "p-value 0". 1e-12 survives the serializer and the mapper shows
  #' any p under 0.0001 as "<0.0001", which is what the reader should see. Table and tile cells only; the answer keeps p.
  p_cell <- function(p) { p <- as.numeric(p); ifelse(is.na(p), NA_real_, ifelse(p < 1e-4, 1e-12, signif(p, 3))) }
  params <- inputs$module_parameters %||% list()
  # THE QUESTION this tool answers: the customer's objective, verbatim, when given.
  question <- (inputs$userContext %||% list())$objective %||%
    "Which drivers predict the outcome class, and how well does the classifier separate the classes on data it has not seen?"

  #' ## Column mapping
  #' The customer maps one binary `outcome` and one to eight `driver_N` columns.
  #' Semantic names are used inside; the customer's own headers are carried in
  #' `col_map` so every table, term and axis names the columns the reader knows.
  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)
  }

  #' ## The outcome
  #' Two classes. A numeric 0/1 column is read as is (1 = positive); a two-level
  #' text column takes the LESS frequent level as positive unless
  #' `module_parameters$positive` names one; yes/no, true/false, y/n are read.
  #' Rows without an outcome are dropped. More than two levels is a refusal.
  n_in <- nrow(df)
  if (!"outcome" %in% names(df)) stop("column_mapping must map an 'outcome' column (the two-class outcome to predict)")
  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)")
  raw <- df$outcome
  raw_chr <- trimws(tolower(as.character(raw)))
  raw_chr[is.na(raw) | raw_chr == "" | raw_chr == "na"] <- NA
  df <- df[!is.na(raw_chr), , drop = FALSE]; raw_chr <- raw_chr[!is.na(raw_chr)]
  lv <- sort(unique(raw_chr))
  #' TWO CLASSES OR MANY (LAT-3077). A price tier, a grade, a segment: three or more
  #' classes is an ordinary classification question and refusing it sent a correctly
  #' mapped dataset away. Above two we fit ONE-VS-REST: a logistic model per class,
  #' an AUC per class with a macro average, K curves on one ROC and a K by K confusion
  #' table at the arg-max class. The binary path below is unchanged.
  MAX_CLASSES <- 10
  if (length(lv) < 2) stop(sprintf("The outcome '%s' has one value; a classifier needs at least two classes.", human("outcome")))
  if (length(lv) > MAX_CLASSES) stop(sprintf("The outcome '%s' has %d distinct values; a classifier handles up to %d classes. If this is a number to predict rather than a class, use the regression tool; if it is an identifier, map a different column.", human("outcome"), length(lv), MAX_CLASSES))
  multiclass <- length(lv) > 2
  positive_param <- tolower(trimws(as.character(params$positive %||% "")))
  yes_words <- c("1", "true", "yes", "y", "t")
  if (multiclass) {
    #' The NAMED class leads the report (the parameter, else the rarest, which is the
    #' one a reader usually cares about); every class still gets its own AUC and row.
    positive <- if (nzchar(positive_param) && positive_param %in% lv) positive_param else names(sort(table(raw_chr)))[1]
    negative <- paste0("not ", positive)
    y <- as.integer(raw_chr == positive)
  } else {
    positive <- if (nzchar(positive_param) && positive_param %in% lv) positive_param
                else if (any(lv %in% yes_words)) lv[lv %in% yes_words][1]
                else names(sort(table(raw_chr)))[1]   # the rarer class
    negative <- setdiff(lv, positive)
    y <- as.integer(raw_chr == positive)
  }
  df$outcome <- y
  df$.class <- raw_chr

  #' ## Drivers
  #' Numeric drivers are median-imputed; text drivers become factors with blanks as
  #' "Missing" and levels beyond twelve lumped into "Other"; identifier-like and
  #' constant drivers are excluded.
  dropped <- character(0); why <- character(0)
  for (dc in driver_cols) {
    v <- df[[dc]]
    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) df[[dc]] <- conv
    }
    v <- df[[dc]]
    if (is.numeric(v)) {
      med <- median(v, na.rm = TRUE)
      if (is.na(med)) { dropped <- c(dropped, dc); why <- c(why, "empty"); next }
      v[is.na(v)] <- med; df[[dc]] <- v
      if (isTRUE(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 }
    } else {
      v <- as.character(v); v[is.na(v) | trimws(v) == ""] <- "Missing"
      tab <- sort(table(v), decreasing = TRUE)
      if (length(tab) > 12) v[!(v %in% names(tab)[1:12])] <- "Other"
      if (length(unique(v)) > nrow(df) / 2) { dropped <- c(dropped, dc); why <- c(why, "identifier-like (almost every row its own value)"); next }
      if (length(unique(v)) <= 1) { dropped <- c(dropped, dc); why <- c(why, "constant"); next }
      df[[dc]] <- factor(v)
    }
  }
  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).")
  n_used <- nrow(df)
  n_pos <- sum(y); n_neg <- n_used - n_pos
  if (n_used < 30) stop(sprintf("Only %d rows carry an outcome; at least 30 are required.", n_used))
  if (min(n_pos, n_neg) < 10) stop(sprintf("The rarer class ('%s') has only %d rows; at least 10 are required.", positive, min(n_pos, n_neg)))

  #' ## Held-out split
  #' A stratified 70/30 split (seed 42) when there are at least 100 rows and 20 of the
  #' rarer class; the model is fitted on the training part and EVERY performance
  #' figure below is measured on the test part. Below that size the model is fitted
  #' on all rows and the figures are in-sample, which the method says plainly.
  set.seed(42)
  holdout <- n_used >= 100 && min(n_pos, n_neg) >= 20
  if (holdout) {
    idx_pos <- which(y == 1); idx_neg <- which(y == 0)
    test_idx <- c(sample(idx_pos, round(0.3 * length(idx_pos))), sample(idx_neg, round(0.3 * length(idx_neg))))
    train <- df[-test_idx, c(model_drivers, "outcome"), drop = FALSE]
    test  <- df[ test_idx, c(model_drivers, "outcome"), drop = FALSE]
    # a factor level seen only in the test part cannot be scored: those rows are set aside and counted
    unseen <- rep(FALSE, nrow(test))
    for (dc in model_drivers) if (is.factor(train[[dc]])) {
      train[[dc]] <- droplevels(train[[dc]])
      unseen <- unseen | !(as.character(test[[dc]]) %in% levels(train[[dc]]))
      test[[dc]] <- factor(as.character(test[[dc]]), levels = levels(train[[dc]]))
    }
    n_unseen <- sum(unseen); test <- test[!unseen, , drop = FALSE]
    cls_train <- df$.class[-test_idx]; cls_test <- df$.class[test_idx][!unseen]
  } else {
    train <- df[, c(model_drivers, "outcome"), drop = FALSE]; test <- train; n_unseen <- 0L
    cls_train <- df$.class; cls_test <- df$.class
  }
  n_train <- nrow(train); n_test <- nrow(test)

  #' ## Logistic regression
  #' `glm(outcome ~ ., binomial)` on the training part. Coefficients are reported as
  #' odds ratios with 95% Wald intervals; a driver's importance is the largest |z|
  #' across its terms, scaled so the strongest driver reads 100.
  model <- suppressWarnings(glm(outcome ~ ., data = train, family = binomial()))
  sm <- summary(model); co <- sm$coefficients
  terms_raw <- rownames(co)
  human_term <- function(t) {
    if (t == "(Intercept)") return("(Intercept)")
    for (dc in model_drivers) if (startsWith(t, dc)) {
      lvl <- substring(t, nchar(dc) + 1)
      return(if (nzchar(lvl)) paste0(human(dc), " = ", lvl) else human(dc))
    }
    t
  }
  est <- unname(co[, "Estimate"]); se <- unname(co[, "Std. Error"])
  #' PER-SD EFFECT (LAT-3181, the logistic twin of LAT-3180's beta): the change in log odds for one standard deviation
  #' of the term's model-matrix column, so terms in different units compare on one axis, with 0 as no effect. The
  #' interval chart plotted raw odds ratios: per-month and per-dollar terms shrank to dots beside plan levels, and the
  #' chart's reference line sat at 0 while an odds ratio's no-effect value is 1.
  Xm <- model.matrix(model); col_sd <- apply(Xm[, -1, drop = FALSE], 2, sd)
  term_sd <- vapply(terms_raw, function(t) if (t %in% names(col_sd)) unname(col_sd[t]) else NA_real_, numeric(1))
  coef_df <- data.frame(
    term       = vapply(terms_raw, human_term, character(1), USE.NAMES = FALSE),
    odds_ratio = tidy(exp(est)),
    log_odds   = tidy(est),
    per_sd     = round(est * unname(term_sd), 3),
    std_error  = tidy(se),
    z_value    = round(unname(co[, "z value"]), 2),
    p_value    = p_cell(unname(co[, "Pr(>|z|)"])),
    low        = tidy(exp(est - 1.96 * se)),
    high       = tidy(exp(est + 1.96 * se)),
    stringsAsFactors = FALSE)
  rownames(coef_df) <- NULL
  slope_df <- coef_df[coef_df$term != "(Intercept)" & is.finite(coef_df$odds_ratio) & coef_df$high < 1e6, , drop = FALSE]
  imp_rows <- lapply(model_drivers, function(dc) {
    idx <- which(startsWith(terms_raw, dc) & terms_raw != "(Intercept)")
    idx <- idx[is.finite(co[idx, "z value"])]
    if (length(idx) == 0) return(NULL)
    best <- idx[which.max(abs(co[idx, "z value"]))]
    data.frame(driver = human(dc), abs_z = abs(unname(co[best, "z value"])),
               direction = if (is.numeric(train[[dc]])) { if (co[best, "Estimate"] > 0) "raises the odds" else "lowers the odds" } else "categorical",
               p_value = p_cell(unname(co[best, "Pr(>|z|)"])), stringsAsFactors = FALSE)
  })
  imp_df <- do.call(rbind, imp_rows)
  imp_df <- imp_df[order(-imp_df$abs_z), , drop = FALSE]
  imp_df$importance <- round(100 * imp_df$abs_z / max(imp_df$abs_z), 1)
  imp_df$abs_z <- NULL; rownames(imp_df) <- NULL
  top_driver <- imp_df$driver[1]; top_dir <- imp_df$direction[1]

  #' ## One-vs-rest, when there are more than two classes
  #' The model above is already the NAMED class against the rest; here every other
  #' class gets the same treatment on the same split, so each class has an AUC and a
  #' curve, and a row is assigned to whichever class scores it highest.
  ovr <- NULL
  if (multiclass) {
    ovr <- list()
    for (k in lv) {
      tr_k <- train; tr_k$outcome <- as.integer(cls_train == k)
      if (sum(tr_k$outcome) < 5 || sum(tr_k$outcome) == nrow(tr_k)) next
      m_k <- suppressWarnings(glm(outcome ~ ., data = tr_k, family = binomial()))
      ovr[[k]] <- list(model = m_k,
                       p_test = suppressWarnings(as.numeric(predict(m_k, newdata = test, type = "response"))),
                       y_test = as.integer(cls_test == k))
    }
    if (length(ovr) < 2) stop("Fewer than two classes had enough rows in the training part to fit a model; the outcome may be too sparse to classify.")
  }

  #' ## Performance on the held-out part
  #' AUC by the rank (Mann-Whitney) formula with a Hanley-McNeil 95% interval; the
  #' ROC curve at every distinct score; the threshold is chosen on the TRAINING
  #' part (Youden's J) and applied unchanged to the test part; the confusion
  #' counts, sensitivity, specificity, precision, accuracy, F1 and the Brier score
  #' are all test-part figures; calibration compares mean predicted probability with
  #' the observed rate in ten equal-count bins of the test part.
  p_train <- suppressWarnings(as.numeric(predict(model, newdata = train, type = "response")))
  p_test  <- suppressWarnings(as.numeric(predict(model, newdata = test,  type = "response")))
  y_test  <- test$outcome
  auc_of <- function(p, yy) {
    r <- rank(p); np <- sum(yy == 1); nn <- sum(yy == 0)
    if (np == 0 || nn == 0) return(NA_real_)
    (sum(r[yy == 1]) - np * (np + 1) / 2) / (np * nn)
  }
  auc <- auc_of(p_test, y_test)
  np_t <- sum(y_test == 1); nn_t <- sum(y_test == 0)
  q1 <- auc / (2 - auc); q2 <- 2 * auc^2 / (1 + auc)
  auc_se <- sqrt((auc * (1 - auc) + (np_t - 1) * (q1 - auc^2) + (nn_t - 1) * (q2 - auc^2)) / (np_t * nn_t))
  auc_low <- max(0, auc - 1.96 * auc_se); auc_high <- min(1, auc + 1.96 * auc_se)
  roc_at <- function(p, yy, thr) {
    pred <- p >= thr
    c(tpr = sum(pred & yy == 1) / sum(yy == 1), fpr = sum(pred & yy == 0) / sum(yy == 0))
  }
  thr_grid <- sort(unique(c(0, p_train, 1)))
  if (length(thr_grid) > 400) thr_grid <- unique(quantile(thr_grid, probs = seq(0, 1, length.out = 400), names = FALSE))
  j_vals <- vapply(thr_grid, function(t) { rr <- roc_at(p_train, train$outcome, t); rr[["tpr"]] - rr[["fpr"]] }, numeric(1))
  threshold <- thr_grid[which.max(j_vals)]
  roc_grid <- sort(unique(c(0, p_test, 1)))
  if (length(roc_grid) > 300) roc_grid <- unique(quantile(roc_grid, probs = seq(0, 1, length.out = 300), names = FALSE))
  roc_df <- do.call(rbind, lapply(rev(roc_grid), function(t) { rr <- roc_at(p_test, y_test, t); data.frame(fpr = round(rr[["fpr"]], 3), tpr = round(rr[["tpr"]], 3)) }))
  roc_df <- roc_df[!duplicated(roc_df), , drop = FALSE]; rownames(roc_df) <- NULL
  # LAT-3146: the ROC frame carries `class` on both paths, so its columns do not depend on the number of
  # classes. With two classes it names the positive class and the chart still draws one line (no series).
  roc_df$class <- as.character(positive)
  pred_pos <- p_test >= threshold
  tp <- sum(pred_pos & y_test == 1); fp <- sum(pred_pos & y_test == 0)
  fn <- sum(!pred_pos & y_test == 1); tn <- sum(!pred_pos & y_test == 0)
  sens <- tp / max(1, tp + fn); spec <- tn / max(1, tn + fp); prec <- tp / max(1, tp + fp)
  acc <- (tp + tn) / n_test; f1 <- if (prec + sens > 0) 2 * prec * sens / (prec + sens) else 0
  brier <- mean((p_test - y_test)^2)
  base_rate <- mean(y_test)
  # LAT-3146: one row per actual and predicted class, the same three columns on both paths; with two
  # classes that is caught, missed, wrongly flagged and correctly cleared, in that order.
  lab <- c(paste0(human("outcome"), " = ", positive), paste0(human("outcome"), " = ", negative))
  confusion_df <- data.frame(
    actual = rep(lab, each = 2), predicted = rep(lab, times = 2),
    n = c(tp, fn, fp, tn), stringsAsFactors = FALSE)
  perf_df <- data.frame(
    metric = c("AUC (held-out)", "AUC 95% low", "AUC 95% high", "threshold", "sensitivity (recall)", "specificity", "precision", "accuracy", "F1", "Brier score", "base rate", "training rows", "test rows"),
    value = round(c(auc, auc_low, auc_high, threshold, sens, spec, prec, acc, f1, brier, base_rate, n_train, n_test), 4),
    stringsAsFactors = FALSE)
  bins <- cut(rank(p_test, ties.method = "first"), breaks = 10, labels = FALSE)
  cal_df <- do.call(rbind, lapply(sort(unique(bins)), function(b) data.frame(
    predicted_rate = round(mean(p_test[bins == b]), 3), observed_rate = round(mean(y_test[bins == b]), 3), n = sum(bins == b))))
  rownames(cal_df) <- NULL

  #' ## Multi-class: the same figures, per class and averaged
  #' Every figure above belongs to the named class against the rest. With more than
  #' two classes the report leads with the MACRO average AUC (the mean of the
  #' per-class AUCs, each class weighted the same however rare it is), the ROC carries
  #' one curve per class, the confusion table is K by K at the arg-max class, and
  #' accuracy is the share of held-out rows whose highest-scoring class was the right
  #' one. The interval on the macro AUC treats the per-class intervals as independent,
  #' which the method text says.
  auc_by_class <- NULL; macro_auc <- NA_real_; argmax_acc <- NA_real_
  if (multiclass) {
    cls <- names(ovr)
    per <- lapply(cls, function(k) {
      a <- auc_of(ovr[[k]]$p_test, ovr[[k]]$y_test)
      npk <- sum(ovr[[k]]$y_test == 1); nnk <- sum(ovr[[k]]$y_test == 0)
      if (is.na(a) || npk == 0 || nnk == 0) return(data.frame(class = k, auc = NA_real_, se = NA_real_, n = npk, stringsAsFactors = FALSE))
      qq1 <- a / (2 - a); qq2 <- 2 * a^2 / (1 + a)
      se <- sqrt((a * (1 - a) + (npk - 1) * (qq1 - a^2) + (nnk - 1) * (qq2 - a^2)) / (npk * nnk))
      data.frame(class = k, auc = a, se = se, n = npk, stringsAsFactors = FALSE)
    })
    auc_by_class <- do.call(rbind, per)
    ok <- !is.na(auc_by_class$auc)
    macro_auc <- mean(auc_by_class$auc[ok])
    macro_se <- sqrt(sum(auc_by_class$se[ok]^2, na.rm = TRUE)) / sum(ok)
    auc <- macro_auc
    auc_low <- max(0, macro_auc - 1.96 * macro_se); auc_high <- min(1, macro_auc + 1.96 * macro_se)
    # arg-max assignment across the per-class scores
    P <- do.call(cbind, lapply(cls, function(k) ovr[[k]]$p_test))
    colnames(P) <- cls
    pred_cls <- cls[max.col(P, ties.method = "first")]
    argmax_acc <- mean(pred_cls == cls_test)
    acc <- argmax_acc
    # K by K confusion as one row per actual and predicted class (LAT-3146): the columns are the same as
    # the two-class table's, instead of one column named after each class.
    confusion_df <- do.call(rbind, lapply(cls, function(a_k) do.call(rbind, lapply(cls, function(p_k)
      data.frame(actual = paste0(human("outcome"), " = ", a_k), predicted = paste0(human("outcome"), " = ", p_k),
                 n = sum(cls_test == a_k & pred_cls == p_k), stringsAsFactors = FALSE)))))
    rownames(confusion_df) <- NULL
    # one ROC curve per class, the class as the series
    roc_df <- do.call(rbind, lapply(cls, function(k) {
      pk <- ovr[[k]]$p_test; yk <- ovr[[k]]$y_test
      g <- sort(unique(c(0, pk, 1)))
      if (length(g) > 120) g <- unique(quantile(g, probs = seq(0, 1, length.out = 120), names = FALSE))
      d <- do.call(rbind, lapply(rev(g), function(t) { rr <- roc_at(pk, yk, t); data.frame(fpr = round(rr[["fpr"]], 3), tpr = round(rr[["tpr"]], 3), class = k, stringsAsFactors = FALSE) }))
      d[!duplicated(d[, c("fpr", "tpr")]), , drop = FALSE]
    }))
    rownames(roc_df) <- NULL
    perf_df <- rbind(
      data.frame(metric = "macro-average AUC (held-out)", value = round(macro_auc, 4), stringsAsFactors = FALSE),
      data.frame(metric = "macro AUC 95% low", value = round(auc_low, 4), stringsAsFactors = FALSE),
      data.frame(metric = "macro AUC 95% high", value = round(auc_high, 4), stringsAsFactors = FALSE),
      data.frame(metric = paste0("AUC, ", auc_by_class$class, " vs rest"), value = round(auc_by_class$auc, 4), stringsAsFactors = FALSE),
      data.frame(metric = "accuracy (arg-max class)", value = round(argmax_acc, 4), stringsAsFactors = FALSE),
      data.frame(metric = "majority-class accuracy", value = round(max(table(cls_test)) / length(cls_test), 4), stringsAsFactors = FALSE),
      data.frame(metric = paste0("sensitivity, ", positive, " vs rest"), value = round(sens, 4), stringsAsFactors = FALSE),
      data.frame(metric = paste0("specificity, ", positive, " vs rest"), value = round(spec, 4), stringsAsFactors = FALSE),
      data.frame(metric = paste0("threshold, ", positive, " vs rest"), value = round(threshold, 4), stringsAsFactors = FALSE),
      data.frame(metric = "classes", value = length(cls), stringsAsFactors = FALSE),
      data.frame(metric = "training rows", value = n_train, stringsAsFactors = FALSE),
      data.frame(metric = "test rows", value = n_test, stringsAsFactors = FALSE))
    rownames(perf_df) <- NULL
  }
  perf_df$value <- tidy(perf_df$value)   # LAT-3181: 0.7449 reads 0.745; row counts stay whole
  band <- if (is.na(auc)) "not measurable" else if (auc < 0.6) "little better than chance" else if (auc < 0.7) "weak" else if (auc < 0.8) "fair" else if (auc < 0.9) "good" else "excellent"

  #' ## Method text and the answer
  method <- paste0(
    if (multiclass) paste0("One-vs-rest logistic regression of ", human("outcome"), " across ", length(lv), " classes (",
                           paste(lv, collapse = ", "), "), one model per class, reported for ", positive, " against the rest, on ")
    else paste0("Logistic regression of ", human("outcome"), " (positive class: ", positive, ") on "),
    paste(vapply(model_drivers, human, character(1)), collapse = ", "), "; ", n_used, " rows with an outcome of ", n_in,
    if (length(dropped)) paste0("; drivers excluded: ", paste(paste0(vapply(dropped, human, character(1)), " (", why, ")"), collapse = ", ")) else "",
    ". ", if (holdout) paste0("Stratified 70/30 split (seed 42): fitted on ", n_train, " rows, every performance figure measured on the ", n_test, " held-out rows",
                              if (n_unseen > 0) paste0(" (", n_unseen, " test rows set aside for a category unseen in training)") else "", ".")
          else paste0("Too few rows for a held-out split; fitted on all ", n_train, " rows and the performance figures are IN-SAMPLE, so they flatter the model."),
    if (multiclass) paste0(" AUC by the rank formula per class with a Hanley-McNeil 95% interval each; the headline is the MACRO average (every class weighted the same), its interval formed by treating the per-class intervals as independent; accuracy is the share of held-out rows assigned to the right class by the highest score, against a majority-class baseline in the table; the confusion table and the ", positive, " threshold ", round(threshold, 3), " (Youden's J on the training part) describe the named class; ")
    else paste0(" AUC by the rank formula with a Hanley-McNeil 95% interval; threshold ", round(threshold, 3), " chosen on the training part by Youden's J and applied unchanged; "),
    "odds ratios with 95% Wald intervals; driver importance = max |z| per driver scaled to 100; calibration in ten equal-count bins.")
  statement <- if (multiclass)
    sprintf("Across %d classes the classifier reaches a macro-average AUC of %s on held-out data (%s), assigning %s%% of rows to the right class; %s %s most for %s",
            length(lv), format(round(auc, 3), nsmall = 3), band, format(round(100 * argmax_acc, 1), nsmall = 1), top_driver, top_dir, positive)
  else
    sprintf("The classifier separates %s = %s from %s with AUC %s on held-out data (%s); %s %s most",
            human("outcome"), positive, negative, format(round(auc, 3), nsmall = 3), band, top_driver, top_dir)
  answer <- list(
    auc = round(auc, 4), auc_low = round(auc_low, 4), auc_high = round(auc_high, 4), auc_band = band,
    threshold = round(threshold, 4), sensitivity = round(sens, 4), specificity = round(spec, 4), precision = round(prec, 4),
    accuracy = round(acc, 4), f1 = round(f1, 4), brier = round(brier, 4), base_rate = round(base_rate, 4),
    positive_class = positive, top_driver = top_driver, top_driver_direction = top_dir,
    n_classes = length(lv), multiclass = multiclass,
    n_drivers = length(model_drivers), n_train = n_train, n_test = n_test, held_out = holdout, n = n_used)

  # ── Results: one entry per place. `columns` carries column order past jsonb's key
  #    sort; `value_order` says which number leads a headline or verdict. ──
  #' ## Summary metrics and assumption checks (LAT-3138): the figures that summarise the
  #' whole analysis as one metrics row, and one verdict per assumption of the method,
  #' each with the statistic behind it. holds / strained / violated are rules of thumb a
  #' reader can re-derive from the statistic shown beside them. Base R throughout.
  verdict_p <- function(p, soft = 0.05, hard = 0.001) if (is.na(p)) "unknown" else if (p >= soft) "holds" else if (p >= hard) "strained" else "violated"
  fmt_p <- function(p) if (is.na(p)) "" else if (p < 1e-4) "<0.0001" else as.character(signif(p, 3))
  summary_vals <- list(auc = round(auc, 3), accuracy = round(acc, 3), sensitivity = round(sens, 3), specificity = round(spec, 3),
                       precision = round(prec, 3), f1 = round(f1, 3), brier = round(brier, 3), n_test = n_test, n_train = n_train, n = n_used)
  n_terms_c <- tryCatch(max(1, length(coef(model)) - 1), error = function(e) length(model_drivers))
  epv <- (min(n_pos, n_neg) * (n_train / max(1, n_used))) / n_terms_c
  brier_ref <- base_rate * (1 - base_rate)
  Xc <- tryCatch(model.matrix(model), error = function(e) NULL)
  # VIF PER DRIVER (not per column): a factor's dummy columns are collinear with each
  # other by construction, so testing a dummy against its siblings reads "violated" for
  # every categorical driver. Each column is regressed on the OTHER drivers' columns only.
  max_vif <- if (!is.null(Xc) && ncol(Xc) > 2 && length(model_drivers) > 1) {
    Xd <- Xc[, -1, drop = FALSE]; cn <- colnames(Xd)
    owner <- vapply(cn, function(c) { m <- model_drivers[startsWith(c, model_drivers)]; if (length(m)) m[which.max(nchar(m))] else c }, character(1))
    max(vapply(seq_len(ncol(Xd)), function(j) {
      others <- which(owner != owner[j]); if (!length(others)) return(1)
      r2j <- summary(lm(Xd[, j] ~ Xd[, others, drop = FALSE]))$r.squared
      if (is.na(r2j) || r2j >= 1) Inf else 1 / (1 - r2j) }, numeric(1)))
  } else NA_real_
  checks_df <- data.frame(
    check = c("Measured on held-out rows", "Separation beyond chance", "Calibrated probabilities", "Enough events per term", "Class balance", "No strong collinearity"),
    statistic = c(if (holdout) paste0(n_test, " held-out rows") else "in-sample only",
                  paste0("AUC lower bound = ", round(auc_low, 3)),
                  if (multiclass) "not assessed for three or more classes" else paste0("Brier ", round(brier, 3), " against ", round(brier_ref, 3), " for guessing the base rate"),
                  paste0(round(epv, 1), " events per term"),
                  paste0("positive rate = ", round(base_rate, 3)),
                  if (is.na(max_vif)) "not assessed" else paste0("max VIF = ", round(max_vif, 2))),
    p_value = c("", "", "", "", "", ""),
    verdict = c(if (holdout) "holds" else "violated",
                if (auc_low > 0.5) "holds" else if (auc > 0.5) "strained" else "violated",
                if (multiclass) "unknown" else if (brier < 0.9 * brier_ref) "holds" else if (brier < brier_ref) "strained" else "violated",
                if (epv >= 10) "holds" else if (epv >= 5) "strained" else "violated",
                if (base_rate >= 0.1 && base_rate <= 0.9) "holds" else if (base_rate >= 0.05 && base_rate <= 0.95) "strained" else "violated",
                if (is.na(max_vif)) "unknown" else if (max_vif > 10) "violated" else if (max_vif > 5) "strained" else "holds"),
    note = c("in-sample performance flatters the model",
             "a lower bound at or below one half means the ranking could be chance",
             "a Brier score no better than the base rate means the probabilities carry no information",
             "fewer than ten events per term makes the odds ratios unstable",
             "a rare class makes accuracy misleading and the threshold fragile",
             "a VIF above 5 inflates that driver's standard error"),
    stringsAsFactors = FALSE)

  results <- list()
  #' THE ROW-BEARING PLACES USE THE CONSTRUCTORS (LAT-3093/3095): the frame plus the
  #' ROLE of each column. `columns` is derived from names(df) inside the constructor,
  #' so it cannot drift from the data, and a role naming a column the frame lacks
  #' stops HERE rather than at a customer. The three value-bearing places below stay
  #' long-hand on purpose: their explicit `value_order` is already correct and a
  #' constructor would reorder the metric tiles for no gain.
  #' The verdict and the headline are NOT places of a library tool (LAT-3130): a
  #' library tool supplies data objects; the answer is written by the last mile,
  #' which is the first stage that reads the objects together.
  #' ## summary_metrics: the whole analysis at a glance (LAT-3138)
  results$summary_metrics <- place_metric(summary_vals, lead = "auc", place = "summary_metrics")
  #' ## assumption_checks: a verdict per assumption (a `checks` place, table-shaped)
  results$assumption_checks <- place_table(checks_df, place = "assumption_checks")
  #' ## roc_curve: false positive rate against true positive rate on the held-out part
  #' A multi-class ROC carries K curves, so the class column must be DECLARED or the
  #' object mapper sees two numeric columns, finds no series, and draws every class as
  #' one tangled line (LAT-3077).
  results$roc_curve <- place_relationship(roc_df, x = "fpr", y = "tpr",
    series = if (multiclass) "class" else NULL, place = "roc_curve")
  #' ## confusion_table: counts at the chosen threshold
  results$confusion_table <- place_table(confusion_df, place = "confusion_table")
  #' ## performance_table: every figure with the split it was measured on
  results$performance_table <- place_table(perf_df, place = "performance_table")
  #' ## calibration: predicted against observed per decile
  results$calibration <- place_relationship(cal_df, x = "predicted_rate", y = "observed_rate",
    place = "calibration")
  #' ## odds_ratio_interval: each term's odds ratio with its interval (no intercept)
  sel <- match(slope_df$term, coef_df$term)
  or_int <- data.frame(term = slope_df$term, per_sd = slope_df$per_sd,
    low = round((est[sel] - 1.96 * se[sel]) * unname(term_sd[sel]), 3),
    high = round((est[sel] + 1.96 * se[sel]) * unname(term_sd[sel]), 3), stringsAsFactors = FALSE)
  or_int <- or_int[is.finite(or_int$per_sd), , drop = FALSE]; rownames(or_int) <- NULL
  results$odds_ratio_interval <- place_interval(or_int, term = "term", value = "per_sd", low = "low", high = "high",
    place = "odds_ratio_interval")
  #' ## driver_importance: ranked drivers
  results$driver_importance <- place_comparison(imp_df, category = "driver", value = "importance",
    series = "direction", place = "driver_importance")
  #' ## coefficient_table: every term
  results$coefficient_table <- place_table(coef_df, place = "coefficient_table")
  #' ## fit_method: how it was done
  results$fit_method <- list(kind = "metric", values = list(
    method = method, n_in = n_in, n_used = n_used,
    excluded = if (length(dropped)) as.list(unname(vapply(dropped, human, character(1)))) else list(),
    # LAT-3181: the method card's "rows · x → y" line reads these; unset, it printed a bare arrow
    x_column = paste(vapply(model_drivers, human, character(1)), collapse = ", "), y_column = human("outcome"),
    assumptions = list(
      "Log-odds are linear in the numeric drivers; a curved effect is understated.",
      "Rows are independent; repeated rows per customer or time inflate the AUC.",
      "The held-out figures are one split; a different seed moves AUC by roughly its interval width.",
      if (multiclass) "One model per class against the rest: the per-class scores are not calibrated against each other, so the arg-max assignment is a ranking, not a probability over classes." else NULL,
      "The threshold trades sensitivity for specificity by Youden's J, not by the cost of each error; set it to the business cost before acting.",
      if (!holdout) "Performance is in-sample and flatters the model." else "Performance is measured on rows the model never saw.")))

  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))
}