Feature Engineering · Notebook 8 companion articleFeature Selection: What to Keep and What to Drop

Start hereIntroduction#

Choosing features using the whole dataset and then cross-validating the result is leakage. It's the same rule that applied to imputers, scalers, encoders, target encoders and group aggregates in earlier notebooks. Notebook 7 showed a return-rate feature scoring 0.81 AUC on a target with no signal. Selecting features the wrong way does the same thing: a model reports a respectable score on data that can't be learned at all.

Before that, the notebook covers a quieter problem with no leakage involved. Several popular methods for judging a column's value each have a specific blind spot, and we'll see each one.

Section 1A dataset where we know the answer#

Feature selection needs a known right answer to grade a method. On real data you can see whether selection improved the score, but never whether it kept the right columns

The dataset has 1,000 rows, a balanced yes/no target and 30 features in four groups, named to match: info_0 to info_4 are truly related to the target. corr_0 to corr_4 are each a copy of the matching info_ column with noise added. const_0 to const_4 are constant or nearly constant. And fifteen noise_ columns carry nothing at all.

A good method should keep the five real signals, notice the five copies are redundant, and drop everything else. Every method in this notebook is graded against that standard.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.inspection import permutation_importance
from sklearn.feature_selection import (
    VarianceThreshold, SelectKBest, SelectFromModel, RFE,
    f_classif, mutual_info_classif,
)
from functools import partial

# mutual_info_classif adds tie-breaking noise and defaults to random_state=None,
# so an unpinned SelectKBest gives a different answer on every run.
mi_scorer = partial(mutual_info_classif, random_state=42)

plt.style.use('seaborn-v0_8-whitegrid')
pd.set_option('display.width', 200)
pd.set_option('display.max_columns', 40)

CV5 = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

print('Ready.')

Ready.

data = pd.read_csv('datasets/feature_engineering_data.csv')
X_all = data.drop(columns='target')
y = data['target'].values

GROUPS = {g: [c for c in X_all.columns if c.startswith(g)]
          for g in ['info', 'corr', 'const', 'noise']}

print('shape:', data.shape)
print('target balance:', np.bincount(y), f'({y.mean():.1%} positive)')
print()
for g, cols in GROUPS.items():
    print(f'  {g:6} {len(cols):2d} columns')
print()
print('how each corr_ column relates to its info_ partner:')
for i in range(5):
    r = np.corrcoef(X_all[f'info_{i}'], X_all[f'corr_{i}'])[0, 1]
    print(f'  corr(info_{i}, corr_{i}) = {r:+.4f}')

shape: (1000, 31)
target balance: [500 500] (50.0% positive)

  info    5 columns
  corr    5 columns
  const   5 columns
  noise  15 columns

how each corr_ column relates to its info_ partner:
  corr(info_0, corr_0) = +0.9543
  corr(info_1, corr_1) = +0.9597
  corr(info_2, corr_2) = +0.9569
  corr(info_3, corr_3) = +0.9525
  corr(info_4, corr_4) = +0.9598

The target is split exactly 50/50, so accuracy would be fine here, but ROC AUC is reported to stay consistent with the rest of the series. Each corr_ column has a correlation of about 0.95 to 0.96 with its partner. That's close enough to be redundant but not an exact duplicate, which is what you usually see in real, messy data.

Section 2Does dropping features actually help?#

Feature selection is often justified by "it reduces overfitting." Before learning any method, let's check that claim where we already know which columns are useless

def auc(cols, frame=None):
    frame = X_all if frame is None else frame
    est = Pipeline([('scale', StandardScaler()),
                    ('clf', LogisticRegression(max_iter=3000))])
    return cross_val_score(est, frame[cols], y, cv=CV5,
                           scoring='roc_auc', n_jobs=1).mean()


print(f'{"feature set":38}{"columns":>9}{"ROC AUC":>10}')
for label, cols in [
    ('everything',                 list(X_all.columns)),
    ('the 5 informative only',      GROUPS['info']),
    ('informative + their copies',  GROUPS['info'] + GROUPS['corr']),
    ('the 15 noise columns only',   GROUPS['noise']),
]:
    print(f'{label:38}{len(cols):9d}{auc(cols):10.4f}')

feature set                             columns   ROC AUC
everything                                   30    0.9959
the 5 informative only                        5    0.9999
informative + their copies                   10    0.9985
the 15 noise columns only                    15    0.5038

Three things to read from that table. The fifteen noise columns alone score 0.5038, which is chance, so the data is what it claims to be. Five features beat thirty, 0.9999 against 0.9959, so dropping twenty-five columns didn't just cost nothing, it made the model better. And the gain, +0.0040, is small, because logistic regression with L2 regularization already handles useless inputs fairly well. With less regularization, fewer rows or a harder problem, that gap gets much bigger.

So selection is worth doing here, and the goal is clear: find those five columns. But keep the real benefit in mind. In an actual project, the main value of feature selection usually isn't the fourth decimal place of AUC. It's a model that trains faster, is smaller to ship, has fewer things to monitor and fewer upstream dependencies to break, and is easier for a colleague to understand.

Section 3The free wins, and why the usual tool misses most of them#

The simplest filter needs no target and no model: a column that barely changes can't separate anything. The usual tool for this gets most of it wrong on this dataset

selector = VarianceThreshold(threshold=0.01).fit(X_all)
dropped = X_all.columns[~selector.get_support()].tolist()

print(f'VarianceThreshold(0.01) drops {len(dropped)}: {dropped}')
print(f'const_ columns that survive: '
      f'{[c for c in GROUPS["const"] if c not in dropped]}')
print()
print(f'{"column":10}{"variance":>12}{"distinct":>10}{"most common value":>20}{"its share":>11}')
for col in GROUPS['const']:
    s = X_all[col]
    top = s.value_counts()
    print(f'{col:10}{s.var():12.5f}{s.nunique():10d}{top.index[0]:20.3f}{top.iloc[0] / len(s):11.1%}')

VarianceThreshold(0.01) drops 3: ['const_0', 'const_1', 'const_2']
const_ columns that survive: ['const_3', 'const_4']

