Whittaker fits GAMs by penalized iteratively reweighted least squares (P-IRLS), a generalization of weighted least squares that handles non-Gaussian responses by iterating over a sequence of penalized weighted regression problems. This page describes the fitting algorithm, smoothness selection criteria, and the practical controls you have over the fitting process.
The P-IRLS algorithm
For a given set of smoothing parameters \boldsymbol{\lambda} = (\lambda_1, \ldots, \lambda_J), fitting proceeds in five steps:
- Initialize: set $= $
family.initialize(y), \eta = g(\mu). The initialization is family-specific: for Gaussian data \mu = y, for Poisson \mu = y + 0.1, for binomial \mu = (y + 0.5)/2.
- Pseudo-data: compute the working response z = \eta + (y - \mu)\, g'(\mu) This is a first-order Taylor expansion of the link function around the current \mu, converting the non-Gaussian problem into a weighted least-squares problem.
- Working weights: compute W = \frac{1}{\bigl[g'(\mu)\bigr]^2 \, V(\mu)} where V(\mu) is the variance function of the family. These weights account for both the link function and the heteroscedasticity of the response.
- Penalized WLS: solve the penalized weighted least-squares system \bigl(X^\top W X + \textstyle\sum_{j=1}^{J} \lambda_j S_j\bigr)\,\hat\beta = X^\top W z where S_j is the penalty matrix for the j-th smooth term. This is the core computational step. The penalty matrices S_j are positive semi-definite and penalize the roughness of each smooth, pulling f_j toward the penalty null space (typically low-order polynomials).
- Update: compute \eta = X\hat\beta and \mu = g^{-1}(\eta).
Repeat steps 2–5 until the relative change in penalized deviance falls below 10^{-7}.
The Gaussian identity-link special case
For the Gaussian family with an identity link, P-IRLS reduces to a single penalized least-squares solve. The link function is the identity (g(\mu) = \mu), so g'(\mu) = 1, and the variance function is constant (V(\mu) = 1). The working weights become W = I and the pseudo-data become z = y, so the system collapses to:
\bigl(X^\top X + \textstyle\sum_j \lambda_j S_j\bigr)\,\hat\beta = X^\top y
No iteration is required. This makes Gaussian GAMs substantially faster to fit than their non-Gaussian counterparts.
Smoothness selection
The roughness of each smooth is controlled by a non-negative smoothing parameter \lambda_j. Larger \lambda_j produces a smoother f_j while \lambda_j = 0 removes the penalty entirely (interpolating spline). Rather than setting \lambda_j by hand, Whittaker selects it automatically by optimizing one of three criteria.
GCV (generalized cross-validation)
GCV selects \boldsymbol\lambda by minimizing:
\text{GCV}(\boldsymbol\lambda) = \frac{n\, D(y,\hat\mu)}{[n - \text{tr}(\mathbf{A})]^2}
where D is the deviance and \mathbf{A} = X(X^\top W X + \sum_j \lambda_j S_j)^{-1}X^\top W is the influence (hat) matrix. The trace \text{tr}(\mathbf{A}) is the effective degrees of freedom of the model.
GCV is an approximation to leave-one-out cross-validation that avoids refitting the model n times. It is minimized over \log\boldsymbol\lambda using a Newton method with analytical gradient and Hessian.
Properties:
- GCV is well-understood theoretically and widely used.
- It tends to undersmooth in finite samples — occasional fits with too many effective degrees of freedom.
- It can be sensitive to influential observations because it weights all residuals equally.
model = wk.GAM("y ~ s(x)")
model.fit(data, method="GCV")
REML (restricted maximum likelihood)
REML treats the smooth coefficients as random effects and maximizes the restricted log-likelihood:
\ell_{\text{REML}}(\boldsymbol\lambda, \phi) = -\frac{1}{2}\Bigl[
n \log(2\pi\phi) + \frac{D}{\phi} +
\log\bigl|X^\top W X + \textstyle\sum_j \lambda_j S_j\bigr| -
\log\bigl|\textstyle\sum_j \lambda_j S_j^{+}\bigr|
\Bigr]
where S_j^{+} denotes the pseudoinverse restricted to the range space of S_j, and \phi is the scale parameter.
Advantages over GCV:
- REML has better finite-sample properties: it is less prone to undersmoothing.
- It is more stable when the true function is close to the penalty null space.
- REML accounts for the uncertainty in \hat\beta when estimating \phi, producing better variance estimates.
- It is the recommended default in both Whittaker and
mgcv.
model = wk.GAM("y ~ s(x)")
model.fit(data, method="REML") # recommended default
ML (maximum likelihood)
ML maximizes the full marginal log-likelihood rather than the restricted version:
\ell_{\text{ML}}(\boldsymbol\lambda, \phi) = -\frac{1}{2}\Bigl[
n \log(2\pi\phi) + \frac{D}{\phi} +
\log\bigl|X^\top W X + \textstyle\sum_j \lambda_j S_j\bigr|
\Bigr]
The key difference from REML is the absence of the \log|\sum_j \lambda_j S_j^{+}| term. This means ML does not adjust for the degrees of freedom consumed by the fixed effects, which can lead to slight undersmoothing relative to REML — the same bias that makes ML variance estimates biased downward in classical linear models.
When to use ML: ML is appropriate when you need to compare models with different fixed-effect structures using likelihood ratio tests, because REML likelihoods are not comparable across models with different fixed effects.
model = wk.GAM("y ~ s(x)")
model.fit(data, method="ML")
Comparing GCV, REML, and ML
import numpy as np
import whittaker as wk
# Generate noisy data with a smooth underlying function
rng = np.random.default_rng(23)
n = 150
x = np.linspace(0, 2 * np.pi, n)
y_true = np.sin(x) + 0.3 * np.cos(3 * x)
y = y_true + rng.normal(0, 0.4, n)
data = {"x": x, "y": y}
# Fit with each method
model_gcv = wk.GAM("y ~ s(x, k=20)")
model_gcv.fit(data, method="GCV")
model_reml = wk.GAM("y ~ s(x, k=20)")
model_reml.fit(data, method="REML")
model_ml = wk.GAM("y ~ s(x, k=20)")
model_ml.fit(data, method="ML")
# Compare EDF and smoothing parameters
print("Method | EDF total | Smoothing param")
print("--------|-----------|----------------")
print(f"GCV | {model_gcv.edf_total:9.3f} | {model_gcv.smoothing_params[0]:.4f}")
print(f"REML | {model_reml.edf_total:9.3f} | {model_reml.smoothing_params[0]:.4f}")
print(f"ML | {model_ml.edf_total:9.3f} | {model_ml.smoothing_params[0]:.4f}")
Method | EDF total | Smoothing param
--------|-----------|----------------
GCV | 11.054 | 0.6030
REML | 9.853 | 1.0530
ML | 9.853 | 1.0530
GCV typically selects the smallest smoothing parameter (largest EDF), REML the largest (fewest EDF), and ML falls between the two.
import altair as alt
# Build predictions for each method
x_plot = np.linspace(0, 2 * np.pi, 200)
preds_gcv = model_gcv.predict({"x": x_plot})
preds_reml = model_reml.predict({"x": x_plot})
preds_ml = model_ml.predict({"x": x_plot})
# Assemble plot data
plot_records = []
for i in range(len(x_plot)):
plot_records.append({"x": float(x_plot[i]), "fit": float(preds_gcv.values[i]), "Method": "GCV"})
plot_records.append({"x": float(x_plot[i]), "fit": float(preds_reml.values[i]), "Method": "REML"})
plot_records.append({"x": float(x_plot[i]), "fit": float(preds_ml.values[i]), "Method": "ML"})
# Observed data
obs_records = [{"x": float(x[i]), "y": float(y[i])} for i in range(n)]
# True function (recomputed on the plot grid)
y_true_plot = np.sin(x_plot) + 0.3 * np.cos(3 * x_plot)
true_records = [{"x": float(x_plot[i]), "y": float(y_true_plot[i])} for i in range(len(x_plot))]
points = alt.Chart({"values": obs_records}).mark_circle(
size=12, opacity=0.25, color="gray"
).encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="y"),
)
lines = alt.Chart({"values": plot_records}).mark_line(strokeWidth=2).encode(
x="x:Q",
y="fit:Q",
color=alt.Color("Method:N", scale=alt.Scale(
domain=["GCV", "REML", "ML"],
range=["#e45756", "#4c78a8", "#72b7b2"]
)),
strokeDash=alt.StrokeDash("Method:N", scale=alt.Scale(
domain=["GCV", "REML", "ML"],
range=[[0], [0], [4, 4]]
)),
)
true_line = alt.Chart({"values": true_records}).mark_line(
color="black", strokeDash=[2, 2], strokeWidth=1, opacity=0.5
).encode(x="x:Q", y="y:Q")
(points + true_line + lines).properties(
width="container", height=300,
title="Smoothness selection: GCV vs REML vs ML"
)
The black dashed line is the true function. REML (blue) tends to produce the smoothest fit, GCV (red) the most flexible, and ML (teal, dashed) sits in between.
REML is the recommended default for nearly all applications. It is less prone to overfitting than GCV and produces more reliable confidence intervals. Switch to GCV only if you have a specific reason (e.g., reproducing a legacy analysis) or to ML when you need comparable likelihoods across models with different fixed effects.
Fixed smoothing parameters
Sometimes you want to bypass automatic smoothness selection and fix \lambda_j at known values. This is useful for:
- Reproducing results from another analysis where \lambda was determined externally.
- Sensitivity analysis: checking how the fit changes as you vary \lambda.
- Simulation studies: fitting at a known truth.
Pass a list of smoothing parameters to fit() in the same order as the smooth terms appear in the formula:
# Fix the smoothing parameter for s(x)
model_fixed = wk.GAM("y ~ s(x, k=20)")
model_fixed.fit(data, smoothing_params=[1.0])
print(f"Fixed lambda: {model_fixed.smoothing_params}")
print(f"EDF: {model_fixed.edf_total:.3f}")
Fixed lambda: [1.0]
EDF: 9.960
When smoothing_params is provided, the method argument is ignored — no optimization over \lambda takes place.
# Multiple smooths: one lambda per smooth
rng = np.random.default_rng(0)
n = 200
x1 = np.linspace(0, 1, n)
x2 = rng.uniform(0, 1, n)
y2 = np.sin(4 * x1) + 0.5 * x2**2 + rng.normal(0, 0.3, n)
data2 = {"x1": x1, "x2": x2, "y": y2}
model_two = wk.GAM("y ~ s(x1) + s(x2)")
model_two.fit(data2, smoothing_params=[1.0, 0.5])
print(f"Lambda for s(x1): {model_two.smoothing_params[0]:.2f}")
print(f"Lambda for s(x2): {model_two.smoothing_params[1]:.2f}")
Lambda for s(x1): 1.00
Lambda for s(x2): 0.50
Fixing smoothing parameters at inappropriate values can produce severely under- or overfit models. Use automatic selection (REML) unless you have a compelling reason to fix \lambda.
Variable selection with select=True
Standard penalized regression cannot shrink a smooth term to zero — the penalty null space (the space of functions not penalized, usually polynomials up to a given order) always survives. The double penalty approach adds an extra penalty on the null space, allowing the entire smooth to be penalized toward zero:
\text{penalty}_j = \lambda_j \boldsymbol\beta^\top S_j \boldsymbol\beta +
\lambda_j^{*} \boldsymbol\beta^\top S_j^{*} \boldsymbol\beta
where S_j^{*} penalizes the null space of S_j and \lambda_j^{*} is an additional smoothing parameter estimated alongside \lambda_j.
When select=True, Whittaker adds this extra penalty for every smooth term. If a predictor has no effect, both \lambda_j and \lambda_j^{*} can grow large enough to drive the smooth’s EDF to approximately zero, effectively removing it from the model.
# Variable selection example: x3 is pure noise
rng = np.random.default_rng(23)
n = 300
x1 = np.linspace(0, 2 * np.pi, n)
x2 = rng.uniform(0, 1, n)
x3 = rng.normal(0, 1, n) # noise predictor
y_sel = np.sin(x1) + x2**2 + rng.normal(0, 0.3, n)
data_sel = {"x1": x1, "x2": x2, "x3": x3, "y": y_sel}
# Without variable selection
model_nosel = wk.GAM("y ~ s(x1) + s(x2) + s(x3)")
model_nosel.fit(data_sel, method="REML")
# With variable selection
model_sel = wk.GAM("y ~ s(x1) + s(x2) + s(x3)")
model_sel.fit(data_sel, method="REML", select=True)
print("Without select=True:")
print(model_nosel.summary())
print("\nWith select=True:")
print(model_sel.summary())
Without select=True:
GAM fit summary
============================================================
Formula: y ~ s(x1) + s(x2) + s(x3)
Family: Gaussian(link='identity')
Inference: REML
Observations: 300
Coefficients: 28
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.3331 0.0165 20.128 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x1) 7.57 8 1924.348 < 1e-16
s(x2) 2.52 3 253.698 < 1e-16
s(x3) 1.90 2 2.021 0.3641
Total EDF: 12.99
Scale est: 0.082157
Deviance: 23.5797
Null dev: 201.2325
Dev. expl: 88.3%
GCV score: 0.085877
AIC: 114.62
BIC: 162.74
With select=True:
GAM fit summary
============================================================
Formula: y ~ s(x1) + s(x2) + s(x3)
Family: Gaussian(link='identity')
Inference: REML
Observations: 300
Coefficients: 28
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.3331 0.0166 20.086 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x1) 7.55 8 1914.498 < 1e-16
s(x2) 2.46 3 250.370 < 1e-16
s(x3) 0.00 1 0.000 0.9969
Total EDF: 11.01
Scale est: 0.082505
Deviance: 23.8434
Null dev: 201.2325
Dev. expl: 88.2%
GCV score: 0.085647
AIC: 113.90
BIC: 154.67
# Grouped bar chart comparing per-term EDF with and without select=True
term_labels = ["s(x1)", "s(x2)", "s(x3)"]
edf_nosel = model_nosel.edf
edf_sel = model_sel.edf
bar_records = []
for label, edf_std, edf_shrink in zip(term_labels, edf_nosel, edf_sel):
bar_records.append({"Term": label, "Model": "Standard", "EDF": float(edf_std)})
bar_records.append({"Term": label, "Model": "select=True", "EDF": float(edf_shrink)})
alt.Chart({"values": bar_records}).mark_bar().encode(
x=alt.X("Term:N", title="Smooth term", axis=alt.Axis(labelAngle=0)),
y=alt.Y("EDF:Q", title="Effective degrees of freedom"),
color=alt.Color("Model:N", scale=alt.Scale(
domain=["Standard", "select=True"],
range=["#4c78a8", "#e45756"]
)),
xOffset="Model:N",
).properties(width="container", height=300, title="Variable selection: per-term EDF")
With select=True, the EDF for s(x3) should be driven close to zero, confirming that Whittaker identified the noise predictor. The EDF values for s(x1) and s(x2) remain largely unchanged.
Use select=True when you have many candidate predictors and want the model to automatically determine which ones have a nonlinear effect. It is especially useful in exploratory analysis where you are unsure which predictors are relevant. For confirmatory analysis where the model structure is known, standard fitting (without select=True) is appropriate.
Observation weights
Observation weights enter the penalized WLS system by modifying the weight matrix. If you supply a weight vector w_i, the working weight matrix becomes \tilde{W} = \text{diag}(w_i) \cdot W, and the penalized WLS system is:
\bigl(X^\top \tilde{W} X + \textstyle\sum_j \lambda_j S_j\bigr)\,\hat\beta =
X^\top \tilde{W} z
This allows you to:
- Downweight outliers: set w_i < 1 for observations you suspect are contaminated.
- Account for known precision: when observations have different known variances, set w_i = 1 / \sigma_i^2.
- Handle aggregated data: when each row represents n_i observations, set w_i = n_i.
# Weights example: downweight outlier observations
rng = np.random.default_rng(23)
n = 100
x = np.linspace(0, 2 * np.pi, n)
y_w = np.sin(x) + rng.normal(0, 0.3, n)
# Add some outliers
outlier_idx = [10, 30, 50, 70, 90]
y_w[outlier_idx] += rng.choice([-3, 3], size=len(outlier_idx))
# Create weights: 1.0 for normal observations, 0.1 for outliers
w = np.ones(n)
w[outlier_idx] = 0.1
data_w = {"x": x, "y": y_w}
# Fit without weights
model_unweighted = wk.GAM("y ~ s(x)")
model_unweighted.fit(data_w, method="REML")
# Fit with weights
model_weighted = wk.GAM("y ~ s(x)")
model_weighted.fit(data_w, method="REML", weights=w)
print(f"Unweighted scale: {model_unweighted.scale:.4f}")
print(f"Weighted scale: {model_weighted.scale:.4f}")
print(f"Unweighted EDF: {model_unweighted.edf_total:.3f}")
print(f"Weighted EDF: {model_weighted.edf_total:.3f}")
Unweighted scale: 0.5376
Weighted scale: 0.1516
Unweighted EDF: 5.587
Weighted EDF: 7.000
# Scatter with outliers highlighted plus weighted vs unweighted fitted curves
x_fit = np.linspace(0, 2 * np.pi, 200)
preds_unw = model_unweighted.predict({"x": x_fit})
preds_w = model_weighted.predict({"x": x_fit})
is_outlier = np.zeros(n, dtype=bool)
is_outlier[outlier_idx] = True
scatter_records = [
{"x": float(x[i]), "y": float(y_w[i]), "Type": "Outlier" if is_outlier[i] else "Normal"}
for i in range(n)
]
fit_records = []
for i in range(len(x_fit)):
fit_records.append({"x": float(x_fit[i]), "fit": float(preds_unw.values[i]), "Model": "Unweighted"})
fit_records.append({"x": float(x_fit[i]), "fit": float(preds_w.values[i]), "Model": "Weighted"})
points = alt.Chart({"values": scatter_records}).mark_circle().encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="y"),
color=alt.Color("Type:N", scale=alt.Scale(
domain=["Normal", "Outlier"],
range=["gray", "#e45756"]
)),
size=alt.Size("Type:N", scale=alt.Scale(
domain=["Normal", "Outlier"],
range=[30, 80]
), legend=None),
opacity=alt.condition(
alt.datum.Type == "Outlier",
alt.value(0.8),
alt.value(0.3),
),
)
lines = alt.Chart({"values": fit_records}).mark_line(strokeWidth=2).encode(
x="x:Q",
y="fit:Q",
color=alt.Color("Model:N", scale=alt.Scale(
domain=["Unweighted", "Weighted"],
range=["#e45756", "#4c78a8"]
)),
strokeDash=alt.StrokeDash("Model:N", scale=alt.Scale(
domain=["Unweighted", "Weighted"],
range=[[4, 4], [0]]
)),
)
(points + lines).properties(
width="container", height=300,
title="Effect of observation weights on outlier robustness"
)
The weighted fit should produce a lower scale estimate and fewer effective degrees of freedom because it is not distorted by the outliers.
Convergence
What convergence means
The P-IRLS algorithm converges when the relative change in penalized deviance between successive iterations falls below the tolerance threshold (10^{-7} by default). Specifically, the algorithm stops when:
\frac{|D^{(k)} - D^{(k-1)}|}{|D^{(k-1)}| + 0.1} < \text{tol}
The 0.1 in the denominator prevents numerical issues when the deviance is near zero.
Checking the fit
After fitting, inspect the model to verify the fit is reasonable:
model_check = wk.GAM("y ~ s(x)")
model_check.fit({"x": x, "y": np.sin(x) + rng.normal(0, 0.3, n)}, method="REML")
print(f"EDF total: {model_check.edf_total:.1f}")
print(f"Scale: {model_check.scale:.4f}")
print(f"Deviance: {model_check.deviance:.2f}")
EDF total: 7.6
Scale: 0.0930
Deviance: 8.59
For Gaussian models with an identity link, the P-IRLS algorithm converges in a single step. For non-Gaussian models, convergence typically occurs within 5–15 iterations.
- Check the data: fitting problems are often caused by complete separation (in logistic regression), extreme outliers, or a poor choice of family.
- Reduce
k: an overly flexible model can be unstable. Reducing the basis dimension can help.
- Try a different method: GCV and REML can behave differently on ill-conditioned problems.
- Run
model.check(): the basis dimension adequacy test can reveal underfitting (see Diagnostics).
Effective degrees of freedom (EDF)
The effective degrees of freedom measures the complexity of each smooth term. It is defined as:
\text{EDF}_j = \text{tr}(\mathbf{A}_j)
where \mathbf{A}_j is the block of the hat matrix corresponding to the j-th smooth. The total model EDF is \text{EDF}_{\text{total}} = \sum_j \text{EDF}_j + p, where p is the number of parametric (unpenalized) coefficients including the intercept.
Interpreting EDF:
| \approx 1 |
The smooth is approximately linear |
| \approx 2 |
The smooth is approximately quadratic |
| 3–5 |
Moderate curvature |
| > 5 |
High complexity: consider checking the data for artifacts |
| \approx k - 1 |
The smooth is using nearly all available basis functions: increase k |
# EDF interpretation
model_edf = wk.GAM("y ~ s(x, k=20)")
model_edf.fit({"x": x, "y": np.sin(x) + rng.normal(0, 0.3, n)}, method="REML")
print(f"Total EDF: {model_edf.edf_total:.3f}")
print(f"Scale estimate: {model_edf.scale:.4f}")
Total EDF: 10.293
Scale estimate: 0.0745
EDF is not the same as the basis dimension k. The parameter k sets an upper bound on complexity, while the smoothing parameter \lambda (selected by REML or GCV) determines how much of that capacity is actually used. Setting k=20 does not mean the smooth uses 20 degrees of freedom — it may use only 4 or 5 if the data support a simpler shape. See the smooth terms page for guidance on choosing k.
Scale estimation
The scale parameter \hat\phi measures the residual variability after accounting for the smooth effects. For the Gaussian family, \hat\phi = \hat\sigma^2 is the residual variance. It is estimated by:
\hat\phi = \frac{D(y, \hat\mu)}{n - \text{EDF}_{\text{total}}}
where D is the deviance and \text{EDF}_{\text{total}} is the total effective degrees of freedom. This is the analogue of s^2 = \text{RSS} / (n - p) in ordinary least squares, replacing p with the (fractional) EDF.
For families with a known scale (Poisson, binomial), \phi is fixed at 1 and is not estimated.
# Scale estimation
print(f"Scale (phi-hat): {model_edf.scale:.4f}")
print(f"Deviance: {model_edf.deviance:.4f}")
print(f"GCV score: {model_edf.gcv_score:.4f}")
Scale (phi-hat): 0.0745
Deviance: 6.6825
GCV score: 0.0830
The scale estimate directly affects confidence intervals and p-values: standard errors are proportional to \sqrt{\hat\phi}, so an overestimated scale produces wider confidence bands.
Practical guidance
When to use REML vs GCV
| Default, general-purpose fitting |
REML |
| Reproducing a legacy analysis that used GCV |
GCV |
| Comparing models with different fixed effects via LRT |
ML |
| Very large n (> 50,000) |
GCV may be faster (both give similar results) |
| Sparse data or small n |
REML (GCV tends to undersmooth) |
When to increase k
The basis dimension k should be large enough that the smooth can capture the true function shape. After fitting, run model.check() (see Model diagnostics) to test whether k is adequate:
If the k-index is below 1 with a significant p-value, double k and re-fit:
model = wk.GAM("y ~ s(x, k=20)")
model.fit(data, method="REML")
model.check()
- Start with the defaults:
method="REML", k=10.
- Fit the model and run
model.check().
- If any smooth fails the k-index test, increase
k for that term and re-fit.
- Review the summary: are the EDF values plausible for the data?
- If you suspect irrelevant predictors, re-fit with
select=True.
- Generate predictions and diagnostic plots.
Summary of fit attributes
After calling model.fit(), the following attributes are available:
model.smoothing_params |
list[float] |
Estimated (or fixed) \lambda_j values |
model.edf_total |
float |
Total effective degrees of freedom |
model.scale |
float |
Estimated scale parameter \hat\phi |
model.deviance |
float |
Model deviance |
model.gcv_score |
float |
GCV score (computed regardless of method) |