13
$\begingroup$

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?

$\endgroup$
1
  • 1
    $\begingroup$ What you're missing is the ridge penalty applied to the spline. Choice of basis function is totally second-order. $\endgroup$ Dec 7, 2016 at 2:03

2 Answers 2

21
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ I thought that thin plate spline basis would be used for 2D or higher-dimensional data, but I do notice that it does work for 1D when I set bs to 'tp' in gam(). What am I missing? $\endgroup$
    – bluepole
    Jun 7, 2020 at 20:57
  • $\begingroup$ That a TPRS can't be used for univariate smooths...? shrug Not sure what your question is; your assumption/thought is just incorrect. TPRS can be used for 2+d but there's a strong assumption of isopropyl, hence the prominence of tensor product smoothers in mgcv. $\endgroup$ Jun 7, 2020 at 21:15
  • 2
    $\begingroup$ The default in mgcv's s() is bs = "tp", so any smooth where bs is not specified is a TPRS. How are these different? In the same way that cubic regression splines are different to P-splines, to, B-splines, .... More specifically the advantage of TPRS is there is no knot selection problem; we have knots at every data point. But now we have the practical problem of computing with so many basis functions. Hence Simon Wood came up with low-rank version TPRS, which uses an eigen decomposition to concentrate the signal in the TPRS basis into the first $k$ eigenvectors = basis functions. $\endgroup$ Jun 8, 2020 at 13:59
  • 1
    $\begingroup$ @bluepole, that seems like a good new question :-) $\endgroup$ Dec 7, 2022 at 13:43
  • 1
    $\begingroup$ @bluepole I think you can use fit = FALSE in the gam() call - I can't recall what this sets up, but it returns an object that you can pass to the G argument of gam() when you need to do fitting. The question then is whether you can update the G object efficiently? If the basis construction is that expensive; I'd just use a CRS as you're unlikely to notice the difference in fit but it will be much more efficient. $\endgroup$ Feb 21 at 7:25
4
$\begingroup$

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.

$\endgroup$

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.