4
$\begingroup$

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)

enter image description here

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

$\endgroup$
2
  • $\begingroup$ (Excuse the corrections made on the fly soon after posting.) $\endgroup$ Jun 21, 2021 at 15:02
  • 1
    $\begingroup$ 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). $\endgroup$ Jun 23, 2021 at 11:39

1 Answer 1

4
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$ Jun 22, 2021 at 10:03
  • 2
    $\begingroup$ 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. $\endgroup$ Jun 23, 2021 at 11:46
  • 2
    $\begingroup$ @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. $\endgroup$ Jun 23, 2021 at 11:52
  • $\begingroup$ @GavinSimpson thanks for the details! $\endgroup$
    – jld
    Jun 23, 2021 at 15:29
  • 1
    $\begingroup$ 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. $\endgroup$ Jun 23, 2021 at 17:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.