How can basis functions in a GAM be described by a single parameter each?

Question

One way of fitting Generalised Additive Models (GAM) involves using cubic splines as basis functions. Cubic splines are splines constructed of piecewise third-order polynomials (source). If I understand correctly, given n knots, the data will be subset in n+1 intervals, and the data in each interval will be modelled by a third-order polynomial.

Given a set of data subdivided by n knots, I would expect a cubic spline to be described by at least n + 1 + (n+1)*4 parameters, namely:

• the coordinates of the n knots,
• four parameters for each cubic spline used to model each of the n-1 subsets of data (a, b, c, and d in y = a + bx + cx^2 + dx^3*),
• at least one parameter for the function of the sum of the basis functions that produces the model estimates.

To test this I generated data according to a sinusoid pattern using R (after this example):

set.seed(1)
x <- seq(0, pi * 2, 0.1)
sin_x <- sin(x)
y <- sin_x + rnorm(n = length(x), mean = 0, sd = sd(sin_x / 2))
Sample_data <- data.frame(y,x) 

And I modelled them as:

gam_y <- mgcv::gam(y ~ s(x, bs="cr"), knots=list(1,2,3,4,5), method = "REML")

# plot data
plot(y~x)
# plot GAM predictions
gam_pred <- predict(gam_y, newdata = data.frame(x = Sample_data[,2]))
lines(x, gam_pred, col="red", lwd=2)
# display basis functions
model_matrix <- predict(gam_y, type = "lpmatrix")
matplot(x, model_matrix[,-1], type = "l", lty = 2, add = T)

(GAM fit as a solid line, basis functions as dashed lines)

If I were right, each basis function should be described by 30 parameters (5 + 1 + (5+1)*4). Instead, each basis function is defined by a single parameter:

gam_y$coefficients
(Intercept)      s(x).1      s(x).2      s(x).3      s(x).4      s(x).5      s(x).6      s(x).7      s(x).8 
 0.05361059  0.67279976  0.99708672  0.88434236  0.28890966 -0.35159312 -0.82555194 -0.96022992 -0.61303024 
     s(x).9 
 0.05663602 

I am clearly confused. Can anyone help me understand how each basis function can be defined by a single parameter?

Answer 1

If we had $n$ knots $\xi_1 < \dots < \xi_n$ and wanted to fit a cubic polynomial on each of the $n+1$ intervals thus defined, then we could do something like $$ f(x) = \sum_{i=0}^n (\beta_{i0} + \beta_{i1} x + \beta_{i2} x^2 + \beta_{i3} x^3) \mathbf 1_{\xi_i < x \leq \xi_{i+1}} $$ (adding knots $\xi_0 = -\infty$ and $\xi_{n+1} = +\infty$ for convenience) and here there are $4 \cdot(n+1)$ parameters. We don't count the knot locations as parameters because we want this to be a linear combination of known basis functions. Specifying the knots in advance like you are doing is treating them like a hyperparameter (we can also dodge the knot issue by using a smoothing spline which shifts the issue to the regularization hyperparameter). Many algorithms set the knots by using either quantiles of x or just an evenly spaced sequence between the smallest and largest values of $x$ in the data, and then we effectively condition on the knots.

The issue with this model is that there's nothing that makes $f$ continuous or smooth. We not only want cubic polynomials between each knots but we want $f$ to be twice continuously differentiable everywhere. We know $f$ already is smooth everywhere that's not at a knot, since it's just a polynomial there, so we only need to make sure things work out on the knots.

Here's how this goes: consider a knot $\xi_i$. As $x$ approaches $\xi_i$ from below, we have $$ f(x) = \beta_{i-1,0 } + \beta_{i-1, 1}x + \beta_{i-1, 2}x^2 + \beta_{i-1, 3}x^3 $$ while as $x$ approaches $\xi_i$ from above we'll have $$ f(x) = \beta_{i,0 } + \beta_{i, 1}x + \beta_{i, 2}x^2 + \beta_{i, 3}x^3. $$ In order for $f$ to be continuous as $\xi_i$ we need these two to agree as $x\to\xi_i$, so this means $$ \beta_{i-1,0 } + \beta_{i-1, 1}\xi_i + \beta_{i-1, 2}\xi_i^2 + \beta_{i-1, 3}\xi_i^3 = \beta_{i,0 } + \beta_{i, 1}\xi_i + \beta_{i, 2}\xi_i^2 + \beta_{i, 3}\xi_i^3. $$ $\xi_i$ is fixed so this is a single linear equation in $8$ unknowns and we eliminate one parameter. Next, requiring being once continuously differentiable at $\xi_i$ gives us $$ \beta_{i-1, 1} + 2\beta_{i-1, 2}\xi_i + 3\beta_{i-1, 3}\xi_i^2 = \beta_{i, 1} + 2\beta_{i, 2}\xi_i + 3\beta_{i, 3}\xi_i^2 $$ which is another linear constraint and we eliminate another parameter. Finally, being twice continuously differentiable means $$ 2\beta_{i-1, 2} + 6\beta_{i-1, 3}\xi_i = 2\beta_{i, 2} + 6\beta_{i, 3}\xi_i $$ and a third parameter is lost. Thus once we have our initial cubic, taking 4 parameters, each additional cubic only requires a single parameter because a cubic that is required to smoothly join a previous one only has so much freedom to wiggle.

All together this means we have $4 + n$ parameters required.

Another way to see this is to use a basis that directly incorporates the constraints. One such basis is to take $h_i(x) = x^{i-1}$ for $i=1,2,3,4$ and then $h_i(x) = (x - \xi_{i-4})^3_+$ for $i = 5, \dots, n+4$. This is called the truncated power basis and I discuss why it works in my answer here.

We then have $f(x) = \sum_{i=1}^{n+4}\beta_i h_i(x)$ and there are $n+4$ parameters.

---

There's one extra complexity in your example. Some of the arguments in s are using default values which change behavior. In particular, you aren't setting k so gam uses choose.k to determine the number of basis functions rather than using the knots. It also penalizes by default so you don't have a standard regression spline as I wrote it above.