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?
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 bygam()
to set up the basis, hence you got 9 basis functions (the default is to setk = 10
and then we lose one for identifiability constraints). If you had doneknots = list(x = c(1:5))
and useds(x, bs = 'cr', k = 5)
you would have 4 basis functions (5 minus 1 lost for identifiability constraints). $\endgroup$