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() 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():
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 containing a list of ComparisonRow objects. You can index into it or use the .best property:
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 (or None for Bayesian fits), scale, and n_obs.
Comparing Bayesian fits
compare() works with VI-fitted models too. Since GCV is not available for Bayesian fits, that column is omitted from the table:
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:
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() 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, which performs sequential deviance-difference tests between nested models.
Use 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.