Take a simple analysis using generalized additive models in mgcv:

 x <- seq(0, 1, length=500)
 y <- sin(2*pi*x)^2 + rnorm(500, sd=0.2)
 m <- gam(y ~ s(x))

I know how to see the final tuning parameter selected (type m$sp). Is there a way to plot the GCV vs. smoothing parameter to graphically see how it was selected? I know you can do stuff like this in glmnet when you choose the elastic net penalty size. Anything like this for mgcv? Thanks

  • $\begingroup$ What would you put on the y axis? OOS mse? $\endgroup$ – generic_user Apr 3 '18 at 15:36
  • $\begingroup$ the smoothing parameter $\endgroup$ – within_person Apr 3 '18 at 16:06
  • $\begingroup$ You'd put the GCV score on the y-axis and the log smoothness parameter on the x-axis. $\endgroup$ – Reinstate Monica - G. Simpson Apr 3 '18 at 16:08
  • $\begingroup$ How do you extract those values from the fitted object? $\endgroup$ – within_person Apr 3 '18 at 16:18
  • $\begingroup$ You can't; the GAM isn't fitted as a path of models in the sense of lassso penalties and similar. It is optimised using very efficient algorithms and we don't need actual crossvalidation to choose the optimal lambda, hence models for a sequence of lambdas are not estimated during fitting. You need to optimise the model for a fixed value of lambda and extract the score you want as I show in the answer below. $\endgroup$ – Reinstate Monica - G. Simpson Apr 3 '18 at 16:40

The way to do this is to fit the GAM for known values of the smoothness parameter $\lambda$, and set the smoothness parameter for each fit. You can then extract the GCV score from each model and plot them against the sequence of $\lambda$ values.

To fix the smoothness parameter you use the sp argument. For your example data this means setting sp to be a length one vector. If you are doing this for multiple parameters you'll need to fit models over a grid of 2 or more smoothness parameters.

Here is an example using your dummy data set:

## packages

## data
x <- seq(0, 1, length=500)
y <- sin(2*pi*x)^2 + rnorm(500, sd=0.2)
df <- data.frame(x = x, y = y)

## Generate GCV traces
lambda <- 1e-6
gcv <- numeric(length = 100)

## loop over values of lambda, fixing sp
for (i in seq_along(gcv)) {
    m <- gam(y ~ s(x, k = 20), data = df, method = "GCV.Cp", sp = lambda)
    gcv[i] <- m$gcv.ubre
    lambda <- lambda * 1.1

## gather results
res <- data.frame(lambda = 1e-6 * 1.1^{0:99}, gcv = gcv)

## plot
ggplot(res, aes(x = lambda, y = gcv)) +
    geom_line() +
    scale_x_log10() +
    xlab(expression(lambda)) +

This produces

enter image description here

You can do this for the other fitting methods available in mgcv also; the score is always in $gcv.ubre, despite the name of that component.

You'll notice that it take a lot longer to fit those models than it does to find the optimum value of $\lambda$ using GCV if you fitted without setting sp. mgcv has very efficient algorithms for selecting the optimal lambda and the steps it takes along the GCV surface (curev) will not produce a nice curve and typically aren't worth saving in the model object as a result. Hence the need to do the fitting iteratively, as shown above, to yield a nice curve of values over $\lambda$.

The above methods requires a little tuning to get right; you need to set lambda small enough initially and not increase it took quickly to home in on the minimum of the curve.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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