# Cross-validation

Smoothing parameter selection (GCV, REML) optimises within a single fit, but cross-validation answers a different question: **how well does the model generalise to unseen data?** Use it to compare formulas, families, or basis configurations on the same dataset.

Whittaker provides [cross_validate()](../reference/cross_validate.md#whittaker.cross_validate), which performs K-fold cross-validation on a GAM specification and returns a [CVResult](../reference/CVResult.md#whittaker.CVResult) with the mean score, per-fold scores, and a standard error.


# Basic usage


``` python
import numpy as np
import whittaker as wk

# Generate data with a known signal
rng = np.random.default_rng(23)
n = 400
x = np.linspace(0, 2 * np.pi, n)
y = np.sin(x) + 0.5 * np.sin(3 * x) + rng.normal(0, 0.3, n)

data = {"x": x, "y": y}

# 10-fold cross-validation (the default)
result = wk.cross_validate("y ~ s(x)", data, method="REML")

print(f"CV score (deviance): {result.cv_score:.4f}")
print(f"SE of CV score:      {result.cv_se:.4f}")
print(f"Number of folds:     {result.n_folds}")
```


    CV score (deviance): 0.0984
    SE of CV score:      0.0060
    Number of folds:     10


The `cv_score` is the mean out-of-sample loss across folds. Lower is better. The `cv_se` gives the standard error of this mean, which is useful for comparing models: two models whose scores differ by less than one standard error are essentially equivalent.


# Comparing models

Cross-validation is most useful for comparing competing specifications. For example, how many basis functions does the smooth need?


``` python
# Compare different basis dimensions
results = {}
for k in [5, 10, 15, 20, 30]:
    cv = wk.cross_validate(f"y ~ s(x, k={k})", data, method="REML", seed=23)
    results[k] = cv

for k, cv in results.items():
    print(f"k={k:2d}: CV = {cv.cv_score:.4f} (SE = {cv.cv_se:.4f})")
```


    k= 5: CV = 0.2044 (SE = 0.0113)
    k=10: CV = 0.0989 (SE = 0.0067)
    k=15: CV = 0.1001 (SE = 0.0066)
    k=20: CV = 0.1004 (SE = 0.0066)
    k=30: CV = 0.1006 (SE = 0.0067)


``` python
import altair as alt

# Plot CV scores with error bars
plot_data = [
    {"k": k, "cv_score": float(cv.cv_score),
     "lower": float(cv.cv_score - cv.cv_se),
     "upper": float(cv.cv_score + cv.cv_se)}
    for k, cv in results.items()
]

points = alt.Chart({"values": plot_data}).mark_point(size=60, color="steelblue").encode(
    x=alt.X("k:Q", title="Basis dimension k", scale=alt.Scale(domain=[3, 32])),
    y=alt.Y("cv_score:Q", title="CV score (deviance)"),
)

errorbars = alt.Chart({"values": plot_data}).mark_rule(color="steelblue").encode(
    x="k:Q", y="lower:Q", y2="upper:Q",
)

(points + errorbars).properties(
    width="container", height=280,
    title="Cross-validation score by basis dimension"
)
```


<style>
  #altair-viz-5ddcdce33b93468cb6e8c5d884e33ac6.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-5ddcdce33b93468cb6e8c5d884e33ac6.vega-embed details,
  #altair-viz-5ddcdce33b93468cb6e8c5d884e33ac6.vega-embed details summary {
    position: relative;
  }
</style>


The CV score drops as `k` increases from 5 to around 15, then levels off. Since REML penalises away unnecessary complexity, larger `k` values do not overfit, but the CV score confirms that 15 basis functions are sufficient for this signal.


# Comparing families

Cross-validation also helps choose between response distributions. Here we simulate count data and compare Poisson and Gaussian fits:


``` python
# Count data
rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2, n)
mu = np.exp(0.5 + 0.8 * np.sin(2 * np.pi * x))
y_counts = rng.poisson(mu).astype(float)
count_data = {"x": x, "y": y_counts}

# Compare families
cv_gauss = wk.cross_validate("y ~ s(x)", count_data, method="REML", seed=23)
cv_pois = wk.cross_validate("y ~ s(x)", count_data, family=wk.Poisson(), method="REML", seed=23)

print(f"Gaussian CV: {cv_gauss.cv_score:.4f} (SE = {cv_gauss.cv_se:.4f})")
print(f"Poisson  CV: {cv_pois.cv_score:.4f} (SE = {cv_pois.cv_se:.4f})")
```


    Gaussian CV: 1.7966 (SE = 0.1420)
    Poisson  CV: 1.0957 (SE = 0.0939)