column        variance  distinct   most common value  its share
const_0        0.00000         1               1.000     100.0%
const_1        0.00000         1               1.000     100.0%
const_2        0.00991         2               0.000      99.0%
const_3        0.70468         2               5.000      97.1%
const_4        0.09760         2              -1.000      97.5%

Only three of the five near-constant columns get caught. const_3 is the interesting one that slips through: 97.1% of its rows have the same value, yet its variance is 0.70, seventy times the 0.01 cutoff. That's because variance measures spread (in the column's units, squared), not how concentrated the values are.

const_3 is 5.0 for 971 rows and far away for the other 29, and squaring those distances does the rest. A column can be 97% one value and still have a big variance if the other 3% are far enough out. That points to the real fix: if you mean "almost always the same value," measure that directly instead of using variance as a stand-in.

def dominant_share(frame):
    # Share of rows taken by each column's single most common value.
    return frame.apply(lambda s: s.value_counts(normalize=True).iloc[0])


shares = dominant_share(X_all)
quasi_constant = shares[shares > 0.95].sort_values(ascending=False)

print('columns where one value covers more than 95% of rows:')
print(quasi_constant.round(4).to_string())
print()
print(f'that finds {len(quasi_constant)} of the {len(GROUPS["const"])} const_ columns, '
      f'and nothing else')

columns where one value covers more than 95% of rows:
const_0    1.000
const_1    1.000
const_2    0.990
const_4    0.975
const_3    0.971

that finds 5 of the 5 const_ columns, and nothing else

All five, with no mistakes, from two lines of code. There's still a threshold (95% here), but it means the same thing on every column, whatever its units. VarianceThreshold has a second problem that matters more, because it fails silently: variance depends on units. A distance in meters has a variance a million times bigger than the same distance in kilometers (1,000 squared), and your threshold has no way of knowing which one it's looking at.

rescaled = X_all.copy()
rescaled['info_0'] = rescaled['info_0'] / 1000       # same column, different units

still_kept = VarianceThreshold(threshold=0.01).fit(rescaled)
now_dropped = rescaled.columns[~still_kept.get_support()].tolist()

print('after expressing info_0 in different units:')
print(f'  dropped: {now_dropped}')
print()
print(f'info_0 is the single most predictive column in the dataset '
      f'(|corr| with target = {abs(np.corrcoef(X_all["info_0"], y)[0, 1]):.4f})')

after expressing info_0 in different units:
  dropped: ['info_0', 'const_0', 'const_1', 'const_2']

info_0 is the single most predictive column in the dataset (|corr| with target = 0.5529)

info_0 is the strongest predictor in the whole file, and dividing it by 1,000 is enough to get it silently dropped. Its relationship with the target didn't change at all. Only its units did.

So the practical advice is narrower than usual. Use VarianceThreshold(0) to remove truly constant columns, since a column that never changes has zero variance in any units. For nearly-constant columns, measure the share of the most common value instead. And if you must use a variance threshold, scale first, and know that you're then thresholding on a number with no meaningful unit.

zero_variance = X_all.columns[X_all.var() == 0].tolist()
X = X_all.drop(columns=zero_variance)
print(f'dropped {zero_variance}, leaving {X.shape[1]} columns')

dropped ['const_0', 'const_1'], leaving 28 columns

Section 4Redundancy, and a heuristic that keeps the wrong column#

Two columns with the same information aren't twice as useful. For a linear model they actually hurt, and the standard fix has a flaw

When columns are highly correlated (multicollinearity), a linear model's coefficients become unstable and hard to interpret, because it can shift weight freely between two columns that say almost the same thing. The standard fix is to find pairs above a correlation cutoff and drop the weaker one in each pair, where "weaker" means less correlated with the target. That sounds obviously right.

corr_matrix = X.corr().abs()
target_corr = X.apply(lambda c: abs(np.corrcoef(c, y)[0, 1]))

to_drop = set()
cols = corr_matrix.columns.tolist()
for i in range(len(cols)):
    for j in range(i + 1, len(cols)):
        if corr_matrix.iloc[i, j] > 0.9:
            weaker = cols[i] if target_corr[cols[i]] < target_corr[cols[j]] else cols[j]
            to_drop.add(weaker)

print(f'pairs above |r| = 0.9, dropping the weaker of each: {sorted(to_drop)}')
print()
print('the two decisions worth looking at:')
for i in [3, 4]:
    print(f'  info_{i}: |corr with target| = {target_corr[f"info_{i}"]:.4f}')
    print(f'  corr_{i}: |corr with target| = {target_corr[f"corr_{i}"]:.4f}   <- kept')

pairs above |r| = 0.9, dropping the weaker of each: ['corr_0', 'corr_1', 'corr_2', 'info_3', 'info_4']

the two decisions worth looking at:
  info_3: |corr with target| = 0.2001
  corr_3: |corr with target| = 0.2052   <- kept
  info_4: |corr with target| = 0.1414
  corr_4: |corr with target| = 0.1415   <- kept

The rule drops info_3 and info_4, the real features, and keeps corr_3 and corr_4, which are those same features with noise added. It's not a coding bug. In this sample, corr_3 really does correlate slightly more with the target, 0.2052 against 0.2001, and corr_4 wins by even less, 0.1415 against 0.1414.

Those gaps are random noise, not signal, and the rule can't tell a real difference from a coin flip. For info_4, the whole decision comes down to the fourth decimal place.

kept = [c for c in X.columns if c not in to_drop]
print(f'before removal, {len(X.columns)} columns: {auc(list(X.columns), X):.4f}')
print(f'after removal,  {len(kept)} columns: {auc(kept, X):.4f}')

before removal, 28 columns: 0.9959
after removal,  23 columns: 0.9959

The score is identical to four decimal places, which makes sense: each column was swapped for a near-copy of itself, so the model loses almost nothing measurable. But it matters for everything else. The model now runs on corr_3 and corr_4, which in a real project would usually be derived fields rather than the actual measurement driving the outcome. Every explanation, every monitoring alert and every discussion about why the model behaves as it does now points at the wrong column.

Two habits help. When a correlated pair is basically tied, choose the column you'd rather depend on: the one that's cheaper to compute, more reliable upstream or easier to explain to whoever maintains it next. And unless you have a good reason, keep the original measurement over anything derived from it.

