How does random variable nesting in GAMs work (mgcv)?

Question

Consider me very new to the world of GAMs, I am actually using it out of necessity rather than by choice but I thought it could be a chance to learn something new anyway.

Originally I had hoped to model my data with a GLMM, which feels fairly natural to me in the way you assign nested random variables. In the case of my study, using lme4 I would write (1|SURVEYGR/YEAR/PlotID). This implies PlotID is nested in YEAR which is nested in SURVEYGR. But when it comes to the mgcv package I haven't managed to find a clear explanation for the equivalent expression. I have seen it is possible to use the bs=re (re=random effect) argument for example s(YEAR, bs="re") while in other cases such as this blog I have seen a number of random variables together as such s(SURVEYGR,YEAR,PlotID, bs="re") but there is a lot going on and it is hard to tease out what this code implies. The problem is, I am not sure how nested random variables should be structured

So ultimately I would like to know what the mgcv code equivalent for a GAM model would be for the following GLMM model written for lme4:

glmer.nb(count ~ vegetation.cover + temperature + 
                (1|SURVEYGR/YEAR/PlotID)
    , na.action = "na.fail"
    , data = dd, verbose=T)

At a guess it might be something like

gam(count ~ s(vegetation.cover) + s(temperature) + 
            s(SURV.GR,YEAR,PlotID, bs="re") 
           , family = nb()
           , data=dd, method="REML")

EDIT

a bit of snooping around turned up this post on SO where translating the comments to my question suggest including the terms individually and then the interaction for nested random effect terms like so:

+ s(SURV.GR, bs="re") + s(SURV.GR,YEAR, bs="re") + s(SURV.GR,YEAR,PlotID, bs="re")

Any thoughts?

Answer 1

There are several ways to specify this with mgcv, each yielding a different way to estimate the smoothing parameter and variance component. For the predicted random intercepts: these will not be identical between different specifications, but I do not expect substantial differences. If you are interested in interpretation of the variance components, the difference between specifications may be quite important.

With s(SURV.GR, bs="re") + s(SURV.GR, YEAR, bs="re") + s(SURV.GR, YEAR, PlotID, bs="re"), three separate smoothing parameters and variance components will be estimated, one for each smooth term specified.

But perhaps you want to allow the amount of smoothing to vary between different levels of YEAR. Or you might be interested in differences of the variance explained by SURV.GR, for each of the different levels of YEAR.

Then you could replace s(SURV.GR, YEAR, bs="re") with s(SURV.GR, by=YEAR, bs="re"); this would yield a separate smoothing parameter and variance component being estimated for the random intercept w.r.t. SURV.GR, for each level of YEAR.

I would recommend after fitting your model(s), to inspect gamfit$sp and varcomp.gam(gamfit) to see whether the smoothing parameters and variance components are estimated as you expected. Also, you may want to compare against VarCorr(glmerfit). (where gamfit refers to the object returned by gam() and glmerfit refers to the object returned by glmer())