# Model diagnostics
Always check your model after fitting. GAMs can fail silently when a basis dimension is too small, when the data contain structure the smoother cannot capture, or when the distributional assumptions are wrong. This page covers the diagnostic tools Whittaker provides.
# A model to diagnose
We start by fitting a model to data with a known structure, so we can verify that the diagnostics behave sensibly.
``` python
import numpy as np
import whittaker as wk
# Generate data: sin(x) + noise
rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2 * np.pi, n)
y = np.sin(x) + rng.normal(0, 0.3, n)
# Fit with default settings
model = wk.GAM("y ~ s(x)")
model.fit({"x": x, "y": y}, method="REML")
print(model.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Gaussian(link='identity')
Inference: REML
Observations: 300
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) -0.0014 0.0180 -0.076 0.9395
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 7.47 8 1457.725 < 1e-16
Total EDF: 8.47
Scale est: 0.097016
Deviance: 28.2828
Null dev: 170.3338
Dev. expl: 83.4%
GCV score: 0.099835
AIC: 159.97
BIC: 191.35
# Model summary
The `summary()` method is the first thing to inspect. It reports:
- **Effective degrees of freedom (EDF)** for each smooth term. An EDF near 1 means the smooth is approximately linear. An EDF close to k - 1 (the maximum) suggests the basis may be too small.
- **Deviance explained**: the proportion of null deviance accounted for by the model, analogous to R^2 in linear regression.
- **Scale estimate** \hat\phi: for Gaussian models, this is the estimated residual variance \hat\sigma^2.
``` python
# Access individual fit statistics
print(f"EDF total: {model.edf_total:.1f}")
print(f"Scale (sigma^2): {model.scale:.4f}")
print(f"Deviance: {model.deviance:.2f}")
```
EDF total: 8.5
Scale (sigma^2): 0.0970
Deviance: 28.28
# Goodness of fit
Rather than accessing each fit statistic individually, the [goodness_of_fit()](../reference/GAM.md#whittaker.GAM.goodness_of_fit) method returns a single [GoodnessOfFit](../reference/GoodnessOfFit.md#whittaker.GoodnessOfFit) object that bundles together the most commonly used measures of model quality. This is especially convenient when you want a quick snapshot of how well the model fits the data, or when comparing several models side by side.
``` python
gof = model.goodness_of_fit()
print(gof)
```
GoodnessOfFit
Deviance: 28.2828
Null deviance: 170.3338
Deviance explained: 83.4%
Adj. R-squared: 0.8291
AIC: 159.97
BIC: 191.35
GCV: 0.099835
Scale: 0.097016
EDF total: 8.47
Observations: 300
The printed summary includes deviance explained, adjusted R^2, AIC, BIC, GCV (when available), the scale estimate, and the total effective degrees of freedom. Adjusted R^2 penalizes for model complexity using the effective degrees of freedom rather than the raw parameter count, so it gives a fairer comparison between models of different flexibility:
R^2\_{\text{adj}} = 1 - (1 - R^2) \cdot \frac{n - 1}{n - \text{EDF} - 1}
Each field is accessible as a plain attribute, which makes it straightforward to build comparison tables or apply decision rules programmatically:
``` python
print(f"Deviance explained: {gof.deviance_explained:.1%}")
print(f"Adjusted R-squared: {gof.r_squared_adj:.4f}")
print(f"AIC: {gof.aic:.2f}")
print(f"BIC: {gof.bic:.2f}")
print(f"GCV: {gof.gcv_score:.6f}")
print(f"Observations: {gof.n_obs}")
```
Deviance explained: 83.4%
Adjusted R-squared: 0.8291
AIC: 159.97
BIC: 191.35
GCV: 0.099835
Observations: 300
For Bayesian fits (VI or MCMC), the GCV score is not applicable and is reported as `None`. All other fields remain available.
> **Tip: Comparing models with goodness of fit**
>
> When comparing two or more models, collect each one's [GoodnessOfFit](../reference/GoodnessOfFit.md#whittaker.GoodnessOfFit) and compare the metrics that matter for your goal. Lower AIC or BIC favors predictive accuracy with a complexity penalty. Higher deviance explained and adjusted R^2 indicate better in-sample fit.
# Basis dimension adequacy
The basis dimension k sets the maximum complexity of each smooth. If k is too small, the model cannot capture the true function shape. The [check()](../reference/check.md#whittaker.check) method runs a basis dimension adequacy test (the k-index test) for each smooth term.
``` python
# Run the check
wk.check(model)
```
The k-index is based on the ratio of the residual variance estimated from neighboring residuals to the overall residual variance. A k-index below 1 with a significant p-value is a warning that the basis dimension may be too small.
> **Note: When to increase k**
>
> If [check()](../reference/check.md#whittaker.check) reports a k-index below 1 with a significant p-value, re-fit with a larger `k`:
>
> ``` python
> model2 = wk.GAM("y ~ s(x, k=20)")
> model2.fit(data, method="REML")
> wk.check(model2)
> ```
>
> Keep increasing `k` until the k-index test is no longer significant. The smoothing parameter selection (REML) will prevent overfitting even with a large `k`. The penalty shrinks away unnecessary complexity.
# Residual analysis
Residuals are the primary tool for checking distributional assumptions. For Gaussian models, well- behaved residuals should be approximately normal with constant variance.
``` python
# Deviance residuals (default)
resids = model.residuals
print(f"Residual shape: {resids.shape}")
print(f"Mean: {resids.mean():.4f}")
print(f"Std: {resids.std():.4f}")
```
Residual shape: (300,)
Mean: 0.0000
Std: 0.3070
## Residuals vs. fitted values
Plotting residuals against fitted values checks the constant-variance assumption. The plot should show no systematic pattern, just a random scatter around zero.
``` python
import altair as alt
# Get fitted values and residuals
fitted = model.fitted_values
resids = model.residuals
# Build the plot data
resid_data = [
{"fitted": float(fitted[i]), "residual": float(resids[i])}
for i in range(len(fitted))
]
alt.Chart({"values": resid_data}).mark_circle(
size=15, opacity=0.4, color="steelblue"
).encode(
x=alt.X("fitted:Q", title="Fitted values"),
y=alt.Y("residual:Q", title="Deviance residuals"),
).properties(
width="container", height=300,
title="Residuals vs. fitted values"
) + alt.Chart({"values": [{"y": 0}]}).mark_rule(
color="firebrick", strokeDash=[4, 4]
).encode(y="y:Q")
```
A funnel shape (variance increasing with fitted values) suggests a non-constant variance and may indicate that a different family or link function is needed.
## QQ plot of residuals
A quantile-quantile plot compares the distribution of residuals against a theoretical normal distribution. Points should fall close to the diagonal line.
``` python
# Compute theoretical and sample quantiles for QQ plot
sorted_resids = np.sort(resids)
n_resids = len(sorted_resids)
theoretical = np.array([
float(x) for x in np.quantile(
rng.normal(0, 1, 10000),
np.linspace(0.5 / n_resids, 1 - 0.5 / n_resids, n_resids)
)
])
qq_data = [
{"theoretical": float(theoretical[i]), "sample": float(sorted_resids[i])}
for i in range(n_resids)
]
# QQ points
qq_points = alt.Chart({"values": qq_data}).mark_circle(
size=15, opacity=0.4, color="steelblue"
).encode(
x=alt.X("theoretical:Q", title="Theoretical quantiles"),
y=alt.Y("sample:Q", title="Sample quantiles"),
)
# Reference line
ref_min = min(theoretical.min(), sorted_resids.min())
ref_max = max(theoretical.max(), sorted_resids.max())
ref_data = [{"x": float(ref_min), "y": float(ref_min)}, {"x": float(ref_max), "y": float(ref_max)}]
ref_line = alt.Chart({"values": ref_data}).mark_line(
color="firebrick", strokeDash=[4, 4]
).encode(x="x:Q", y="y:Q")
(qq_points + ref_line).properties(
width="container", height=400,
title="Normal QQ plot of residuals"
)
```
Systematic departures from the line indicate non-normality. Heavy tails (S-shaped departures) suggest a heavier-tailed family might be more appropriate.
## Histogram of residuals
``` python
# Histogram of deviance residuals
hist_data = [{"residual": float(r)} for r in resids]
alt.Chart({"values": hist_data}).mark_bar(
opacity=0.7, color="steelblue"
).encode(
x=alt.X("residual:Q", bin=alt.Bin(maxbins=30), title="Deviance residuals"),
y=alt.Y("count():Q", title="Frequency"),
).properties(
width="container", height=250,
title="Distribution of residuals"
)
```
# Diagnosing an inadequate model
To illustrate what diagnostics look like when something is wrong, let's deliberately under-fit by using too few basis functions for a complex signal.
``` python
# Generate data with high-frequency oscillation
rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 4 * np.pi, n)
y = np.sin(x) + 0.5 * np.sin(3 * x) + rng.normal(0, 0.3, n)
# Fit with k=5 (too few basis functions for this signal)
model_bad = wk.GAM("y ~ s(x, k=5)")
model_bad.fit({"x": x, "y": y}, method="REML")
# Check: the k-index should flag the problem
wk.check(model_bad)
```
The k-index test flags the model as inadequate. The fix is to increase the basis dimension.
``` python
# Re-fit with more basis functions
model_good = wk.GAM("y ~ s(x, k=20)")
model_good.fit({"x": x, "y": y}, method="REML")
# The check should now pass
wk.check(model_good)
```
With k = 20, the k-index test passes. Now compare the two fits visually to see the difference.
``` python
# Compare fits visually
x_plot = np.linspace(0, 4 * np.pi, 300)
pred_bad = model_bad.predict({"x": x_plot})
pred_good = model_good.predict({"x": x_plot})
true_vals = np.sin(x_plot) + 0.5 * np.sin(3 * x_plot)
# Build plot data
plot_data = []
for i in range(len(x_plot)):
plot_data.append({"x": float(x_plot[i]), "y": float(pred_bad.values[i]), "model": "k=5 (underfit)"})
plot_data.append({"x": float(x_plot[i]), "y": float(pred_good.values[i]), "model": "k=20 (adequate)"})
plot_data.append({"x": float(x_plot[i]), "y": float(true_vals[i]), "model": "Truth"})
alt.Chart({"values": plot_data}).mark_line().encode(
x=alt.X("x:Q"),
y=alt.Y("y:Q"),
color=alt.Color("model:N", title="Model"),
strokeDash=alt.condition(
alt.datum.model == "Truth",
alt.value([4, 4]),
alt.value([0])
),
).properties(width="container", height=300, title="Effect of basis dimension on fit quality")
```
The underfit model (k=5) misses the high-frequency component entirely, while k=20 captures both sine components. The [check()](../reference/check.md#whittaker.check) method correctly flagged the k=5 model.
# Practical diagnostic workflow
A good diagnostic workflow after fitting any GAM:
1. **`model.summary()`**: check that EDF values make sense, deviance explained is reasonable
2. **`model.goodness_of_fit()`**: get a compact snapshot of AIC, BIC, adjusted R^2, and other quality metrics in one call
3. **`wk.check(model)`**: verify basis dimensions are adequate (k-index test)
4. **Residuals vs. fitted**: check for patterns indicating wrong family or missing terms
5. **QQ plot**: check the distributional assumption
6. **If any diagnostic fails**: consider increasing `k`, changing the family, adding terms, or restructuring the model
> **Tip: The most common fix**
>
> The most common diagnostic issue is an inadequate basis dimension. REML will never overfit even with a generous `k`, so it is always safe to increase `k`. Start with the default (10), check, and double if the k-index test is significant. Repeat until the test passes.
You can now inspect a fitted GAM for basis adequacy, residual patterns, and distributional assumptions using the diagnostic tools on this page.
# Where to go next
- **[Model fitting](fitting.md)**: smoothness selection methods (REML, GCV, ML) that control how flexible the smooth terms are.
- **[Response families](families.md)**: choosing a family that matches the data-generating process, which is the most common fix when diagnostics reveal distributional problems.
- **[Posterior predictive checks](ppc.md)**: a complementary diagnostic that tests whether the model generates realistic data.
- **[Advanced diagnostics](advanced-diagnostics.md)**: influence, concurvity, dispersion tests, and quantile residuals for deeper model checking.
- **[Smoothing parameter sensitivity](sensitivity.md)**: check how much predictions change as smoothing parameters vary, to assess robustness of conclusions.
- **[Prediction and inference](prediction.md)**: confidence intervals and term-level predictions from a model that passes diagnostics.