Generalized Additive Models

A regression can miss a pattern even when every predictor has a clear relationship with the response: a straight-line term cannot follow a sine-shaped association or a U-shaped one. Spline Regression gave us curves for an individual predictor. Combining several curves raises two further questions: how do we keep their contributions identifiable, and how do we choose their flexibility without judging them on the same observations used for fitting?

Generalized additive models (GAMs) fill that gap. They keep the additive structure of a linear model, replace each linear term with a smooth function, and let us pair that structure with the link functions familiar from generalized linear models. Formally, a GAM has the form

where:

• $g$: the link function connecting the conditional mean to the additive predictor
• $\mathbb{E}[y \mid \mathbf{x}]$: the conditional mean of the response given the predictors
• $\beta_0$: the intercept
• $f_j$: a smooth function of predictor $x_j$
• $x_j$: the $j$-th predictor value, with $j=1,\ldots,p$ and $p$ predictors

Two questions guide the implementation:

• how to estimate the amount of flexibility (the smoothing parameter $\lambda$) without leaking information from validation into training,
• how to interpret a fitted additive model, including what it cannot capture (interactions).

This chapter develops those ideas from first principles. We will build a GAM by hand using spline bases and a difference penalty, show how the fit reduces to a closed-form ridge-style solve, extend the framework to non-Gaussian responses via a link function, and demonstrate the whole pipeline on a synthetic dataset where the true component functions are known. Implementation is done in NumPy for the mechanism and compared against statsmodels' GLMGam for a reference. We will close with practical guidance on when GAMs make sense and what to watch out for.

The additive restriction limits the number of functions that need estimating. For example, a Cartesian grid with $m$ positions per predictor contains $m^p$ grid points, while $p$ separate one-dimensional grids contain $pm$ positions. A GAM imposes that lower-dimensional additive structure rather than estimating an unrestricted surface. The comparison describes the representation, not a universal sample-size bound. The additive restriction represents a sum of $p$ one-dimensional functions instead of an unrestricted surface. Their coefficients still have to be fitted jointly, and computation grows with the total basis dimension.

Second, a GAM remains interpretable. Because the components enter additively, we can plot each $f_j$ as a curve and read off how the predicted value (or the transformed mean) changes as a function of $x_j$, holding the others fixed. These are direct plots of model terms on the additive-predictor scale. With a nonlinear link, they differ from response-scale partial-dependence curves, which average predictions over other inputs; see Part XI: Model Interpretation. That makes GAMs attractive in domains where the answer "the model predicts this" is not enough and stakeholders want to know why.

Third, a GAM extends the machinery of generalized linear models to non-Gaussian responses without abandoning interpretability. We can use the same spline bases we built for continuous responses inside a logistic, Poisson, or gamma regression by adding a link function. This gives us nonlinear logistic and count regression through penalized likelihood. As a consequence, when a linear GLM is clearly misspecified, such as a binary outcome whose logit is obviously curved in a predictor, a GAM is a natural upgrade.

Finally, the smoothing penalty plays the role of a continuous regularizer. Rather than choosing the number of knots arbitrarily, we let the data guide how wiggly each curve should be via cross-validation or marginal likelihood. This is the same idea ridge regularization will formalize in the next chapter Ridge Regularization (L2 Regularization); in fact, a GAM with a difference penalty is a generalized quadratic-penalty regression on an expanded feature space; the penalty matrix differs from ordinary ridge.

GAMs are not free. Their main limitation is that the additive structure itself is a strong assumption. If the true relationship involves an interaction between two predictors, say the effect of $x_1$ flips sign depending on $x_2$, an additive model cannot represent it. We can add explicit tensor-product smooths $f_{12}(x_1, x_2)$ to handle interactions, but each interaction term costs us both flexibility and sample size, and the resulting components are much harder to plot and interpret. In high dimensions the number of potential pairwise interactions explodes.

Second, GAMs depend on several tuning choices: the basis dimension and knot arrangement, the penalty order, and the smoothing parameter $\lambda$ per component. Choosing those badly produces either an overfit curve that wiggles through noise or an over-smoothed curve that misses real structure. We can automate the choice with cross-validation or generalized cross-validation, but their estimates have sampling variability. Neither criterion guarantees that the selected model will generalize well; strong overlap between smooth terms can also complicate interpretation.

