cross_validate()

K-fold cross-validation for a GAM specification.

Usage

cross_validate(
    formula,
    data,
    *,
    family=None,
    n_folds=10,
    method="GCV",
    metric="deviance",
    select=False,
    seed=None,
)

Estimates the out-of-sample predictive performance of a GAM by repeatedly refitting it on n_folds - 1 folds of the data and scoring the fit on the remaining held-out fold, then aggregating the resulting losses across folds. Because each fold is scored on data that was not used to fit that particular model, cross_validate() gives a more honest estimate of generalization error than simply scoring a single fit on the data it was trained on, which is optimistic (the fit has already adapted to that data’s noise). Use it to compare candidate formulas, families, fitting methods, or basis choices on equal footing, or as a sanity check that a model selected by GCV/REML/ML also performs well out of sample.

Two loss metrics are available via metric:

"deviance" (default): for each fold, the family’s deviance() between the held-out responses and the predictions from a model fit on the training folds, divided by the number of test observations in that fold. This matches the deviance-based loss the family itself uses during fitting, so it is comparable across method choices for the same family.
"mse": the mean squared error, mean((y_test - pred) ** 2), on the response scale. This is family-agnostic and directly interpretable in the response’s original units, but does not account for family-specific variance structure the way deviance does.

Folds are constructed by randomly permuting the row indices with rng.permutation(n) (where rng is seeded from seed) and then assigning fold id i * n_folds // n to the row that ends up in position i of the permutation. This produces a non-stratified partition into n_folds contiguous-in-permutation-order, roughly (but not exactly, when n is not a multiple of n_folds) equal-size groups; no attempt is made to balance the distribution of the response or any covariate across folds.

Parameters

formula: str: GAM formula string, e.g. "y ~ s(x1) + s(x2) + x3". The response named on the left-hand side is looked up in data to build the fold assignment and to compute the loss; the right-hand side is passed unchanged to ~whittaker.gam.GAM for every fold.
data: dict[str, numpy.ndarray] or InputData: Column-oriented data as {name: 1-D array} (or any InputData-compatible object, such as a pandas.DataFrame or polars.DataFrame). Must contain every column referenced by formula, all of equal length.
family: Family = None: Response distribution family, e.g. Gaussian(), Binomial(), Poisson(), Gamma(), or Tweedie(). Used both to fit each fold’s ~whittaker.gam.GAM and, when metric="deviance", to compute each fold’s loss via family.deviance(). Defaults to Gaussian().
n_folds: int = 10: Number of folds to split the data into. Must be at least 2 and, for every fold to receive at least one test observation, should not exceed the number of rows in data. Defaults to 10. See the Notes section below for guidance on choosing this value.
method: str = "GCV": Smoothing-parameter selection method passed through to GAM.fit() for every fold. One of "GCV" (default), "REML", or "ML"; see GAM.fit() for what each criterion optimizes.
metric: str = "deviance": Loss metric to compute on each held-out fold: "deviance" (default) or "mse". See the discussion above for exactly how each is computed.
select: bool = False: Whether to add shrinkage penalties for automatic smooth-term selection, forwarded to GAM.fit(select=...) for every fold. Defaults to False.
seed: int = None: Seed for the numpy.random.default_rng() random number generator used to build the fold assignment. Pass a fixed integer to make the fold split (and hence the resulting CVResult) reproducible across calls; None (the default) uses a fresh, non-reproducible seed.

Returns

CVResult: Cross-validation result holding the mean out-of-sample loss (cv_score), the per-fold losses (cv_scores), their standard error (cv_se), and the number of folds used (n_folds).

Notes

The number of folds controls a bias-variance tradeoff in the CV estimate itself. With a small n_folds (e.g. 3-5), each training fold omits a large fraction of the data, so the fitted model is somewhat different from (typically smoother/less flexible than) a model fit on the full dataset; the resulting cv_score tends to be pessimistically biased, but because there are only a few, relatively large folds, cv_scores tends to have lower variance across repeated runs. With a large n_folds (up to the leave-one-out limit, n_folds = n), each training fold is nearly the full dataset, so bias shrinks toward the true generalization error of the full-data fit — but the individual test folds are tiny (a single point at n_folds = n), so cv_scores becomes noisier (higher variance), and fitting cost grows linearly with n_folds since a full GAM.fit() is performed once per fold. In practice, n_folds = 5 or n_folds = 10 are common compromises between these effects. Leave-one-out cross-validation is rarely used directly for GAMs because of its cost; method="GCV" in GAM.fit() already computes an efficient analytical approximation to the leave-one-out error from a single fit, without refitting the model n times.

Examples

import numpy as np
import whittaker as wt

rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 1, 200))
y = np.sin(2 * np.pi * x) + rng.normal(scale=0.2, size=200)

result = wt.cross_validate("y ~ s(x)", {"x": x, "y": y}, n_folds=5, seed=0)
print(result.cv_score, result.cv_se)

0.044761796080802364 0.0034832382870951125

# Compare two candidate formulas on the same folds using the seed.
linear_result = wt.cross_validate("y ~ x", {"x": x, "y": y}, n_folds=5, seed=0)
smooth_result = wt.cross_validate("y ~ s(x)", {"x": x, "y": y}, n_folds=5, seed=0)
linear_result.cv_score, smooth_result.cv_score

(0.2517163130471863, 0.044761796080802364)