6.3 Resampling and cross-validation

A model must be judged on data it has not seen. Cross-validation recycles the data we have: every part takes a turn as the unseen data.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from matplotlib.axes import Axes
from sklearn import datasets
from sklearn.base import BaseEstimator
from sklearn.dummy import DummyRegressor
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import KFold, TimeSeriesSplit, cross_val_score, train_test_split
from sklearn.neighbors import KNeighborsRegressor
from sklearn.pipeline import Pipeline, make_pipeline
from sklearn.preprocessing import StandardScaler

from statly import style

style.use()
rng = np.random.default_rng(1)

diabetes = datasets.load_diabetes(as_frame=True, scaled=False)
X, y = diabetes.data, diabetes.target.to_numpy()  # 442 patients: 10 measurements, progression
NEIGHBOURS = np.arange(1, 61)


def knn(k: int) -> Pipeline:
    """K nearest neighbours on standardised measurements."""
    return make_pipeline(StandardScaler(), KNeighborsRegressor(k))


def rmse(
    model: BaseEstimator,
    X_fit: pd.DataFrame,
    y_fit: np.ndarray,
    X_eval: pd.DataFrame,
    y_eval: np.ndarray,
) -> float:
    """Fit on one set, return the root mean squared error on another."""
    fitted = model.fit(X_fit, y_fit)
    return float(np.sqrt(mean_squared_error(y_eval, fitted.predict(X_eval))))

6.3.1 Training error lies

The 442 diabetes patients from Chapter 4.2 have ten measurements each, and a number for how much their disease progressed within a year. We predict that number with k nearest neighbours, which averages the k most similar patients. Small k is flexible; large k is rigid. Which k is best?

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=3)
train_error = [rmse(knn(k), X_train, y_train, X_train, y_train) for k in NEIGHBOURS]
test_error = [rmse(knn(k), X_train, y_train, X_test, y_test) for k in NEIGHBOURS]


def plot_errors(ax: Axes) -> None:
    """Training and test error against flexibility (1/k)."""
    ax.plot(1 / NEIGHBOURS, train_error, color=style.GREY, label="training")
    ax.plot(1 / NEIGHBOURS, test_error, color=style.CORAL, label="test")
    best = NEIGHBOURS[int(np.argmin(test_error))]
    ax.plot(1 / best, min(test_error), "o", ms=7, color=style.CORAL)
    ax.annotate(
        f"best: k = {best}",
        (1 / best, min(test_error)),
        xytext=(0, -18),
        textcoords="offset points",
        ha="center",
        fontsize=9,
        color=style.CORAL,
    )
    ax.set_xscale("log")
    ax.set_xticks(1 / np.array([50, 20, 10, 5, 2, 1]), ["50", "20", "10", "5", "2", "1"])
    ax.minorticks_off()
    ax.legend(fontsize=8)
    ax.set(xlabel="neighbours k (flexible to the right)", ylabel="prediction error (RMSE)")


fig, ax = plt.subplots(figsize=(8, 3.3))
plot_errors(ax)
plt.show()

Figure 6.3.1: Prediction error of k nearest neighbours, on the patients used for fitting (training) and on patients held back (test), as k goes from 60 (rigid) to 1 (flexible). Training error falls all the way to zero; test error has a minimum in between.

With k = 1, each training patient is its own nearest neighbour, so the training error is zero. On new patients it is the worst of all. This is the bias and variance trade-off of Chapter 1.3, in real data: the error we care about is the test error, and the training error cannot see it.

6.3.2 One split is not enough

The obvious fix is to hold back a validation set, as above. But which patients end up in it is luck, and with 442 patients that luck matters:

def plot_splits(ax: Axes, splits: int = 10) -> None:
    """One validation curve per random split, with its minimum marked."""
    for seed in range(splits):
        parts = train_test_split(X, y, test_size=0.3, random_state=100 + seed)
        errors = [rmse(knn(k), parts[0], parts[2], parts[1], parts[3]) for k in NEIGHBOURS]
        ax.plot(NEIGHBOURS, errors, color=style.BLUE, lw=1, alpha=0.6)
        best = int(np.argmin(errors))
        ax.plot(NEIGHBOURS[best], errors[best], "o", ms=5, color=style.CORAL)
    ax.set(xlabel="neighbours k", ylabel="test error (RMSE)", ylim=(48, 80))


