Feature importance: MDI, MDA, SFI and SHAP.

What importance is for (and what it is not for)

A backtest says whether a strategy made money in a sample; it does not say why. Feature importance tries to answer the second question: which variables a model uses and how much it loses without them. López de Prado (2018, chapter 8) sums it up in a working rule: the backtest is not a research tool, feature importance is. You research with importance (which phenomenon lies behind, which features are superfluous, which hypothesis deserves a new experiment) and use the backtest only once, at the end, to validate (lessons 2.5.4 and 5.4.4).

There are four families of measures, and each one answers a different question:

flowchart TB
  F[Features] --> M[Model]
  M --> I[Importance]
  I --> H[Hypothesis]
  H --> F
  I -.-> B[Final backtest]

The problem with judging these measures on real data is that you never know what the answer is. The level 5 panel does have it.

The test bench: known truth and two targets

The panel has 400 fictitious stocks, 729 labeled weeks and 22 weekly features (without revision_analistas, the leak from lesson 5.2.4). The level's README publishes what each one does in the generating model:

To those 22, a control feature built in the cell itself is added here, ruido_binario: a coin toss, with only two values. The features enter as each week's normal rank (lesson 5.2.4) and there are two targets:

• The real label: next week's return minus the weekly mean, in basis points.
• The true alpha minus its mean: the predictable part of the label, without the noise. No real model has this target; it serves to see what the methods would say if the signal were not buried.

import numpy as np
import pandas as pd
from scipy import stats

f = pd.read_parquet("panel5_features.parquet")
e = pd.read_parquet("panel5_etiquetas.parquet")
v = pd.read_parquet("panel5_verdad.parquet")
d = (f.merge(e[["fecha", "ticker", "ret_sig_semana",
               "fin_sig_semana"]], on=["fecha", "ticker"])
      .merge(v[["fecha", "ticker", "alfa_verdadero"]],
             on=["fecha", "ticker"])
      .dropna(subset=["ret_sig_semana"])
      .reset_index(drop=True))
FEATS = [c for c in f.columns
         if c not in ("fecha", "ticker", "sector",
                      "revision_analistas")]

sem = d.groupby("fecha")
def gauss(col):
    n = sem[col].transform("count")
    return stats.norm.ppf((sem[col].rank() - 0.5) / n)

X = pd.DataFrame({c: gauss(c) for c in FEATS}).fillna(0.0)
rng = np.random.default_rng(0)
X["ruido_binario"] = rng.integers(0, 2, len(X)).astype(float)
cols = list(X.columns)

def menos_media(col):   # in basis points
    return ((d[col] - sem[col].transform("mean")) * 1e4).to_numpy()

y = menos_media("ret_sig_semana")      # the real label
a = menos_media("alfa_verdadero")      # the noise-free target
print(X.shape)
print(f"label std {y.std():.0f} bp · "
      f"alpha std {a.std():.1f} bp")

(291600, 22)
label std 366 bp · alpha std 12.7 bp

The label has a standard deviation of 366 bp and the alpha, of 12.7: the noise is 29 times the signal. That ratio decides almost everything that follows.

The model is a random forest of 50 trees with López de Prado's choices for measuring importance: max_features=1 (at each split one feature is chosen at random, so that none hides its correlated peers, the masking), large leaves (min_samples_leaf=1000) and 20% subsamples per tree. The out-of-sample score is the weekly IC (lesson 5.1.3), and the validation, a purged k-fold of five blocks of weeks with an embargo of four (lesson 5.4.2): with a one-week label the purge barely removes rows; the embargo prevents persistent features from connecting the end of one block with the start of the next.

from sklearn.ensemble import RandomForestRegressor

fechas = np.sort(d["fecha"].unique())
semana = np.searchsorted(fechas, d["fecha"].to_numpy())
fin = d["fin_sig_semana"].to_numpy()
inicio = d["fecha"].to_numpy()

def pliegues(k=5, embargo=4):
    for b in np.array_split(np.arange(len(fechas)), k):
        te = np.isin(semana, b)
        ini_te, fin_te = fechas[b[0]], fin[te].max()
        solapa = (fin > ini_te) & (inicio <= fin_te)   # purge
        emb = (semana > b[-1]) & (semana <= b[-1] + embargo)
        yield ~te & ~solapa & ~emb, te

