Executive Summary
What the model leans on most when predicting SalePrice, with the uncertainty attached.
Gr Liv Area is what the model leans on most: shuffling it costs 0.274 R-squared (95% interval 0.138 to 0.410) against a baseline of 0.795. However, this lead is not reliable—Gr Liv Area outranked Overall Qual in only 60.8 of every 100 resampling repeats, so the top two should be read as tied. Three predictors have intervals that include zero (Total Bsmt SF, TotRms AbvGrd, Full Bath), meaning no measurable contribution in this fitted model. Gr Liv Area and TotRms AbvGrd, correlated at 0.808, together account for 0.042 R-squared; the joint total is identified but the split between them is not. Of 8 ranking steps, only 1 is genuinely separated by resampling, so the list is reliable only where it explicitly says so. Importance describes what this model uses, not what would happen if you changed a predictor.
Analysis Overview
How importance is measured for 9 predictors of SalePrice.
The analysis fits a linear model to predict SalePrice using 9 predictors across 2,930 rows, then measures importance by shuffling each column and observing the fall in R-squared. The baseline model achieves R-squared 0.795 on held-out predictions. Permutation importance is measured on cross-fitted predictions—each row is scored by a model that never saw it during training—to avoid inflated numbers. The interval on each importance comes from resampling the entire procedure 102 times with folds redrawn and models refitted, which captures both data and split uncertainty. Drop-column importance is reported alongside for contrast: it refits the model without the predictor entirely, answering a different question about what is lost if the column is never collected. Importance measures association, not causation; a high rank can reflect that a column drives the outcome, shares a cause with it, or is a consequence of it.
Data Quality
Rows used, columns dropped, imputation, and the evaluation split.
All 2,930 rows loaded had usable SalePrice values and all mapped predictor columns were usable; no rows were removed. The one categorical predictor had levels beyond the eight most common pooled into 'Other' and blanks into 'Missing', then dummy-coded. Missing values in numeric predictors were filled with each column's median. The 9 surviving predictors were split into 5 folds. The full model reaches R-squared 0.795 on held-out predictions, which is the baseline from which every importance is measured as a fall.
What the Model Uses
Permutation importance with 95% intervals, ranked.
Gr Liv Area leads at 0.274 R-squared, ahead of Overall Qual at 0.245; the smallest is Full Bath at 0.003. The top two predictors together account for the largest share of measured importance. Bars whose whiskers overlap are not distinguishable from one another—the ordering between them is an artifact of where the sample happened to land. The baseline from which all falls are measured is R-squared 0.795. Gr Liv Area's interval spans 0.138 to 0.410, while Overall Qual's spans 0.1844 to 0.3059; these ranges overlap substantially, consistent with their near-tie in the resampling checks.
Importance Detail
Permutation and drop-column importance, direction, and standalone performance for every predictor of SalePrice.
| Predictor | Importance | CI Low | CI High | Share Of Measured | Drop Column Importance | Drop CI Low | Drop CI High | Alone Performance | Direction | Marginal Correlation | Verdict |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Gr Liv Area | 0.2742 | 0.138 | 0.4103 | 41 | 0.0282 | 0.0108 | 0.0456 | 0.497 | higher goes with higher outcome | 0.707 | measurable |
| Overall Qual | 0.2451 | 0.1844 | 0.3059 | 36.7 | 0.0449 | 0.033 | 0.0568 | 0.6382 | higher goes with higher outcome | 0.799 | measurable |
| Year Built | 0.0386 | 0.0253 | 0.0518 | 5.8 | 0.0061 | 0.0039 | 0.0083 | 0.3112 | higher goes with higher outcome | 0.558 | measurable |
| Total Bsmt SF | 0.0336 | -0.0023 | 0.0696 | 5 | 0.0083 | -0.004 | 0.0207 | 0.3949 | higher goes with higher outcome | 0.632 | not distinguishable from zero |
| Neighborhood | 0.0331 | 0.0163 | 0.0499 | 4.9 | 0.0126 | 0.006 | 0.0191 | 0.3196 | highest level: NridgHt | — | measurable |
| Garage Area | 0.0218 | 0.0104 | 0.0332 | 3.3 | 0.0062 | 0.0032 | 0.0093 | 0.4093 | higher goes with higher outcome | 0.64 | measurable |
| Lot Area | 0.0128 | 0.0044 | 0.0212 | 1.9 | 0.0045 | 0.0011 | 0.008 | 0.0567 | higher goes with higher outcome | 0.267 | measurable |
| TotRms AbvGrd | 0.006 | -0.0038 | 0.0158 | 0.9 | 0.0003 | -0.0017 | 0.0023 | 0.2447 | higher goes with lower outcome | 0.495 | not distinguishable from zero |
| Full Bath | 0.0034 | -0.002 | 0.0089 | 0.5 | 0.0001 | -0.0017 | 0.0019 | 0.2965 | higher goes with lower outcome | 0.546 | not distinguishable from zero |
Gr Liv Area, Overall Qual, and Year Built lose far less under drop-column importance than under permutation importance—0.0282, 0.0449, and 0.0061 R-squared respectively versus their permutation values of 0.2742, 0.2451, and 0.0386—meaning the model can rebuild most of what they carry from other columns once refitted. Gr Liv Area achieves 0.497 R-squared on its own and Overall Qual 0.6382, yet both are heavily discounted when other predictors are present. TotRms AbvGrd and Full Bath perform well alone (0.245 and 0.296 R-squared) but show near-zero permutation importance, the signature of redundancy. Total Bsmt SF, TotRms AbvGrd, and Full Bath have intervals that include zero, marking them as not distinguishable from zero contribution. Of the 6 measurable predictors, Gr Liv Area and Overall Qual together hold 77.7% of the positive importance (41% + 36.7%).
Which Ranks Are Real
Adjacent-rank comparisons and whether each step in the order survives resampling.
| Rank | Predictor | Importance | Compared With | Beats Share | Verdict |
|---|---|---|---|---|---|
| 1 | Gr Liv Area | 0.2742 | Overall Qual | 60.8 | not distinguishable |
| 2 | Overall Qual | 0.2451 | Year Built | 100 | separated |
| 3 | Year Built | 0.0386 | Total Bsmt SF | 62.7 | not distinguishable |
| 4 | Total Bsmt SF | 0.0336 | Neighborhood | 52 | neither is measurable |
| 5 | Neighborhood | 0.0331 | Garage Area | 88.2 | not distinguishable |
| 6 | Garage Area | 0.0218 | Lot Area | 89.2 | not distinguishable |
| 7 | Lot Area | 0.0128 | TotRms AbvGrd | 87.3 | not distinguishable |
| 8 | TotRms AbvGrd | 0.006 | Full Bath | 70.6 | neither is measurable |
| 9 | Full Bath | 0.0034 | — | — | last place |
Only 1 of 8 adjacent ranking steps survives the 95% resampling threshold. The first step—Gr Liv Area over Overall Qual—held in only 60.8 of 100 repeats, so these two are not distinguishable. The only separated step is Overall Qual over Year Built, which held in 100% of repeats. All other steps (Year Built over Total Bsmt SF at 62.7%, Neighborhood over Garage Area at 88.2%, Garage Area over Lot Area at 89.2%, Lot Area over TotRms AbvGrd at 87.3%) fall below the 95% bar. Below the Overall Qual–Year Built boundary, the ordering carries no information you can act on.
Overlap and Redundancy
Correlated predictor groups and the difference between redundant and irrelevant.
| Item | Predictors | Max Correlation | Max Vif | Group Importance | Reading |
|---|---|---|---|---|---|
| Correlated group 1 | Gr Liv Area, TotRms AbvGrd | 0.808 | 4.7 | 0.0424 | Together they are worth 0.042 R-squared, and the largest single share credited to any one of them is 0.274. With a variance inflation of 4.70 that split is unstable: a different sample would move it, so read the group total rather than the individual ranks. |
| Redundant: Total Bsmt SF | Total Bsmt SF | 0.547 | 1.77 | — | On its own this column reaches R-squared 0.395, yet its permutation importance is not distinguishable from zero. That is the signature of redundancy, not irrelevance: the other predictors already carry what it knows. |
| Redundant: TotRms AbvGrd | TotRms AbvGrd | 0.808 | 3.13 | — | On its own this column reaches R-squared 0.245, yet its permutation importance is not distinguishable from zero. That is the signature of redundancy, not irrelevance: the other predictors already carry what it knows. |
| Redundant: Full Bath | Full Bath | 0.63 | 2.19 | — | On its own this column reaches R-squared 0.296, yet its permutation importance is not distinguishable from zero. That is the signature of redundancy, not irrelevance: the other predictors already carry what it knows. |
Gr Liv Area and TotRms AbvGrd, correlated at 0.808 with a variance inflation factor of 4.70, together account for 0.042 R-squared; this joint contribution is identified, but the split of 0.274 to Gr Liv Area versus 0.006 to TotRms AbvGrd is unstable and would shift on a different sample. Three predictors are redundant, not irrelevant: Total Bsmt SF reaches R-squared 0.395 alone but shows near-zero permutation importance; TotRms AbvGrd reaches 0.245 alone; Full Bath reaches 0.296 alone. All three perform well in isolation yet add nothing on top of the other predictors, meaning the model has already captured their information through other columns. Variance inflation factors range up to 4.70, below the usual instability threshold of 5.
Methods and Disclosure
The exact procedure, what it can support, and what it cannot.
| Item | Detail |
|---|---|
| Model | Ordinary least squares linear model with SalePrice as the response. |
| Performance metric | R-squared, computed as one minus the residual sum of squares over the total sum of squares. |
| Evaluation scheme | 5-fold cross-fitting: every row is scored by a model fitted without it, so no importance is measured on rows the model was trained on. |
| Permutation importance | Each predictor's values are shuffled inside every held-out fold and the already-fitted model is re-scored. The importance is the fall in R-squared. No model is refitted, so the number answers 'how much does this fitted model lean on the column'. |
| Uncertainty interval | 102 bootstrap replications of the ENTIRE procedure: the rows are resampled, the folds are redrawn, every model is refitted and every predictor is reshuffled. The interval is the mean across replications plus or minus 1.96 standard deviations, so it covers the sampling of the data and of the split, not only the shuffle. Intervals built from shuffles alone make every ranking look far more certain than it is. |
| What the interval cannot settle | Permutation importance cannot fall below zero for a model that leans on a column at all, so its interval is a statement about precision, not proof of a real relationship. A lower bound at or below zero is read here as 'not distinguishable from a column the model does not use'. The drop-column column is the signed check: it can come out negative, and does when a predictor was only adding noise. |
| Drop-column importance | The model is refitted from scratch without the predictor and re-scored under the same cross-fitting, inside every bootstrap replication. It answers a different question from permutation: what is lost if the column is never collected. |
| Standalone performance | The model is refitted using ONLY that predictor. A high standalone score with near-zero permutation importance means the column is redundant rather than uninformative. |
| Collinearity diagnostic | Variance inflation factors from the inverse correlation matrix of the model terms, plus connected groups of predictors correlated at 0.8 or above. Largest variance inflation observed: 4.70. |
| Shapley / SHAP values | Not computed. This analysis reports GLOBAL importance for the whole model; SHAP values are per-observation attributions and are a different quantity. Nothing here is a Shapley value. |
| Causal status | Association only. A predictor can rank high because it causes the outcome, because it shares a cause with it, because it is a consequence of it, or because it leaked from it. This analysis cannot tell those apart. |
The analysis fits an ordinary least squares linear model of SalePrice, scored by R-squared on 5-fold cross-fitted predictions so no importance is measured on training rows. Permutation importance shuffles each predictor inside held-out folds and re-scores the already-fitted model; the importance is the fall in R-squared. The procedure is rerun 102 times on resampled rows with folds redrawn and models refitted, producing the interval from sampling both data and split, not shuffling alone. Drop-column importance refits the model without each predictor and answers what is lost if the column is never collected. No Shapley or SHAP values are computed; this is global permutation importance, a different quantity from per-row attribution. The result is conditional on the 9 mapped predictors; unmapped columns cannot appear and their influence is absorbed by whichever mapped column they correlate with. Importance is association under this model and carries no claim about what would follow from changing any predictor in the real world.
Feature Importance — What the Model Uses
Given an outcome and a set of candidate predictors, this ranks the predictors by how much the model's out-of-sample performance actually depends on each one — and attaches an uncertainty interval to every rank, because a ranking without intervals is a story, not a measurement.
Why This Method?
Permutation importance is model-agnostic and refit-free: shuffle one predictor, re-score the already-fitted model, and measure how much performance falls. Repeating the shuffle many times — each time on a resample of the evaluation rows — turns a single number into an interval, which is what decides whether one predictor really outranks another. A drop-column importance (refit without the predictor) is reported beside it: the two answer different questions and disagreeing is informative.
What This Analysis Covers
- Permutation importance per predictor with a 95% interval
- Drop-column importance for comparison
- Direction of effect and the predictor's standalone performance
- Which adjacent ranks are genuinely separated and which are not
- A collinearity diagnostic that names correlated groups when credit splits
Standard Library
Platform standard-library module (LAT-1441): runs on ANY dataset via the semantic mapping {outcome, predictor_1..predictor_N}. All narrative is derived from the user's own column names and computed values. Importance is an association with the outcome, never a proven cause.
suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))Step 1: Row accounting + semantic column discovery
initial_rows <- nrow(df)
if (!"outcome" %in% names(df)) {
stop(paste0("column_mapping must map an 'outcome' column — the numeric or ",
"two-level outcome whose predictors you want ranked."))
}
pred_cols <- grep("^predictor_[0-9]+$", names(df), value = TRUE)
pred_cols <- pred_cols[order(as.integer(sub("^predictor_", "", pred_cols)))]
outcome_name <- humanize_semantic("outcome", col_map)
if (length(pred_cols) < 2) {
stop(paste0("Feature importance compares predictors against each other, so ",
"column_mapping must map at least two candidate predictor ",
"columns(predictor_1, predictor_2) alongside the outcome '",
outcome_name, "'."))
}
pred_names <- setNames(humanize_semantic(pred_cols, col_map), pred_cols)Step 2: Type the outcome — two distinct values means a logistic model
y_chr <- trimws(as.character(df$outcome))
y_chr[y_chr == ""] <- NA_character_
lev <- sort(unique(y_chr[!is.na(y_chr)]))
if (length(lev) < 2) {
stop(sprintf(paste0("The outcome column '%s' has %d distinct value(s) once ",
"blanks are removed. An outcome that never varies cannot ",
"be predicted, so no predictor can be ranked."),
outcome_name, length(lev)))
}
is_binary <- length(lev) == 2
positive_class <- NA_character_
if (is_binary) {
positive_class <- lev[2]
keep_rows <- !is.na(y_chr)
df <- df[keep_rows, , drop = FALSE]
y <- as.numeric(y_chr[keep_rows] == positive_class)
if (min(table(y)) < 10) {
stop(sprintf(paste0("The outcome column '%s' has only %d row(s) in its ",
"smaller class. At least 10 of each are required before ",
"predictor importance can be estimated."),
outcome_name, min(table(y))))
}
} else {
y_num <- suppressWarnings(as.numeric(y_chr))
n_nonblank <- sum(!is.na(y_chr))
if (sum(!is.na(y_num)) < 0.95 * n_nonblank) {
stop(sprintf(paste0("The outcome column '%s' is neither numeric nor a ",
"two-level yes/no column — only %s of its non-blank ",
"values could be read as numbers."),
outcome_name,
fmt_pct(sum(!is.na(y_num)) / max(1, n_nonblank), 0)))
}
keep_rows <- !is.na(y_num)
df <- df[keep_rows, , drop = FALSE]
y <- y_num[keep_rows]
}
metric_label <- if (is_binary) "AUC" else "R-squared"Step 3: Build a numeric design block per predictor
Numeric columns pass the 95% coercion rule and are median-imputed. Categorical columns are lumped to the 8 most common levels plus "Other" ("Missing" for blanks) and dummy-coded against the most common level. Encoding everything as numbers up front means the shuffle and the refits never hit an unseen-factor-level error.
design_list <- list()
ptype <- character(0)
cat_ref <- list()
drop_col <- character(0)
drop_reason <- character(0)
n_rows_now <- nrow(df)
for (pc in pred_cols) {
v_chr <- trimws(as.character(df[[pc]]))
v_chr[v_chr == ""] <- NA_character_
n_nonblank <- sum(!is.na(v_chr))
if (n_nonblank == 0) {
drop_col <- c(drop_col, pc); drop_reason <- c(drop_reason, "empty"); next
}
v_num <- suppressWarnings(as.numeric(v_chr))
if (sum(!is.na(v_num)) >= 0.95 * n_nonblank) {
med <- median(v_num, na.rm = TRUE)
if (!is.finite(med)) {
drop_col <- c(drop_col, pc); drop_reason <- c(drop_reason, "empty"); next
}
v_num[is.na(v_num)] <- med
vv <- var(v_num)
if (!is.finite(vv) || vv == 0) {
drop_col <- c(drop_col, pc)
drop_reason <- c(drop_reason, "constant"); next
}
uniq_n <- length(unique(v_num))
if (uniq_n == n_rows_now && n_rows_now > 20 &&
all(v_num == round(v_num)) &&
(max(v_num) - min(v_num) + 1) == n_rows_now) {
drop_col <- c(drop_col, pc)
drop_reason <- c(drop_reason, "identifier-like"); next
}
mm <- matrix(v_num, ncol = 1)
colnames(mm) <- pc
design_list[[pc]] <- mm
ptype[pc] <- "numeric"
} else {
v_chr[is.na(v_chr)] <- "Missing"
uniq_n <- length(unique(v_chr))
if (uniq_n < 2) {
drop_col <- c(drop_col, pc)
drop_reason <- c(drop_reason, "constant"); next
}
if (uniq_n > 20 && uniq_n > 0.5 * n_rows_now) {
drop_col <- c(drop_col, pc)
drop_reason <- c(drop_reason, "identifier-like"); next
}
tab <- sort(table(v_chr), decreasing = TRUE)
keep_lev <- names(tab)[seq_len(min(8L, length(tab)))]
v_lump <- ifelse(v_chr %in% keep_lev, v_chr, "Other")
lv <- names(sort(table(v_lump), decreasing = TRUE))
if (length(lv) < 2) {
drop_col <- c(drop_col, pc)
drop_reason <- c(drop_reason, "constant"); next
}
dummies <- vapply(lv[-1], function(L) as.numeric(v_lump == L),
numeric(length(v_lump)))
mm <- matrix(as.numeric(dummies), nrow = length(v_lump))
colnames(mm) <- paste0(pc, "::", lv[-1])
design_list[[pc]] <- mm
ptype[pc] <- "categorical"
cat_ref[[pc]] <- lv
}
}
used <- names(design_list)
if (length(used) < 2) {
survivors <- if (length(used) == 0) "none" else
paste(unname(pred_names[used]), collapse = ", ")
stop(sprintf(paste0("Feature importance needs at least two usable predictor ",
"columns; only %d of the %d mapped columns survived ",
"cleaning(usable: %s). Columns excluded as constant, ",
"empty, or identifier-like: %s."),
length(used), length(pred_cols), survivors,
paste(unname(pred_names[drop_col]), collapse = ", ")))
}
X <- do.call(cbind, design_list)
owner <- rep(names(design_list), vapply(design_list, ncol, integer(1)))Step 4: Remove exactly-redundant design columns (QR pivot)
A predictor that loses ALL of its columns is a perfect linear combination of the others — it is reported, never silently kept.
qr_fit <- qr(cbind(1, X))
piv_keep <- sort(qr_fit$pivot[seq_len(qr_fit$rank)])
keep_idx <- piv_keep[piv_keep > 1] - 1L
if (length(keep_idx) < ncol(X)) {
X <- X[, keep_idx, drop = FALSE]
owner <- owner[keep_idx]
}
lost <- setdiff(used, unique(owner))
if (length(lost) > 0) {
drop_col <- c(drop_col, lost)
drop_reason <- c(drop_reason, rep("perfectly redundant", length(lost)))
used <- intersect(used, unique(owner))
}
if (length(used) < 2) {
stop(sprintf(paste0("Only %d predictor column carries information that is ",
"not already a copy of another mapped column, so there ",
"is nothing to rank. Usable: %s."),
length(used), paste(unname(pred_names[used]), collapse = ", ")))
}
dropped_info <- data.frame(
column = unname(pred_names[drop_col]),
reason = drop_reason,
stringsAsFactors = FALSE
)Cap the modelling rows so the bootstrap below stays inside a sane runtime.
n_all <- nrow(X)
subsampled <- 0L
if (n_all > 4000) {
set.seed(4242)
keep_s <- sort(sample.int(n_all, 4000))
X <- X[keep_s, , drop = FALSE]
y <- y[keep_s]
subsampled <- n_all - 4000L
}
n <- nrow(X); m <- ncol(X); p <- length(used)
final_rows <- n
rows_removed <- initial_rows - final_rowsStep 5: Sample-size guards, stated in the user's own column names
if (n < 60) {
stop(sprintf(paste0("Only %d row(s) have a usable '%s' value. Feature ",
"importance splits the data to score predictors on rows ",
"the model never saw, so at least 60 rows are required."),
n, outcome_name))
}
if (n < 5 * (m + 1)) {
stop(sprintf(paste0("%s rows is too few for the %d model terms produced by ",
"the %d mapped predictor column(s) (%s). Provide at least ",
"%s rows or map fewer predictors."),
format(n, big.mark = ","), m + 1L, p,
paste(unname(pred_names[used]), collapse = ", "),
format(5 * (m + 1), big.mark = ",")))
}Step 6: Cross-fitting — every row is scored by a model that never
saw it, so the importances are not inflated by in-sample fit.
n_folds <- max(3L, min(5L, as.integer(floor(n / 20))))
set.seed(20260728)
fold <- sample(rep(seq_len(n_folds), length.out = n))
if (is_binary) {
ok_folds <- vapply(seq_len(n_folds), function(k) {
tr <- y[fold != k]
length(unique(tr)) == 2 && min(table(tr)) >= 3
}, logical(1))
if (!all(ok_folds)) {
stop(sprintf(paste0("The two levels of '%s' are too unevenly spread to ",
"cross-fit reliably — at least one split left a ",
"training set with almost no examples of one level."),
outcome_name))
}
}
fit_coef <- function(Xtr, ytr) {
A <- cbind(Intercept = rep(1, nrow(Xtr)), Xtr)
b <- if (is_binary) {
f <- tryCatch(
suppressWarnings(glm.fit(A, ytr, family = binomial(),
control = glm.control(maxit = 50))),
error = function(e) NULL)
if (is.null(f)) lm.fit(A, ytr)$coefficients else f$coefficients
} else {
lm.fit(A, ytr)$coefficients
}
b <- as.numeric(b)
b[!is.finite(b)] <- 0
b
}
link_out <- function(b, Xte) {
lp <- as.numeric(cbind(1, Xte) %*% b)
if (is_binary) 1 / (1 + exp(-lp)) else lp
}
metric_fn <- function(yy, pp) {
ok <- is.finite(yy) & is.finite(pp)
if (sum(ok) < 3) return(NA_real_)
yy <- yy[ok]; pp <- pp[ok]
if (is_binary) {
n1 <- sum(yy == 1); n0 <- sum(yy == 0)
if (n1 == 0 || n0 == 0) return(NA_real_)
rk <- rank(pp, ties.method = "average")
(sum(rk[yy == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
} else {
sst <- sum((yy - mean(yy))^2)
if (!is.finite(sst) || sst <= 0) return(NA_real_)
1 - sum((yy - pp)^2) / sst
}
}Out-of-fold predictions from a chosen subset of design columns.
oof_pred <- function(cols) {
out <- numeric(n)
for (k in seq_len(n_folds)) {
tr <- fold != k
b <- fit_coef(X[tr, cols, drop = FALSE], y[tr])
out[!tr] <- link_out(b, X[!tr, cols, drop = FALSE])
}
out
}
all_cols <- seq_len(m)
pred_base <- oof_pred(all_cols)
base_perf <- metric_fn(y, pred_base)Standalone performance: refit using ONLY that predictor. A column that scores well alone but adds nothing to the full model is redundant, which is a different finding from irrelevant.
pred_cols_of <- lapply(seq_len(p), function(jj) which(owner == used[jj]))
alone_perf <- vapply(seq_len(p),
function(jj) metric_fn(y, oof_pred(pred_cols_of[[jj]])),
numeric(1))Step 7: Uncertainty by bootstrapping the WHOLE procedure
Each replication resamples the rows, re-draws the fold split, refits every model, and re-shuffles. The spread across replications therefore carries the sampling of the data, the split, and the shuffle — not just the shuffle, which on its own would make every ranking look certain. The number of replications is set from a fit budget so the runtime does not explode with the predictor count or the row count.
fit_scale <- min(1, 1500 / n)
fit_budget <- (if (is_binary) 3000 else 10000) * fit_scale
n_repeats <- as.integer(floor(fit_budget / (n_folds * (1 + p))))
n_repeats <- max(40L, min(if (n > 2000) 120L else 200L, n_repeats))
perm_imp <- matrix(NA_real_, nrow = p, ncol = n_repeats,
dimnames = list(used, NULL))
drop_imp <- matrix(NA_real_, nrow = p, ncol = n_repeats,
dimnames = list(used, NULL))
set.seed(998877)
for (r in seq_len(n_repeats)) {
bi <- sample.int(n, n, replace = TRUE)
Xb <- X[bi, , drop = FALSE]
yb <- y[bi]
if (is_binary && (length(unique(yb)) < 2 || min(table(yb)) < 5)) next
fr <- sample(rep(seq_len(n_folds), length.out = n))
models <- vector("list", n_folds)
pb <- numeric(n)
for (k in seq_len(n_folds)) {
tr <- fr != k
models[[k]] <- fit_coef(Xb[tr, , drop = FALSE], yb[tr])
pb[!tr] <- link_out(models[[k]], Xb[!tr, , drop = FALSE])
}
base_r <- metric_fn(yb, pb)
if (!is.finite(base_r)) next
for (jj in seq_len(p)) {
jc <- pred_cols_of[[jj]]
pp <- numeric(n); pd <- numeric(n)
for (k in seq_len(n_folds)) {
te <- which(fr == k); tr <- fr != k
Xte <- Xb[te, , drop = FALSE]
if (length(te) > 1) {
pm <- sample.int(length(te))
Xte[, jc] <- Xte[pm, jc, drop = FALSE]
}
pp[te] <- link_out(models[[k]], Xte)
b_drop <- fit_coef(Xb[tr, -jc, drop = FALSE], yb[tr])
pd[te] <- link_out(b_drop, Xb[te, -jc, drop = FALSE])
}
perm_imp[jj, r] <- base_r - metric_fn(yb, pp)
drop_imp[jj, r] <- base_r - metric_fn(yb, pd)
}
}
n_effective <- sum(is.finite(perm_imp[1, ]))
if (n_effective < 20) {
stop(paste0("The resampling needed to put an uncertainty interval on each ",
"importance did not produce enough usable replications for '",
outcome_name, "'. This usually means the outcome is too ",
"imbalanced or the predictor set is too large for the number ",
"of rows supplied."))
}Step 8: Direction of effect on the full-data fit
b_full <- fit_coef(X, y)
direction <- character(p)
marginal_r <- rep(NA_real_, p)
sign_flip <- rep(FALSE, p)
for (jj in seq_len(p)) {
pc <- used[jj]
jc <- which(owner == pc)
if (identical(ptype[[pc]], "numeric")) {
coef_j <- b_full[jc[1] + 1L]
mr <- suppressWarnings(cor(X[, jc[1]], y))
if (!is.finite(mr)) mr <- NA_real_
marginal_r[jj] <- mr
direction[jj] <- if (!is.finite(coef_j) || coef_j == 0) "flat" else
if (coef_j > 0) "higher goes with higher outcome" else
"higher goes with lower outcome"
if (is.finite(mr) && is.finite(coef_j) && mr != 0 && coef_j != 0 &&
sign(mr) != sign(coef_j)) sign_flip[jj] <- TRUE
} else {
coefs <- b_full[jc + 1L]
lv <- cat_ref[[pc]]
best <- if (all(!is.finite(coefs)) || max(coefs, na.rm = TRUE) <= 0) lv[1] else
lv[-1][which.max(coefs)]
direction[jj] <- paste0("highest level: ", best)
}
}Step 9: Collinearity — variance inflation plus correlated groups
Cm <- suppressWarnings(cor(X))
Cm[!is.finite(Cm)] <- 0
vif_col <- tryCatch({
d <- diag(solve(Cm))
d[!is.finite(d)] <- Inf
d
}, error = function(e) rep(Inf, m))
vif_pred <- vapply(seq_len(p), function(jj)
max(vif_col[which(owner == used[jj])]), numeric(1))
max_vif <- suppressWarnings(max(vif_pred[is.finite(vif_pred)]))
if (!is.finite(max_vif)) max_vif <- Inf
Coff <- abs(Cm); diag(Coff) <- 0
cross <- outer(owner, owner, FUN = "!=")
max_cross_cor <- if (m > 1) max(Coff * cross) else 0Connected components over design columns correlated above 0.8.
adj <- (Coff >= 0.8) & cross
comp <- rep(NA_integer_, m); cid <- 0L
for (i in seq_len(m)) {
if (!is.na(comp[i])) next
cid <- cid + 1L
stack <- i
while (length(stack) > 0) {
cur <- stack[1]; stack <- stack[-1]
if (!is.na(comp[cur])) next
comp[cur] <- cid
nb <- which(adj[cur, ])
stack <- c(stack, nb[is.na(comp[nb])])
}
}
groups <- list()
for (g in unique(comp)) {
members <- unique(owner[comp == g])
if (length(members) >= 2) groups[[length(groups) + 1]] <- members
}
group_names <- lapply(groups, function(g) unname(pred_names[g]))The one number a correlated group DOES identify: drop every member at once and refit. The split between members is arbitrary; the joint contribution is not.
group_joint <- if (length(groups) > 0) {
vapply(groups, function(g) {
gc <- which(owner %in% g)
base_perf - metric_fn(y, oof_pred(setdiff(all_cols, gc)))
}, numeric(1))
} else numeric(0)
collinear <- length(groups) > 0 || (is.finite(max_vif) && max_vif >= 5) ||
!is.finite(max_vif)Step 10: Rank, summarise, and judge separation
The interval is the bootstrap mean plus or minus 1.96 standard deviations across replications. A normal interval rather than a percentile one, because the replication count is set by a fit budget and percentile tails are unstable when that count is modest.
band <- function(M) {
mu <- rowMeans(M, na.rm = TRUE)
sdv <- apply(M, 1, function(z) sd(z[is.finite(z)]))
sdv[!is.finite(sdv)] <- 0
list(mu = as.numeric(mu), lo = as.numeric(mu - 1.96 * sdv),
hi = as.numeric(mu + 1.96 * sdv))
}
pb_band <- band(perm_imp); db_band <- band(drop_imp)
imp_mean <- pb_band$mu; imp_lo <- pb_band$lo; imp_hi <- pb_band$hi
drop_mean <- db_band$mu; drop_lo <- db_band$lo; drop_hi <- db_band$hi
verdict <- ifelse(imp_lo > 0, "measurable", "not distinguishable from zero")
pos_total <- sum(pmax(imp_mean, 0))
share <- if (pos_total > 0) 100 * pmax(imp_mean, 0) / pos_total else rep(0, p)A predictor with no measurable permutation importance but real standalone performance is REDUNDANT, not irrelevant — the distinction matters.
alone_floor <- if (is_binary) 0.55 else 0.05
redundant <- verdict == "not distinguishable from zero" &
is.finite(alone_perf) & alone_perf >= alone_floor
imp_df <- data.frame(
semantic = used,
predictor = unname(pred_names[used]),
importance = round(as.numeric(imp_mean), 4),
ci_low = round(as.numeric(imp_lo), 4),
ci_high = round(as.numeric(imp_hi), 4),
share = round(as.numeric(share), 1),
drop_importance = round(as.numeric(drop_mean), 4),
drop_ci_low = round(as.numeric(drop_lo), 4),
drop_ci_high = round(as.numeric(drop_hi), 4),
alone = round(as.numeric(alone_perf), 4),
direction = direction,
marginal_r = round(as.numeric(marginal_r), 3),
vif = round(as.numeric(ifelse(is.finite(vif_pred), vif_pred, 999)), 2),
verdict = verdict,
redundant = redundant,
sign_flip = sign_flip,
stringsAsFactors = FALSE
)
ord <- order(-imp_df$importance)
imp_df <- imp_df[ord, , drop = FALSE]
rownames(imp_df) <- NULL
perm_ord <- perm_imp[ord, , drop = FALSE]
n_measurable <- sum(imp_df$verdict == "measurable")
n_redundant <- sum(imp_df$redundant)
beats <- rep(NA_real_, p - 1)
for (i in seq_len(p - 1)) {
d <- perm_ord[i, ] - perm_ord[i + 1, ]
beats[i] <- mean(d > 0, na.rm = TRUE)
}A step is only a real ranking claim when the HIGHER predictor cleared zero. Two predictors that both fail that test can still order consistently, but the ordering says nothing worth acting on.
sep_verdict <- ifelse(imp_df$verdict[seq_len(p - 1)] != "measurable",
"neither is measurable",
ifelse(beats >= 0.95, "separated", "not distinguishable"))
rank_separation_df <- data.frame(
rank = seq_len(p),
predictor = imp_df$predictor,
importance = imp_df$importance,
compared_with = c(imp_df$predictor[-1], "—"),
beats_share = c(round(100 * beats, 1), NA_real_),
verdict = c(sep_verdict, "last place"),
stringsAsFactors = FALSE
)
n_separated <- sum(sep_verdict == "separated")
top_name <- imp_df$predictor[1]
top_imp <- imp_df$importance[1]
top_lo <- imp_df$ci_low[1]
top_hi <- imp_df$ci_high[1]
top_separated <- p > 1 && sep_verdict[1] == "separated"Step 11: Leakage screen — a suspiciously perfect model, computed
leak_perf <- if (is_binary) 0.99 else 0.98
best_alone <- suppressWarnings(max(imp_df$alone[is.finite(imp_df$alone)]))
if (!is.finite(best_alone)) best_alone <- NA_real_
leakage_flag <- is.finite(base_perf) && base_perf >= leak_perf &&
is.finite(best_alone) && best_alone >= 0.95 * base_perf
leak_name <- if (isTRUE(leakage_flag))
imp_df$predictor[which.max(replace(imp_df$alone, !is.finite(imp_df$alone), -Inf))] else NA_character_Which failure mode actually happened here is read off the numbers, not assumed: either each member is individually droppable while the group is not, or the fit divided the credit unevenly.
shape <- if (is.finite(joint) && is.finite(max_drop) &&
max_drop < 0.4 * joint) {
paste0("Removing any single one of them costs at most ", fmt(max_drop),
" ", metric_label, ", because whichever is left rebuilds the ",
"signal. Removing all ", length(members), " at once costs ",
fmt(joint), ". The group's contribution is real; the split of ",
"it between the members is not identified, so do not rank them ",
"against each other.")
} else {
paste0("Together they are worth ", fmt(joint), " ", metric_label,
", and the largest single share credited to any one of them is ",
fmt(max_perm),
". With a variance inflation of ",
fmt(min(999, max(vif_col[mcols])), 2),
" that split is unstable: a different sample would move it, so ",
"read the group total rather than the individual ranks.")
}
group_readings <- c(group_readings,
paste0(paste(hn, collapse = " and "), " (correlated at ",
fmt(max(sub), 3), ") together account for ",
fmt(joint), " ", metric_label))
red_rows[[length(red_rows) + 1]] <- data.frame(
item = paste0("Correlated group ", gi),
predictors = paste(hn, collapse = ", "),
max_correlation = round(max(sub), 3),
max_vif = round(min(999, max(vif_col[mcols])), 2),
group_importance = round(joint, 4),
reading = shape,
stringsAsFactors = FALSE
)
}
}
if (n_redundant > 0) {
for (rn in which(imp_df$redundant)) {
red_rows[[length(red_rows) + 1]] <- data.frame(
item = paste0("Redundant: ", imp_df$predictor[rn]),
predictors = imp_df$predictor[rn],
max_correlation = round(max(abs(Cm[which(owner == imp_df$semantic[rn]),
which(owner != imp_df$semantic[rn]),
drop = FALSE])), 3),
max_vif = imp_df$vif[rn],
group_importance = NA_real_,
reading = paste0("On its own this column reaches ",
metric_label, " ", fmt(imp_df$alone[rn]),
", yet its permutation importance is not ",
"distinguishable from zero. That is the signature of ",
"redundancy, not irrelevance: the other predictors ",
"already carry what it knows."),
stringsAsFactors = FALSE
)
}
}
if (length(red_rows) == 0) {
red_rows[[1]] <- data.frame(
item = "No overlapping group detected",
predictors = paste(imp_df$predictor, collapse = ", "),
max_correlation = round(max_cross_cor, 3),
max_vif = round(min(999, max_vif), 2),
group_importance = NA_real_,
reading = paste0("The strongest correlation between any two mapped ",
"predictors is ", fmt(max_cross_cor, 3),
" and the largest variance inflation is ",
fmt(min(999, max_vif), 2),
", both below the level at which credit starts to split, ",
"so each importance can be read on its own."),
stringsAsFactors = FALSE
)
}
redundancy_df <- do.call(rbind, red_rows)
rownames(redundancy_df) <- NULL
methods_df <- data.frame(
item = c(
"Model",
"Performance metric",
"Evaluation scheme",
"Permutation importance",
"Uncertainty interval",
"What the interval cannot settle",
"Drop-column importance",
"Standalone performance",
"Collinearity diagnostic",
"Shapley / SHAP values",
"Causal status"
),
detail = c(
if (is_binary)
paste0("Logistic regression(binomial GLM) on the two levels of ",
outcome_name, "; '", positive_class,
"' is scored as the event.")
else
paste0("Ordinary least squares linear model with ", outcome_name,
" as the response."),
if (is_binary)
"Area under the ROC curve, computed from the rank statistic of the predicted probabilities."
else
"R-squared, computed as one minus the residual sum of squares over the total sum of squares.",
paste0(n_folds, "-fold cross-fitting: every row is scored by a model ",
"fitted without it, so no importance is measured on rows the ",
"model was trained on."),
paste0("Each predictor's values are shuffled inside every held-out fold ",
"and the already-fitted model is re-scored. The importance is the ",
"fall in ", metric_label,
". No model is refitted, so the number answers 'how much does this ",
"fitted model lean on the column'."),
paste0(n_effective, " bootstrap replications of the ENTIRE procedure: the ",
"rows are resampled, the folds are redrawn, every model is refitted ",
"and every predictor is reshuffled. The interval is the mean across ",
"replications plus or minus 1.96 standard deviations, so it covers ",
"the sampling of the data and of the split, not only the shuffle. ",
"Intervals built from shuffles alone make every ranking look far ",
"more certain than it is."),
paste0("Permutation importance cannot fall below zero for a model that ",
"leans on a column at all, so its interval is a statement about ",
"precision, not proof of a real relationship. A lower bound at or ",
"below zero is read here as 'not distinguishable from a column the ",
"model does not use'. The drop-column column is the signed check: ",
"it can come out negative, and does when a predictor was only ",
"adding noise."),
paste0("The model is refitted from scratch without the predictor and ",
"re-scored under the same cross-fitting, inside every bootstrap ",
"replication. It answers a different question from permutation: ",
"what is lost if the column is never collected."),
paste0("The model is refitted using ONLY that predictor. A high standalone ",
"score with near-zero permutation importance means the column is ",
"redundant rather than uninformative."),
paste0("Variance inflation factors from the inverse correlation matrix of ",
"the model terms, plus connected groups of predictors correlated at ",
"0.8 or above. Largest variance inflation observed: ",
fmt(min(999, max_vif), 2), "."),
paste0("Not computed. This analysis reports GLOBAL importance for the ",
"whole model; SHAP values are per-observation attributions and are ",
"a different quantity. Nothing here is a Shapley value."),
paste0("Association only. A predictor can rank high because it causes the ",
"outcome, because it shares a cause with it, because it is a ",
"consequence of it, or because it leaked from it. This analysis ",
"cannot tell those apart.")
),
stringsAsFactors = FALSE
)Step 13: KPI metrics
metrics <- list(
`Rows Used` = final_rows,
`Predictors Ranked` = p,
`Model` = if (is_binary) "Logistic regression" else "Linear regression",
`Performance Metric` = metric_label,
`Baseline Performance` = round(base_perf, 3),
`Top Predictor` = top_name,
`Top Importance` = round(top_imp, 3),
`Measurable Predictors` = as.integer(n_measurable),
`Separated Rank Steps` = as.integer(n_separated)
)Step 14: json_output machine channel
sep_clause <- if (p > 1) {
if (top_separated)
paste0(" It is separated from ", imp_df$predictor[2], ", which comes second at ",
fmt(imp_df$importance[2]), ": ", top_name, " ranked ahead in ",
fmt(100 * beats[1], 1), "% of the repeats.")
else
paste0(" It is NOT reliably separated from ", imp_df$predictor[2],
" (second at ", fmt(imp_df$importance[2]), "); ", top_name,
" ranked ahead in only ", fmt(100 * beats[1], 1),
"% of the repeats, so treat the top two as tied.")
} else ""
collinear_clause <- if (collinear && length(groups) > 0) {
paste0(" Collinearity warning: ", paste(group_readings, collapse = "; "),
", so the individual importances inside that overlap are not ",
"separately identified and should not be ranked against each other.")
} else if (collinear) {
paste0(" Collinearity warning: the largest variance inflation across the ",
"mapped predictors is ", fmt(min(999, max_vif), 2),
", high enough that individual importances are unstable.")
} else ""
json_output <- list(
answer = paste0(
"Permutation importance for ", outcome_name, " across ",
n_things(final_rows, "row"), " and ", n_things(p, "candidate predictor"),
", measured on ", n_folds, "-fold cross-fitted predictions(",
metric_label, " of the full model: ", fmt(base_perf), "). ",
top_name, " is the most important, costing ", fmt(top_imp), " ",
metric_label, " when shuffled(95% interval ", fmt(top_lo), " to ",
fmt(top_hi), ").", sep_clause, " ",
n_things(n_measurable, "predictor"),
" of ", p, vb(n_measurable, " has", " have"),
" an interval entirely above zero.", collinear_clause,
" Importance is association with ", outcome_name,
", not evidence that changing a predictor would change it."
),
cards = lapply(
c("tldr", "overview", "preprocessing", "importance_chart",
"importance_table", "rank_reliability", "redundancy", "methods"),
function(cid) list(id = cid, metrics = metrics)
)
)
list(
initial_rows = initial_rows, final_rows = final_rows, rows_removed = rows_removed,
outcome_name = outcome_name, pred_names = pred_names, used = used,
dropped_info = dropped_info, ptype = ptype,
is_binary = is_binary, positive_class = positive_class,
metric_label = metric_label, base_perf = base_perf,
n_folds = n_folds, n_repeats = n_effective, subsampled = subsampled,
imp_df = imp_df, importance_df = importance_df,
importance_detail_df = importance_detail_df,
rank_separation_df = rank_separation_df,
redundancy_df = redundancy_df, methods_df = methods_df,
groups = groups, group_names = group_names,
group_readings = group_readings, group_joint = group_joint,
collinear = collinear, max_vif = max_vif, max_cross_cor = max_cross_cor,
beats = beats, sep_verdict = sep_verdict, n_separated = n_separated,
top_name = top_name, top_imp = top_imp, top_lo = top_lo, top_hi = top_hi,
top_separated = top_separated,
n_measurable = n_measurable, n_redundant = n_redundant,
leakage_flag = leakage_flag, leak_name = leak_name, best_alone = best_alone,
metrics = metrics, json_output = json_output
)
}