Section 5Univariate filters, and why the top five are not the best five#

The next group of methods scores each feature against the target on its own. They're fast and need no model, but they have one built-in blind spot

f_classif runs a one-way ANOVA on each column. For a yes/no target, that basically asks how far apart the two class averages are, compared with the spread inside each class.

f_scores = pd.Series(f_classif(X, y)[0], index=X.columns).sort_values(ascending=False)
print('top 10 by univariate F score:')
print(f_scores.head(10).round(1).to_string())
print()
print(f'best noise column: {f_scores[GROUPS["noise"]].max():.1f} '
      f'({f_scores[GROUPS["noise"]].idxmax()})')

top 10 by univariate F score:
info_0    439.5
corr_0    380.3
info_1    205.5
corr_1    193.8
info_2     80.4
corr_2     73.7
corr_3     43.8
info_3     41.6
corr_4     20.4
info_4     20.4

best noise column: 3.3 (noise_10)

The ranking gets one thing right: every real and near-real feature ranks above every noise column. The best noise column scores 3.3, against info_0's 439.5. Now look at the order: info_0, corr_0, info_1, corr_1, info_2. Take the top five and count how many different signals you actually get.

top5 = f_scores.head(5).index.tolist()
print(f'top 5 by F score : {top5}')
print(f'distinct underlying signals: {len({c.split("_")[1] for c in top5})}')
print()
print(f'{"feature set":34}{"columns":>9}{"ROC AUC":>10}')
print(f'{"top 5 by F score":34}{5:9d}{auc(top5, X):10.4f}')
print(f'{"the true info_ five":34}{5:9d}{auc(GROUPS["info"], X):10.4f}')
print(f'{"top 10 by F score":34}{10:9d}{auc(f_scores.head(10).index.tolist(), X):10.4f}')

top 5 by F score : ['info_0', 'corr_0', 'info_1', 'corr_1', 'info_2']
distinct underlying signals: 3

feature set                         columns   ROC AUC
top 5 by F score                          5    0.9674
the true info_ five                       5    0.9999
top 10 by F score                        10    0.9985

Three signals, not five, because two of the five spots went to near-copies of columns already picked. The score, 0.9674, is well behind the 0.9999 from the correct five, and even behind the top ten. This is the main weakness of scoring features one at a time: each one is scored alone, so it can't know that corr_0 says almost the same thing as info_0. Each looks great on its own, but together the pair adds little.

Two things follow. One-at-a-time filters work badly with redundant features, which describes most feature sets that come out of a notebook like the last one. And if you do use one, either remove redundancy first (as section 4 tried, not perfectly), or use it as a quick first cut before a method that looks at features together, not as the final decision.

Section 6Mutual information, and its blind spot#

Correlation and the F test both look for one kind of relationship. Mutual information asks a broader question, and that comes at a cost

A feature could be strongly related to the target in a way correlation and the F test can't see, like a U-shape where both very high and very low values point to the same class. Mutual information asks a broader question: how much does knowing this feature reduce uncertainty about the target? It assumes nothing about the shape of the relationship, and a true value of zero means the two are completely independent.

mi = pd.Series(mutual_info_classif(X, y, random_state=42),
               index=X.columns).sort_values(ascending=False)

comparison = pd.DataFrame({
    'mutual_info': mi.round(4),
    'F_score': f_scores.round(1),
    'abs_corr': target_corr.round(4),
}).loc[GROUPS['info'] + [mi[GROUPS['noise']].idxmax()]]
print(comparison.to_string())
print()
print(f'best noise column by MI : {mi[GROUPS["noise"]].max():.4f} '
      f'({mi[GROUPS["noise"]].idxmax()})')
print(f'worst info_ column by MI: {mi[GROUPS["info"]].min():.4f} '
      f'({mi[GROUPS["info"]].idxmin()})')

          mutual_info  F_score  abs_corr
info_0         0.1934    439.5    0.5529
info_1         0.0643    205.5    0.4132
info_2         0.0398     80.4    0.2731
info_3         0.0343     41.6    0.2001
info_4         0.0000     20.4    0.1414
noise_11       0.0222      0.2    0.0128


best noise column by MI : 0.0222 (noise_11)
worst info_ column by MI: 0.0000 (info_4)

info_4 scores exactly 0.0000, and a noise column scores 0.0222. So mutual information ranks a real feature below pure noise. info_4 is real but weak, with a correlation of only 0.1414 with the target.

Scikit-learn estimates mutual information using a nearest-neighbor method, which has its own random error. With 1,000 rows, a relationship this weak falls inside that error. The raw estimate comes out at zero or slightly below, and scikit-learn reports anything below zero as 0. noise_11 gets a small positive score for the same reason, random error in the other direction.

So mutual information catches relationships that correlation and F tests miss, which is a real advantage when relationships are curved. But it's an estimate, not an exact measurement, and its error is biggest where it matters most: on weak signals.

A zero from it doesn't prove independence. It means the estimator didn't find anything, which is a weaker claim. On the same weak column, the F test scored info_4 at 20.4, against the best noise column's 3.3. When a relationship is roughly a straight line, the simpler tool is more reliable.

One more practical detail matters more than it looks. mutual_info_classif adds a little random noise internally to break ties, and its random_state defaults to None. If you don't set it, you get a different answer every time you run it.

def mi_pipeline(scorer):
    return Pipeline([('scale', StandardScaler()),
                     ('select', SelectKBest(scorer, k=10)),
                     ('clf', LogisticRegression(max_iter=3000))])


unpinned = [cross_val_score(mi_pipeline(mutual_info_classif), X, y, cv=CV5,
                            scoring='roc_auc', n_jobs=1).mean() for _ in range(5)]
pinned = [cross_val_score(mi_pipeline(mi_scorer), X, y, cv=CV5,
                          scoring='roc_auc', n_jobs=1).mean() for _ in range(3)]

print('same code, five runs, random_state left at its default:')
print('  ', [f'{v:.4f}' for v in unpinned])
print(f'   spread: {max(unpinned) - min(unpinned):.4f}')
print()
print('three runs with random_state=42 pinned via mi_scorer:')
print('  ', [f'{v:.4f}' for v in pinned])
print(f'   spread: {max(pinned) - min(pinned):.4f}')

