Generalized Additive Models › Implementation in NumPy and Statsmodels
Machine Learning from ScratchAll chapters Generalized Additive Models › Implementation in NumPy and StatsmodelsSection view · Sections in this chapter7 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
We will carry out the fit in two stages. First we implement the GAM from primitives using NumPy, building the design matrix with the B-spline transformation in patsy and constructing the penalty matrix by hand. Then we compare against statsmodels.gam.api.GLMGam, which provides the same kind of model with additional diagnostics. The choice of statsmodels as reference is deliberate: it exposes both the fitted coefficients and the effective degrees of freedom for each smooth, both of which we want for interpretation.
Building the smooth design matrices
The fit proceeds by constructing a B-spline basis for each predictor, stacking them into a single design matrix (with an intercept), and solving the penalized normal equations. We start with 12 basis functions per smooth. Each training-sample centering constraint leaves 11 free smooth coefficients, so the two-feature design has 23 columns including the intercept. The transformation must be applied to the penalty as well as the basis. The bounds are known from our synthetic design before seeing the data; interior knots are learned from each training fold. These example helpers deliberately use that fixed domain. They assume finite complete inputs with enough distinct predictor values to construct the basis and an identifiable penalized system. They are a teaching implementation, rather than a general-purpose missing-data or degenerate-design handler.
PythonIn [13]Expand
def b_spline_basis(x, df=12, degree=3):
"""Fit a full basis on training values within the known synthetic domain."""
design = patsy.dmatrix(
f"0 + bs(x, df={df}, degree={degree}, include_intercept=True, "
"lower_bound=-3, upper_bound=3)",
{"x": x},
)
return np.asarray(design), design.design_info
def difference_penalty(k, order=2):
"""Quadratic penalty on differences of the full coefficient sequence."""
D = np.diff(np.eye(k), n=order, axis=0)
return D.T @ D
def fit_additive_gam(x_list, y, df=12, lam=1.0):
"""Fit centered difference-penalized splines with training basis metadata."""
y = np.asarray(y, dtype=float)
if lam < 0 or not np.isfinite(lam) or y.ndim != 1 or len(y) == 0:
raise ValueError(
"Use nonempty one-dimensional responses and finite lam >= 0"
)
if not x_list or any(np.asarray(x).shape != y.shape for x in x_list):
raise ValueError("Each predictor must match the response shape")
blocks = [np.ones((len(y), 1))]
penalty_blocks = [np.zeros((1, 1))]
infos, transforms = [], []
for x in x_list:
raw, info = b_spline_basis(x, df=df)
means = raw.mean(axis=0)
# Write the last coefficient in terms of the others so means @ Q = 0.
Q = np.vstack([np.eye(df - 1), -means[:-1] / means[-1]])
blocks.append(raw @ Q)
penalty_blocks.append(lam * Q.T @ difference_penalty(df) @ Q)
infos.append(info)
transforms.append(Q)
B = np.hstack(blocks)
P = np.zeros((B.shape[1], B.shape[1]))
offset = 0
for block in penalty_blocks:
k = len(block)
P[offset : offset + k, offset : offset + k] = block
offset += k
beta_hat = np.linalg.solve(B.T @ B + P, B.T @ y)
return B, P, beta_hat, infos, transforms
def prediction_design(x_list, infos, transforms):
"""Reuse training knots and training centering on new rows."""
blocks = [np.ones((len(x_list[0]), 1))]
for x, info, Q in zip(x_list, infos, transforms):
raw = np.asarray(patsy.build_design_matrices([info], {"x": x})[0])
blocks.append(raw @ Q)
return np.hstack(blocks)Now we fit to the synthetic data and extract the fitted components. Note that each B_j block is multiplied by its own set of coefficients; to obtain the fitted smooth at training points we slice the relevant block.
PythonIn [14]Expand
B, P, beta_hat, basis_infos, center_transforms = fit_additive_gam(
[x1, x2], y, df=12, lam=1.0
)
b0 = beta_hat[0]
n_basis = 11 # 12 full basis functions minus one centering constraint
b1 = beta_hat[1 : 1 + n_basis]
b2 = beta_hat[1 + n_basis : 1 + 2 * n_basis]
f1_hat_train = B[:, 1 : 1 + n_basis] @ b1
f2_hat_train = B[:, 1 + n_basis : 1 + 2 * n_basis] @ b2
y_hat = b0 + f1_hat_train + f2_hat_train
rmse = np.sqrt(np.mean((y - y_hat) ** 2))We can now compare the fitted smooths against the true curves. Because each smooth is only identified up to a constant (the intercept absorbs it), we center both the fitted and true curves to have zero mean before plotting.
PythonIn [15]Expand
f1_hat_centered = f1_hat_train - f1_hat_train.mean()
f2_hat_centered = f2_hat_train - f2_hat_train.mean()
f1_true_centered = f1_true - f1_true.mean()
f2_true_centered = f2_true - f2_true.mean()
order1 = np.argsort(x1)
order2 = np.argsort(x2)VisualizationsOut [16]Expand
The fitted smooths can be compared directly with the known generating curves. Deviations reflect sampling noise, the finite basis, and the smoothing penalty; they are not a guarantee of any particular boundary behavior. The relevant residual statistic is computed below.
PythonIn [17]Expand
noise_sd = 0.4ConsoleOut [18]Expand
RMSE on training data: 0.366 Noise standard deviation: 0.400
Training RMSE describes the fitted residuals, but cannot establish whether the model overfits or generalizes well. We compare fitted components with the known synthetic truth here and evaluate a selected model on independent test observations below.
Choosing the smoothing parameter
The choice of matters a great deal. We can visualize the trade-off by refitting at several values and plotting the smooth for each. We compute the fits first, then render the figure.
PythonIn [19]Expand
lambdas = [0.01, 1.0, 1000.0]
grid1 = np.linspace(x1.min(), x1.max(), 300)
smooth_paths = {}
for lam in lambdas:
_, _, beta_lam, infos_lam, Q_lam = fit_additive_gam([x1, x2], y, lam=lam)
grid_design = prediction_design(
[grid1, np.zeros_like(grid1)], infos_lam, Q_lam
)
smooth_paths[lam] = (
grid_design[:, 1 : 1 + n_basis] @ beta_lam[1 : 1 + n_basis]
)VisualizationOut [20]Expand
The plotted paths show how the fixed basis responds to different penalties. The largest penalty removes much of the sinusoidal structure. Whether the smaller penalty is too flexible requires validation rather than a judgment from training RMSE or appearance alone.
Comparing with statsmodels
statsmodels.gam.api.GLMGam implements the same general penalized-GLM framework, but its BSplines smoother uses an integrated squared-second-derivative penalty, rather than our coefficient-difference penalty. Equal numeric values of alpha and do not specify equal smoothing. We fit a centered statsmodels model with an explicit intercept, then independently solve its own Gaussian penalized equations in NumPy. This comparison tests the fitting mechanism on identical inputs; it does not claim our earlier P-spline is the same model.
PythonIn [21]Expand
from statsmodels.gam.api import BSplines, GLMGam
n_splines = 12
smoothers = BSplines(
np.column_stack([x1, x2]),
df=[n_splines, n_splines],
degree=[3, 3],
constraints="center",
)
gam_sm = GLMGam(y, exog=np.ones((n, 1)), smoother=smoothers, alpha=[0.1, 0.1])
gam_res = gam_sm.fit()
sm_fitted = gam_res.predict()
sm_rmse = np.sqrt(np.mean((y - sm_fitted) ** 2))
# statsmodels 0.14.6 PIRLS doubles the smooth penalty in its WLS solve.
B_reference = gam_sm.exog
P_reference = 2 * gam_sm.penal.penalty_matrix()
beta_reference = np.linalg.solve(
B_reference.T @ B_reference + P_reference, B_reference.T @ y
)
np.testing.assert_allclose(beta_reference, gam_res.params, atol=1e-9, rtol=1e-9)
reference_difference = np.max(np.abs(B_reference @ beta_reference - sm_fitted))We now overlay the NumPy solution of this reference system with the statsmodels result. Both curves use the fitted statsmodels basis and its training centering transform, evaluated on the same grid.
PythonIn [22]Expand
reference_grid = smoothers.smoothers[0].transform(grid1)
reference_width = smoothers.smoothers[0].basis.shape[1]
smooth_scratch = reference_grid @ beta_reference[1 : 1 + reference_width]
smooth_sm = reference_grid @ gam_res.params[1 : 1 + reference_width]VisualizationOut [23]Expand
The assertion above checks agreement of all reference coefficients. The numerical maximum prediction difference is printed below. Agreement verifies the Gaussian solve for this basis and penalty; it does not verify non-Gaussian inference or the choice of smoothing strength.
ConsoleOut [24]Expand
From-scratch training RMSE: 0.3659 statsmodels training RMSE: 0.3653 Reference maximum prediction difference: 6.661e-15
The earlier P-spline RMSE and the statsmodels RMSE need not match: they use different bases, penalties, and smoothing parameters. The reference assertion compares the library to an independent solve of its own system, which is the controlled comparison.
Key parameters
- Number of basis functions per smooth (
df): controls maximum flexibility. A full cubic basis withdf = 4spans cubic polynomials; increasing the dimension permits more detail. The penalty cannot recover shapes excluded by a basis that is too small, so check sensitivity to basis size as well as smoothing. - Smoothing parameter : the actual knob that trades off flexibility against overfitting. Choose it by cross-validation, as we will see in the next section.
- Penalty order: 2 for the standard second-difference penalty, 1 for a roughness penalty on first differences (encourages flatness rather than linearity), 3 for third coefficient differences. The order describes the coefficient sequence and is distinct from the spline degree.
- Link function: identity for Gaussian, logit for binary, log for counts, inverse for gamma.
Key methods
GLMGam.fit()estimates coefficients and standard errors in a single call (using PIRLS by default in statsmodels 0.14.6).gam_res.predict(exog=linear_values, exog_smooth=smooth_values)evaluates the fitted reference model at new inputs within its supported knot domain.smoothers.transform(X)reuses the fitted spline basis for within its knot bounds, which is what we need to evaluate the fitted smooth on a grid.gam_res.edfexposes effective degrees of freedom by coefficient; block sums describe components.gam_res.summary()reports the fitted model summary.
Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!