Generalized Additive Models › Formula

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$:

$$\eta_i = \beta_0 + f_1(x_{i1}) + f_2(x_{i2}) + \cdots + f_p(x_{ip}),$$

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:

$$g(\mu_i) = \eta_i, \qquad \mu_i = \mathbb{E}[y_i \mid \mathbf{x}_i].$$

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:

$$f_j(x) = \sum_{k=1}^{K_j} \beta_{jk} \, b_{jk}(x),$$

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

$$\eta_i = \beta_0 + \sum_{j=1}^p \sum_{k=1}^{K_j} \beta_{jk} \, b_{jk}(x_{ij}).$$

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:

$$\text{Pen}(\mathbf{f}_j) = \lambda_j \sum_{k=3}^{K_j} \left( \beta_{jk} - 2\beta_{j,k-1} + \beta_{j,k-2} \right)^2.$$

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

$$\text{Pen}(\mathbf{f}_j) = \lambda_j \, \boldsymbol\beta_j^\top \mathbf{D}^\top \mathbf{D} \, \boldsymbol\beta_j = \lambda_j \, \boldsymbol\beta_j^\top \mathbf{S}_j \, \boldsymbol\beta_j,$$

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

$$\mathcal{L}(\boldsymbol\beta) = \sum_{i=1}^n \left( y_i - \beta_0 - \sum_{j=1}^p \sum_{k=1}^{K_j} \beta_{jk} \, b_{jk}(x_{ij}) \right)^2 + \sum_{j=1}^p \lambda_j \, \boldsymbol\beta_j^\top \mathbf{S}_j \boldsymbol\beta_j.$$

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,

$$\mathcal{L}(\boldsymbol\beta) = \| \mathbf{y} - \mathbf{B}\boldsymbol\beta \|_2^2 + \boldsymbol\beta^\top \mathbf{S} \boldsymbol\beta,$$

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

$$\hat{\boldsymbol\beta} = \left( \mathbf{B}^\top \mathbf{B} + \mathbf{S} \right)^{-1} \mathbf{B}^\top \mathbf{y}.$$

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

$$\hat{\boldsymbol\beta}^{(t+1)} = \left( \mathbf{B}^\top \mathbf{W}^{(t)} \mathbf{B} + \mathbf{S} \right)^{-1} \mathbf{B}^\top \mathbf{W}^{(t)} \boldsymbol{z}^{(t)},$$

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.