same code, five runs, random_state left at its default:
   ['0.9916', '0.9917', '0.9915', '0.9917', '0.9916']
   spread: 0.0002

three runs with random_state=42 pinned via mi_scorer:
   ['0.9915', '0.9915', '0.9915']
   spread: 0.0000

Without a fixed seed, the result changes in the fourth decimal place on every run, so nobody else can reproduce your exact numbers. Every use of mutual information in this notebook goes through a wrapper that fixes the seed. This applies beyond mutual information: any scorer, selector or model with a random_state left at its default makes results unrepeatable. In feature selection it hurts twice, because both the score and the list of kept columns can change between runs.

fig, axes = plt.subplots(1, 2, figsize=(15, 4.6))
palette = {'info': '#2b6cb0', 'corr': '#dd6b20', 'const': '#805ad5', 'noise': '#a0aec0'}

def colours_for(index):
    return [palette[c.split('_')[0]] for c in index]

ax = axes[0]
ordered = target_corr.sort_values(ascending=False)
ax.bar(range(len(ordered)), ordered.values, color=colours_for(ordered.index))
ax.set_xticks(range(len(ordered)))
ax.set_xticklabels(ordered.index, rotation=90, fontsize=7)
ax.set_ylabel('|correlation with target|')
ax.set_title('Linear association: the four groups separate cleanly', fontsize=10)

ax = axes[1]
ordered_mi = mi.sort_values(ascending=False)
ax.bar(range(len(ordered_mi)), ordered_mi.values, color=colours_for(ordered_mi.index))
ax.set_xticks(range(len(ordered_mi)))
ax.set_xticklabels(ordered_mi.index, rotation=90, fontsize=7)
ax.set_ylabel('mutual information')
ax.set_title('Mutual information: info_4 falls to zero, below several noise columns',
             fontsize=10)

handles = [plt.Rectangle((0, 0), 1, 1, color=v) for v in palette.values()]
axes[0].legend(handles, palette.keys(), fontsize=8, ncol=4)
plt.tight_layout()
plt.show()

Section 7Tree importance, and why the default one is misleading#

A random forest gives you a feature importance for free. It's the most widely used importance measure, and it has a well-known bias that this dataset shows clearly

Every split on a feature reduces impurity a little. Adding up those reductions across all the trees gives a score for each column, scaled so all the scores add up to 1.

forest = RandomForestClassifier(n_estimators=300, min_samples_leaf=2,
                                random_state=42, n_jobs=1).fit(X, y)
impurity = pd.Series(forest.feature_importances_, index=X.columns)

perm = permutation_importance(forest, X, y, n_repeats=10,
                              random_state=42, n_jobs=1)
permutation = pd.Series(perm.importances_mean, index=X.columns)

side_by_side = pd.DataFrame({'impurity': impurity.round(4),
                             'permutation': permutation.round(4)})
print(side_by_side.sort_values('impurity', ascending=False).head(12).to_string())
print()
print(f'{"":34}{"impurity":>10}{"permutation":>13}')
print(f'{"total over the 15 noise columns":34}'
      f'{impurity[GROUPS["noise"]].sum():10.4f}{permutation[GROUPS["noise"]].sum():13.4f}')
print(f'{"total over the 5 info columns":34}'
      f'{impurity[GROUPS["info"]].sum():10.4f}{permutation[GROUPS["info"]].sum():13.4f}')
print(f'{"grand total":34}{impurity.sum():10.4f}{permutation.sum():13.4f}')
print(f'{"noise as a share of the total":34}'
      f'{impurity[GROUPS["noise"]].sum() / impurity.sum():9.2%}'
      f'{permutation[GROUPS["noise"]].sum() / permutation.sum():13.2%}')

          impurity  permutation
info_0      0.2046       0.1567
corr_0      0.1647       0.0338
info_1      0.1069       0.0437
corr_1      0.1047       0.0357
info_2      0.0595       0.0190
corr_2      0.0555       0.0128
info_3      0.0376       0.0062
corr_3      0.0360       0.0044
info_4      0.0213       0.0055
corr_4      0.0211       0.0018
noise_11    0.0139       0.0000
noise_0     0.0138       0.0000

                                    impurity  permutation
total over the 15 noise columns       0.1872       0.0007
total over the 5 info columns         0.4299       0.2311
grand total                           1.0000       0.3203
noise as a share of the total        18.72%        0.22%

Impurity importance gives 18.7% of its total to columns that are pure noise. Permutation importance gives those same fifteen columns 0.0007 in total, which is 0.22% of its own total of 0.3203. The two totals can't be compared directly: impurity importances are scaled to add up to exactly 1, while permutation importances are raw drops in score that add up to whatever they add up to. Compare the shares, not the raw numbers.

The cause is built in. A continuous noise column offers the tree many possible split points, and with enough of them, some split will reduce impurity on the training rows by pure chance. That chance gets recorded as importance, even though it won't help on any new data. Impurity importance is calculated only on training data, so it rewards the ability to overfit.

Permutation importance asks a different question: shuffle this one column and see how much the score drops. One caution about the code above: it scores the forest on the same rows it was trained on. Here the noise columns still came out near zero, but that isn't guaranteed. A forest that memorized training rows partly through a noise column can lose some training accuracy when that column is shuffled, so training-set permutation importance can still give noise some credit.

To measure what actually helps on new data, compute permutation importance on a held-out validation set. There, a column the model only used for memorizing drops to about zero, or even below zero.

The correlated columns show a second, separate problem with impurity importance.

pairs = pd.DataFrame({
    'impurity': [impurity[f'{g}_{i}'] for g in ['info', 'corr'] for i in range(5)],
    'permutation': [permutation[f'{g}_{i}'] for g in ['info', 'corr'] for i in range(5)],
}, index=[f'{g}_{i}' for g in ['info', 'corr'] for i in range(5)]).round(4)
print(pairs.to_string())
print()
print(f'info_0 vs corr_0 under impurity   : {impurity["info_0"]:.4f} vs {impurity["corr_0"]:.4f}')
print(f'info_0 vs corr_0 under permutation: {permutation["info_0"]:.4f} vs {permutation["corr_0"]:.4f}')

        impurity  permutation