fig, ax = plt.subplots(figsize=(8, 3.3))
plot_splits(ax)
plt.show()

Figure 6.3.2: Test error against k for ten different random splits into 70% for fitting and 30% for testing. The curves have the same shape, but their level and their best k differ from split to split.

A single split wastes data, since 30% of the patients never help fit the model, and its answer depends on the split.

6.3.3 Cross-validation

K-fold cross-validation fixes both problems. Cut the data into k parts, often 5 or 10. Fit on all parts but one, measure the error on the one left out, and repeat until every part has been left out once. The average of the k errors is the estimate.

FOLDS = KFold(10, shuffle=True, random_state=1)
cv = np.array(
    [
        -cross_val_score(knn(k), X, y, cv=FOLDS, scoring="neg_root_mean_squared_error")
        for k in NEIGHBOURS
    ]
)
cv_mean, cv_se = cv.mean(axis=1), cv.std(axis=1, ddof=1) / np.sqrt(cv.shape[1])


def plot_folds(ax: Axes, k: int = 5) -> None:
    """Rows of blocks: which fifth of the data is held out in each round."""
    for round_ in range(k):
        for part in range(k):
            color = style.CORAL if part == round_ else style.BLUE
            ax.add_patch(
                plt.Rectangle((part, k - 1 - round_), 0.92, 0.8, color=color, alpha=0.6)
            )
    ax.set(xlim=(0, k), ylim=(-0.2, k), title="5 folds", xlabel="the data, in five parts")
    ax.set_yticks(np.arange(k) + 0.4, [f"round {k - i}" for i in range(k)])
    ax.set_xticks([])
    ax.grid(False)
    ax.tick_params(length=0)


def plot_cv(ax: Axes) -> None:
    """The cross-validation curve with its standard errors and the one-standard-error line."""
    ax.fill_between(
        NEIGHBOURS, cv_mean - cv_se, cv_mean + cv_se, color=style.BLUE, alpha=0.15, lw=0
    )
    ax.plot(NEIGHBOURS, cv_mean, color=style.BLUE)
    best = int(np.argmin(cv_mean))
    ax.axhline(cv_mean[best] + cv_se[best], color=style.INK, ls="--", lw=1)
    ax.plot(NEIGHBOURS[best], cv_mean[best], "o", ms=7, color=style.CORAL)
    simplest = int(np.max(np.where(cv_mean <= cv_mean[best] + cv_se[best])))
    ax.plot(NEIGHBOURS[simplest], cv_mean[simplest], "o", ms=7, color=style.GOLD)
    style.note(ax, f"lowest: k = {NEIGHBOURS[best]}", x=0.45, y=0.95, color=style.CORAL)
    style.note(
        ax,
        f"simplest within one\nstandard error: k = {NEIGHBOURS[simplest]}",
        x=0.45,
        y=0.8,
        color=style.GOLD,
    )
    ax.set(title="10-fold cross-validation", xlabel="neighbours k", ylabel="error (RMSE)")


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_folds(axes[0])
plot_cv(axes[1])
fig.tight_layout()
plt.show()

Figure 6.3.3: Left: 5-fold cross-validation. Each row is one round; the red part is held out and the blue parts are used for fitting. Right: the 10-fold cross-validation error for every k, with ± one standard error. The dashed line is one standard error above the minimum.

Every patient is used for fitting k − 1 times and for testing once, and the estimate no longer hangs on one lucky split. The curve is flat near its minimum, and the standard errors show that many k are about equally good. A common rule picks the simplest model whose error is within one standard error of the best: a little less accuracy on paper, much less risk of fitting noise.

With k equal to the number of data points, each round leaves out one case: leave-one-out cross-validation. It uses the most data for fitting but needs as many fits as there are cases. For least squares there is a shortcut that needs only one fit, using the leverage $h_i$ of Chapter 4.3:

$$\text{CV}_{\text{LOO}} = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{y_i - \hat y_i}{1 - h_i} \right)^2$$

each residual inflated by how much its own point pulled the fit towards itself. Five or ten folds are the usual compromise: leave-one-out has little bias but high variance, because its n training sets are nearly identical, while 5 or 10 folds have a little more bias (each fit uses less data) and less variance.

