05 - Cross-Validation and Model Selection

You have trained a model and want to know how much it will get wrong on new data. The most immediate way is to set aside 20% of the examples and measure the error on them; the defect is that this number depends on which examples happened to land in the 20%. Cross-validation corrects the defect by repeating the measurement on several partitions, but it has its own traps: the order of the data, related examples, time, and the temptation to use the same score both to choose the model and to declare its quality. This article shows with reproducible experiments where each of these traps bites and which split avoids it. It starts from the ideas of chapters 01 (baseline, test set, leakage) and 04 (bias and variance of a model).

The examples run with scikit-learn 1.9.0 and numpy 2.4.6 on bundled datasets (Breast Cancer, 569 examples; Iris, 150) or on synthetic data with a fixed seed. Every output was obtained twice with the same result.

Why a single partition is not enough

K-fold cross-validation divides the data into K blocks (folds) of similar size; in turn each block acts as the test set and the other K − 1 as the training set; the final estimate is the mean of the K scores. With K equal to the number of examples you get the leave-one-out variant. The textbook by James, Witten, Hastie and Tibshirani (An Introduction to Statistical Learning) describes the method in section 5.1.3 and discusses the choice of K in 5.1.4: with K equal to 5 or 10 the training sets contain about (K − 1)/K of the data, so the estimate is only slightly pessimistic, whereas leave-one-out has minimal bias but higher variance, because averaging models trained on almost identical sets gives highly correlated results. The authors conclude that 5 or 10 are choices that empirically suffer neither from excessive bias nor from high variance.

The listing implements K-fold by hand, compares it with cross_val_score, and measures how much an estimate varies when only the partition changes. The model is the standardised logistic regression of chapter 03, in the Pipeline that recomputes the standardisation inside each fold.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000))

# K-fold by hand: each block is the test set once, the other K-1 are the training set
def manual_kfold(K, seed):
    idx = np.random.RandomState(seed).permutation(len(y))      # same permutation as KFold(shuffle=True)
    blocks = np.array_split(idx, K)
    acc = []
    for k in range(K):
        train = np.concatenate([b for j, b in enumerate(blocks) if j != k])
        acc.append(model.fit(X[train], y[train]).score(X[blocks[k]], y[blocks[k]]))
    return np.mean(acc)

manual = manual_kfold(5, 0)
sk = cross_val_score(model, X, y, cv=KFold(5, shuffle=True, random_state=0)).mean()
print("5-fold by hand:", round(manual, 4), " cross_val_score:", round(sk, 4), " equal:", np.isclose(manual, sk))

# how much does the estimate vary when only the partition changes?
hold = []
for s in range(100):
    Xa, Xb, ya, yb = train_test_split(X, y, test_size=0.2, random_state=s)
    hold.append(model.fit(Xa, ya).score(Xb, yb))
cv5 = [cross_val_score(model, X, y, cv=KFold(5, shuffle=True, random_state=s)).mean() for s in range(100)]
print(f"hold-out 80/20, 100 partitions: min={min(hold):.3f} max={max(hold):.3f} std={np.std(hold):.4f}")
print(f"5-fold, 100 partitions:         min={min(cv5):.3f} max={max(cv5):.3f} std={np.std(cv5):.4f}")

# Output:
# 5-fold by hand: 0.9719  cross_val_score: 0.9719  equal: True
# hold-out 80/20, 100 partitions: min=0.947 max=1.000 std=0.0138
# 5-fold, 100 partitions:         min=0.967 max=0.984 std=0.0037

The hand-made version and the scikit-learn one coincide (0.9719) because they use the same permutation, with blocks of equal size to within one. The result that matters is the last one: repeating an 80/20 split 100 times, the estimated accuracy goes from 0.947 to 1.000, with a standard deviation of 0.0138; repeating a 5-fold 100 times, from 0.967 to 0.984, with a standard deviation of 0.0037. With a single split you could declare "100%" or "94.7%" for the very same model, depending on the seed. The 5-fold uses all 569 examples in turn in both training and test and averaging reduces the standard deviation by about 3.7 times in this experiment; the cost is doing K fits instead of one.