> **Note: Comparing scores across families**
>
> The `"deviance"` metric uses each family's own deviance, so scores from different families are not directly comparable on the same scale. Switch to `metric="mse"` for an apples-to-apples comparison on the response scale:
>
> <div id="33ffbee7" class="cell" data-execution_count="6">
>
> ``` python
> cv_gauss_mse = wk.cross_validate(
>     "y ~ s(x)", count_data, method="REML", metric="mse", seed=23
> )
> cv_pois_mse = wk.cross_validate(
>     "y ~ s(x)", count_data, family=wk.Poisson(), method="REML", metric="mse", seed=23
> )
>
> print(f"Gaussian MSE: {cv_gauss_mse.cv_score:.4f}")
> print(f"Poisson  MSE: {cv_pois_mse.cv_score:.4f}")
> ```
>
> <div class="cell-output cell-output-stdout">
>
>     Gaussian MSE: 1.7966
>     Poisson  MSE: 1.7976
>
> <div id="per-fold-scores" class="section level1">
>
> # Per-fold scores
>
> The `cv_scores` array gives the loss for each fold, which is useful for checking whether a single fold is driving the overall score:
>
> <div id="8e9d6db7" class="cell" data-execution_count="7">
>
> ``` python
> result = wk.cross_validate("y ~ s(x)", data, method="REML", seed=23)
>
> print(f"Per-fold scores: {result.cv_scores.round(4)}")
> print(f"Mean: {result.cv_scores.mean():.4f}")
> print(f"Std:  {result.cv_scores.std():.4f}")
> ```
>
> <div class="cell-output cell-output-stdout">
>
>     Per-fold scores: [0.0936 0.0863 0.0791 0.1154 0.1017 0.0837 0.0924 0.1437 0.1171 0.0762]
>     Mean: 0.0989
>     Std:  0.0200
>
> </div>
>
> </div>
>
> <div id="443d3894" class="cell" data-execution_count="8">
>
> ``` python
> # Plot per-fold scores
> fold_data = [
>     {"fold": i + 1, "score": float(s)}
>     for i, s in enumerate(result.cv_scores)
> ]
>
> bars = alt.Chart({"values": fold_data}).mark_bar(color="steelblue", opacity=0.7).encode(
>     x=alt.X("fold:O", title="Fold"),
>     y=alt.Y("score:Q", title="Fold score (deviance)"),
> )
>
> mean_line = alt.Chart({"values": [{"y": float(result.cv_score)}]}).mark_rule(
>     color="firebrick", strokeDash=[4, 4], strokeWidth=1.5
> ).encode(y="y:Q")
>
> (bars + mean_line).properties(
>     width="container", height=250,
>     title="Per-fold CV scores (red line = mean)"
> )
> ```
>
> <div class="cell-output cell-output-display" data-execution_count="8">
>
> <style>
>   #altair-viz-961d84d191714001b4a8c874210477ed.vega-embed {
>     width: 100%;
>     display: flex;
>   }
>
>   #altair-viz-961d84d191714001b4a8c874210477ed.vega-embed details,
>   #altair-viz-961d84d191714001b4a8c874210477ed.vega-embed details summary {
>     position: relative;
>   }
> </style>
>
> <div id="altair-viz-961d84d191714001b4a8c874210477ed">
>
> </div>
>
> </div>
>
> </div>
>
> <div id="choosing-the-number-of-folds" class="section level1">
>
> # Choosing the number of folds
>
> <div id="a402298a" class="cell" data-execution_count="9">
>
> ``` python
> # Compare 5-fold vs 10-fold vs leave-one-out-ish (n_folds=n)
> for k in [5, 10, 20]:
>     cv = wk.cross_validate("y ~ s(x)", data, n_folds=k, method="REML", seed=23)
>     print(f"{k:2d}-fold: CV = {cv.cv_score:.4f} (SE = {cv.cv_se:.4f})")
> ```
>
> <div class="cell-output cell-output-stdout">
>
>      5-fold: CV = 0.0992 (SE = 0.0052)
>     10-fold: CV = 0.0989 (SE = 0.0067)
>     20-fold: CV = 0.0986 (SE = 0.0064)
>
> </div>
>
> </div>
>
> > **Tip: How many folds?**
> >
> > - **5-fold**: faster, slightly higher bias, lower variance. Good for large datasets.
> > - **10-fold** (default): a good balance for most datasets.
> > - **20-fold or more**: lower bias but higher variance and slower. Useful for small datasets where you want to use as much training data as possible per fold.
>
> <div id="smooth-selection-via-cross-validation" class="section level1">
>
> # Smooth selection via cross-validation
>
> The `select=True` option enables double-penalty smooth selection, which can shrink entire smooth terms to zero. Cross-validation can verify whether this helps:
>
> <div id="e960f79e" class="cell" data-execution_count="10">
>
> ``` python
> # Restore the original data from the first example
> rng = np.random.default_rng(23)
> n = 400
> x = np.linspace(0, 2 * np.pi, n)
> y = np.sin(x) + 0.5 * np.sin(3 * x) + rng.normal(0, 0.3, n)
> data = {"x": x, "y": y}
>
> # Add a noise variable that should be selected out
> noise = rng.uniform(0, 1, n)
> data_noise = {"x": x, "noise": noise, "y": y}
>
> cv_no_select = wk.cross_validate(
>     "y ~ s(x) + s(noise)", data_noise, method="REML", seed=23
> )
> cv_select = wk.cross_validate(
>     "y ~ s(x) + s(noise)", data_noise, method="REML", select=True, seed=23
> )
>
> print(f"Without selection: {cv_no_select.cv_score:.4f}")
> print(f"With selection:    {cv_select.cv_score:.4f}")
> ```
>
> <div class="cell-output cell-output-stdout">
>
>     Without selection: 0.1002
>     With selection:    0.0994
>
> </div>
>
> </div>
>
> <div id="reproducibility" class="section level1">
>
> # Reproducibility
>
> Pass `seed` to get reproducible fold assignments:
>
> <div id="20c1566c" class="cell" data-execution_count="11">
>
> ``` python
> cv1 = wk.cross_validate("y ~ s(x)", data, method="REML", seed=123)
> cv2 = wk.cross_validate("y ~ s(x)", data, method="REML", seed=123)
>
> print(f"Same seed, same result: {cv1.cv_score == cv2.cv_score}")
> ```
>
> <div class="cell-output cell-output-stdout">
>
>     Same seed, same result: True
>
> </div>
>
> </div>
>
> Without a seed, the fold assignment is randomised on each call.
>
> <div id="the-cvresult-object" class="section level1">
>
> # The CVResult object
>
> [CVResult](../reference/CVResult.md#whittaker.CVResult) is a simple dataclass with four fields:
>
> | Field       | Type      | Description                          |
> |-------------|-----------|--------------------------------------|
> | `cv_score`  | `float`   | Mean out-of-sample loss across folds |
> | `cv_scores` | `NDArray` | Per-fold loss values                 |
> | `cv_se`     | `float`   | Standard error of the mean score     |
> | `n_folds`   | `int`     | Number of folds                      |
>
> <div id="practical-workflow" class="section level1">
>
> # Practical workflow
>
> A typical model-selection workflow combines cross-validation with REML fitting:
>
> 1.  **Candidate models**: vary the formula, basis dimension, or family
> 2.  **Cross-validate each**: use the same `seed` so folds are identical
> 3.  **Compare scores**: pick the model with the lowest CV score, or the simplest model within one SE of the best
> 4.  **Final fit**: refit the chosen model on all the data
>
> <div id="c1f8b00e" class="cell" data-execution_count="12">
>
> ``` python
> # Step 1-3: compare candidates
> candidates = {
>     "s(x, k=5)": "y ~ s(x, k=5)",
>     "s(x, k=10)": "y ~ s(x, k=10)",
>     "s(x, k=20)": "y ~ s(x, k=20)",
>     "s(x) + s(noise)": "y ~ s(x) + s(noise)",
> }
>
> best_name, best_score = None, float("inf")
> for name, formula in candidates.items():
>     cv = wk.cross_validate(formula, data_noise, method="REML", seed=23)
>     flag = ""
>     if cv.cv_score < best_score:
>         best_score = cv.cv_score
>         best_name = name
>         flag = " <-- best"
>     print(f"  {name:20s}: {cv.cv_score:.4f} (SE {cv.cv_se:.4f}){flag}")
>
> # Step 4: refit the winner on all data
> print(f"\nBest model: {best_name}")
> final = wk.GAM(candidates[best_name])
> final.fit(data_noise, method="REML")
> print(f"Final EDF: {final.edf_total:.1f}")
> ```
>
> <div class="cell-output cell-output-stdout">
>
>       s(x, k=5)           : 0.2044 (SE 0.0113) <-- best
>       s(x, k=10)          : 0.0989 (SE 0.0067) <-- best
>       s(x, k=20)          : 0.1004 (SE 0.0066)
>       s(x) + s(noise)     : 0.1002 (SE 0.0069)
>
>     Best model: s(x, k=10)
>     Final EDF: 9.9
>
> </div>
>
> </div>
>
> > **Tip: The one-SE rule**
> >
> > A common heuristic: instead of picking the model with the lowest CV score, pick the simplest model whose score is within one standard error of the best. This guards against overfitting to the validation folds and tends to produce more parsimonious models.
>
> You can now use K-fold cross-validation to compare models, select basis dimensions, and apply the one-SE rule for parsimonious model selection.
>
> <div id="where-to-go-next" class="section level1">
>
> # Where to go next
>
> - **[Model comparison with LOO](loo.md)**: PSIS-LOO cross-validation for Bayesian model comparison without refitting.
> - **[Model diagnostics](diagnostics.md)**: residual checks and basis adequacy tests as a complement to cross-validation.
> - **[Model fitting](fitting.md)**: how smoothness selection (REML, GCV) relates to cross-validation.
>
> </div>
>
> </div>
>
> </div>
>
> </div>
>
> </div>
>
> </div>
>
> </div>
>
> </div>
>
> </div>