# Model comparison with compare()

When you have several candidate models for the same data, you need a quick way to line them up and see which one fits best. The [compare()](../reference/compare.md#whittaker.compare) function collects AIC, BIC, deviance explained, adjusted R-squared, EDF, and GCV from each model and presents them in a single sorted table.


# Basic usage

Fit two or more models on the same data, then pass them all to [compare()](../reference/compare.md#whittaker.compare):


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

rng = np.random.default_rng(0)
n = 200
x = np.linspace(0, 2 * np.pi, n)
y = np.sin(x) + rng.normal(0, 0.3, n)
data = {"x": x, "y": y}

m_linear = wk.GAM("y ~ x").fit(data)
m_smooth = wk.GAM("y ~ s(x, k=5)").fit(data)
m_flex   = wk.GAM("y ~ s(x, k=15)").fit(data)

result = wk.compare(m_linear, m_smooth, m_flex)
print(result)
```


    Model Comparison (3 models, 200 observations)

      # Formula                               AIC     ΔAIC        BIC  Dev.Expl.   Adj.R²     EDF        GCV
    --- ------------------------------ ---------- -------- ---------- ---------- -------- ------- ----------
      1 y ~ s(x, k=5)                       78.32    +0.00      94.67     86.7%   0.8636     5.0   0.086642
      2 y ~ s(x, k=15)                      79.21    +0.89     103.14     87.0%   0.8646     7.3   0.087059
      3 y ~ x                              332.62  +254.30     339.22     51.1%   0.5065     2.0   0.308909


Models are sorted by AIC (lowest first). The **ΔAIC** column shows how far each model is from the best: a ΔAIC of 0 marks the winner, and values above ~10 indicate models with essentially no empirical support relative to the best.


# Accessing individual rows

The result is a [ComparisonResult](../reference/ComparisonResult.md#whittaker.ComparisonResult) containing a list of [ComparisonRow](../reference/ComparisonRow.md#whittaker.ComparisonRow) objects. You can index into it or use the `.best` property:


``` python
best = result.best
print(f"Best model: {best.label}")
print(f"  AIC:              {best.aic:.2f}")
print(f"  BIC:              {best.bic:.2f}")
print(f"  Dev. explained:   {best.deviance_explained:.1%}")
print(f"  Adj. R²:          {best.r_squared_adj:.4f}")
print(f"  EDF:              {best.edf_total:.1f}")
```


    Best model: y ~ s(x, k=5)
      AIC:              78.32
      BIC:              94.67
      Dev. explained:   86.7%
      Adj. R²:          0.8636
      EDF:              5.0


Each row also carries [gcv_score](../reference/GAM.md#whittaker.GAM.gcv_score) (or `None` for Bayesian fits), `scale`, and `n_obs`.


# Comparing Bayesian fits

[compare()](../reference/compare.md#whittaker.compare) works with VI-fitted models too. Since GCV is not available for Bayesian fits, that column is omitted from the table:


``` python
m_vi1 = wk.GAM("y ~ s(x, k=5)").fit(data, method="VI")
m_vi2 = wk.GAM("y ~ s(x, k=15)").fit(data, method="VI")

print(wk.compare(m_vi1, m_vi2))
```


    Model Comparison (2 models, 200 observations)

      # Formula                               AIC     ΔAIC        BIC  Dev.Expl.   Adj.R²     EDF
    --- ------------------------------ ---------- -------- ---------- ---------- -------- -------
      1 y ~ s(x, k=5)                       78.32    +0.00      94.72     86.7%   0.8636     5.0
      2 y ~ s(x, k=15)                      80.16    +1.83     110.66     87.2%   0.8653     9.2


# Non-Gaussian models

The same interface works for any response family. Here we compare Poisson models:


``` python
from whittaker.families.poisson import Poisson

x_p = np.linspace(0, 3, 150)
y_p = rng.poisson(np.exp(0.5 * np.sin(x_p))).astype(float)
pois_data = {"x": x_p, "y": y_p}

p1 = wk.GAM("y ~ x", family=Poisson()).fit(pois_data)
p2 = wk.GAM("y ~ s(x)", family=Poisson()).fit(pois_data)

print(wk.compare(p1, p2))
```


    Model Comparison (2 models, 150 observations)

      # Formula                               AIC     ΔAIC        BIC  Dev.Expl.   Adj.R²     EDF        GCV
    --- ------------------------------ ---------- -------- ---------- ---------- -------- ------- ----------
      1 y ~ s(x)                           467.18    +0.00     481.78      5.8%   0.0261     4.8   1.268470
      2 y ~ x                              470.13    +2.95     476.15      1.2%  -0.0014     2.0   1.279296


# When to use compare()

[compare()](../reference/compare.md#whittaker.compare) is most useful when you have a small set of candidate models and want a quick side-by-side summary. It is *not* a formal hypothesis test. For that, see [ANOVA for GAMs](anova.md), which performs sequential deviance-difference tests between nested models.

Use [compare()](../reference/compare.md#whittaker.compare) for:

- **Model selection**: choosing among several candidate formulas or basis dimensions.
- **Quick screening**: narrowing down a large set of models before deeper analysis.
- **Reporting**: producing a tidy summary table for a paper or presentation.

> **Tip: Combining with other tools**
>
> - Use [smoothing parameter sensitivity](sensitivity.md) to check whether the winning model's conclusions are robust to the choice of smoothing parameters.
> - Use [cross-validation](cross-validation.md) for out-of-sample prediction accuracy, which [compare()](../reference/compare.md#whittaker.compare) does not assess.
> - For Bayesian model comparison, see [LOO](loo.md), [WAIC](waic.md), and [stacking](stacking.md).


# Where to go next

- **[Model diagnostics](diagnostics.md)**: residual analysis and goodness-of-fit metrics for individual models.
- **[ANOVA for GAMs](anova.md)**: formal deviance-difference tests for nested model comparison.
- **[Cross-validation](cross-validation.md)**: out-of-sample evaluation with K-fold CV.
- **[Model comparison with LOO](loo.md)**: leave-one-out cross-validation for Bayesian models.