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

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?

• (Excuse the corrections made on the fly soon after posting.) Jun 21 at 15:02
• You made a mistake in coding the knots, you have to pass a list with named components. Because your list didn't name any components the knots you passed weren't used by gam() to set up the basis, hence you got 9 basis functions (the default is to set k = 10 and then we lose one for identifiability constraints). If you had done knots = list(x = c(1:5)) and used s(x, bs = 'cr', k = 5) you would have 4 basis functions (5 minus 1 lost for identifiability constraints). Jun 23 at 11:39

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.

• Many thanks! I have a follow-up question. How comes GAM is not simply using a single cubic spline to describe the data? with knots=list(1,2,3,4,5), the minimum k value allowed by gam() is 4, but using a single cubic spline with lines(smooth.spline(x, y)) gives a very similar fit as the GAM with 5 knots and 4 basis functions. I am happy to ask this as a stand-alone question if you think it's worth it. Jun 22 at 10:03
• Just a note; when k is not set by the user it takes the default value of -1 which forces gam to use the essentially arbitrary default of k = 10 for univariate smooths, this sets the number of knots in the cubic regresison spline basis which results in the same number of basis function; the knots are spread evenly through the data. Identifiability constraints then alter the basis so that we can include a separate intercept, and that forces the loss of one function. Jun 23 at 11:46
• @MarcoPlebani the knots need to be list(x = c(1:5)) as I mentioned above; the minimum isn't 4, it is 3 if you don't specify knots. If you specify knots correctly, then k must equal the number of knots specified. Jun 23 at 11:52
• @GavinSimpson thanks for the details!
– jld
Jun 23 at 15:29
• The red line in your plot (the function fitted in the GAM) is a penalised cubic regression spline. The functions drawn by your matplot() command are the basis functions of the the cubic regression spline. As such I believe your second option is correct. Jun 23 at 17:46