Third, extrapolation is dangerous in a way that is easy to forget. Out-of-domain behavior depends on the basis implementation: it may refuse prediction or impose an extrapolation rule. In particular, the Patsy and statsmodels B-spline transforms used here reject points beyond their boundary knots. A GAM should not be used to predict far outside the training support of any predictor.

A spline expansion also widens the design matrix relative to a linear model of the original predictors. Dense cross-products and coefficient solves then grow with the total basis dimension; specialized implementations can exploit sparsity or low-rank structure. We give the dense costs in Practical Implications.

Component plots show the fitted terms on the additive-predictor scale. They need care when predictors are dependent, because varying one predictor while holding the others fixed can leave the region supported by the data. Overlapping smooth-term spaces can also make individual components unstable.

A GAM is best read as the answer to a design question: how do we let every predictor bend, yet keep the model as tractable and interpretable as a linear one? The formula answers that question in three interlocking layers, and each layer exists to solve a specific problem.

1. An additive predictor built from smooths replaces the single straight line of linear regression with a sum of curves. This buys flexibility while preserving the "one term per predictor" structure that makes the model readable.
2. A link function connects that additive predictor to the mean of the response, so that the same structure can serve a Gaussian, binary, or count outcome.
3. A smoothness penalty discourages rough fitted curves. Identifiability constraints are also needed, and the amount of smoothing still requires selection.

It helps to see each problem before its formula. Picture two failure modes we must avoid. If we let the smooths be arbitrarily wiggly, they will thread through noise; if we over-restrict them, we are back to a straight line and lose the whole point. The penalty is the dial between these extremes. And if we ignore the scale of the response, the same smooths cannot respect a probability or a count. The link is what fixes the scale. With both ideas in place, the additive predictor is the piece that ties the predictors together. We build the layers in that order of need.

Once the basis is fixed, the additive predictor is linear in its unknown coefficients. For a Gaussian response with an identity link, fitting reduces to a quadratic penalized least-squares problem. Other response families generally require iterative likelihood fitting. The following steps distinguish the predictor, link, and fitting objective.

The additive predictor sums an intercept and one component per predictor for observation $i$:

where:

• $i = 1, \dots, n$: index over observations,
• $\eta_i$: the value of the additive predictor for observation $i$,
• $\beta_0$: an intercept shared by all observations,
• $p$: the number of predictors,
• $x_{ij}$: the value of predictor $j$ for observation $i$,
• $f_j$: a smooth function of predictor $j$, not necessarily linear.

If $f_j(x)=\beta_j x$ for every predictor, this becomes the predictor from Multiple Linear Regression. A GAM relaxes that restriction one term at a time: each $f_j$ may curve, but the terms still add, so we can still read the model one predictor at a time.

Why insist on additivity instead of fitting one flexible surface $f(x_1, \dots, x_p)$? Additivity restricts the model to a smaller family of surfaces. The additive decomposition restricts the unknown surface to a sum of one-dimensional curves, whose coefficients are fitted jointly. That trade between flexibility within each predictor and additivity across predictors is the structural bet of a GAM, and it is exactly what we will exploit when we expand each $f_j$ in a basis.

There is a second gap to close. The smooths naturally produce a real number, but the response need not live on the real line: a binary outcome is zero or one, while a finite logit predictor gives a fitted mean in $(0,1)$. A count mean need not be an integer; a log link keeps it positive. The link function bridges that gap:

Here $g$ is a monotone, differentiable link function. Its role is to connect the permissible mean to the additive predictor. The log and logit links map their mean domains onto the real line; an inverse link has a restricted predictor domain and requires corresponding constraints. Common choices are:

• identity link $g(\mu) = \mu$ for Gaussian regression (this is the case we will spend most time on),
• logit link $g(\mu) = \log(\mu / (1 - \mu))$ for binary responses,
• log link $g(\mu) = \log(\mu)$ for counts,
• inverse link $g(\mu) = 1/\mu$ for gamma or inverse-Gaussian responses.

For a Gaussian GAM with identity link, the model predicts the conditional mean directly. For a binary GAM with logit link, the model predicts the log-odds of success, and each $f_j$ describes how the log-odds change as a smooth function of $x_j$. For a count GAM with log link, each $f_j$ describes how the log expected count changes.

Notice the clean division of labor this creates: the smooths describe shape on the scale of $\eta$, and the link describes scale. That separation is what allows one spline-and-penalty mechanism to serve radically different response types without being redesigned for each.

