Penalized Splines as Mixed Models

Why a Penalty Can Be a Variance Component

A flexible spline can follow detailed patterns in longitudinal measurements, but an overly flexible trajectory may also follow noise. Penalized splines address this by allowing a rich basis while discouraging unnecessary variation. Their mixed-model representation gives a useful interpretation of that trade-off: the variance of spline coefficients controls how much the trajectory can depart from its unpenalized part.

Consider measurements \(y_{ij}\) from subject \(i\) at time \(t_{ij}\), where \(i=1,\ldots,N\) and \(j=1,\ldots,n_i\). For the derivation, stack all measurements into a vector \(y\) and write the spline contribution as \(B\gamma\). Here \(B\) is the matrix of evaluated spline basis functions, and \(\gamma\) is the vector of basis coefficients. Other fixed covariates can be included in the unpenalized part of the model.

For a simple Gaussian model with independent errors of common variance, the criterion is

The penalty matrix need not penalize every coefficient. For example, common spline penalties leave a low-dimensional space of simple functions unpenalized. Treating all coefficients as equally penalized would therefore change the model rather than represent the same spline fit.

Separating Penalized and Unpenalized Directions

Because \(S\) is positive semidefinite, it has an eigendecomposition. Let \(U_0\) contain orthonormal eigenvectors for its zero eigenvalues, and let \(U_r\) contain orthonormal eigenvectors for its \(r\) positive eigenvalues. Write the positive eigenvalues in the diagonal matrix \(D\). Then \(S=U_rDU_r^{\mathsf T}\), and the columns of \(U_0\) span the null space of \(S\).

Reparameterize the coefficients as \(\gamma=U_0a+U_rD^{-1/2}b\). The vector \(a\) describes the unpenalized part, while \(b\) describes the penalized part in scaled coordinates. In these coordinates, the penalty becomes \(b^{\mathsf T}b\). Defining \(X=BU_0\) and \(Z=BU_rD^{-1/2}\), the spline fit is therefore expressed as

This transformation is the key to the mixed-model connection: the penalty is now a squared length of the coefficient vector \(b\), the form of a Gaussian quadratic penalty.

Mixed-Model Representation and Its Derivation

Assume for now that the conditional errors are independent and normally distributed, \(\varepsilon\sim N(0,\sigma^2I)\), and that the penalized coefficients follow an independent Gaussian distribution \(b\sim N(0,\tau^2I)\), independently of \(\varepsilon\). Treat \(a\) as a fixed-effect coefficient vector. The conditional model is \(y=Xa+Zb+\varepsilon\). Given the data, the negative log likelihood for \(b\), up to terms that do not depend on \(a\) or \(b\), is proportional to \((y-Xa-Zb)^{\mathsf T}(y-Xa-Zb)/\sigma^2\). The negative log density of \(b\), up to a constant, is \(b^{\mathsf T}b/\tau^2\). Adding these terms and multiplying by \(\sigma^2\) gives the joint criterion

Comparing this with the penalized least-squares criterion establishes the relationship:

Integrating over \(b\) also gives the marginal model \(y\sim N(Xa,\tau^2ZZ^{\mathsf T}+\sigma^2I)\). Thus the same model can be viewed either as a spline fit with a penalty or as a mixed model with a particular covariance contribution. The equivalence concerns the likelihood and coefficient shrinkage under the stated assumptions; it does not imply that spline coefficients are observed subject-level effects.

The variance ratio explains the direction of the smoothing effect. If \(\tau^2\) is small relative to \(\sigma^2\), then \(\lambda\) is large and penalized coefficients are pulled strongly toward zero. As \(\tau^2\) increases relative to \(\sigma^2\), the penalty weakens and the fitted trajectory can express more variation. At the boundary \(\tau^2=0\), the penalized component is absent; this corresponds to the limiting case \(\lambda\to\infty\). A large variance does not guarantee a wiggly fitted curve, but it permits greater variation when supported by the data.

Estimating Smoothness and Checking the Representation

In practice, the variance components \(\sigma^2\) and \(\tau^2\) can be estimated from the data, commonly by maximum likelihood or restricted maximum likelihood under the Gaussian mixed model. Their ratio then determines the corresponding smoothing parameter. This lets the data inform the degree of shrinkage instead of fixing it in advance. However, a variance estimate at or near zero is a boundary case, and uncertainty about it may not be well described by routine large-sample approximations.

The equivalence above assumes a correctly specified Gaussian model, the stated penalty, and a clear separation between penalized and unpenalized directions. In longitudinal data, observations from the same subject may have correlated errors or additional subject-specific effects. If those features are present, an independent-error model is not adequate merely because the spline has been written as a mixed model. A covariance structure for the errors or other random effects must be included as appropriate. When the residual covariance is \(\sigma^2R_0\) for a known matrix \(R_0\), the corresponding fit uses generalized least squares; the penalty-to-variance relationship remains \(\lambda=\sigma^2/\tau^2\) when the criterion is scaled consistently by \(\sigma^2\).

Diagnostics should therefore examine both the trajectory fit and the model assumptions. Residual patterns over time can reveal misspecified mean structure or remaining temporal dependence. A trajectory that changes sharply where the data are sparse may signal weak support for the fitted shape. A variance component near zero indicates that the data provide little evidence for variation in the penalized directions under the chosen basis and penalty; it does not prove that all time trends are linear. Conversely, a large estimated spline variance can permit complexity, but does not itself establish that the curve is well supported. Basis dimension, penalty form, residual behavior, and uncertainty all matter.

A common alternative is to fit a fixed, unpenalized spline with a selected number of basis terms. That approach gives direct fixed coefficients but does not use a variance component to shrink penalized directions. Penalized splines instead allow a relatively rich basis and regulate its effective complexity through shrinkage. The mixed-model representation is especially useful for variance-component estimation and for combining spline terms with other random effects. Neither approach removes the need to choose a basis and assess fit; the representation also does not make the spline automatically appropriate for every longitudinal response.

Common Pitfall and Practical Limits

Wrong reasoning: “The spline coefficients are random effects, so they must represent individual subjects’ deviations from a population trajectory.” Correction: In this representation, \(b\) is a random vector in the penalized spline basis used to express the trajectory. Its distribution supplies the quadratic penalty. It is not automatically a set of subject-specific deviations. Additional subject-level random effects, if required, are distinct model components.

Another mistake is to compare numerical values of \(\lambda\) across bases without checking their scaling. Rescaling a basis or penalty changes the numerical coefficient associated with the same fitted curve. Interpret \(\lambda\) through the specified penalty and its variance ratio, not as a scale-free smoothness score. Finally, a good mixed-model representation does not guarantee that the chosen spline basis captures the trajectory of interest or that Gaussian errors are appropriate.