info_0    0.2046       0.1567
info_1    0.1069       0.0437
info_2    0.0595       0.0190
info_3    0.0376       0.0062
info_4    0.0213       0.0055
corr_0    0.1647       0.0338
corr_1    0.1047       0.0357
corr_2    0.0555       0.0128
corr_3    0.0360       0.0044
corr_4    0.0211       0.0018

info_0 vs corr_0 under impurity   : 0.2046 vs 0.1647
info_0 vs corr_0 under permutation: 0.1567 vs 0.0338

Under impurity, corr_0 scores 0.1647 against info_0's 0.2046, so they look like two separately strong features. The forest had to pick between the two at each relevant split and used both, which splits the credit roughly in proportion. Permutation importance shows the same pair as 0.1567 against 0.0338: shuffling corr_0 barely hurts the score, because info_0 is still there carrying the same information.

Both behaviors are worth knowing, and they point in different directions. Impurity importance overstates noise and splits credit between correlated columns, so don't read it as "how much this feature matters." Permutation importance answers "what would I lose if this column disappeared," which is usually the question you want.

But with correlated columns, it can mark both columns in a pair as unimportant, since each one covers for the other. Removing duplicates first (as in section 4), or dropping correlated columns as a group and re-measuring, fixes this. Permutation importance is also slower, needing a full scoring run per repeat per column. That doesn't matter with thirty columns but adds up with thousands.

fig, axes = plt.subplots(1, 2, figsize=(15, 5.4), sharey=True)
order = impurity.sort_values(ascending=True).index

for ax, series, title in [
    (axes[0], impurity,
     f'Impurity importance\n(noise takes '
     f'{impurity[GROUPS["noise"]].sum() / impurity.sum():.1%} of the total)'),
    (axes[1], permutation,
     f'Permutation importance\n(noise takes '
     f'{permutation[GROUPS["noise"]].sum() / permutation.sum():.2%} of the total)'),
]:
    ax.barh(range(len(order)), series[order].values, color=colours_for(order))
    ax.set_yticks(range(len(order)))
    ax.set_yticklabels(order, fontsize=7)
    ax.set_title(title, fontsize=10)
    ax.set_xlabel('importance')

plt.tight_layout()
plt.show()

Section 8L1, where selection happens during training#

Every method so far decides before or after training. An L1 penalty decides during training, so selecting and fitting become one step

An L1 penalty is based on the sum of the coefficients' absolute sizes. It pushes weights toward zero, and unlike an L2 penalty, it can push them all the way to exactly zero. A feature with a zero coefficient is effectively not in the model. In scikit-learn's logistic regression, C is the inverse of the penalty strength, so a small C means a strong penalty and few surviving features.

A version note: the code uses LogisticRegression(l1_ratio=1, solver='liblinear'), which means pure L1 in scikit-learn 1.8 and later, where the penalty argument is deprecated. On older versions, l1_ratio is ignored unless penalty='elasticnet', so the same line quietly fits an ordinary L2 model and no coefficients reach zero. On those versions, write penalty='l1' instead.

X_scaled = StandardScaler().fit_transform(X)   # L1 is meaningless on unscaled columns

print(f'{"C":>7}{"kept":>7}{"info":>7}{"corr":>7}{"noise":>7}{"const":>7}   surviving info_ columns')
for C in [0.01, 0.03, 0.1, 0.3, 1.0, 3.0]:
    model = LogisticRegression(l1_ratio=1, solver='liblinear', C=C,
                               max_iter=3000, random_state=42).fit(X_scaled, y)
    alive = X.columns[model.coef_[0] != 0]
    counts = {g: sum(1 for c in alive if c.startswith(g))
              for g in ['info', 'corr', 'noise', 'const']}
    print(f'{C:7.2f}{len(alive):7d}{counts["info"]:7d}{counts["corr"]:7d}'
          f'{counts["noise"]:7d}{counts["const"]:7d}   '
          f'{sorted(c for c in alive if c.startswith("info"))}')

      C   kept   info   corr  noise  const   surviving info_ columns
   0.01      5      4      1      0      0   ['info_0', 'info_1', 'info_2', 'info_3']
   0.03      6      5      1      0      0   ['info_0', 'info_1', 'info_2', 'info_3', 'info_4']
   0.10      8      5      2      1      0   ['info_0', 'info_1', 'info_2', 'info_3', 'info_4']
   0.30     12      5      1      5      1   ['info_0', 'info_1', 'info_2', 'info_3', 'info_4']
   1.00     18      5      1     11      1   ['info_0', 'info_1', 'info_2', 'info_3', 'info_4']
   3.00     24      5      3     14      2   ['info_0', 'info_1', 'info_2', 'info_3', 'info_4']

The results are easy to read, and both ends tell you something. At C = 0.03, the selection is exactly right: all five info_ columns, one corr_ column and no noise. Tighten to C = 0.01 and it starts dropping real signal, losing info_4, the weakest of the five. Loosen it and noise creeps back: one noise column at C = 0.1, five at 0.3 and eleven at 1.0.

Notice also that L1 mostly keeps just one column from each correlated pair. At C = 0.03, only one of the five copies survives. Given two nearly identical columns, L1 tends to put the weight on one and push the other to zero, because paying the penalty twice for the same information is wasteful. That's the opposite of impurity importance, which split credit between the pair, and it's usually what you want from a selector.

The catch is that which column of the pair survives is fairly random and unstable. Resample the data and L1 might swap info_3 for corr_3 next time. If you need stable selection across correlated groups, that's what elastic net is for. Also, C controls both the penalty strength and, as a side effect, how many features survive. The honest way to set it is to tune it against the cross-validated score, not to pick a number on the spot.

Section 9Recursive elimination, and the payoff for thinking jointly#

Section 5 left a question open. What does a method that judges features together, not one at a time, do with the same redundancy problem?

Recursive feature elimination (RFE) fits the model, drops the least useful column, and repeats on what's left. Because the model is refit after every removal, each remaining column is judged alongside all the others, not alone. Once info_0 is in the kept set, corr_0 has nothing left to add, so it gets removed in a later round.

