2

Consider the simple GAM fit as below:

library(mgcv)
my.gam <- gam(y~s(x), data=mydata)
  1. Is there anyway to return the estimated smoothing parameter (lambda) so that I can save it? I know that lambda is given in output as 'GCV score', but I need a specific code for returning it.
  2. How can I set lambda to a desired value?
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
M. Er
  • 881
  • 1
  • 7
  • 13

1 Answers1

5

summary() does not return smoothing parameters. You have mixed up GCV score with smoothing parameter. Consult a local statistician if you don't understand those concepts, or raise a question on Cross Validated. I will only show you how to extract and set smoothing parameters.

Consider an example:

library(mgcv)
set.seed(2)
dat <- gamSim(1, n=400, dist="normal", scale=2)
b <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data=dat)

You can get internal smoothing parameters from:

b$sp
#       s(x0)        s(x1)        s(x2)        s(x3) 
#3.648590e+00 3.850127e+00 1.252710e-02 4.986399e+10

But these are not lambda. They differ from lambda by some positive scaling factors. It is usually sufficient to use sp for smoothing parameters. If you want to set it to a fixed value, do for example:

b1 <- gam(y ~ s(x0, sp = 0) + s(x1, sp = 0) + s(x2, sp = 0) + s(x3, sp = 0),
          data = dat)

This essentially disables penalization for all smooth terms. Note that setting sp to a negative value implies auto-selection of sp.

lambda and sp:

enter image description here

sapply(b$smooth, "[[", "S.scale") / b$sp
#       s(x0)        s(x1)        s(x2)        s(x3) 
#6.545005e+00 5.326938e+00 1.490702e+03 4.097379e-10

Sometimes getting lambda is necessary. When considering smooth functions as random effects or random fields, there is

variance_parameter_of_random_effect = scale_parameter / lambda

where the scale parameter can be found in b$scale (for Gaussian model this is also b$sig2). See a related question: GAM with "gp" smoother: how to retrieve the variogram parameters?

Follow-up

Yes, I need the exact value of lambda, so thanks for the neat code. Yet I'm interested in knowing more about the scaling factor. Where can I read more about it in addition to the package manual?

Have a read on ?smoothCon:

smoothCon(object,data,knots=NULL,absorb.cons=FALSE,
          scale.penalty=TRUE,n=nrow(data),dataX=NULL,
          null.space.penalty=FALSE,sparse.cons=0,
          diagonal.penalty=FALSE,apply.by=TRUE,modCon=0)

scale.penalty: should the penalty coefficient matrix be scaled to have
          approximately the same 'size' as the inner product of the
          terms model matrix with itself? ...

In the source code of smoothCon, there is:

if (scale.penalty && length(sm$S) > 0 && is.null(sm$no.rescale)) {
    maXX <- norm(sm$X, type = "I")^2
    for (i in 1:length(sm$S)) {
        maS <- norm(sm$S[[i]])/maXX
        sm$S[[i]] <- sm$S[[i]]/maS
        sm$S.scale[i] <- maS
    }
}

Briefly speaking, for a model matrix X and raw penalty matrix S, scaling factor maS is:

norm(S) / norm(X, type = "I")^2
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • Maybe that makes a difference but I have something like this: `gam.out <-gam(y~s(x1, bs = "tp", k = 5) + (x2, bs = "tp", k = 5)+s(x1,x2, bs = "tp", k = 24)), data=data,select=TRUE, method="REML",family = poisson(link = "log"))` I guess 'select' is doing the 2X3 matrix. – M. Beausoleil Jul 29 '22 at 16:29
  • 1
    @M.Beausoleil Oh, with "select", each `s()` has two smoothing parameters. My code will give you a vector of 6 values. You can of course reshape it to a matrix, using `matrix(results, nrow = 2, ncol = 3)`. So each column contains 2 smoothing parameters for an `s()` term. This is a good point. I should mention this in the answer. Thanks. Also, you are welcome to vote it up, if it has helped you. :D – Zheyuan Li Jul 29 '22 at 16:33
  • What would be the meaning? In the end you'd have 6 different lambdas? – M. Beausoleil Jul 29 '22 at 16:42
  • 1
    @M.Beausoleil For each `s()` under `select = TRUE`, the first smoothing parameter is to penalize curvature, the second is to penalize linear trend. `select = FALSE` only has the first one, which at most shrink the function to a linear line or bilinear plane; `select = TRUE` adds the second one, which can further shrink the line or plane to 0. – Zheyuan Li Jul 29 '22 at 16:44
  • Also, in your example, is lambda computed like this: `scale.par = sapply(b$smooth, "[[", "S.scale") / b$species` and then `lambda = b$sp/scale.par` so the smallest lambda is 0.000008 and the largest is 2228963173967866624.00? (in your example) – M. Beausoleil Aug 02 '22 at 14:56
  • 1
    @M.Beausoleil What is `b$species`? But .... oops, I did something wrong there... I should reverse the numerator and denominator. So (1) the scaling parameter `alpha` is `alpha <- sapply(b$smooth, "[[", "S.scale")`; (2) the correct code for getting `lambda` from `sp` is: `b$sp / alpha` (not `alpha / b$sp` in the answer). It will take bigger effort to revise this answer for better clarity. I should probably add more examples as well. I will definitely improve this answer when I am available. Thank you!! – Zheyuan Li Aug 02 '22 at 19:19
  • Sorry, my computer corrected "sp" to "species"... That was my understanding also, that the denominator and numerator were reversed... also what is the reference for the image under the "lambda and sp:" section? I'm curious to read around this. – M. Beausoleil Aug 02 '22 at 19:33
  • 1
    @M.Beausoleil I wrote the equation in the image, to explain the scaling that **mgcv** does under the hood. It is not to be found anywhere else. This is a low-level technical detail, and is not mentioned in the book *GAMs: and introduction with R*. It is only mentioned in package documentation `?smoothCon`, as I quoted in the answer. The scaling is used to enhance numerical stability during computations. – Zheyuan Li Aug 02 '22 at 19:37