def ic_semanal(pred, obj, sem_idx):
    r = pd.DataFrame({"p": pred, "y": obj}).groupby(sem_idx).rank()
    r = r - r.groupby(sem_idx).transform("mean")
    num = (r["p"] * r["y"]).groupby(sem_idx).sum()
    den = np.sqrt((r["p"] ** 2).groupby(sem_idx).sum()
                  * (r["y"] ** 2).groupby(sem_idx).sum())
    return (num / den).mean()

def bosque(n=50, muestra=0.2):
    return RandomForestRegressor(
        n_estimators=n, max_features=1, min_samples_leaf=1000,
        max_samples=muestra, n_jobs=4, random_state=0)

for i, (tr, te) in enumerate(pliegues()):
    print(i, f"test {te.sum():6d}  train {tr.sum():6d}  "
          f"out {(~tr & ~te).sum():5d}")

0 test  58400  train 231600  out  1600
1 test  58400  train 231600  out  1600
2 test  58400  train 231600  out  1600
3 test  58400  train 231600  out  1600
4 test  58000  train 233600  out     0

MDI: in-sample importance

Each split of a tree reduces the impurity (in regression, the squared error) of its node. A feature's MDI is the sum of those reductions over all the splits that use it, weighted by the fraction of samples that pass through the node, averaged over the trees and normalized so that the total is 100%. It is free (it comes out of training) and it is what feature_importances_ returns in scikit-learn.

mdi = pd.DataFrame({
    "alpha": bosque().fit(X, a).feature_importances_,
    "label": bosque().fit(X, y).feature_importances_,
}, index=cols) * 100
mdi.sort_values("alpha", ascending=False).round(2)

                    alpha  label
ret_1s              21.60   6.47
vol_63d             12.75   6.64
valor_bp             8.17   4.76
ey                   8.12   5.59
rsi_14               7.98   5.30
sue                  7.60   4.63
mom_12_1             6.11   5.47
sentimiento          5.66   3.05
vol_21d              5.46   5.41
ret_1m               4.78   4.84
roe                  3.98   5.06
dist_max_52s         3.06   4.36
ruido_4              1.90   4.13
log_capitalizacion   0.92   4.22
ret_3m               0.86   4.57
volumen_rel          0.53   3.87
dias_desde_informe   0.13   4.85
ruido_3              0.10   3.67
ruido_1              0.08   4.15
ruido_2              0.08   3.63
ruido_5              0.08   4.22
ruido_binario        0.02   1.13

With the noise-free target, MDI resembles the truth: ret_1s (reversal, the strongest effect) takes 21.6%, the pure noise features 0.1% or less, and the coin 0.02%. But it already shows its first bias: it shares the effect among correlated features. ey (8.1%) and rsi_14 (8.0%), with no effect of their own, beat mom_12_1 (6.1%), which does have one, and ret_1m (4.8%) comes close to it: with max_features=1, when the tree cannot pick the true one, it picks its neighbor.

With the real label, MDI separates nothing. All continuous features sit between 3.6% and 6.6%, around the 4.5% of a uniform split; ruido_1 (4.2%) and ruido_5 (4.2%) are within tenths of sue (4.6%). The reason is that MDI is measured in sample: if 99.9% of the label's variance is noise, the splits that reduce it most are the ones that fit noise, and any continuous feature will do. The two exceptions are sentimiento (3.1%), which has 32% of empty values filled with 0, and the coin (1.1%): with fewer distinct values there are fewer possible splits. It is MDI's cardinality bias, documented in scikit-learn's own help: it favors features with many values, even if they are noise.

MDA: what is lost out of sample

Permutation importance (MDA) is measured on data the model has not seen: it is trained on each fold, scored on the test set, one feature is shuffled in the test set (breaking its relationship with the label while keeping its distribution) and it is scored again. The drop in score is the importance; it is averaged over the folds and its standard error comes from the dispersion between them. Here the drop is measured in thousandths of IC.