RFE decides which column is "least useful" from the size of each coefficient. That's only a fair comparison when the columns are on the same scale, which is why the code runs it on X_scaled.

rfe = RFE(LogisticRegression(max_iter=3000), n_features_to_select=5).fit(X_scaled, y)
rfe_choice = X.columns[rfe.support_].tolist()

print(f'RFE keeping 5      : {sorted(rfe_choice)}')
print(f'distinct signals   : {len({c.split("_")[1] for c in rfe_choice})}')
print()
print(f'univariate top 5   : {sorted(top5)}')
print(f'distinct signals   : {len({c.split("_")[1] for c in top5})}')
print()
print(f'{"selection":26}{"ROC AUC":>10}')
print(f'{"univariate top 5":26}{auc(top5, X):10.4f}')
print(f'{"RFE top 5":26}{auc(rfe_choice, X):10.4f}')
print(f'{"the true answer":26}{auc(GROUPS["info"], X):10.4f}')

RFE keeping 5      : ['info_0', 'info_1', 'info_2', 'info_3', 'info_4']
distinct signals   : 5

univariate top 5   : ['corr_0', 'corr_1', 'info_0', 'info_1', 'info_2']
distinct signals   : 3

selection                    ROC AUC
univariate top 5              0.9674
RFE top 5                     0.9999
the true answer               0.9999

On this dataset, RFE finds info_0 to info_4 exactly, all five real signals, and scores 0.9999, matching the right answer. That's the case for methods that look at features together: the one-at-a-time filter used the same five slots and reached only 0.9674, because it couldn't account for what it had already picked.

But don't over-generalize from one dataset. Because RFE judges each column alongside the others, it's much more likely to find the truly separate signals. It's not guaranteed on every dataset: how much it helps depends on how redundant the columns are and how the elimination plays out. What is guaranteed is that RFE can't fall into section 5's specific trap, because judging features together is the whole point of refitting after each removal.

The cost is computing time. RFE refits the model once per removal, so starting from 1,000 features means about 1,000 fits, and RFECV multiplies that by the number of folds. The usual shortcut is to filter cheaply down to a smaller set first, then run RFE on that. Setting step to a fraction, like 0.1 to drop 10% of the remaining columns each round instead of one at a time, also helps a lot.

RFE doesn't tell you how many features to keep. Here we told it five because we already knew the answer. RFECV picks the number using cross-validation, which is the right approach when you don't know. The next section explains why that step needs real care.

Section 10The mistake that makes all of this worthless#

Most methods in this notebook look at the target. That one difference decides which of them can leak

Section 3's two filters, VarianceThreshold and the most-common-value share, only use the feature columns and never touch y. A selector that never sees the target can't leak it, so those two are safe to run before splitting. Everything from section 4 on, including the correlation tie-break, learns from the labels. So notebook 3's rule applies: anything that learns from data must learn only from training rows.

If you select features using the whole dataset and then cross-validate, the test folds have helped choose their own features. The cross-validation score no longer measures anything real. The clearest way to see how bad this gets is to run the same steps on data built to have no signal at all.

rng = np.random.default_rng(0)
n_rows, n_features, keep = 200, 2000, 20

X_noise = rng.normal(size=(n_rows, n_features))
y_noise = rng.integers(0, 2, n_rows)

print(f'{n_rows} rows, {n_features} pure-noise features, a coin-flip target.')
print('The correct answer is 0.5000 and nothing can beat it.')
print()

# The wrong way: rank all features against the target, then cross-validate the winners.
scores = np.abs([np.corrcoef(X_noise[:, j], y_noise)[0, 1] for j in range(n_features)])
winners = np.argsort(scores)[::-1][:keep]
leaky = cross_val_score(LogisticRegression(max_iter=3000), X_noise[:, winners],
                        y_noise, cv=CV5, scoring='accuracy', n_jobs=1).mean()

# The right way: selection is a pipeline step, refitted inside every fold.
honest_pipe = Pipeline([('select', SelectKBest(f_classif, k=keep)),
                        ('clf', LogisticRegression(max_iter=3000))])
honest = cross_val_score(honest_pipe, X_noise, y_noise, cv=CV5,
                         scoring='accuracy', n_jobs=1).mean()

print(f'select on all data, then cross-validate : {leaky:.4f}')
print(f'select inside each fold                 : {honest:.4f}')
print()
print(f'strongest correlation found among {n_features} noise columns: {scores.max():.4f}')

200 rows, 2000 pure-noise features, a coin-flip target.
The correct answer is 0.5000 and nothing can beat it.

select on all data, then cross-validate : 0.7450
select inside each fold                 : 0.4750

strongest correlation found among 2000 noise columns: 0.2343

74.5% accuracy on a coin flip. There's nothing to learn in that data, yet the procedure reports a model that looks respectable. Done properly, the same pipeline returns 0.4750, which is chance plus normal fold variation. That's the truth.

Here's how it happens. With 2,000 columns and only 200 rows, some columns will correlate with the target by pure luck. Here the best one reaches 0.2343. Ranking all 2,000 against the target picks exactly those lucky columns. Their luck was measured on all 200 rows, so it carries into every test fold. The model isn't learning a pattern. It's re-reading an accident it already saw during selection.

That also predicts how the problem grows: the more columns you search, the luckier the best one looks.

print(f'{"features searched":>19}{"leaky":>9}{"honest":>9}{"best |corr|":>13}')
sizes = [50, 200, 1000, 5000]
leaky_curve, honest_curve = [], []
for p in sizes:
    Xp = rng.normal(size=(n_rows, p))
    yp = rng.integers(0, 2, n_rows)
    s = np.abs([np.corrcoef(Xp[:, j], yp)[0, 1] for j in range(p)])
    w = np.argsort(s)[::-1][:keep]
    a = cross_val_score(LogisticRegression(max_iter=3000), Xp[:, w], yp,
                        cv=CV5, scoring='accuracy', n_jobs=1).mean()
    b = cross_val_score(Pipeline([('select', SelectKBest(f_classif, k=keep)),
                                  ('clf', LogisticRegression(max_iter=3000))]),
                        Xp, yp, cv=CV5, scoring='accuracy', n_jobs=1).mean()
    leaky_curve.append(a); honest_curve.append(b)
    print(f'{p:>19}{a:9.4f}{b:9.4f}{s.max():13.4f}')

  features searched    leaky   honest  best |corr|
                 50   0.5900   0.5300       0.1422
                200   0.6550   0.4950       0.2139
               1000   0.7650   0.4850       0.2223
               5000   0.8200   0.5150       0.2968