6.3.4 Nested cross-validation

Chapter 6.2 kept three parts apart: training to fit, validation to choose, test to report. With little data, a separate test set wastes cases and measures little. Nested cross-validation does both jobs with folds. An outer loop of folds does the reporting; inside each outer training part, an inner cross-validation does the choosing. The outer test fold never takes part in any choice, so the average over outer folds is an honest estimate of the whole procedure, choosing included.

outer = KFold(5, shuffle=True, random_state=0)
inner_choice, outer_error = [], []
for train_rows, test_rows in outer.split(X):
    inner_scores = [
        -cross_val_score(
            knn(k),
            X.iloc[train_rows],
            y[train_rows],
            cv=4,
            scoring="neg_root_mean_squared_error",
        ).mean()
        for k in NEIGHBOURS
    ]
    best = int(NEIGHBOURS[np.argmin(inner_scores)])
    inner_choice.append(best)
    outer_error.append(
        rmse(knn(best), X.iloc[train_rows], y[train_rows], X.iloc[test_rows], y[test_rows])
    )


def plot_nested(ax: Axes) -> None:
    """Outer folds as rows of five cells; under the first, its inner folds as smaller rows."""
    y_pos = 0.0
    for r in range(5):
        for c in range(5):
            color = style.CORAL if c == r else style.BLUE
            ax.add_patch(plt.Rectangle((c, y_pos), 0.95, 0.8, color=color, alpha=0.55))
        ax.text(-0.15, y_pos + 0.4, f"outer {r + 1}", ha="right", va="center", fontsize=8.5)
        y_pos -= 1
        if r == 0:
            for i in range(4):
                for c in range(4):
                    color = style.GOLD if c == i else style.BLUE
                    ax.add_patch(
                        plt.Rectangle((1 + c, y_pos + 0.2), 0.95, 0.5, color=color, alpha=0.4)
                    )
                ax.text(
                    0.85, y_pos + 0.45, f"inner {i + 1}", ha="right", va="center", fontsize=7.5
                )
                y_pos -= 0.65
            y_pos -= 0.3
    ax.set(xlim=(-1.3, 5.1), ylim=(y_pos + 0.6, 1))
    ax.set_axis_off()


fig = plt.figure(figsize=(10.5, 4))
grid = fig.add_gridspec(1, 2, width_ratios=[1.3, 1])
plot_nested(fig.add_subplot(grid[0]))
ax = fig.add_subplot(grid[1])
ax.bar(range(1, 6), outer_error, color=style.CORAL, alpha=0.7)
for i, (err, k) in enumerate(zip(outer_error, inner_choice), start=1):
    ax.text(i, err + 0.8, f"k = {k}", ha="center", fontsize=8.5)
ax.axhline(np.mean(outer_error), color=style.INK, ls="--", lw=1)
style.note(ax, f"average {np.mean(outer_error):.1f}", x=0.03, y=0.97, color=style.INK)
ax.set(xlabel="outer fold", ylabel="test error", ylim=(0, max(outer_error) * 1.2))
fig.tight_layout()
plt.show()

Figure 6.3.4: Nested cross-validation with 5 outer and 4 inner folds. Each outer round (one row) holds out one fold (red) for testing. Its remaining data are cross-validated again (the small rows below it) to choose k, and the chosen model is refitted on the whole outer training part and tested once on the red fold. Right: the 5 outer test errors for k nearest neighbours on the diabetes patients, with the k chosen in each round; their average is the honest estimate.

6.3.5 Folds that respect the structure

Random folds assume the cases are independent. Two common structures break that, and then random folds leak information from the held-out cases into the training data.

• Groups. Several rows can come from the same patient, household or firm (Chapter 4.7). If one patient’s rows land in both the training and the held-out part, the model is tested on someone it has partly seen. Grouped cross-validation keeps all rows of a group in the same fold. Stratified cross-validation, a different idea, keeps the share of each class the same in every fold, which matters for rare classes.
• Time. A model for the future must be tested on the future. Shuffling time points lets the model peek at days after the ones it predicts, and neighbouring days are nearly identical. A rolling origin (forward-chaining) scheme trains on everything up to a point and tests on what comes right after, then moves the point forward.

