Generalized additive models: What exactly is being penalized when using a P-spline smoother with $\texttt{mgcv}$?

Question

Generalized additive models (GAMs) avoid overfitting by introducing a $\color{#D55E00}{\text{penalty term}}$ in the loss function: $$ ||\mathbf{y}-\mathbf{X\beta}||^{2} + \color{#D55E00}{\lambda\mathbf{\beta}^{\intercal}\mathbf{S\beta}} $$ where $\lambda$ controls the wiggliness of the smoother or basis. The optimal $\lambda$ can be estimated from the data using several methods, such as REML (restricted maximum likelihood), ML (maximum likelihood) or on the basis of the AIC (see the wonderful introductiory video on GAMs by @GavinSimpson).

The function s() in Rs mgcv package let's you choose the type of smoother or basis (thin plate splines being the default). Among those smoothers are P-splines (Penalized B-Splines), introduced by Eilers & Marx (1996). I must admit that I don't fully understand what they're doing, but the documentation of mgcv states:

They [P-splines] combine a B-spline basis, with a discrete penalty on the basis coefficients, and any sane combination of penalty and basis order is allowed.

So it seems to me that we're penalizing twice: Once by using P-splines and once when gam is fitting the model using the wiggliness-penalty $\lambda$.

My questions are: What exactly is penalized when we use P-splines in a GAM that already penalizes the wiggliness of the smooth? And what's the benefit of penalizing twice (if that's what's happening)?

Answer 1

P splines in mgcv are not penalised twice, they just use a different form of penalty matrix where we penalize some particular order of differences between adjacent $\beta_i$.

It's important to note that GCV, REML, etc are algorithms for choosing $\boldsymbol{\lambda}$, the smoothness parameters; because of the way the model complexity is defined as a function of the model coefficients, fitting a penalised GAM basically involves finding $\boldsymbol{\beta}$ and $\boldsymbol{\lambda}$ to minimise the negative penalised log-likelihood. How we measure wiggliness of the fitted functions is determined by the form of the penalty matrix $\mathbf{S}$ but in defining these matrices we're not penalising anything with them in and by themselves. The penalization or shrinkage comes in when we choose values for $\boldsymbol{\lambda}$ because as we increase $\boldsymbol{\lambda}$ we have to shrink the values of $\boldsymbol{\beta}$.

The wigglines measure in a penalized GAM is

$$\boldsymbol{\beta}^{\mathsf{T}} \mathbf{S} \boldsymbol{\beta}$$

where $\mathbf{S}$ is the penalty matrix. For P splines, instead of a directly defining penalty matrix $\mathbf{S}$ that measures the wiggliness of the basis functions, as we do with thin plate splines, we define a matrix $\mathbf{P}$ such that our wiggliness measure is now, following the argument of Wood (2017, pp. 206; slightly modified to fix a typo)

$$\boldsymbol{\beta}^{\mathsf{T}} \mathbf{P}^{\mathsf{T}}\mathbf{P} \boldsymbol{\beta}$$

and hence $\mathbf{P}^{\mathsf{T}}\mathbf{P}$ is now our penalty matrix, and $\mathbf{P}$ the square root of the penalty matrix.

In a P spline we want to penalise squared differences between adjacent coefficients

$$\sum_{i=1}^{k-1}(\beta_{i+1} - \beta_i)^2 = \boldsymbol{\beta}^{\mathsf{T}} \mathbf{P}^{\mathsf{T}}\mathbf{P} \boldsymbol{\beta}$$

we define $\mathbf{P}$ as

$$\mathbf{P} = \begin{bmatrix} -1 & 1 & 0 & \cdot & \cdot \\ 0 & -1 & 1 & 0 & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix}$$

so that we have

$$\begin{bmatrix} \beta_2 - \beta_1 \\ \beta_3 - \beta_2 \\ \cdot \\ \cdot \end{bmatrix} = \mathbf{P}\boldsymbol{\beta}$$

$\mathbf{P}$ itself is the square root of the penalty matrix, hence we end up with

$$\boldsymbol{\beta}^{\mathsf{T}} \mathbf{P}^{\mathsf{T}}\mathbf{P} \boldsymbol{\beta} = \boldsymbol{\beta}^{\mathsf{T}} \begin{bmatrix} 1 & -1 & 0 & \cdot & \cdot & \cdot \\ -1 & 2 & -1 & 0 & \cdot & \cdot\\ 0 & -1 & 2 & -1 & 0 & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix} \boldsymbol{\beta}$$

To be clear; all of the above is cribbed/summarised from Wood (2017, pp 206).

Hopefully it's clear that the penalty is defined directly in terms of differences of adjacent groups (the difference can be higher order than the 1st-order differences shown here) of coefficients. In other words more complex models are those where there are large differences between adjacent groups of spline/model coefficients.

In the thin plate splines and cubic regression splines by way of contrast, the penalty matrix $\mathbf{S}$ measures the wiggliness of each basis function and how the wiggliness of one function modifies the wiggliness of the other, i.e. how the wigglinesses of basis functions covary. These penalty matrices are defined in terms of the wiggliness of the basis functions themselves, but in a P spline the penalty is applied to the coefficients and doesn't measure the wiggliness of the basis functions themselves.

This is to say that P splines aren't typically defined in terms of derivative-base penalties, while the penalties for thin plate splines and cubic regression splines are derivative-based. I say typically, because Simon Wood has a section (5.3.4) in his GAM book (Wood 2017) showing how you can create P splines with derivative-based penalties, but you can read that for yourself if you're interested).

If you want to play with this in R (or translate to the language of your choice), Simon includes the following simple snippet to illustrate how to create $\mathbf{P}$ for a given size of $k$, and hence what the resulting penalty matrix $\mathbf{S}$ looks like:

k <- 6
P <- diff(diag(k), differences = 1)
S <- crossprod(P)

with S now being

> S
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1   -1    0    0    0    0
[2,]   -1    2   -1    0    0    0
[3,]    0   -1    2   -1    0    0
[4,]    0    0   -1    2   -1    0
[5,]    0    0    0   -1    2   -1
[6,]    0    0    0    0   -1    1

To get higher-order differences, change the differences argument

> P <- diff(diag(k), differences = 2)
> S <- crossprod(P)
> S
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1   -2    1    0    0    0
[2,]   -2    5   -4    1    0    0
[3,]    1   -4    6   -4    1    0
[4,]    0    1   -4    6   -4    1
[5,]    0    0    1   -4    5   -2
[6,]    0    0    0    1   -2    1

As you increase the order of differences in the penalty, the penalty will act to pull, shrink, larger groups of consecutive $\beta_i$ towards one another as the non-zero elements in $\mathbf{S}$ cover more of the $\beta_i$.

References

Wood, S.N., 2017. Generalized Additive Models: An Introduction with R, Second Edition. CRC Press.