When naive K-fold misleads

K-fold gives an honest estimate only if each block resembles "new data" as much as future data will. Three situations violate the assumption, and each has a splitting scheme designed for it.

Sorted data and rare classes: stratify

The Iris dataset has its 150 rows sorted by species: the first 50 are all of species 0, the last 50 of species 2. A 3-block K-fold without shuffling takes each block of 50 rows as the test set: the test block contains a species that the model, trained on the other two, has never seen.

import numpy as np
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, StratifiedKFold, cross_val_score

X, y = load_iris(return_X_y=True)                  # the 150 rows are sorted by species
m = LogisticRegression(max_iter=1000)
print("labels in the first/last rows:", y[:3].tolist(), y[-3:].tolist())
print("KFold(3) without shuffling:   ", cross_val_score(m, X, y, cv=KFold(3)).round(3).tolist())
print("KFold(3) shuffling:           ", cross_val_score(m, X, y, cv=KFold(3, shuffle=True, random_state=0)).round(3).tolist())
print("StratifiedKFold(3):           ", cross_val_score(m, X, y, cv=StratifiedKFold(3)).round(3).tolist())
print("cv=3 (default for classifiers):", cross_val_score(m, X, y, cv=3).round(3).tolist())

# Output:
# labels in the first/last rows: [0, 0, 0] [2, 2, 2]
# KFold(3) without shuffling:    [0.0, 0.0, 0.0]
# KFold(3) shuffling:            [0.98, 0.96, 0.96]
# StratifiedKFold(3):            [0.98, 0.96, 0.98]
# cv=3 (default for classifiers): [0.98, 0.96, 0.98]

The result is 0.0 on all three folds: the model gets every example wrong because it does not know the class it has to predict. By shuffling (shuffle=True) or using StratifiedKFold, which keeps in every fold the same proportion of classes as the dataset, the accuracy goes back to 0.96-0.98. The scikit-learn documentation specifies that when cv is an integer, cross_val_score already uses StratifiedKFold for classifiers (last line), whereas a KFold passed explicitly overrides this choice: the error in the first line comes from the KFold without shuffling. Stratification matters above all with rare classes, where a fold might not contain a single example of one.

Related examples: group

If the dataset contains several rows for the same person, machine or customer, a random K-fold puts rows of the same person in both training and test, and the model can recognise the person instead of learning the phenomenon. The listing builds the extreme case: 40 patients with 10 visits each, a personal trait that repeats in every visit, and an outcome chosen at random per patient, so with no link to the measurements. The honest value is 0.5.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GroupKFold, KFold, cross_val_score

rng = np.random.default_rng(0)
n_patients, n_visits = 40, 10
fingerprint = rng.standard_normal((n_patients, 8))                 # trait specific to each patient
X = np.repeat(fingerprint, n_visits, axis=0) + 0.3 * rng.standard_normal((n_patients * n_visits, 8))
outcome = rng.integers(0, 2, n_patients)                           # outcome per patient, chosen at random: no link with the measurements
y = np.repeat(outcome, n_visits)
groups = np.repeat(np.arange(n_patients), n_visits)

rf = RandomForestClassifier(n_estimators=100, random_state=0)
random_cv = cross_val_score(rf, X, y, cv=KFold(5, shuffle=True, random_state=0)).mean()
by_patient = cross_val_score(rf, X, y, cv=GroupKFold(5), groups=groups).mean()
print("400 visits of 40 patients, random outcome (honest limit: 0.5)")
print("random KFold:", round(random_cv, 3))
print("GroupKFold by patient:", round(by_patient, 3))

# Output:
# 400 visits of 40 patients, random outcome (honest limit: 0.5)
# random KFold: 0.95
# GroupKFold by patient: 0.508

With the random K-fold the estimated accuracy is 0.95 on a problem in which nothing can be predicted: the model recognised the patients by the personal trait and remembered their outcome. With GroupKFold, which keeps all the visits of a patient in the same fold and receives the identifier through groups=, the estimate goes back to 0.508. The scikit-learn documentation describes the same scenario, subjects with several samples and models flexible enough to learn personal traits. The defect is the same as the leakage of chapter 01: information that you would not have in production, because new patients have no visits in the training set.