def mda(obj, grupos=None, n=50):
    grupos = grupos or [[c] for c in cols]
    rng = np.random.default_rng(0)
    caidas, base = {" + ".join(g): [] for g in grupos}, []
    for tr, te in pliegues():
        m = bosque(n).fit(X.values[tr], obj[tr])
        Xte = X.values[te]
        s0 = ic_semanal(m.predict(Xte), obj[te], semana[te])
        base.append(s0)
        for g in grupos:
            Xp = Xte.copy()
            j = [cols.index(c) for c in g]
            Xp[:, j] = Xp[rng.permutation(len(Xp))][:, j]
            s = ic_semanal(m.predict(Xp), obj[te], semana[te])
            caidas[" + ".join(g)].append(s0 - s)
    c = pd.DataFrame(caidas).T * 1000      # thousandths of IC
    return pd.DataFrame({"mda": c.mean(axis=1),
                         "se": c.std(axis=1) / np.sqrt(c.shape[1])}
                        ), np.mean(base)

mda_a, base_a = mda(a)
mda_y, base_y = mda(y)
print(f"out-of-sample IC: alpha {base_a:.3f} · "
      f"label {base_y:.4f}")
tabla = pd.DataFrame({"mda alpha": mda_a["mda"],
                      "mda label": mda_y["mda"],
                      "se label": mda_y["se"]})
tabla.sort_values("mda label", ascending=False).round(2)

out-of-sample IC: alpha 0.652 · label 0.0148
                    mda alpha  mda label  se label
ret_1s                 156.82       4.04      1.09
valor_bp                25.71       3.27      0.89
ey                      22.81       1.84      0.54
mom_12_1                26.92       1.67      0.58
ruido_1                  0.01       1.20      0.50
ruido_binario           -0.08       1.02      0.66
sentimiento             34.83       0.82      0.47
ret_1m                  -0.11       0.80      0.66
sue                     39.69       0.70      0.60
dias_desde_informe      -0.11       0.58      0.59
ruido_4                  5.47       0.47      0.24
roe                      2.17       0.29      0.50
dist_max_52s             8.70       0.23      0.23
volumen_rel              2.28       0.22      0.21
rsi_14                  14.20       0.19      0.81
ruido_2                  0.18       0.09      0.52
vol_63d                 57.46      -0.04      0.49
log_capitalizacion       1.37      -0.08      0.41
ruido_5                  0.29      -0.23      0.82
vol_21d                  9.30      -0.41      1.02
ruido_3                  0.14      -0.41      0.30
ret_3m                  -1.19      -0.74      0.54

Out of sample, the forest achieves an IC of 0.652 with the noise-free target and of 0.0148 with the real label; the oracle's ceiling with the real label is 0.029 (lesson 5.2.4). With the noise-free target, MDA is clean: all noise features stay at zero and those with their own effect at the top. But it has the opposite bias to MDI's: ret_1m (−0.1) and roe (2.2) fall to almost nothing and rsi_14 (14.2) sits well below sue or mom_12_1, because permuting a feature when its neighbor is still in the model barely costs anything: the model leans on the other one. It is the substitution effect.

