Standard Feature Importance
Executive Summary

Executive Summary

What the model leans on most when predicting SalePrice, with the uncertainty attached.

Rows Used
2930
Predictors Ranked
9
Model
Linear regression
Performance Metric
R-squared
Baseline Performance
0.795
Top Predictor
Gr Liv Area
Top Importance
0.274
Measurable Predictors
6
Separated Rank Steps
1
Across 2,930 rows, 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. That lead is not reliable: Gr Liv Area outranked Overall Qual in only 60.8 of every 100 repeats, so the top two should be read as tied. 3 predictors have an interval that includes zero (Total Bsmt SF, TotRms AbvGrd, Full Bath), which means no measurable contribution here — not proof of irrelevance, since a column can be redundant with one that is already in. Read the individual ranks with care: Gr Liv Area and TotRms AbvGrd (correlated at 0.808) together account for 0.042 R-squared. That joint total is identified; the split of it between those columns is not identified, so their relative order carries no information. Of the 8 steps down the ranking, 1 is genuinely separated and 7 are not, so the list is a ranking only where it says so. Importance describes what this model uses, not what would happen if you changed a predictor.
What this means

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.

Overview

Analysis Overview

How importance is measured for 9 predictors of SalePrice.

N Observations2930
N Predictors9
Baseline Performance0.795
N Folds5
What this means

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 Preparation

Data Quality

Rows used, columns dropped, imputation, and the evaluation split.

Initial Rows2930
Final Rows2930
Rows Removed0
N Predictors9
Baseline Performance0.795
What this means

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.

Visualization

What the Model Uses

Permutation importance with 95% intervals, ranked.

What this means

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.

Data Table

Importance Detail

Permutation and drop-column importance, direction, and standalone performance for every predictor of SalePrice.

PredictorImportanceCI LowCI HighShare Of MeasuredDrop Column ImportanceDrop CI LowDrop CI HighAlone PerformanceDirectionMarginal CorrelationVerdict
Gr Liv Area0.27420.1380.4103410.02820.01080.04560.497higher goes with higher outcome0.707measurable
Overall Qual0.24510.18440.305936.70.04490.0330.05680.6382higher goes with higher outcome0.799measurable
Year Built0.03860.02530.05185.80.00610.00390.00830.3112higher goes with higher outcome0.558measurable
Total Bsmt SF0.0336-0.00230.069650.0083-0.0040.02070.3949higher goes with higher outcome0.632not distinguishable from zero
Neighborhood0.03310.01630.04994.90.01260.0060.01910.3196highest level: NridgHtmeasurable
Garage Area0.02180.01040.03323.30.00620.00320.00930.4093higher goes with higher outcome0.64measurable
Lot Area0.01280.00440.02121.90.00450.00110.0080.0567higher goes with higher outcome0.267measurable
TotRms AbvGrd0.006-0.00380.01580.90.0003-0.00170.00230.2447higher goes with lower outcome0.495not distinguishable from zero
Full Bath0.0034-0.0020.00890.50.0001-0.00170.00190.2965higher goes with lower outcome0.546not distinguishable from zero
What this means

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%).

Data Table

Which Ranks Are Real

Adjacent-rank comparisons and whether each step in the order survives resampling.

RankPredictorImportanceCompared WithBeats ShareVerdict
1Gr Liv Area0.2742Overall Qual60.8not distinguishable
2Overall Qual0.2451Year Built100separated
3Year Built0.0386Total Bsmt SF62.7not distinguishable
4Total Bsmt SF0.0336Neighborhood52neither is measurable
5Neighborhood0.0331Garage Area88.2not distinguishable
6Garage Area0.0218Lot Area89.2not distinguishable
7Lot Area0.0128TotRms AbvGrd87.3not distinguishable
8TotRms AbvGrd0.006Full Bath70.6neither is measurable
9Full Bath0.0034last place
What this means

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.

Data Table

Overlap and Redundancy

Correlated predictor groups and the difference between redundant and irrelevant.

ItemPredictorsMax CorrelationMax VifGroup ImportanceReading
Correlated group 1Gr Liv Area, TotRms AbvGrd0.8084.70.0424Together 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 SFTotal Bsmt SF0.5471.77On 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 AbvGrdTotRms AbvGrd0.8083.13On 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 BathFull Bath0.632.19On 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.
What this means

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.

Data Table

Methods and Disclosure

The exact procedure, what it can support, and what it cannot.

ItemDetail
ModelOrdinary least squares linear model with SalePrice as the response.
Performance metricR-squared, computed as one minus the residual sum of squares over the total sum of squares.
Evaluation scheme5-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 importanceEach 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 interval102 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 settlePermutation 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 importanceThe 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 performanceThe 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 diagnosticVariance 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 valuesNot 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 statusAssociation 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.
What this means

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.

Rate this report Was this the answer you needed?
The exact source that produced this report — yours to keep, read, and re-run.
Download PDF
How this was computed method · R source · citation
The code that did it

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 &#x27;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 &#x27;",
                outcome_name, "&#x27;."))
  }
  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 &#x27;%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 &#x27;%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 &#x27;%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_rows

Step 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 &#x27;%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 &#x27;%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 &#x27;",
                outcome_name, "&#x27;. 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 0

Connected 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&#x27;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, "; &#x27;", positive_class,
               "&#x27; 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&#x27;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 &#x27;how much does this ",
             "fitted model lean on the column&#x27;."),
      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 &#x27;not distinguishable from a column the ",
             "model does not use&#x27;. 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
  )
}
Your data has more stories to tell.Run any analysis on your own data — validated R modules, interactive reports, AI insights, and PDF export. 500 free credits on signup.
Try Free — No SignupSign Up Free

Cite this analysis

Report an Issue

Tell us what's wrong. You'll get a free re-run of this analysis so you can try again with different parameters. If the re-run still doesn't meet your expectations, we'll refund your credits.

Want to run this analysis on your own data? Upload CSV — Free Analysis See Pricing