Smoothing methods for gam in mgcv package?

Question

I am currently working with gam models in the mgcv package and for me the smoothing methods are a bit confusing and I hope that you guys can help me to understand that better. So here is what I've understood so far:

gam models have the advantage that I can model complex functions, meaning I can model continuous variables as functions $f(z)$. But the challenge is to find those $f$ functions.

Some basic approaches are B-Splines, P-Splines, Cubic Splines and Thin Plate Splines.

The default setting in gam is the thin plate spline. So I was wondering what it basically does. From my understanding so far, the advantage in that method is, that you don't have to specify the number of knots k. You start with the maximal number of knots, then gam chooses via GCV which k suits best for the function.

Did I get this right?

Answer 1

mgcv uses a thin plate spline basis as the default basis for it's smooth terms. To be honest it likely makes little difference in many applications which of these you choose, though in some situations or with very large data set sizes, other basis types might be used to good effect. Thin plate splines tend to have better RMSE performance than the other three you mention but are more computationally expensive to set up. Unless you have a reason to use the P or B spline bases, use thin plate splines unless you have a lot of data and if you have a lot of data consider the cubic spline option.

k doesn't set the number of knots, at least not in the default thin plate spline basis. What k does is to set the dimensionality of the basis expansion; you'll end up with k - 1 basis functions. In mgcv Simon Wood does a trick to reduce the rank of basis dimension. IIRC, in the usual thin plate spline basis there is a knot at each data location, but this is wasteful as once you've set up this large basis you end up using far fewer degrees of freedom in the fitted function. What Simon does is to eigen decompose the matrix of basis functions and choose the eigenvectors of the decomposition corresponding to the k - 1 largest eigenvalues. This has the effect of concentrating the main wiggliness "information" of the full basis in a reduced rank form.

The choice of k is important and the default is arbitrary and something you want to check (see gam.check()), but the critical observation is that you want to set k to be large enough to contain the envisioned dimensionality of the underlying function you are trying to recover from the data. In practice, one tends to fit with a modest k given the data set size and then use gam.check() on the resulting model to check if k was large enough. If it wasn't, increase k and refit. Rinse and repeat...

You are most likely going to want to fit the model using REML (or ML) smoothness selection via method = "REML" or method = "ML": this treats the model as a mixed effects one with the wiggly parts of the spline bases being treated as special random effects terms. Simon Wood has shown that REML (or ML) selection performs better than GCV, which can undersmooth in situations where the objective function is flat around the optimal smoothness parameter value.

The ridge penalty mentioned by @generic_user is taken care of for you, so you can ignore this part of setting up the model.

Answer 2

Try this:

#Fake data
library(mgcv)
x <- runif(10000, 0, 10)
y <- log(x) + sin((x^2)/10)+rnorm(10000)
plot(x,y)

Try three different basis functions with default "knot locations"

plot(gam(y~s(x, bs = 'ps')))
plot(gam(y~s(x, bs = 'cr')))
plot(gam(y~s(x, bs = 'ts')))

They're basically identical. Now the same ones with 40 knots each.

plot(gam(y~s(x, bs = 'ps', k = 40)))
plot(gam(y~s(x, bs = 'cr', k = 40)))
plot(gam(y~s(x, bs = 'ts', k = 40)))

...a little over-wiggly, but really not that different.

Now futz with the smoothing parameter, and you see where the action is:

plot(gam(y~s(x, sp = 0)))
plot(gam(y~s(x, sp = 1e6)))
plot(gam(y~s(x))) #default smoothing parameter, estimated by generalized cross-validation

So in sum, knot location doesn't matter much when you specify more knots than you probably need, and then penalize them down via the ridge penalty.

If you don't know what I'm talking about with the ridge penalty, pick up Simon Wood's book on GAMs. Note the author is the same guy who wrote the R package.