Model Evaluation II — Interpretability: impurity vs permutation vs SHAP¶
Which features matter, how much to trust each measure, and what the model actually does¶
A model that predicts well is only half of what a quant or a risk committee needs; they also need to know why. But "feature importance" is not one thing — different measures can disagree sharply, and one of the most common (a tree's built-in importance) is biased. This notebook lines up the three standard tools, shows where each misleads, and settles on the principled one:
- Impurity / Gini importance (MDI) — a tree's built-in measure: total impurity reduction from splits on each feature. Free, but biased toward high-cardinality / continuous features (Strobl et al., 2007) — it can rank pure noise above real signal.
- Permutation importance — model-agnostic: shuffle a feature and measure the drop in test performance. Much less biased, though it can be misled when features are strongly correlated.
- SHAP (SHapley Additive exPlanations) — game-theoretic per-prediction attributions that are signed, locally exact (TreeSHAP), and consistent; the modern standard, and uniquely it explains individual predictions, not just global averages. Note what it is faithful to: SHAP explains the model. If the model has learned something spurious, a correct SHAP attribution reports that it did.
We reproduce the impurity bias, use SHAP for global and local explanations, and read effect shapes with partial dependence and ICE. Data: the Taiwan credit-default problem. Python-lead.
1. Three measures, three answers¶
To test the measures we inject two pure-noise features into the credit data — one high-cardinality (random integers 0–9999), one binary — that have no relationship to default, then fit a random forest and score every feature three ways. A trustworthy measure should place both at the bottom.
The results need to be read as magnitudes, not ranks, and the distinction matters more than it sounds. A rank of 17 out of 25 seems damning until you notice that most of those 25 features have importances statistically indistinguishable from zero, so their ordering is close to arbitrary — which is why a rank can move ten places when the random seed changes while the underlying number does not move at all. Ranks are what most write-ups report; they are the least stable thing in the table.
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
import shap
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("credit_default.csv"); feat=[c for c in d.columns if c!="default"]
X=d[feat].copy(); y=d["default"].values; rng=np.random.default_rng(0)
X["NOISE_hi"]=rng.integers(0,10000,len(X)); X["NOISE_bin"]=rng.integers(0,2,len(X)); cols=list(X.columns)
Xtr,Xte,ytr,yte=train_test_split(X.values,y,test_size=0.3,random_state=0,stratify=y)
sub=rng.choice(len(Xtr),8000,replace=False)
rf=RandomForestClassifier(150,min_samples_leaf=20,random_state=0,n_jobs=-1).fit(Xtr[sub],ytr[sub])
imp=rf.feature_importances_
pi=permutation_importance(rf,Xte[:3000],yte[:3000],n_repeats=20,random_state=0,scoring="roc_auc")
perm=pi.importances_mean; perm_sd=pi.importances_std # 20 repeats, so the noise is measurable
sv=shap.TreeExplainer(rf).shap_values(Xte[:800]); sv=sv[...,1] if np.ndim(sv)==3 else sv; shp=np.abs(sv).mean(0)
R=pd.DataFrame({"impurity (MDI)":imp/imp.sum(),"permutation":perm/np.abs(perm).sum(),"SHAP":shp/shp.sum()},index=cols)
print("What each measure says about the two features known to carry no information:")
print(f" {'measure':16s} {'NOISE_hi':>22} {'NOISE_bin':>22}")
for meas in R.columns:
r=R[meas].rank(ascending=False)
print(f" {meas:16s} {R.loc['NOISE_hi',meas]:>+9.4f} (rank {int(r['NOISE_hi']):>2}/{len(cols)})"
f" {R.loc['NOISE_bin',meas]:>+9.4f} (rank {int(r['NOISE_bin']):>2}/{len(cols)})")
_top=R['impurity (MDI)'].idxmax()
print(f"\n for scale, the strongest real feature is {_top}: MDI {R.loc[_top,'impurity (MDI)']:.4f}, SHAP {R.loc[_top,'SHAP']:.4f}.")
print(f" So MDI awards the high-cardinality noise {R.loc['NOISE_hi','impurity (MDI)']/R['impurity (MDI)'].max():.0%} of the strongest real feature's importance, and SHAP "
f"{R.loc['NOISE_hi','SHAP']/R['SHAP'].max():.0%}.")
print(f" Permutation gives it {pi.importances_mean[cols.index('NOISE_hi')]:+.5f} AUC -- NEGATIVE. Shuffling it slightly IMPROVES")
print(" held-out performance, which is the right answer for a column that contains nothing.")
print("\n Read as ranks the three measures look similar and the story is muddled: permutation puts NOISE_bin at rank")
print(f" {int(R['permutation'].rank(ascending=False)['NOISE_bin'])} while MDI puts it at {int(R['impurity (MDI)'].rank(ascending=False)['NOISE_bin'])}, which would suggest MDI is the better measure. Read as magnitudes")
print(" the ordering is unambiguous: permutation is right, SHAP is nearly right, MDI is wrong by a wide margin.")
print("\n Why the ranks are the wrong summary -- permutation importance with its standard error:")
print(f" {'feature':12s} {'importance':>12} {'sd':>9} {'mean/sd':>9}")
_ord=np.argsort(perm)[::-1]
for j in list(_ord[:4])+[cols.index('NOISE_hi'),cols.index('NOISE_bin')]:
print(f" {cols[j]:12s} {perm[j]:>12.5f} {perm_sd[j]:>9.5f} {(perm[j]/perm_sd[j] if perm_sd[j]>0 else 0):>9.2f}")
_ind=sum(1 for j in range(len(cols)) if perm[j] < 2*perm_sd[j])
print(f" {_ind} of {len(cols)} features have importance under two standard errors of zero. Ordering THOSE is arbitrary,")
print(" and that is exactly the region both noise features live in.")
show=["PAY_1","PAY_2","LIMIT_BAL","AGE","BILL_AMT1","NOISE_hi","NOISE_bin"]
Rs=R.loc[show]; fig,ax=plt.subplots(figsize=(11,4.4)); w=0.27; xp=np.arange(len(show))
for k,(m,c) in enumerate(zip(R.columns,[GREY,BLUE,GREEN])): ax.bar(xp+(k-1)*w,Rs[m],w,color=c,label=m)
ax.set_xticks(xp); ax.set_xticklabels(show,rotation=30,ha="right"); ax.set_ylabel("normalized importance"); ax.set_title("Three importance measures — note the two NOISE features")
ax.legend(); plt.tight_layout(); plt.show()
print("Impurity (MDI) hands the high-cardinality NOISE feature a tenth of the strongest real feature's importance. SHAP")
print("gives it less but still not zero. Only permutation, measured on held-out data, calls it worthless outright --")
print("and the reason is the next section.")
What each measure says about the two features known to carry no information: measure NOISE_hi NOISE_bin impurity (MDI) +0.0253 (rank 13/25) +0.0043 (rank 25/25) permutation -0.0072 (rank 24/25) +0.0017 (rank 15/25) SHAP +0.0118 (rank 18/25) +0.0035 (rank 25/25) for scale, the strongest real feature is PAY_1: MDI 0.2516, SHAP 0.2303. So MDI awards the high-cardinality noise 10% of the strongest real feature's importance, and SHAP 5%. Permutation gives it -0.00114 AUC -- NEGATIVE. Shuffling it slightly IMPROVES held-out performance, which is the right answer for a column that contains nothing. Read as ranks the three measures look similar and the story is muddled: permutation puts NOISE_bin at rank 15 while MDI puts it at 25, which would suggest MDI is the better measure. Read as magnitudes the ordering is unambiguous: permutation is right, SHAP is nearly right, MDI is wrong by a wide margin. Why the ranks are the wrong summary -- permutation importance with its standard error: feature importance sd mean/sd PAY_1 0.07605 0.00393 19.37 LIMIT_BAL 0.01965 0.00233 8.43 PAY_2 0.01682 0.00297 5.66 PAY_AMT1 0.00695 0.00156 4.46 NOISE_hi -0.00114 0.00070 -1.62 NOISE_bin 0.00028 0.00025 1.09 18 of 25 features have importance under two standard errors of zero. Ordering THOSE is arbitrary, and that is exactly the region both noise features live in.
Impurity (MDI) hands the high-cardinality NOISE feature a tenth of the strongest real feature's importance. SHAP gives it less but still not zero. Only permutation, measured on held-out data, calls it worthless outright -- and the reason is the next section.
2. Why they disagree — cardinality and correlation¶
Two mechanisms drive the whole disagreement, and both can be demonstrated rather than asserted.
Cardinality. A tree measures a feature's importance by the impurity it removes, and a feature offering more distinct split points gets more chances to remove impurity by luck. Strobl et al. (2007) identified this as the source of MDI's bias. The test below is direct: five columns of pure noise, identical in every respect except how many distinct values they take — 2, 10, 100, 1,000 and 10,000. If cardinality is the mechanism, MDI should rise along that sequence while the held-out measure stays flat at zero.
Correlation. Permutation importance asks what happens when one column is scrambled. If a second column carries nearly the same information, the model simply leans on the substitute and the measured drop is small — so two genuinely useful correlated features can both look unimportant. The credit data has a natural test case: six monthly bill-amount columns. Permuting them one at a time and permuting the whole block together answer different questions, and the gap between the two is the size of the problem.
# --- cardinality: five noise columns differing only in how many distinct values they take -----------------
lv={"N_2":2,"N_10":10,"N_100":100,"N_1000":1000,"N_10000":10000}
X2=d[feat].copy(); r2_=np.random.default_rng(1)
for k,v in lv.items(): X2[k]=r2_.integers(0,v,len(X2))
c2=list(X2.columns)
X2tr,X2te,y2tr,y2te=train_test_split(X2.values,y,test_size=0.3,random_state=0,stratify=y)
s2=r2_.choice(len(X2tr),8000,replace=False)
rf2=RandomForestClassifier(150,min_samples_leaf=20,random_state=0,n_jobs=-1).fit(X2tr[s2],y2tr[s2])
i2=rf2.feature_importances_; i2=i2/i2.sum()
p2=permutation_importance(rf2,X2te[:3000],y2te[:3000],n_repeats=10,random_state=0,scoring="roc_auc").importances_mean
v2=shap.TreeExplainer(rf2).shap_values(X2te[:800]); v2=v2[...,1] if np.ndim(v2)==3 else v2
v2=np.abs(v2).mean(0); v2=v2/v2.sum()
print("Five columns of pure noise, differing ONLY in cardinality:")
print(f" {'column':10s} {'distinct values':>16} {'MDI':>9} {'SHAP':>9} {'permutation':>13}")
for k,v in lv.items():
j=c2.index(k); print(f" {k:10s} {v:>16,} {i2[j]:>9.4f} {v2[j]:>9.4f} {p2[j]:>+13.5f}")
print(f" {'(PAY_1)':10s} {'real feature':>16} {i2[c2.index('PAY_1')]:>9.4f} {v2[c2.index('PAY_1')]:>9.4f} {p2[c2.index('PAY_1')]:>+13.5f}")
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ks=list(lv.values())
ax[0].semilogx(ks,[i2[c2.index(k)] for k in lv],"o-",color=GREY,lw=2,label="impurity (MDI)")
ax[0].semilogx(ks,[v2[c2.index(k)] for k in lv],"s-",color=GREEN,lw=2,label="SHAP")
ax[0].semilogx(ks,[max(p2[c2.index(k)],0) for k in lv],"^-",color=BLUE,lw=2,label="permutation (held out)")
ax[0].axhline(0,color="k",lw=.6); ax[0].set_xlabel("distinct values in the noise column")
ax[0].set_ylabel("normalized importance"); ax[0].set_title("Importance awarded to pure noise, by cardinality"); ax[0].legend(fontsize=8)
# --- correlation: the six bill-amount columns ------------------------------------------------------------
from sklearn.metrics import roc_auc_score
bill=[c for c in cols if c.startswith("BILL_AMT")]
cb=np.corrcoef(d[bill].values.T)[np.triu_indices(len(bill),1)].mean()
base=roc_auc_score(yte[:3000],rf.predict_proba(Xte[:3000])[:,1])
rp=np.random.default_rng(0); joint=[]
for _ in range(10):
Xp=Xte[:3000].copy(); pidx=rp.permutation(len(Xp))
for b in bill: Xp[:,cols.index(b)]=Xp[pidx,cols.index(b)]
joint.append(base-roc_auc_score(yte[:3000],rf.predict_proba(Xp)[:,1]))
indiv=[perm[cols.index(b)] for b in bill]
ax[1].bar(range(len(bill)),indiv,color=BLUE,label="permuted individually")
ax[1].axhline(np.mean(joint),color=RED,ls="--",lw=2,label=f"whole block permuted together ({np.mean(joint):.4f})")
ax[1].axhline(sum(indiv),color=ORANGE,ls=":",lw=2,label=f"sum of individual ({sum(indiv):.4f})")
ax[1].set_xticks(range(len(bill))); ax[1].set_xticklabels(bill,rotation=30,ha="right")
ax[1].set_ylabel("drop in test AUC"); ax[1].set_title(f"Correlated block (mean corr {cb:.2f})"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"\nMDI rises from {i2[c2.index('N_2')]:.4f} to {max(i2[c2.index(k)] for k in lv):.4f} across that sequence -- a {max(i2[c2.index(k)] for k in lv)/i2[c2.index('N_2')]:.0f}-fold increase driven by nothing but the")
print("number of distinct values, and it plateaus once the column offers enough split points to be effectively continuous.")
print(f"Permutation stays at zero throughout (all five between {min(p2[c2.index(k)] for k in lv):+.5f} and {max(p2[c2.index(k)] for k in lv):+.5f}). SHAP sits in between,")
print("and for a reason worth being precise about: SHAP is faithful to THE MODEL, and this forest really did split on")
print("those columns. A correct attribution has to report that. It is not measuring whether the feature is useful --")
print("it is measuring what the model does with it, and those are different questions whenever the model is imperfect.")
print(f"\nOn correlation: the {len(bill)} bill-amount columns correlate at {cb:.2f}. Permuted one at a time they cost {sum(indiv):.4f} AUC in")
print(f"total; permuted together they cost {np.mean(joint):.4f}, about {np.mean(joint)/sum(indiv):.1f} times as much. Tested individually each column looks")
print("dispensable because the other five still carry its information -- so a low permutation importance means 'this")
print("feature adds nothing GIVEN THE OTHERS', which is not the same as 'this feature does not matter'. Group correlated")
print("columns and permute them together when that distinction matters.")
Five columns of pure noise, differing ONLY in cardinality: column distinct values MDI SHAP permutation N_2 2 0.0047 0.0039 -0.00024 N_10 10 0.0116 0.0062 -0.00045 N_100 100 0.0222 0.0100 +0.00014 N_1000 1,000 0.0258 0.0117 -0.00165 N_10000 10,000 0.0223 0.0085 -0.00004 (PAY_1) real feature 0.2369 0.2149 +0.07564
MDI rises from 0.0047 to 0.0258 across that sequence -- a 5-fold increase driven by nothing but the number of distinct values, and it plateaus once the column offers enough split points to be effectively continuous. Permutation stays at zero throughout (all five between -0.00165 and +0.00014). SHAP sits in between, and for a reason worth being precise about: SHAP is faithful to THE MODEL, and this forest really did split on those columns. A correct attribution has to report that. It is not measuring whether the feature is useful -- it is measuring what the model does with it, and those are different questions whenever the model is imperfect. On correlation: the 6 bill-amount columns correlate at 0.89. Permuted one at a time they cost 0.0090 AUC in total; permuted together they cost 0.0137, about 1.5 times as much. Tested individually each column looks dispensable because the other five still carry its information -- so a low permutation importance means 'this feature adds nothing GIVEN THE OTHERS', which is not the same as 'this feature does not matter'. Group correlated columns and permute them together when that distinction matters.
3. SHAP — global and local explanations¶
SHAP assigns each feature, for each prediction, a signed contribution such that the contributions sum to the prediction (minus a baseline) — Shapley values from cooperative game theory, computed exactly and fast for trees (TreeSHAP). This gives two things no importance ranking can:
- a global view (the beeswarm) — every feature's impact across all clients, coloured by feature value, showing direction and magnitude at once;
- a local view — the exact attribution for one client's prediction, the explanation a credit decision or a regulator actually requires.
The beeswarm confirms the whole tree family's finding — recent repayment status PAY_1 dominates and being in arrears pushes default risk up — and the local plot decomposes a single high-risk client's score into its drivers.
Xte_df=pd.DataFrame(Xte[:800],columns=cols)
shap.summary_plot(sv,Xte_df,show=False,max_display=10,plot_size=(9,5)); plt.title("SHAP beeswarm — global feature impact on default risk"); plt.tight_layout(); plt.show()
# local explanation for the highest-risk client in the sample
proba=rf.predict_proba(Xte[:800])[:,1]; i=int(np.argmax(proba))
contrib=pd.Series(sv[i],index=cols).sort_values(key=np.abs,ascending=False).head(8)[::-1]
fig,ax=plt.subplots(figsize=(8,4)); ax.barh(contrib.index,contrib.values,color=[RED if v>0 else BLUE for v in contrib.values])
ax.axvline(0,color="k",lw=.6); ax.set_xlabel("SHAP value (push toward default →)"); ax.set_title(f"Local explanation: client #{i} (predicted P(default)={proba[i]:.2f})")
plt.tight_layout(); plt.show()
print(f"Global: PAY_* (recent repayment) dominates, high values (red) push risk up. Local: client #{i}'s {proba[i]:.0%} risk")
print("decomposes into specific drivers -- the per-decision explanation a risk committee or regulator needs, which no global ranking provides.")
Global: PAY_* (recent repayment) dominates, high values (red) push risk up. Local: client #491's 83% risk decomposes into specific drivers -- the per-decision explanation a risk committee or regulator needs, which no global ranking provides.
4. Partial dependence and ICE — the shape of an effect¶
Importance says how much a feature matters; partial dependence (PDP) says in which direction and shape. It plots the model's average predicted probability as one feature varies, marginalising over the others. ICE (Individual Conditional Expectation) draws one line per client instead of the average, revealing heterogeneity and interactions the PDP's average can hide. Below, default risk vs the credit limit and vs recent repayment status: risk falls with a higher limit and rises sharply once a client is months behind — the model's learned relationships, made legible.
from sklearn.inspection import PartialDependenceDisplay
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
PartialDependenceDisplay.from_estimator(rf,Xte[:2000],[cols.index("LIMIT_BAL")],feature_names=cols,kind="both",ax=ax[0],ice_lines_kw={"alpha":.1,"color":GREY},pd_line_kw={"color":RED,"lw":2})
ax[0].set_title("PDP + ICE: default risk vs credit limit"); ax[0].set_xlabel("credit limit (NT dollars)")
PartialDependenceDisplay.from_estimator(rf,Xte[:2000],[cols.index("PAY_1")],feature_names=cols,kind="both",ax=ax[1],ice_lines_kw={"alpha":.1,"color":GREY},pd_line_kw={"color":RED,"lw":2})
ax[1].set_title("PDP + ICE: default risk vs recent repayment (PAY_1)"); ax[1].set_xlabel("PAY_1 (months delay)")
plt.tight_layout(); plt.show()
print("PDP (red) shows the average effect: risk decreasing in the credit limit, jumping up once PAY_1 signals arrears.")
print("The ICE lines (grey) show per-client curves -- their spread flags interactions the single average PDP line hides.")
PDP (red) shows the average effect: risk decreasing in the credit limit, jumping up once PAY_1 signals arrears. The ICE lines (grey) show per-client curves -- their spread flags interactions the single average PDP line hides.
5. Summary¶
Feature importance is not one number, and the three standard ones disagree in a way that is systematic rather than random. Given a pure-noise column, impurity (MDI) awarded it 10% of the strongest real feature's importance and SHAP 5%, while permutation returned a negative value — shuffling it slightly improves held-out AUC, the right answer for a column containing nothing.
The mechanism behind the disagreement is worth more than the ranking. MDI's bias tracks cardinality: five noise columns identical except in how many distinct values they take receive steadily more MDI as that number rises, saturating once there are enough split points to matter. SHAP's residual mass is not a bug at all — SHAP explains the model, and the forest genuinely did split on the noise, so a faithful attribution must say so. Permutation is the only one of the three scored on held-out data, which is why it is the only one that answers "does this feature help the model generalise?"
A second caveat is usually stated and rarely shown: with correlated features, individual permutation importances understate a group. The six bill-amount columns correlate at 0.89, and permuting the block jointly costs about 1.5× the sum of permuting them one at a time — credit that is shared between correlated columns goes missing when each is tested alone. SHAP then went further than any ranking: a global beeswarm showing each feature's signed impact, and local attributions explaining a single prediction — the decision-level transparency finance actually requires. Partial dependence and ICE supplied the effect shapes (risk falls with the credit limit, spikes with arrears) and flagged heterogeneity.
A third, structural point: importances should be read as magnitudes with uncertainty, not as ranks. Most of the features here have permutation importance within two standard errors of zero, so their ordering is close to arbitrary and moves ten places between random seeds while the underlying quantity does not move at all. A published table of ranks conveys a precision that the numbers do not have.
Guidance: never rely on impurity importance across mixed-cardinality features; permutation importance on held-out data is the measure to trust for "does this feature help," tested in groups when features are correlated and reported with a standard error; SHAP when you need signed, consistent, per-prediction explanations — remembering that it faithfully describes the model rather than the world; PDP/ICE for the shape of an effect. Cross-links: this makes explicit the impurity-bias caveat from the Random Forests notebook and generalises the SHAP used in the XGBoost and capstone notebooks. Next in the subsection: calibration — are the predicted probabilities themselves correct, and how to fix them when they are not.