We are foreshadowing here because the later chapter Generalized Linear Models develops the link-function framework in full. For this chapter the important point is that the same additive machinery, namely smooths, penalties, and penalized least-squares solves, extends to suitable GLM families and links through penalized likelihood, with additional iterative fitting (iterative reweighting, which we will sketch but not build from scratch).

To compute the fit using finitely many coefficients, we represent each smooth $f_j$ with a spline basis:

where:

• $K_j$ is the number of basis functions we allocate to predictor $j$,
• $b_{jk}(x)$ is the $k$-th basis function, a known deterministic function of $x$,
• $\beta_{jk}$ are unknown coefficients we will estimate.

We use B-spline basis functions: low-degree piecewise polynomials with local support. Natural cubic spline bases are another option, but their individual basis functions need not have local support. Taken together, they form a flexible family of curves. The dimension $K_j$ controls the maximum possible wiggliness: with more basis functions we can represent more complex shapes.

A mental picture makes the next step natural. Each $b_{jk}$ is a localized bump, and the fitted smooth is a weighted sum of those bumps, with the weights given by the coefficients $\beta_{jk}$. Increasing $K_j$ adds more bumps and therefore more potential detail, but it also adds more coefficients that can conspire to fit noise.

Substituting the basis expansion into the additive predictor gives

This is now linear in the coefficients $\beta_{jk}$. If we stack all coefficients into a single vector $\boldsymbol\beta$ and expand each predictor into its basis, the model becomes $\boldsymbol\eta = \mathbf{B}\boldsymbol\beta$ for an augmented design matrix $\mathbf{B}$ whose columns are the basis functions evaluated at each observation.

We now reach the second problem hinted at earlier. With enough independent basis columns, a spline can interpolate a finite sample at distinct inputs. Fitting noise that closely can hurt predictions on new observations. Too small a basis can instead exclude real structure. Rather than solve the flexibility problem by choosing the basis size, we keep the basis generous and control curvature directly. We penalize it with a quadratic penalty on the (second) differences of neighboring coefficients:

This is the coefficient-difference penalty used in P-splines. Classical P-splines pair it with regularly spaced B-splines; our later Patsy example uses quantile knots with the same coefficient-difference regularizer, so its smoothing is tied to that parameterization. Each difference measures curvature in coefficient-index space. For smooth coefficient sequences on a fine, uniform knot grid with compatible boundary treatment, appropriate knot-spacing scaling makes this a discrete approximation to a derivative penalty. It is not generally identical to integrating the squared second derivative of the curve. Squaring and summing accumulates a measure of total curvature, so the penalty grows when the fitted coefficients oscillate from one basis function to the next. The parameter $\lambda_j \ge 0$ controls how strongly that penalty acts:

• $\lambda_j = 0$ gives an unpenalized spline, which can be as wiggly as the basis allows;
• $\lambda_j \to \infty$ forces the coefficient sequence to be affine in its index; the resulting curve depends on the knot arrangement;
• moderate $\lambda_j$ interpolates between the two, giving the data room to bend but discouraging noise.

The penalty is deliberately a local statement about neighboring coefficients. It does not insist that the whole curve be flat; it only makes sharp reversals from one coefficient to the next expensive. That is why a moderate $\lambda_j$ still permits broad bends (like a sinusoid) while a huge one suppresses the penalized directions. An affine coefficient sequence need not produce a globally straight line with quantile knots and repeated boundary knots.

In matrix form the penalty for smooth $j$ can be written as

where:

• $\boldsymbol\beta_j = (\beta_{j1}, \dots, \beta_{jK_j})^\top$ collects the smooth coefficients,
• $\mathbf{D}$ is the $(K_j - 2) \times K_j$ second-difference matrix, each row of the form $(0, \dots, 0, 1, -2, 1, 0, \dots, 0)$,
• $\mathbf{S}_j = \mathbf{D}^\top \mathbf{D}$ is a $K_j \times K_j$ symmetric positive semidefinite matrix that encodes the penalty.

Putting the pieces together, the Gaussian GAM is fit by minimizing

Read the two terms as a negotiation. The first term asks the curve to explain the data; the second asks it to stay smooth. The smoothing parameters $\lambda_j$ set the exchange rate between those demands.

In matrix notation, augmenting $\mathbf{B}$ with a leading column of ones for the intercept and stacking the $p$ coefficient blocks,

