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

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)?

• If you want to have a deeper dive into P-splines, you could give a look to the following website psplines.bitbucket.io/index.html. There is also a R-package and a book (I think it is in press now) on the topic Commented Mar 6, 2021 at 12:19
• @Gi_F. Thanks, thats very interesting. I will take a look. Commented Mar 6, 2021 at 12:58

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.

• Great answer Gavin, thank you! Just to be clear: The penalty matrix $\mathbf{S}$ differs for different type of smoothers/bases, correct? And it basically encodes what we mean by "wiggliness"? Then, an algorithm is used to find the best $\lambda$ to penalize the kind of wiggliness that is defined in terms of $\mathbf{S}$? Commented Mar 3, 2021 at 18:54
• Yes, the form of $\mathbf{S}$ varies btw different types of spline; if you want to look at them graphically, try out my gratia package & the penalty() function therein. Yes, P splines define wiggliness or complexity in terms of differences between adjacent groups of model coefficients. Thin plate splines define wiggliness in terms of some integrated squared derivative. Random effects are just encoded as an identity matrix for $\mathbf{S}$. However, the wiggliness of the function(s) is actually defined by $\boldsymbol{\beta}^{\mathsf{T}} \mathbf{S} \boldsymbol{\beta}$ not just $\mathbf{S}$ Commented Mar 3, 2021 at 19:03
• Got it, thank you. I thought from your video that it is always the integrated second derivative that is being penalized but that seems to be only for thin plate splines. I will definitely look at gratia, thanks. Commented Mar 3, 2021 at 19:05
• Yeah, I could have been clearer about that in the video that wiggliness was not derivative-based in P splines, and some of the other specialist ones. Commented Mar 3, 2021 at 19:06
• It's fine. I could have read Simon Wood's book more carefully. Commented Mar 3, 2021 at 19:08