From 0.59 to 0.82 as the search gets wider, while the honest score stays at chance. The size of the lie grows with how many features were searched. That's why this failure is becoming more common: automated feature construction, covered at the end of notebook 7, makes very wide feature sets normal.

Now the part that makes this more than a simple warning. Run the same comparison on the real dataset used throughout this notebook.

f_all = pd.Series(f_classif(X, y)[0], index=X.columns)
top10_all = f_all.nlargest(10).index.tolist()

leaky_real = auc(top10_all, X)
honest_real = cross_val_score(
    Pipeline([('scale', StandardScaler()),
              ('select', SelectKBest(f_classif, k=10)),
              ('clf', LogisticRegression(max_iter=3000))]),
    X, y, cv=CV5, scoring='roc_auc', n_jobs=1).mean()

print(f'select on all data, then cross-validate : {leaky_real:.4f}')
print(f'select inside each fold                 : {honest_real:.4f}')
print(f'difference                              : {leaky_real - honest_real:+.4f}')

select on all data, then cross-validate : 0.9985
select inside each fold                 : 0.9985
difference                              : +0.0000

No difference at all. On this dataset, the leak costs nothing. That's not a contradiction, and the reason is the useful part. Look at what each fold actually selects.

def selections_per_fold(frame, labels, k, splitter):
    picks = []
    for train_idx, _ in splitter.split(frame, labels):
        sel = SelectKBest(f_classif, k=k).fit(frame.iloc[train_idx], labels[train_idx])
        picks.append(set(frame.columns[sel.get_support()]))
    return picks


real_picks = selections_per_fold(X, y, 10, CV5)
print('real dataset, 10 features chosen per fold:')
print(f'  identical in all 5 folds : {len(set.intersection(*real_picks))}')
print(f'  ever chosen at all       : {len(set.union(*real_picks))}')

noise_frame = pd.DataFrame(X_noise, columns=[f'n{i}' for i in range(n_features)])
noise_picks = selections_per_fold(noise_frame, y_noise, keep, CV5)
print()
print('pure-noise problem, 20 features chosen per fold:')
print(f'  identical in all 5 folds : {len(set.intersection(*noise_picks))}')
print(f'  ever chosen at all       : {len(set.union(*noise_picks))}')

real dataset, 10 features chosen per fold:
  identical in all 5 folds : 10
  ever chosen at all       : 10

pure-noise problem, 20 features chosen per fold:
  identical in all 5 folds : 1
  ever chosen at all       : 68

Four numbers tell the whole story. On the real data, all five folds pick the same ten columns. The signal is strong enough that the choice was never in doubt, so selecting before the split changed nothing. On the noise data, only one column is picked in all five folds, and 68 different columns are picked by at least one. The selection there is basically random, which is exactly when seeing the answer in advance helps a lot.

fig, axes = plt.subplots(1, 2, figsize=(14, 4.4))

ax = axes[0]
ax.plot(sizes, leaky_curve, 'o-', color='#c53030', lw=2, label='selected on all data')
ax.plot(sizes, honest_curve, 's-', color='#2b6cb0', lw=2, label='selected inside each fold')
ax.axhline(0.5, color='#666666', ls='--', lw=1.4, label='the truth (0.5)')
ax.set_xscale('log')
ax.set_xlabel('number of noise features searched')
ax.set_ylabel('cross-validated accuracy')
ax.set_title('The wider you search, the bigger the illusion', fontsize=10)
ax.legend(fontsize=8)

ax = axes[1]
labels = ['real data\n(stable selection)', 'pure noise\n(unstable selection)']
stable = [len(set.intersection(*real_picks)) / len(set.union(*real_picks)),
          len(set.intersection(*noise_picks)) / len(set.union(*noise_picks))]
gaps = [leaky_real - honest_real, leaky - honest]
x = np.arange(2)
ax.bar(x - 0.2, stable, 0.4, color='#2b6cb0', label='selection agreement across folds')
ax.bar(x + 0.2, gaps, 0.4, color='#c53030', label='score inflation from leaking')
ax.set_xticks(x); ax.set_xticklabels(labels, fontsize=9)
ax.axhline(0, color='#666666', lw=1)
ax.set_title('Stable selection, small leak. Unstable selection, large leak.', fontsize=10)
ax.legend(fontsize=8)

plt.tight_layout()
plt.show()

The more unstable the selection, the bigger the leak. We only have two cases here, so this isn't an exact formula, but it's a useful check. If refitting the selector on each fold gives the same columns every time, the selection is stable and the leak is small, even if you broke the rule. If the chosen columns change a lot between folds, you're in the danger zone, and the reported score is closer to fiction than fact.

This also answers the fair objection that doing selection properly is expensive. Sometimes it is, and the fold-agreement check above tells you whether skipping the proper method is a safe risk.

But there's a catch: we could only show the leak was harmless on the real data by doing the selection properly first and comparing. You can't reach that conclusion from inside the leaky version alone, which echoes notebook 6's point about vectorizer vocabularies: a leak's size depends on the setup, so the proper version is the only reliable measure. So the rule stands: put the selector inside the pipeline. It's one line of code.

Section 11An honest comparison#

With the rule in place, here's every method from this notebook tested the same way: the selector sits inside the pipeline and is refit on each fold's training rows

def evaluate(estimator, label):
    score = cross_val_score(estimator, X, y, cv=CV5, scoring='roc_auc', n_jobs=1).mean()
    estimator.fit(X, y)
    n_kept = estimator[:-1].transform(X).shape[1]
    return {'method': label, 'ROC AUC': round(score, 4), 'features kept': n_kept}


S = StandardScaler
LR = lambda: LogisticRegression(max_iter=3000)