where $\mathbf{S}$ is the block-diagonal matrix with zero in the intercept block and $\lambda_j \mathbf{S}_j$ in the block for smooth $j$. After imposing centering, the following matrix formulas use the free coefficients and the correspondingly transformed design and penalty; the symbol $\boldsymbol\beta$ now denotes those identifiable coordinates. Setting the gradient of $\mathcal{L}$ with respect to $\boldsymbol\beta$ equal to zero and rearranging (exactly as we did for OLS in Ordinary Least Squares (OLS)) gives the closed-form solution

This is identical in form to ridge regression, with $\mathbf{S}$ playing the role of a feature-specific ridge penalty. The inverse exists only when the design identifies every direction left unpenalized: $\ker(\mathbf{B}) \cap \ker(\mathbf{S}) = \{0\}$. Indeed, for any coefficient direction $\mathbf{v}$, $\mathbf{v}^\top(\mathbf{B}^\top\mathbf{B}+\mathbf{S})\mathbf{v}=\|\mathbf{B}\mathbf{v}\|^2+\mathbf{v}^\top\mathbf{S}\mathbf{v}$. A positive smoothing parameter alone cannot remove an unpenalized constant ambiguity. The inverse formula therefore applies after identifiability constraints; in code we solve the system rather than form an inverse. The fitted curve for predictor $j$ is $\hat f_j(x) = \sum_k \hat\beta_{jk} \, b_{jk}(x)$, and the full fitted predictor is $\hat\eta_i = \beta_0 + \sum_j \hat f_j(x_{ij})$.

For many family and link choices, such as Bernoulli responses with a logit link or Poisson counts with a log link, the penalized likelihood is non-quadratic. These models are commonly fitted by iterative reweighted least squares (IRLS), which replaces one least-squares solve with a sequence of penalized weighted least-squares solves. The likelihood and link jointly determine the objective; a nonlinear inverse link alone does not determine whether it is quadratic. A common Fisher-scoring update is

where:

• $t$: the IRLS iteration index
• $\mathbf{W}^{(t)}$: the diagonal matrix of GLM working weights at iteration $t$, with entries depending on the variance function and link derivative
• $\boldsymbol{z}^{(t)}$: the working response at iteration $t$, formed from the current linear predictor and the link
• $\mathbf{S}$: the same block-diagonal penalty matrix as in the Gaussian case

