Tree-Based Feature Importance

A stakeholder asks which of your 40 input columns drive fraud, and you show model.feature_importances_ as a bar chart. Half the columns are integer IDs, transaction counts, and raw amounts; one column is a per-row random number you added as a sanity check. If that random column ranks in the top five, the chart you just presented is lying — and it frequently will, because mean decrease in impurity (MDI) importance is biased toward features with many distinct split points. You need to know the formula before you show the chart.

For a tree with N training rows, MDI importance of feature j sums the impurity reduction of every node that splits on j, weighted by the share of rows that reach that node:

importance(j) = Σ over nodes t splitting on j of  (nₜ / N) · Δimpurity(t)

where Δimpurity(t) = impurity(parent) − [n_left·impurity(left) + n_right·impurity(right)] / nₜ

The values are then normalized across features to sum to 1. For Gini impurity, impurity(node) = 1 − Σₖ pₖ²; for entropy, −Σₖ pₖ log₂ pₖ. In scikit-learn this is the same number your gradient boosting model reports:

import numpy as np
from sklearn.ensemble import GradientBoostingClassifier

gbc = GradientBoostingClassifier(n_estimators=300, max_depth=3, random_state=0).fit(X, y)

imp = gbc.feature_importances_          # sums to 1.0
order = np.argsort(imp)[::-1]
for i in order[:8]:
    print(f"{names[i]:<28} {imp[i]:.4f}")

Two properties fall straight out of the formula. It is computed on the training data — no holdout is involved — and it rewards impurity reduction wherever it happens, greedily, wherever the tree found it. That is where the bias enters: a continuous or high-cardinality feature (an ID, a timestamp, a count with thousands of values) offers the split search dozens or hundreds of candidate cut points, so it gets dozens of chances to discover a lucky reduction. A binary feature offers exactly one candidate split. Over hundreds of nodes, the high-cardinality column accumulates importance from noise-driven splits it would never have won on fresh data.

The bias is not hypothetical: the classic simulation in the random-forest literature adds pure-noise columns to real data and watches MDI rank them among the top predictors. Gradient-boosted libraries have their own flavor of the same quantity — LightGBM reports importance by "split" (count of splits) or "gain" (sum of loss reduction, the closer analogue to MDI), and the two rankings often disagree, especially when a feature gets one huge early split versus many small later ones.

So treat MDI as a fast, in-sample screening tool: it tells you which columns the model used, computed at zero extra cost during training. It cannot tell you whether a column earns its keep on unseen data, and it cannot adjudicate causality. When a number from feature_importances_ is going to influence a business decision, you validate it with a measurement that respects the holdout — shuffling one column at a time and watching what breaks, which is permutation importance.

Permutation importance answers a counterfactual question instead of an accounting question: if this feature's relationship to the target were destroyed, how much would my score on held-out data drop? Because it measures the drop on data the model never trained on, it sidesteps the in-sample optimism that undermines impurity importance.

