ANOVA for GAMs

When you have two or more nested models (a simpler model that is a special case of a more complex one) you can test whether the additional complexity is justified by comparing their deviances. The anova() method performs sequential deviance-difference tests.

This is different from information-criterion comparisons (AIC, BIC) or predictive criteria (LOO, WAIC). Those methods do not require the models to be nested and estimate out-of-sample predictive accuracy. ANOVA tests a null hypothesis (“does the more complex model explain significantly more deviance?”) using the sampling distribution of the deviance difference.

When to use ANOVA

Use model.anova() when:

  • you have a clear nesting hierarchy (e.g., a linear model vs. a smooth model, or a model with fewer terms vs. one with more)
  • you want a formal p-value for whether the extra terms are needed
  • the models use the same family and are fitted to the same data

For non-nested models, use loo_compare(), waic_compare(), or stacking() instead.

Basic usage

Fit two or more models and call anova() on any one of them, passing the others as arguments. The models are automatically sorted from simplest to most complex.

import numpy as np
import whittaker as wk

rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2 * np.pi, n)
y = np.sin(x) + 0.3 * np.cos(3 * x) + rng.normal(0, 0.3, n)

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

# Model 1: intercept only (null model)
m_null = wk.GAM("y ~ 1").fit(data)

# Model 2: linear effect
m_linear = wk.GAM("y ~ x").fit(data)

# Model 3: smooth effect with moderate flexibility
m_smooth = wk.GAM("y ~ s(x, k=10)").fit(data)

# Model 4: smooth effect with high flexibility
m_flex = wk.GAM("y ~ s(x, k=20)").fit(data)
result = m_null.anova(m_linear, m_smooth, m_flex)
print(result)
Analysis of Deviance Table

  Model   Resid.Df    Resid.Dev       Df     Deviance          F         Pr(>F)
  ----- ---------- ------------ -------- ------------ ---------- --------------
      1     299.00     182.9233                                                
      2     298.00     101.3513     1.00      81.5720   837.2556    1.51127e-88
      3     290.55      28.2360     7.45      73.1153   100.7087    2.37349e-76
      4     288.69      28.1266     1.86       0.1094     0.6050       0.534843

Reading the ANOVA table

The table has one row per model, sorted from simplest (fewest EDF) to most complex. For each successive pair of models:

  • Df: the difference in effective degrees of freedom between the two models. This measures the additional complexity.
  • Deviance: the reduction in deviance from the simpler to the more complex model.
  • Statistic: an F-statistic (for unknown-scale families like Gaussian and Gamma) or chi-squared statistic (for known-scale families like Poisson and Binomial).
  • p-value: the probability of seeing this large a deviance reduction by chance, under the null that the simpler model is adequate.

A small p-value means the more complex model is significantly better.

# Access individual rows
for i, row in enumerate(result.rows):
    p_str = f"{row.p_value:.4g}" if row.p_value is not None else "---"
    print(f"Model {i+1}: Resid.Df={row.resid_df:.1f}, "
          f"Resid.Dev={row.resid_dev:.2f}, p={p_str}")
Model 1: Resid.Df=299.0, Resid.Dev=182.92, p=---
Model 2: Resid.Df=298.0, Resid.Dev=101.35, p=1.511e-88
Model 3: Resid.Df=290.5, Resid.Dev=28.24, p=2.373e-76
Model 4: Resid.Df=288.7, Resid.Dev=28.13, p=0.5348

Interpreting the results

In the example above, each step adds complexity:

  1. Null → Linear: tests whether x has any linear effect on y. A significant p-value means a linear trend is present.
  2. Linear → Smooth (k=10): tests whether the relationship is nonlinear. A significant p-value means the smooth captures structure that a straight line misses.
  3. Smooth (k=10) → Smooth (k=20): tests whether extra flexibility is needed. If the p-value is not significant, the simpler smooth is adequate (the additional basis functions are being penalized away).

This sequential testing workflow is the standard approach for building up a GAM: start simple, add complexity, and stop when additional terms no longer yield significant deviance reductions.

Poisson example with chi-squared tests

For known-scale families, the test uses a chi-squared distribution instead of an F distribution.

from whittaker.families.poisson import Poisson

rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2 * np.pi, n)
lam = np.exp(0.8 * np.sin(x))
y_pois = rng.poisson(lam).astype(float)

pois_data = {"x": x, "y": y_pois}

m1 = wk.GAM("y ~ 1", family=Poisson()).fit(pois_data)
m2 = wk.GAM("y ~ x", family=Poisson()).fit(pois_data)
m3 = wk.GAM("y ~ s(x)", family=Poisson()).fit(pois_data)

result_pois = m1.anova(m2, m3)
print(result_pois)
Analysis of Deviance Table

  Model   Resid.Df    Resid.Dev       Df     Deviance      Chisq     Pr(>Chisq)
  ----- ---------- ------------ -------- ------------ ---------- --------------
      1     299.00     443.7721                                                
      2     298.00     374.4750     1.00      69.2970    69.2970    8.46962e-17
      3     294.27     314.9843     3.73      59.4907    59.4907     2.4478e-12

The test type is shown in the result:

print(f"Test type: {result_pois.test}")
print(f"Scale used: {result_pois.scale}")
Test type: Chisq
Scale used: 1.0

Two-model comparison

anova() works with just two models for a simple A vs. B test:

result_ab = m_linear.anova(m_smooth)
print(result_ab)
Analysis of Deviance Table

  Model   Resid.Df    Resid.Dev       Df     Deviance          F         Pr(>F)
  ----- ---------- ------------ -------- ------------ ---------- --------------
      1     298.00     101.3513                                                
      2     290.55      28.2360     7.45      73.1153   100.9635    1.82492e-76

This is the most common use case: you have a model and want to know whether adding a smooth term (or extra covariates) is justified.

Requirements and limitations

  • Same family: all models must use the same response family. Comparing a Gaussian model against a Poisson model is not meaningful.
  • Same data: all models must be fitted to the same observations. The method checks this and raises an error if the observation counts differ.
  • Nesting: the test is valid when models are nested (each simpler model is a special case of the more complex one). For non-nested models, the p-values are approximate at best.
  • Frequentist fits only: anova() requires frequentist fits (GCV, REML, or ML). For Bayesian model comparison, use loo_compare(), waic_compare(), or stacking().

Where to go next