table = [
    evaluate(Pipeline([('s', S()), ('m', LR())]),
             'no selection'),
    evaluate(Pipeline([('s', S()), ('sel', SelectKBest(f_classif, k=5)), ('m', LR())]),
             'univariate F, k=5'),
    evaluate(Pipeline([('s', S()), ('sel', SelectKBest(f_classif, k=10)), ('m', LR())]),
             'univariate F, k=10'),
    evaluate(Pipeline([('s', S()), ('sel', SelectKBest(mi_scorer, k=10)), ('m', LR())]),
             'mutual information, k=10'),
    evaluate(Pipeline([('sel', SelectFromModel(
                RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=1),
                max_features=10, threshold=-np.inf)), ('s', S()), ('m', LR())]),
             'random forest importance, top 10'),
    evaluate(Pipeline([('s', S()), ('sel', SelectFromModel(
                LogisticRegression(l1_ratio=1, solver='liblinear', C=0.03,
                                   max_iter=3000, random_state=42))), ('m', LR())]),
             'L1, C=0.03'),
    evaluate(Pipeline([('s', S()), ('sel', RFE(LR(), n_features_to_select=5)), ('m', LR())]),
             'RFE, k=5'),
]
print(pd.DataFrame(table).to_string(index=False))
print()
print(f'for reference, the true five columns score {auc(GROUPS["info"], X):.4f}')

                          method  ROC AUC  features kept
                    no selection   0.9959             28
               univariate F, k=5   0.9674              5
              univariate F, k=10   0.9985             10
        mutual information, k=10   0.9915             10
random forest importance, top 10   0.9985             10
                      L1, C=0.03   0.9996              6
                        RFE, k=5   0.9999              5

for reference, the true five columns score 0.9999

One thing to own up to first. The comparisons in sections 4, 5 and 9 picked their columns using all 1,000 rows and then cross-validated them, which is exactly the leak from section 10. In the table above, selection happens inside each fold, and univariate F with k=5 and RFE still score 0.9674 and 0.9999, the same as before. Those earlier conclusions hold because selection on this data is stable, as section 10's fold check showed.

Ranked honestly, the results match what earlier sections predicted. RFE hits the right answer exactly: 0.9999 with the right five columns. L1 is a close second at 0.9996 with six columns, and needs just one selection fit instead of RFE's 23 rounds of refitting.

These are the only two methods on the table that are pushed to avoid redundancy. RFE refits after each removal, and L1 makes every coefficient compete for one shared penalty budget. That's why they're also the only two that avoided collecting redundant copies.

The one-at-a-time filters below them score each column alone and can't do that. The random forest does use features together, but its impurity scores split credit between copies instead of punishing the second one, so on the full data its top ten by impurity is the five real columns plus all five copies.

Univariate F with k=5 is the worst method on the table, at 0.9674, worse even than no selection, for the redundancy reason from section 5. The same filter with k=10 is fine, because ten slots are enough to cover all five real signals even with half of them wasted on copies. Mutual information scores 0.9915 against the F test's 0.9985 at the same k, most likely because of the random error from section 6.

Those are the only two methods that score below no selection (0.9959). Every other method beats it while keeping at most ten of the twenty-eight columns, and the best two use only five and six.

Section 12Two good methods, two different answers#

A final warning. It's tempting to treat whatever a selector returns as the definitive list of important features. Run several at once and that idea falls apart

l1_model = LogisticRegression(l1_ratio=1, solver='liblinear', C=0.03,
                              max_iter=3000, random_state=42).fit(X_scaled, y)
picks = {
    'univariate F': set(f_all.nlargest(6).index),
    'mutual information': set(mi.nlargest(6).index),
    'forest impurity': set(impurity.nlargest(6).index),
    'permutation': set(permutation.nlargest(6).index),
    'L1 (C=0.03)': set(X.columns[l1_model.coef_[0] != 0]),
}
for name, chosen in picks.items():
    hit = sorted(c for c in chosen if c.startswith('info'))
    print(f'  {name:20} {len(hit)}/5 info columns   {sorted(chosen)}')

agreed = set.intersection(*picks.values())
print()
print(f'chosen by every method : {sorted(agreed)}')
print(f'chosen by at least one : {len(set.union(*picks.values()))} columns')

  univariate F         3/5 info columns   ['corr_0', 'corr_1', 'corr_2', 'info_0', 'info_1', 'info_2']
  mutual information   4/5 info columns   ['corr_0', 'corr_1', 'info_0', 'info_1', 'info_2', 'info_3']
  forest impurity      3/5 info columns   ['corr_0', 'corr_1', 'corr_2', 'info_0', 'info_1', 'info_2']
  permutation          3/5 info columns   ['corr_0', 'corr_1', 'corr_2', 'info_0', 'info_1', 'info_2']
  L1 (C=0.03)          5/5 info columns   ['corr_3', 'info_0', 'info_1', 'info_2', 'info_3', 'info_4']

chosen by every method : ['info_0', 'info_1', 'info_2']
chosen by at least one : 9 columns

Five methods, and the pattern is more interesting than simple disagreement. Three of them, univariate F, forest impurity and permutation, pick exactly the same six columns: the three strongest signals plus all three of their noisy copies. They rank those six very differently, as section 7 showed, but the set is the same.

Mutual information swaps one copy for info_3. Only L1 finds all five real columns, because it's the only one of these five methods that's penalized for keeping a redundant column. (RFE, left out of this comparison, found them too in section 9.)

So five methods give three different answers, and only info_0, info_1 and info_2, the three strongest signals, are picked by all of them.

None of the five methods is broken. They're answering different questions, and a column can be important by one definition and useless by another. The correlated pairs in this dataset are the clearest example of that.

In practice, agreement between methods means something. A column every method keeps is safe to keep, and one every method drops is safe to drop. Disagreement is a prompt to look closer, not something to settle by picking a favorite method. It usually means the column is redundant with something already chosen, genuinely weak, or related to the target in a curved way that one method can't see.

In the end, the only real judge is the cross-validated score on the metric that matters for your problem, with the selector inside the pipeline. Everything else in this notebook is a way to generate good candidates for that final test, not a replacement for it.

One last caution. If you try many selectors and keep whichever scores best in cross-validation, that best score is itself a little optimistic, for the same reason as section 10: you chose it by looking at those folds. A final check on data you set aside at the start, or nested cross-validation, gives the honest number.