The procedure, with fixed notation: fit the model, score the validation set once to get baseline metric M₀. Then take feature j, randomly shuffle its values across the validation rows (this preserves the column's marginal distribution but severs its link to y), score again to get M_j, and record the drop. Repeat n_repeats times with different shuffles and report the mean ± standard deviation:

from sklearn.inspection import permutation_importance

result = permutation_importance(
    model, X_val, y_val,
    n_repeats=10,
    random_state=42,
    scoring="roc_auc",       # match your business metric, not a default
    n_jobs=-1,
)

order = np.argsort(result.importances_mean)[::-1]
for i in order[:10]:
    mean, std = result.importances_mean[i], result.importances_std[i]
    print(f"{names[i]:<28} {mean:+.4f} ± {std:.4f}")

Reading the output requires three rules of thumb. A mean drop within about one standard deviation of zero means the model does not reliably use that column — rank it as noise, and expect negative values for genuinely useless columns (shuffling noise can accidentally help; do not "promote" those either). A large mean with a large std means the importance is real but unstable across permutations, which itself is a signal worth reporting rather than hiding. And a positive drop across all repeats is the claim you can defend: this feature carries signal the ensemble depends on.

A worked reading: baseline ROC-AUC 0.912; shuffling merchant_id across 10 repeats gives 0.885 ± 0.004 — a drop of 0.027, almost seven times the repeat noise, so the model genuinely leans on it. Shuffling region_code gives 0.911 ± 0.006 — the drop sits inside one standard deviation, indistinguishable from a column the model ignores. Report the mean and the std together; a bare ranking collapses the one distinction that matters.

Two design choices matter more than the parameters. First, run it on held-out data. On the training set, tree predictions are partly memorized, so shuffling may barely dent the score even for a feature that matters at deployment, or dent it in ways that reflect memorization structure rather than predictive value. Second, match scoring to the decision you care about — accuracy on a 97% negative dataset can show near-zero importance for the feature that actually finds your positives; use ROC-AUC, average precision, or a cost-weighted metric instead.

The cost is real — 20 features × 10 repeats = 200 full scoring passes — so on large validation sets, subsample a few thousand rows or raise n_repeats only for the top candidates from a first MDI pass. A common pipeline: screen 40 columns with MDI for free, then spend the permutation budget on the top 15 plus a few suspects.

One failure mode survives both methods untouched, and it is the one that fools the most careful analysts: correlated features, where shuffling one column changes almost nothing because its twin still carries the signal.

Build this deliberately: take age, and add age_copy — the same column plus ±0.1 noise. Train a gradient boosting model, run permutation importance on validation data, and both columns will land near zero importance. Now drop both columns and retrain: accuracy falls off a cliff. The model genuinely depends on them, yet permutation reports them as useless because each one's twin keeps the signal alive when it is shuffled. The opposite bias appears in MDI: impurity importance splits the credit between the two copies arbitrarily, sometimes giving 80% to one and 20% to the other based on tiny, luck-driven differences in which copy a node happened to split on first.

The general shape of the problem: importance methods assign credit per column, but predictive information lives in sets of columns. With correlated features, three distinct phenomena show up in practice:

That last row is the subtle one. Permuting latitude while keeping longitude fixed produces houses at impossible coordinates; the model's confusion on those phantom rows can exceed its real-world sensitivity to latitude, inflating the recorded importance. Shuffling breaks joint structure — usually harmless, occasionally misleading.

Practical remedies, in order of cost:

1. Group permutation — shuffle a block of correlated columns with one shared permutation and measure the joint drop. This restores the signal that individual shuffling hid:

def group_permutation_importance(model, X, y, groups, scoring, n_repeats=5):
    base = scoring(y, model.predict_proba(X)[:, 1])
    out = {}
    for name, cols in groups.items():
        drops = []
        for _ in range(n_repeats):
            Xp = X.copy()
            perm = np.random.permutation(len(Xp))
            Xp[cols] = Xp[cols].to_numpy()[perm]      # one permutation for the whole block
            drops.append(base - scoring(y, model.predict_proba(Xp)[:, 1]))
        out[name] = np.mean(drops)
    return out

You do not need to know the duplicates in advance: cluster columns by pairwise |Spearman ρ| above a threshold (say 0.9) and permute each cluster as a unit. That catches the pairs you created yourself — one-hot columns of the same category, ratios derived from a shared denominator, a raw value and its log.

2. Drop-column importance — retrain without the group and compare. This is the gold standard (it measures what happens in the world you actually deploy), and the reason nobody uses it by default: it costs one full retrain per column.
3. SHAP (TreeSHAP) — allocates credit per prediction and distributes it among correlated features instead of double counting, which is a reasonable tie-breaker for credit assignment — but it does not make a correlational ranking causal.

The interview question hiding here is always the same: this feature ranks first — does that mean it causes the outcome? No. Importance measures dependence of this model, on this data, under this perturbation — not causality, not stability under distribution shift. Say that plainly and you have said more than the person asking.

You have now spent four topics on features that came pre-embedded as vectors. Distance and Similarity Metrics starts from the more primitive question those embeddings assume you can answer: how far apart are two vectors, and which notion of "far" does your problem actually need?