For independent GLM observations, a common working response is $z_i=\eta_i+(y_i-\mu_i)g'(\mu_i)$, where $g'(\mu_i)$ is the derivative of the link at the current fitted mean. The corresponding working weight is $w_i=1/[\operatorname{Var}(y_i\mid\mathbf{x}_i)\{g'(\mu_i)\}^2]$. The response linearizes the mean relationship locally, while the weights account for conditional variance and the link's slope. These expressions require valid interior means, finite positive variances, and nonzero link derivatives. A common dispersion factor can instead be absorbed into the penalty convention. The resulting Fisher-scoring iteration may need damping and convergence checks; the displayed update is not a guarantee of convergence.

We can retain the same chosen coefficient penalty when changing the response family and link, while changing the working weights and response. Each IRLS step couples the additive coefficients through $\mathbf{W}^{(t)}$, and a well-defined solution requires the smooth coefficients to satisfy the identifiability conditions described next.

The difference-penalized Gaussian model has several properties worth stating:

• The second-difference matrix has a two-dimensional null space: constant and affine coefficient sequences have zero second differences. With arbitrary knot spacing this statement concerns coefficients, not a guarantee of exact linearity in the predictor.
• A constant can be transferred between a smooth and the intercept without changing predictions. We impose $\sum_i f_j(x_{ij})=0$ on the training observations for each component. This is a constraint during fitting, rather than merely subtracting a mean from a plotted curve afterward.
• With basis and smoothing parameters fixed, the Gaussian fit is linear in the responses. Its total effective degrees of freedom are $\operatorname{tr}(\mathbf{H})$, where $\mathbf{H}=\mathbf{B}(\mathbf{B}^{\top}\mathbf{B}+\mathbf{S})^{-1}\mathbf{B}^{\top}$ for the identifiable design. This counts the intercept too; a centered component does not generally have degrees of freedom ranging from $K_j$ to $2$.

The matrix $\mathbf{H}$ maps observed responses to fitted values: $\hat{\mathbf{y}}=\mathbf{H}\mathbf{y}$. With the basis and smoothing parameters fixed, its diagonal entry $H_{ii}$ measures how sensitive fitted value $i$ is to observed response $i$. The trace, $\operatorname{tr}(\mathbf{H})$, sums these diagonal entries to give an effective parameter count. For an unpenalized, full-column-rank least-squares fit, it equals the number of fitted coefficients. Smoothing shrinks the penalized directions, so the effective count can be fractional.

The next figure isolates one unconstrained, full-rank, 12-column spline basis with no separate intercept. Its total effective degrees of freedom approach the two unpenalized coefficient directions as smoothing grows. This is a separate demonstration from the centered two-feature model, where the intercept and each constraint affect the count.

For independent errors of variance $\sigma^2$ and fixed smoothing parameters, write $\mathbf{A}=\mathbf{B}^{\top}\mathbf{B}+\mathbf{S}$. Then $\operatorname{Cov}(\hat{\boldsymbol\beta})=\sigma^2\mathbf{A}^{-1}\mathbf{B}^{\top}\mathbf{B}\mathbf{A}^{-1}$. Here $\operatorname{Cov}$ collects the sampling covariances of the estimated free coefficients. Its diagonal entries are their variances; their square roots give standard errors conditional on the fixed design and smoothing choices. The factor $\mathbf{A}^{-1}\mathbf{B}^{\top}$ maps response noise to coefficient noise. These conditional sampling variances do not include smoothing-selection uncertainty or remove smoothing bias.

Before we fit anything, we need to see the shape of the problem. Suppose we generate a dataset whose two predictors each have a known non-linear effect on a Gaussian response:

with $x_{i1}, x_{i2}$ drawn uniformly from $[-3, 3]$ and $\sigma = 0.4$. The first component is a sinusoid; the second is a convex quadratic (upward-opening parabola). Both are clearly not linear, and the additive structure holds exactly because the two components enter as a sum.

Each basis curve contributes only to a short segment of the $x$ range, and they overlap smoothly so that a weighted sum produces a smooth curve within the knot domain. Increasing the number of basis functions makes the space of possible curves richer. The job of $\lambda$ is to determine how aggressively we smooth the resulting weighted sum.

The formula above is compact, but the algebra inside the penalized solve deserves to be seen once by hand. We now work through a tiny example where every matrix can be written out in full, using a deliberately small basis so the arithmetic stays readable. The lesson generalizes directly to cubic B-splines with many more knots; only the numbers change.

Suppose we have three points $x = (0.1, 0.5, 0.9)$ and responses $y = (1.2, 2.0, 1.5)$. On the domain $0 \le x \le 1$, place knots at $0$, $0.5$, and $1.0$, giving three linear hat functions:

Evaluating these at our three $x$ values gives the $3 \times 3$ basis matrix

We will add an intercept column $\mathbf{1}$ as the first column, so the full design becomes

The coefficient vector is $\boldsymbol\beta = (\beta_0, \beta_1, \beta_2, \beta_3)^\top$ with $\beta_0$ the intercept and $(\beta_1, \beta_2, \beta_3)$ giving the smooth.

The unpenalized normal-equation quantities are

These come from the same rules as OLS: each entry of $\mathbf{B}^\top \mathbf{B}$ is a dot product of two columns of $\mathbf{B}$, and each entry of $\mathbf{B}^\top \mathbf{y}$ is a dot product of a column of $\mathbf{B}$ with the response. If we stopped here we would just be doing OLS on the basis expansion, flexible in principle, but with no safeguard against a wiggly fit. The penalty is what adds that safeguard.

The penalty comes next. Our smooth has three coefficients, so the second-difference matrix $\mathbf{D}$ has $3 - 2 = 1$ row, $\mathbf{D} = (1, -2, 1)$. The full penalty matrix has a zero row and column for the intercept (the intercept is not penalized) and $\mathbf{D}^\top \mathbf{D}$ in the smooth block:

Take $\lambda = 0.1$. The penalized system is

with

Now compare the matrix above with the unpenalized $\mathbf{B}^\top \mathbf{B}$. The off-diagonal entries (the $-0.04$ and $0.1$ entries) are the change produced by the penalty. Those couplings come from second coefficient differences; they are the mechanism by which curvature is made expensive.

The displayed four-parameter system is still singular: the three hat functions sum to one, so their constant direction duplicates the intercept. Complete the fit by requiring the smooth to sum to zero at the three observations. Its constraint row is $\mathbf{c}^{\top}=(0,0.8,1.4,0.8)$. We solve the normal equations and that constraint together, using a Lagrange multiplier (an extra unknown that enforces the constraint).

The intercept is the sample response mean because the fitted smooth is centered. The fitted responses no longer interpolate all observations at this positive penalty.

If we now think about the extreme choices:

• With $\lambda = 0$ we recover the unpenalized fit, which on these three points fits every observation exactly (because the $3 \times 4$ design, once the smooth is constrained by identifiability, leaves enough freedom to interpolate all three points), including any noise.
• With $\lambda \to \infty$ the penalty forces $\beta_1 - 2\beta_2 + \beta_3 = 0$, meaning the fitted curve is exactly a straight line (a linear trend in the hat basis).
• For intermediate $\lambda$ (here $0.1$) the fit is a compromise that follows the data closely without allowing spurious curvature.

Both fits use a quadratic penalty, but a difference penalty couples neighboring coefficients, whereas ordinary ridge penalizes coefficient magnitudes. Their common normal-equation structure explains the connection. The next chapter Ridge Regularization (L2 Regularization) develops the ridge view of the penalty in more depth.

The worked example shows the arithmetic. The more important question is whether we can recover the sinusoid and the quadratic from the earlier synthetic data. That is what we implement next.

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.

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 $[-3,3]$ 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.

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.

The choice of $\lambda$ 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.

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.

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 $\lambda$ 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.

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.

• Number of basis functions per smooth (df): controls maximum flexibility. A full cubic basis with df = 4 spans 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 $\lambda$: 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.

• 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 $X$ within its knot bounds, which is what we need to evaluate the fitted smooth on a grid.
• gam_res.edf exposes effective degrees of freedom by coefficient; block sums describe components. gam_res.summary() reports the fitted model summary.

The smoothing parameter $\lambda$ is a hyperparameter. Choose it by scoring fits on observations excluded from fitting those coefficients. An untouched test set then assesses the selected procedure. For dependent observations, choose folds that match the generalization question. To evaluate transfer to unseen subjects, keep each subject in one fold. Forecasting requires a temporal split that matches deployment, including forecasting within already observed subjects when that is the target. Random folds below are appropriate for our independently generated observations.

The basis belongs inside that procedure too. Patsy's bs and statsmodels' BSplines can learn interior knot positions from training-predictor quantiles. Constructing that basis before splitting gives the fit information about validation inputs. That is a leakage path, although the size and direction of its effect on a score depend on the problem. Keep the fitted knot specification and training centering transformation, then reuse them on validation, test, and production rows. Our synthetic bounds $[-3,3]$ were fixed by the data-generation design; these are not estimated from held-out inputs. Both libraries' B-spline transforms reject points outside their boundary knots rather than silently performing polynomial extrapolation.

We use five folds and a predeclared grid of nine smoothing parameters. Each fold constructs its own basis and centering constraint using only that fold's training rows. The fold errors are averaged to select a shared $\lambda$ for the two components. A separate synthetic test sample is generated only after selection, then scored once. The earlier full-sample component plots illustrate fitting; their training RMSE does not select the winning parameter.

The reported test RMSE describes this one synthetic sample from the stated generating model. It is not evidence that GAMs will attain that error on another dataset. A separate $\lambda_j$ per smooth allows different flexibility across predictors, at the cost of a larger tuning problem. Sharing a parameter is a simplifying restriction here, not a universal GAM default.

For a Gaussian linear smoother, generalized cross-validation (GCV) offers another selection criterion: $\mathrm{GCV}=n\|\mathbf{y}-\hat{\mathbf{y}}\|^2/(n-\operatorname{tr}\mathbf{H})^2$, provided the denominator is nonzero. Here $n$ is the training-sample size and $\operatorname{tr}\mathbf{H}$ is the total effective degrees of freedom for that fixed basis and penalty. GCV approximates leave-one-out assessment by replacing individual leverage values with their average; it is not identical to every held-out split. Packages also provide likelihood-based smoothness selection methods. Choose a criterion that respects the data structure and evaluation question.

For a dense implementation with $q$ identifiable coefficients and $n$ rows, forming the cross-product matrix costs $O(nq^2)$ and a direct solve costs $O(q^3)$, with $O(nq+q^2)$ storage when the design is retained. Within a fixed fold, the basis and its cross-products can be reused across smoothing candidates. For non-Gaussian models fitted by IRLS, each iteration solves a weighted penalized system; runtime and convergence depend on the data and link. Effective degrees of freedom describe model flexibility; neither their size nor training RMSE alone diagnoses generalization.

Our constraint makes each fitted component average zero over its training observations, so the intercept and smooths have a specified reference. At two values of $x_j$, the difference in $\hat f_j$ gives the change in the fitted additive predictor while the other inputs are held fixed. With an identity link, a rise of $0.5$ corresponds to a $0.5$ increase in the fitted conditional mean. With a logit link it multiplies fitted odds by $e^{0.5}\approx1.65$; with a log link it multiplies the fitted mean by that factor. It does not imply an additive change of $0.5$ in probability for a logit GAM. On the response scale, a nonlinear inverse link also makes the size of a mean change depend on the other predictors' values.

Component curves describe fitted associations. Interpreting them as causal effects needs additional identification assumptions and a design that supports those assumptions. Correlation between inputs can limit which combinations the data support. Concurvity is the more specific problem in which a smooth term can be approximated by other model terms: it can make the allocation among fitted components unstable even when total predictions look adequate. Exact overlap in unpenalized directions can leave a model unidentifiable despite centering. Check supported predictor combinations and basis conditioning before interpreting a component in isolation.

Inspect residuals for patterns involving pairs of predictors, then compare a plausible interaction model using the same validation protocol. Residual plots can suggest a failure of additivity; absence of a visible pattern does not prove additivity. A tensor-product term $f_{jk}(x_j,x_k)$ starts from the row-wise products of two marginal bases. Before constraints, it has $K_jK_k$ coefficients: two 12-function bases give 144 products. Identifiability constraints and removal of overlapping main-effect directions can reduce that count. The interaction needs its own penalties and evidence that the added flexibility helps. There are $p(p-1)/2$ possible pairs, so domain knowledge can help limit the candidates.

Boundary regions often contain less information than the interior. The displayed example has known finite support; production data need explicit handling of out-of-domain inputs and monitoring for distribution shifts. Natural cubic splines impose linear tails, whereas the B-spline transforms used here refuse points beyond their boundary knots. Neither policy establishes that predictions beyond observed support are reliable. Rebuilding a basis on a new training window changes the fitted model and requires validation again.

For uncertainty, component intervals require a covariance calculation matched to the fitted likelihood, penalty, and centering convention. The fixed-smoothing Gaussian covariance given earlier is one such conditional calculation. It omits smoothing-selection uncertainty and does not turn an association into a causal effect. These limits matter when a curve is used to explain a prediction as well as when it is used for inference.

A generalized additive model is the natural upgrade of linear and spline regression when multiple predictors each have non-linear effects but no interaction structure. The additive predictor $\eta = \beta_0 + \sum_j f_j(x_j)$ is connected to the response mean via a link function $g$, allowing the same machinery to handle Gaussian, binary, count, and other response types. Each smooth $f_j$ is expanded in a spline basis, and a difference penalty regularizes how wiggly the resulting curve can be, turning the Gaussian identity-link fit into a penalized least-squares problem. After identifiability constraints, an invertible penalized system has the solution

which is structurally identical to ridge regression on an expanded feature space. For non-identity links, a weighted penalized system reappears inside an IRLS loop for suitable non-Gaussian GLMs.

The worked procedure selects a shared smoothing parameter by cross-validation within training data, fitting the basis and centering inside every fold, then scores once on independent synthetic test observations. Other smoothness-selection criteria answer related questions under their own assumptions. Interpretation is done by plotting each component separately after centering it to have zero mean, and reading off the effect on the scale of the link function. The additive restriction excludes interactions unless they are included explicitly. When interactions are suspected, tensor-product smooths can be added at the cost of sample size and interpretability.

A GAM extends a GLM by allowing smooth terms in its predictor. It contains a linear GLM when those terms are linear, and its component curves expose the fitted additive associations. It is a good first model to try whenever a linear regression or a linear GLM leaves obvious non-linear residuals. Its Gaussian fitting mechanism, quadratic penalized least squares on an expanded basis, is closely related to ridge regression. The next chapter Ridge Regularization (L2 Regularization) develops that connection further. Lasso instead uses an absolute-value penalty and requires a different optimization problem.

Ready to test your understanding? Take this quick quiz to reinforce what you've learned about generalized additive models.