# 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()](../reference/GAM.md#whittaker.GAM.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()](../reference/loo_compare.md#whittaker.loo_compare), [waic_compare()](../reference/waic_compare.md#whittaker.waic_compare), or [stacking()](../reference/stacking.md#whittaker.stacking) instead. # Basic usage Fit two or more models and call [anova()](../reference/GAM.md#whittaker.GAM.anova) on any one of them, passing the others as arguments. The models are automatically sorted from simplest to most complex. ``` python 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) ``` ``` python 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. ``` python # 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. ``` python 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: ``` python 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()](../reference/GAM.md#whittaker.GAM.anova) works with just two models for a simple A vs. B test: ``` python 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()](../reference/GAM.md#whittaker.GAM.anova) requires frequentist fits (GCV, REML, or ML). For Bayesian model comparison, use [loo_compare()](../reference/loo_compare.md#whittaker.loo_compare), [waic_compare()](../reference/waic_compare.md#whittaker.waic_compare), or [stacking()](../reference/stacking.md#whittaker.stacking). # Where to go next - **[Model comparison with LOO](loo.md)**: PSIS-LOO for non-nested model comparison. - **[Model comparison with WAIC](waic.md)**: WAIC as an alternative to LOO. - **[Model averaging with stacking](stacking.md)**: combine multiple models instead of choosing one. - **[Model fitting](fitting.md)**: how smoothness selection criteria (REML, GCV, ML) affect the models being compared. - **[Derivatives and marginal effects](derivatives.md)**: interpret the smooth effects once you have settled on a model.