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.
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.
# 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() method returns a single 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.
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:
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.
When comparing two or more models, collect each one’s 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() method runs a basis dimension adequacy test (the k-index test) for each smooth term.
# 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.
If check() reports a k-index below 1 with a significant p-value, re-fit with a larger k:
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.
# 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.
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.
# 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
# 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.
# 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.
# 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.
# 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() method correctly flagged the k=5 model.
Practical diagnostic workflow
A good diagnostic workflow after fitting any GAM:
model.summary(): check that EDF values make sense, deviance explained is reasonablemodel.goodness_of_fit(): get a compact snapshot of AIC, BIC, adjusted R^2, and other quality metrics in one callwk.check(model): verify basis dimensions are adequate (k-index test)- Residuals vs. fitted: check for patterns indicating wrong family or missing terms
- QQ plot: check the distributional assumption
- If any diagnostic fails: consider increasing
k, changing the family, adding terms, or restructuring the model
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: smoothness selection methods (REML, GCV, ML) that control how flexible the smooth terms are.
- Response families: choosing a family that matches the data-generating process, which is the most common fix when diagnostics reveal distributional problems.
- Posterior predictive checks: a complementary diagnostic that tests whether the model generates realistic data.
- Advanced diagnostics: influence, concurvity, dispersion tests, and quantile residuals for deeper model checking.
- Smoothing parameter sensitivity: check how much predictions change as smoothing parameters vary, to assess robustness of conclusions.
- Prediction and inference: confidence intervals and term-level predictions from a model that passes diagnostics.