Cross-validation Revisited#
The previous chapter left us with a dilemma. Feature selection improves our models substantially, but it looks at the target, so it cannot be done on the data we cross-validate on. Our workaround was to spend a second set of 100 participants on it — which is precisely the kind of extravagance cross-validation was supposed to save us from.
The way out is to apply the same trick again: if a step must not see the test data, then perform it inside the cross-validation loop, separately in every round.
Feature selection inside the loop#
We revisit the SelectKBest example, but this time we unfold the cross-validation by hand, as we
did in chapter 3, and put the selection where it belongs.
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.feature_selection import SelectKBest, f_regression
from sklearn.model_selection import cross_val_predict, KFold
from sklearn.metrics import mean_absolute_error
DATA_URL = "https://raw.githubusercontent.com/pni-lab/predmod_lecture/master/ex_data/IXI/ixi.csv"
ixi = pd.read_csv(DATA_URL).sample(frac=1, random_state=42).reset_index(drop=True)
TARGET = 'Age'
FEATURES = [c for c in ixi.columns if c.endswith('_volume')]
df_cv = ixi.iloc[:100].reset_index(drop=True) # no external dataset this time
cv = KFold(n_splits=5)
cv_predictions = np.zeros(len(df_cv))
for train_index, test_index in cv.split(df_cv):
# 1. choose the best 5 features using THIS FOLD'S TRAINING DATA ONLY
selector = SelectKBest(score_func=f_regression, k=5).fit(
X=df_cv.loc[train_index, FEATURES], y=df_cv.loc[train_index, TARGET])
# 2. fit the model on those features
model = LinearRegression().fit(X=selector.transform(df_cv.loc[train_index, FEATURES]),
y=df_cv.loc[train_index, TARGET])
# 3. apply the SAME selection to the held-out fold (even if these are not its own top 5)
# and predict
cv_predictions[test_index] = model.predict(selector.transform(df_cv.loc[test_index, FEATURES]))
manual_mae = mean_absolute_error(df_cv[TARGET], cv_predictions)
print(f'feature selection inside the loop: MAE {manual_mae:.2f} years')
feature selection inside the loop: MAE 11.02 years
Note step 3 carefully. The selector was fitted on the training data, and we transform the test
data with it — we do not re-select on the test fold. The five features may well not be the best
five for those participants, and that is exactly the point: the model must be applied to new data
as it is, not re-tuned to it.
For comparison, in chapter 4 the same analysis with a separate 100-participant selection set gave 11.37 years. We have just matched it — slightly bettered it, in fact — with half as many participants, and the 100 people we saved can now go into training instead.
Pipelines#
The code above works but is verbose, and every line of hand-written fold bookkeeping is a chance
to make a mistake. Scikit-learn’s
Pipeline chains steps into a
single estimator that behaves like a model: when it is cross-validated, every step is refitted
on each training fold automatically.
from sklearn.pipeline import Pipeline
multistage_model = Pipeline(steps=[
('feature_selection', SelectKBest(score_func=f_regression, k=5)),
('regression_model', LinearRegression()),
])
cv_predictions = cross_val_predict(estimator=multistage_model, X=df_cv[FEATURES], y=df_cv[TARGET],
cv=KFold(n_splits=5))
print(f'the same thing, as a pipeline : '
f'{mean_absolute_error(df_cv[TARGET], cv_predictions):.2f} years')
print(f'identical to the manual version : '
f'{np.allclose(manual_mae, mean_absolute_error(df_cv[TARGET], cv_predictions))}')
the same thing, as a pipeline : 11.02 years
identical to the manual version : True
The same answer, in a handful of lines instead of a dozen.
Anything with fit and transform methods can be a pipeline step, and pipelines can be
arbitrarily long. A realistic one might look like this:
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.decomposition import PCA
from sklearn.linear_model import Ridge
realistic_model = Pipeline(steps=[
('imputation', SimpleImputer(strategy='median')), # fill in missing values
('standardization', StandardScaler()), # centre and scale
('dimensionality_reduction', PCA(n_components=10)), # compress correlated features
('regression_model', Ridge(alpha=10)), # regularized regression
])
predictions = cross_val_predict(realistic_model, X=df_cv[FEATURES], y=df_cv[TARGET], cv=KFold(5))
print(f'four-step pipeline: MAE {mean_absolute_error(df_cv[TARGET], predictions):.2f} years')
four-step pipeline: MAE 10.93 years
Every one of those four steps learns something from data — a median, a mean and standard deviation, a rotation, a set of coefficients — and every one of them would leak if it were fitted outside the loop. Inside the pipeline, all four are refitted on each fold’s training data and then merely applied to the held-out participants.
Tip
Put everything that learns from data inside the pipeline. This one habit prevents the majority
of leakage in practice. If you find yourself calling fit or fit_transform on your full
dataset, stop and ask what would happen to that step inside a fold.
What pipelines do not solve#
Look again at realistic_model. It contains four numbers we simply invented: the imputation
strategy, the number of components, the alpha, and implicitly the choice of Ridge over LASSO.
We know from the leakage page what happens if we try several and keep the best: the winner’s score is optimistic, because the held-out data helped choose it. The pipeline has protected every step within a round of cross-validation, but the choice of pipeline is made across rounds, using their pooled result.
So we need to apply the trick one more time — and that is nested cross-validation.
Exercise 5.1
Deliberately break the rule, in the way it usually gets broken: fit the selector once, on all
100 participants, before the loop, and then cross-validate a plain LinearRegression on the
five features it chose. How much does the MAE improve, and what have you actually measured?
selector = SelectKBest(score_func=f_regression, k=5).fit(X=df_cv[FEATURES], y=df_cv[TARGET])
leaky_predictions = cross_val_predict(LinearRegression(), X=selector.transform(df_cv[FEATURES]),
y=df_cv[TARGET], cv=KFold(5))
mean_absolute_error(df_cv[TARGET], leaky_predictions)
Solution
The error falls from 11.02 to 10.70 years. Every fold’s model has been handed features chosen with the help of the very participants it is about to predict, so what you have measured is how the procedure performs when it is allowed to peek — a quantity with no counterpart in the real world, where the next patient’s age is exactly what you do not know.
Note how undramatic the improvement is: a third of a year, on real data with real signal. That is what makes this error so durable in the literature. It does not produce an implausible result that someone would question — it produces a slightly better one.
Exercise 5.2
Build a pipeline that standardizes, applies PLS with 5 components, and cross-validate it. Then
move the StandardScaler outside: standardize df_cv[FEATURES] once, before cross-validation,
and compare. Is the difference large?
Solution
The difference is usually tiny — with 100 participants, the mean and standard deviation barely change when 20 of them are left out. This is worth knowing for two reasons. First, not all leakage is equally damaging, and scaling is at the mild end. Second, “the effect is small here” is something you can only say after measuring it, and the cost of doing it correctly is one line of code. Do it correctly.