With the real label, MDA separates somewhat: ret_1s (4.0 ± 1.1), valor_bp (3.3 ± 0.9), ey (1.8 ± 0.5) and mom_12_1 (1.7 ± 0.6) stay at the top. But the list continues with ruido_1 (1.2 ± 0.5, $t \approx 2.4$) and ruido_binario (1.0 ± 0.7): two noise features ahead of sue and sentimiento, which do have an effect. And vol_63d, with a planted effect, comes out at zero: the label only subtracts the weekly mean, and in the panel volatile stocks have more beta, so the low-volatility premium is masked by beta (the level's README measures it).

The figure at the top puts the two real-label columns side by side. With 22 features, some noise feature passing $t = 2$ is to be expected (lesson 2.1.8). That is why it pays to always put several known noise controls in the model, like ruido_binario here: their importance gives the yardstick of what chance produces, and a real feature has to sit clearly above all of them.

SFI: each feature separately

SFI removes substitution by construction: it trains a model with a single feature and measures its out-of-sample score, with the same validation. It is expensive (one model per feature and fold), so here each one has 10 trees and 5% subsamples.

def sfi(obj, n=10):
    res = {}
    for j, c in enumerate(cols):
        s = [ic_semanal(bosque(n, 0.05)
                        .fit(X.values[tr][:, [j]], obj[tr])
                        .predict(X.values[te][:, [j]]),
                        obj[te], semana[te])
             for tr, te in pliegues()]
        res[c] = {"sfi": np.mean(s) * 1000,
                  "se": np.std(s, ddof=1) / np.sqrt(len(s)) * 1000}
    return pd.DataFrame(res).T

sfi_a, sfi_y = sfi(a), sfi(y)
pd.DataFrame({"sfi alpha": sfi_a["sfi"], "sfi label": sfi_y["sfi"],
              "se label": sfi_y["se"]}
             ).sort_values("sfi alpha", ascending=False).round(1)

                    sfi alpha  sfi label  se label
ret_1s                  463.8       12.5       3.0
vol_63d                 300.3        2.7       1.9
sue                     283.5        2.1       1.4
ey                      276.4        1.9       0.9
mom_12_1                269.2        5.2       1.4
valor_bp                245.8        3.8       2.6
vol_21d                 232.0        1.1       2.3
rsi_14                  223.6        3.5       1.0
sentimiento             211.6        1.4       0.8
roe                     185.9        2.7       1.0
ret_1m                  183.1        3.3       2.2
dist_max_52s            173.9        2.6       1.4
ruido_4                 159.4       -2.5       3.4
ret_3m                   63.2        0.4       2.1
log_capitalizacion       51.2        0.8       1.8
ruido_3                   8.7       -1.0       2.1
volumen_rel               5.5        0.7       1.7
dias_desde_informe        1.2        0.2       1.3
ruido_binario             1.1       -2.5       1.3
ruido_1                   1.0        0.3       1.9
ruido_5                  -1.2        1.5       1.7
ruido_2                  -2.0       -1.3       2.1

With the noise-free target, SFI shows the other face of substitution: ret_1m (183), rsi_14 (224), ey (276) and roe (186) look almost as predictive as the true ones, because on their own they carry their neighbor's effect. And it does not see interactions: volumen_rel (5.5) and dias_desde_informe (1.2) predict nothing on their own, even though the full model uses them to modulate reversal and SUE.

With the real label, SFI singles out ret_1s (12.5 ± 3.0) and, at a distance, mom_12_1 (5.2 ± 1.4) and rsi_14 (3.5 ± 1.0), which on its own carries part of the reversal. Features with no effect of their own such as roe (2.7 ± 1.0) sit at the level of valor_bp (3.8 ± 2.6) and sue (2.1 ± 1.4), and ruido_4 and ruido_binario come out negative (−2.5): a single-feature model with no effect learns noise from the training set and pays for it in the test.

Importance by groups of correlated features

López de Prado's solution (2020, chapter 6) to substitution is to stop asking about each feature and ask about groups: cluster the features by their correlation and measure the MDA by permuting the whole group, so that the model cannot take refuge in a neighbor. The usual distance is $\sqrt{(1 - \rho)/2}$, which is 0 for two identical features and 1 for two opposite ones; with average linkage and a cut at 0.5 the groups come out:

from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import squareform

C = X[FEATS].corr()
dist = np.sqrt(0.5 * (1 - C.to_numpy()))
Z = linkage(squareform(dist, checks=False), method="average")
etiq = fcluster(Z, t=0.5, criterion="distance")
grupos = {}
for c, k in zip(FEATS, etiq):
    grupos.setdefault(k, []).append(c)
grupos = list(grupos.values()) + [["ruido_binario"]]
for g in grupos:
    if len(g) > 1:
        print(" + ".join(g))
print(len(grupos), "groups")

ret_1s + ret_1m + rsi_14
ret_3m + dist_max_52s
mom_12_1 + ruido_4
vol_21d + vol_63d
valor_bp + ey + roe
15 groups

The five groups with more than one feature make economic sense: short-term reversal (ret_1s, ret_1m, rsi_14), medium-term return (ret_3m, dist_max_52s), momentum (mom_12_1 and its noisy copy ruido_4), volatility (vol_21d, vol_63d) and value-quality (valor_bp, ey, roe). López de Prado proposes choosing the number of groups with a clustering algorithm using the silhouette as the criterion (ONC); a fixed threshold on a dendrogram, as here, is the simple version and is enough when the groups are this clear.

cmda_a, _ = mda(a, grupos)
cmda_y, _ = mda(y, grupos)
pd.DataFrame({"alpha": cmda_a["mda"], "label": cmda_y["mda"],
              "se label": cmda_y["se"]}
             ).sort_values("label", ascending=False).round(2)

                           alpha  label  se label
ret_1s + ret_1m + rsi_14  258.38   6.31      2.27
valor_bp + ey + roe       124.07   4.45      1.90
vol_21d + vol_63d         112.37   2.44      1.29
mom_12_1 + ruido_4         41.51   2.40      0.71
ruido_binario              -0.01   1.27      0.82
log_capitalizacion          1.29   0.81      0.46
ruido_1                     0.10   0.68      0.43
sue                        40.23   0.58      0.63
sentimiento                34.80   0.53      0.31
dias_desde_informe          0.15   0.50      0.75
volumen_rel                 2.28   0.46      0.20
ret_3m + dist_max_52s       8.43   0.28      0.66
ruido_5                     0.20  -0.09      0.71
ruido_3                     0.10  -0.18      0.20
ruido_2                     0.08  -0.19      0.23

With the noise-free target, the reversal group (258) collects what the individual MDA shared out badly, and the value-quality group (124) no longer depends on which of the three features the model looks at. With the real label, the four groups with a planted effect are the top four: reversal (6.3 ± 2.3), value-quality (4.5 ± 1.9), volatility (2.4 ± 1.3) and momentum (2.4 ± 0.7). Behind them, the coin (1.3 ± 0.8) is a reminder of the size of the measurement's noise.

SHAP: the contribution to each prediction

Shapley values come from cooperative game theory: they split the “payoff” of a coalition among its players according to their average contribution over all possible coalitions. SHAP (Lundberg and Lee, 2017) applies them to a prediction: the difference between a row's prediction and the mean is split among the features, and the sum of the values is exactly that difference. For trees, shap.TreeExplainer computes them exactly and fast. The usual global importance is the mean of the absolute values:

import shap

tr, te = list(pliegues())[-1]
m_y = bosque().fit(X.values[tr], y[tr])
muestra = np.random.default_rng(1).choice(np.flatnonzero(te), 3000,
                                          replace=False)
sv = shap.TreeExplainer(m_y).shap_values(X.values[muestra])
shap_y = pd.Series(np.abs(sv).mean(axis=0), index=cols)
shap_y.sort_values(ascending=False).round(2).head(10)

ret_1s      1.02
rsi_14      0.94
ret_1m      0.83
ruido_1     0.71
mom_12_1    0.67
ruido_5     0.66
valor_bp    0.64
ey          0.59
roe         0.53
vol_63d     0.51
dtype: float64

Computed on test rows of the last fold with the forest trained on the real label, the mean SHAP ranks ret_1s first, but puts ruido_1 (0.71) and ruido_5 (0.66) in the top 10 ahead of valor_bp: SHAP explains what the model does, and this model, trained on noise, uses noise. It is not a measure of how much a feature predicts, but of how much it moves the predictions.

SHAP's advantage is in the detail: it says how a feature's contribution depends on its value and on those of the others. The generator plants an interaction: reversal is multiplied by $(1 - 0.5\cdot g(\text{volumen\_rel}))$, stronger when volume is low. For a term $c\cdot x\cdot(1 - 0.5z)$, the Shapley value of $x$ is $c\cdot x\cdot(1 - 0.25z)$: half of the cross term is attributed to $z$. The slope of the SHAP of ret_1s in each tercile of volumen_rel, relative to the middle tercile, can be compared with that exact answer:

import lightgbm as lgb

m_a = bosque().fit(X.values[tr], a[tr])
m_l = lgb.LGBMRegressor(
    n_estimators=300, learning_rate=0.05, num_leaves=15,
    min_child_samples=500, subsample=0.5, subsample_freq=1,
    random_state=0, n_jobs=4, deterministic=True,
    force_row_wise=True, verbose=-1).fit(X.values[tr], a[tr])
xs = X.iloc[muestra]
j = cols.index("ret_1s")
tercio = pd.cut(xs["volumen_rel"], [-np.inf, -0.43, 0.43, np.inf],
                labels=["low", "medium", "high"])

def pendientes(modelo):
    sv = shap.TreeExplainer(modelo).shap_values(X.values[muestra])
    p = pd.Series({k: np.polyfit(xs["ret_1s"][tercio == k],
                                 sv[(tercio == k).to_numpy(), j], 1)[0]
                   for k in ["low", "medium", "high"]})
    return p / p["medium"]

media_g = xs.groupby(tercio, observed=True)["volumen_rel"].mean()
t = pd.DataFrame({"forest": pendientes(m_a),
                  "LightGBM": pendientes(m_l),
                  "exact SHAP": 1 - 0.25 * media_g,
                  "true effect": 1 - 0.5 * media_g})
t.round(2)

        forest  LightGBM  exact SHAP  true effect
low       1.05      1.25        1.27         1.55
medium    1.00      1.00        1.00         1.00
high      0.91      0.72        0.72         0.45

The LightGBM (trees of up to 15 leaves, which can combine features along the same path) gives 1.25 and 0.72, right on the exact Shapley value (1.27 and 0.72). The forest of the previous sections gives 1.05 and 0.91: with one random feature per split and leaves of 1,000 rows, it barely learned the interaction, and SHAP reflects that faithfully. The explanation is of the model, not of the data: if the model did not learn something, SHAP is not going to show it, and if it learned noise, it will show it as if it were signal.

How to use all this

flowchart TB
  A[Candidates] --> G[Groups]
  G --> P[MDA with purging]
  P --> C{Above controls?}
  C -- no --> D[Discard]
  C -- yes --> S[SHAP: shape]
  S --> V[Final validation]

• Never MDI alone. With financial labels it is almost uniform and biased by cardinality. If it is used, it should be with the target of a problem that has signal or as a diagnostic of the model, not of the data.
• MDA with the model's same honest validation (purging, embargo) and with the economic score (IC, Sharpe of the signal), not with raw accuracy.
• Groups before individual features when there are strong correlations, and always with noise controls.
• Standard error and multiple testing. An MDA with $t = 2$ among 22 features is not a finding.
• Selection overfits too. Choosing features by their importance over the whole sample and then validating on the same sample is a selection bias: the choice is made within each training window or a period that nobody has looked at is set aside (lessons 5.4.3 and 5.4.4).

Exercises

1. Repeat the MDI with the noise-free target and max_features=None (all features at each split). What happens to the features correlated with ret_1s?

todas = RandomForestRegressor(
    n_estimators=20, max_features=None, min_samples_leaf=1000,
    max_samples=0.2, n_jobs=4, random_state=0).fit(X, a)
comp = pd.DataFrame({
    "max_features=1": mdi["alpha"],
    "max_features=None": todas.feature_importances_ * 100})
comp.loc[["ret_1s", "ret_1m", "rsi_14", "valor_bp", "ey", "roe",
          "mom_12_1", "ruido_4"]].round(2)

          max_features=1  max_features=None
ret_1s             21.60              52.35
ret_1m              4.78               0.00
rsi_14              7.98               0.00
valor_bp            8.17               2.03
ey                  8.12               8.60
roe                 3.98               0.00
mom_12_1            6.11               1.84
ruido_4             1.90               0.00

2. Add revision_analistas (the leak from lesson 5.2.4) to the model and compute the MDA of four features with the real label. Does importance detect the leak?

X["revision_analistas"] = gauss("revision_analistas")
X["revision_analistas"] = X["revision_analistas"].fillna(0.0)
cols = list(X.columns)
mda_fuga, base_fuga = mda(y, [["revision_analistas"], ["ret_1s"],
                              ["valor_bp"], ["ruido_1"]])
print(f"out-of-sample IC with the leak: {base_fuga:.3f}")
mda_fuga.round(2)

out-of-sample IC with the leak: 0.157
                       mda    se
revision_analistas  155.61  2.89
ret_1s                0.68  0.53
valor_bp              0.53  0.24
ruido_1              -0.15  0.13

3. You have 300 candidate features for a weekly stock selection model and ten years of data. Describe a protocol for keeping a few of them without fooling yourself.