Data over time: do not look at the future

When the rows are ordered in time and you want to predict the future, a test block in the past with respect to the training set is a rigged exam. TimeSeriesSplit makes every test set always come after its training set.

import numpy as np
from sklearn.model_selection import TimeSeriesSplit

for k, (tr, te) in enumerate(TimeSeriesSplit(n_splits=4).split(np.arange(12)), 1):
    print(f"split {k}: train = {tr.tolist()}  test = {te.tolist()}")

# Output:
# split 1: train = [0, 1, 2, 3]  test = [4, 5]
# split 2: train = [0, 1, 2, 3, 4, 5]  test = [6, 7]
# split 3: train = [0, 1, 2, 3, 4, 5, 6, 7]  test = [8, 9]
# split 4: train = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]  test = [10, 11]

The training sets grow and always contain everything that precedes the test set: in split 3 the model is trained on rows 0 to 7 and evaluated on 8 and 9. As the documentation says, unlike the other schemes the successive training sets are supersets of the previous ones. The price is that the first splits train on very little data, and that you cannot shuffle: if the data have seasonality or a change of regime, the estimate is still that of the past.

Choosing the hyperparameters

A hyperparameter is a choice that is not learned from the data but fixed beforehand: for logistic regression, the strength of the penalty, C (smaller means more penalty and a simpler model). You try a grid of values, compute for each the cross-validation score and take the best: this is what GridSearchCV does, which at the end refits (refit) the best model on the whole training set. Before the search you set aside a test set, which is used only once.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV, StratifiedKFold, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=0, stratify=y)
pipe = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000))
Cs = np.logspace(-3, 2, 11)
gs = GridSearchCV(pipe, {"logisticregression__C": Cs}, cv=StratifiedKFold(5, shuffle=True, random_state=0)).fit(Xtr, ytr)

means, stds = gs.cv_results_["mean_test_score"], gs.cv_results_["std_test_score"]
i = int(np.argmax(means))
se = stds[i] / np.sqrt(5)                              # standard error of the best estimate
# one-standard-error rule: the simplest model (smallest C) within 1 SE of the best
j = min(k for k in range(len(Cs)) if means[k] >= means[i] - se)
print(f"best C: {Cs[i]:.3g}  CV accuracy = {means[i]:.4f} (SE {se:.4f})")
print(f"C with the 1-SE rule: {Cs[j]:.3g}  CV accuracy = {means[j]:.4f}")
print(f"test, once only, with the best C: {gs.score(Xte, yte):.4f}")

# Output:
# best C: 0.316  CV accuracy = 0.9812 (SE 0.0054)
# C with the 1-SE rule: 0.1  CV accuracy = 0.9765
# test, once only, with the best C: 0.9650

The best value is C = 0.316 with a validation accuracy of 0.9812, and its standard error is 0.0054 (standard deviation of the five scores divided by the square root of 5). The one-standard-error rule instead chooses the simplest model within one standard error of the best: here C = 0.1 with 0.9765, almost indistinguishable from the maximum but with more penalty. The cited textbook describes it in section 6.1.3, with the motivation that among more or less equivalent models the simplest is preferable, and it is preferable when the difference between two models is below the noise threshold of the estimate. The accuracy on the test set, computed only once, is 0.9650: lower than the validation one, as expected from a value that was chosen precisely because it was the maximum on validation, and with 143 examples it would be a mistake to consider it significantly different.

Nested validation

The best score of a grid search is the maximum of many noisy estimates, so it is systematically optimistic: the winning configuration got part of its luck from the folds on which it was evaluated. To measure the optimism we build a case in which the honest answer is known: 60 examples and 40 features of pure noise, labels with no link to the measurements. We repeat on 10 different datasets, comparing the best score of the search with nested validation: an outer loop splits the data, and inside each outer training set an independent grid search chooses the hyperparameters; the outer fold, never seen by the search, gives the honest estimate.

import numpy as np
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV, StratifiedKFold, cross_val_score
from sklearn.pipeline import Pipeline

