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.