sun = sm.datasets.sunspots.load_pandas().data.SUNACTIVITY.to_numpy()
LAGS = 3
lagged = np.column_stack([sun[i : len(sun) - LAGS + i] for i in range(LAGS)])
target = sun[LAGS:]


def plot_schemes(ax: Axes, rounds: int = 5, length: int = 30) -> None:
    """Training and test blocks for three cross-validation schemes."""
    y = 0
    schemes = [
        ("shuffled k-fold", None),
        ("expanding window", "expanding"),
        ("fixed window", "fixed"),
    ]
    for name, kind in schemes:
        for r in range(rounds):
            if kind is None:
                test = set(
                    np.random.default_rng(r).choice(length, length // rounds, replace=False)
                )
                train = set(range(length)) - test
            else:
                end = 10 + r * 4
                train = set(range(end - 10 if kind == "fixed" else 0, end))
                test = set(range(end, end + 4))
            for t in range(length):
                color = style.CORAL if t in test else style.BLUE if t in train else "none"
                ax.add_patch(
                    plt.Rectangle(
                        (t, -y), 0.9, 0.8, color=color, alpha=0.6 if color != "none" else 0
                    )
                )
            y += 1
        ax.text(-1, -y + rounds / 2 + 0.4, name, ha="right", va="center", fontsize=9)
        y += 1
    ax.set(xlim=(-14, length), ylim=(-y + 1, 1), xticks=[], yticks=[], xlabel="time →")
    ax.grid(False)
    for side in ("left", "bottom"):
        ax.spines[side].set_visible(False)


def plot_honesty(ax: Axes) -> None:
    """Error estimates from shuffled and time-respecting folds."""
    model = make_pipeline(StandardScaler(), KNeighborsRegressor(5))
    shuffled = -cross_val_score(
        model,
        lagged,
        target,
        cv=KFold(5, shuffle=True, random_state=0),
        scoring="neg_root_mean_squared_error",
    ).mean()
    rolling = -cross_val_score(
        model, lagged, target, cv=TimeSeriesSplit(5), scoring="neg_root_mean_squared_error"
    ).mean()
    names, values = ["shuffled k-fold", "rolling origin"], [shuffled, rolling]
    ax.barh(names, values, color=[style.GREY, style.BLUE], alpha=0.7)
    for i, v in enumerate(values):
        ax.text(v + 0.5, i, f"{v:.1f}", va="center", fontsize=9)
    ax.invert_yaxis()
    ax.set(
        title="Estimated prediction error",
        xlabel="error (RMSE, sunspots)",
        xlim=(0, max(values) * 1.2),
    )


fig = plt.figure(figsize=(10.5, 3.6))
grid = fig.add_gridspec(1, 2, width_ratios=[1.3, 1])
plot_schemes(fig.add_subplot(grid[0]))
plot_honesty(fig.add_subplot(grid[1]))
fig.tight_layout()
plt.show()

Figure 6.3.5: Left: three ways to cut a series into training (blue) and test (red) parts, one row per round: shuffled k-fold, which mixes past and future; rolling origin with an expanding window; and rolling origin with a fixed window. Right: predicting next year’s sunspots from the previous three years with 5 nearest neighbours. Shuffled folds promise a smaller error (about 15) than the model achieves on years it has really not seen (about 19); the rolling origin gives the honest answer.

Chapter 9.2 uses the rolling origin to evaluate forecasts.

6.3.6 Comparing models

Cross-validation compares any models on equal terms, including very different ones:

MODELS = [
    ("the average", DummyRegressor()),
    (
        f"neighbours, k = {NEIGHBOURS[int(np.argmin(cv_mean))]}",
        knn(int(NEIGHBOURS[int(np.argmin(cv_mean))])),
    ),
    ("linear regression", LinearRegression()),
]
fold_errors = [
    -cross_val_score(model, X, y, cv=FOLDS, scoring="neg_root_mean_squared_error")
    for _, model in MODELS
]


def plot_folds_by_model(ax: Axes) -> None:
    """Each fold's error per model, with lines joining the same fold."""
    for fold in range(10):
        ax.plot(
            range(len(MODELS)),
            [errors[fold] for errors in fold_errors],
            color=style.GREY,
            lw=0.7,
            alpha=0.6,
        )
    for i, errors in enumerate(fold_errors):
        ax.plot(np.full(10, i), errors, "o", ms=5, color=style.BLUE, alpha=0.7)
        ax.plot([i - 0.2, i + 0.2], [errors.mean()] * 2, color=style.CORAL, lw=2.5)
    ax.set_xticks(range(len(MODELS)), [name for name, _ in MODELS])
    ax.set(ylabel="error in the fold (RMSE)", xlim=(-0.5, len(MODELS) - 0.5))


fig, ax = plt.subplots(figsize=(8, 3.3))
plot_folds_by_model(ax)
plt.show()

Figure 6.3.6: Error in each of the ten folds for three models of disease progression: the overall average (no model), the best k nearest neighbours, and a plain linear regression. Linear regression is best in most folds.

The lines join the same fold across models: some folds are harder for every model, so compare models fold by fold. Here the simple line beats the flexible neighbours, because with 442 patients and ten measurements there is too little data for flexibility to pay off.

6.3.7 The classic mistake

Cross-validation only works if everything learned from the data happens inside the folds. A common error is to pick the best predictors first, on all the data, and cross-validate afterwards. The held-out cases then helped choose the model, and the estimate is too good.

CASES, NOISE_FEATURES, KEEP = 50, 5_000, 20


def leaked_and_honest(seed: int) -> tuple[float, float]:
    """Cross-validated accuracy with selection outside the folds, and inside them."""
    noise_rng = np.random.default_rng(seed)
    features = noise_rng.normal(size=(CASES, NOISE_FEATURES))
    outcome = noise_rng.integers(0, 2, CASES)
    folds = KFold(5, shuffle=True, random_state=seed)
    chosen = SelectKBest(f_classif, k=KEEP).fit(features, outcome).transform(features)
    leaked = cross_val_score(LogisticRegression(), chosen, outcome, cv=folds).mean()
    pipeline = make_pipeline(SelectKBest(f_classif, k=KEEP), LogisticRegression())
    honest = cross_val_score(pipeline, features, outcome, cv=folds).mean()
    return leaked, honest


results = np.array([leaked_and_honest(seed) for seed in range(30)])


def plot_leakage(ax: Axes) -> None:
    """Accuracy from 30 noise data sets, with leaked and honest selection."""
    for i, (name, color) in enumerate(
        (("selected on all data", style.CORAL), ("selected inside folds", style.BLUE))
    ):
        ax.scatter(
            results[:, i],
            np.full(len(results), i) + rng.uniform(-0.12, 0.12, len(results)),
            s=18,
            color=color,
            alpha=0.7,
        )
        ax.plot([results[:, i].mean()] * 2, [i - 0.3, i + 0.3], color=color, lw=2.5)
    ax.axvline(0.5, color=style.INK, ls="--", lw=1)
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.set_yticks([0, 1], ["selected on all data", "selected inside folds"])
    ax.set(xlabel="cross-validated accuracy", ylim=(-0.6, 1.6), xlim=(0.2, 1.02))


fig, ax = plt.subplots(figsize=(8, 2.6))
plot_leakage(ax)
plt.show()

Figure 6.3.7: Pure noise: 50 cases, 5,000 random measurements, and a random yes/no outcome; no classifier can beat 50% here. Choosing the 20 measurements most related to the outcome on all the data, then cross-validating, reports far better accuracy (red). Choosing them inside each fold gives the honest answer (blue).

The leaked version reports about 95% accuracy on data where 50% is the truth. Selecting predictors, scaling, filling in missing values, tuning: put it all inside the cross-validation loop, for example as one pipeline that is fitted fresh in every fold.

6.3.8 In one picture

• Training error always flatters a model; judge it on unseen data.
• A single validation split wastes data and depends on luck.
• K-fold cross-validation lets every case be unseen once. Use 5 or 10 folds.
• Among models with similar error, prefer the simplest (the one-standard-error rule).
• Choosing and reporting need separate data (Chapter 6.2); with little data, use nested cross-validation: inner folds choose, outer folds report.
• Keep groups together in folds; for time, train on the past and test on the future (rolling origin).
• Everything the model learns, including choosing predictors, must happen inside the folds.

Next: models with many predictors, and how to keep them from overfitting.