grid = {"sel__k": [1, 2, 3, 5, 8, 12, 20, 40], "clf__C": [0.01, 0.1, 1, 10, 100]}   # 40 configurations
inner = StratifiedKFold(5, shuffle=True, random_state=0)
outer = StratifiedKFold(5, shuffle=True, random_state=1)
best, nested = [], []
for seed in range(10):
    rng = np.random.default_rng(seed)
    X = rng.standard_normal((60, 40))                  # 60 examples, 40 features of pure noise
    y = np.repeat([0, 1], 30)                          # labels with no link to X: honest limit = 0.5
    pipe = Pipeline([("sel", SelectKBest(f_classif)), ("clf", LogisticRegression(max_iter=1000))])
    gs = GridSearchCV(pipe, grid, cv=inner).fit(X, y)
    best.append(gs.best_score_)
    nested.append(cross_val_score(GridSearchCV(pipe, grid, cv=inner), X, y, cv=outer).mean())
print("mean over 10 noise datasets (true value: 0.5)")
print("best score of the search (best_score_):", round(np.mean(best), 3))
print("nested validation:                     ", round(np.mean(nested), 3))

# Output:
# mean over 10 noise datasets (true value: 0.5)
# best score of the search (best_score_): 0.588
# nested validation:                      0.49

On the grid of 40 configurations (8 values of k for SelectKBest and 5 values of C) the best mean score is 0.588, but the true value is 0.5 because the data are noise: almost nine percentage points above the true value, which look like a signal. Nested validation gives 0.490, that is, chance. The computation is costly: with 5 outer folds and a grid of 11 values evaluated on 5 inner folds it is 280 fits, against 56 for the simple search (counts in the exercise below). The practical conclusion is to use simple validation to choose the hyperparameters, but not to report its score as an estimate of performance: either you use the nested one, or you use a test set that has never been touched.

Worked exercise: how many fits, and on how many examples

With the 569 examples of the dataset, a 5-fold and a grid of 11 values of C, compute the sizes of the folds and the number of fits of the simple search and of the nested one.

import numpy as np
from sklearn.model_selection import KFold

n, K = 569, 5
sizes = [(len(tr), len(te)) for tr, te in KFold(K).split(np.zeros(n))]
print("(train, test) for each fold:", sizes)
configs = 11                                          # values of C in the grid
print("fits of a grid search with 5-fold CV:", configs * K, "+ 1 refit =", configs * K + 1)
print("fits of a nested validation, 5 outer x search:", K * (configs * K + 1))

# Output:
# (train, test) for each fold: [(455, 114), (455, 114), (455, 114), (455, 114), (456, 113)]
# fits of a grid search with 5-fold CV: 55 + 1 refit = 56
# fits of a nested validation, 5 outer x search: 280

569 is not divisible by 5: four folds have 114 test examples and one has 113, and as a consequence the training sets have 455 and 456 examples. The simple search costs 11 × 5 = 55 fits plus one refit; the nested one repeats the whole search in each of the 5 outer folds, hence 5 × 56 = 280. When the model is slow, the factor five of the nested search decides whether it is feasible, and the realistic answer is often a single test set of adequate size.

When cross-validation is not enough

• Future data will have a different distribution from the collected data: the estimate concerns the past, and monitoring is needed (chapter 12).
• The examples are not independent in ways you do not know: you need to know the structure (groups, time) to choose the scheme.
• The metric does not reflect the cost of errors: a precise estimate of the wrong number is still wrong (chapter 11).

What the three items have in common is that cross-validation only guarantees an estimate of the error on data similar to those already seen; the similarity is an assumption of the designer, not a property of the method.

What you can do now

You can implement a K-fold and explain why a single partition gives estimates that jump by a few points, choose between stratified, grouped and time-based schemes by looking at the structure of the data, choose a hyperparameter with grid search and the one-standard-error rule, and recognise when the score of a search is inflated and nested validation is needed. The next step is to make use of trees, which have many hyperparameters to validate: 06 - Decision Trees and Random Forest.

The same series also includes: 02 - Linear Regression from Scratch in Python and 07 - Support Vector Machines: Kernel Trick and Non-Linear SVM.