Determining derivative penalization for hierarchical GAM

Question

I am trying to fit 4 hierarchical GAMs to my dataset and compare them to determine the best-fit model. I have been relying on Pedersen et al 2019 to help me specify these models, but I don't think I fully grasp setting m for each smooth.

My dataset has the following variables: year, region, site, and mass. Sites are nested under region.

In the first two models, I'm interested in comparing the trend in mass across sites with the same smoothness/wiggliness and the with the same smoothness, but different wiggliness. In these models, the global trend is across all sites. These are basically taken exactly out of Pedersen et al 2019, so I'm pretty confident that they are correct.

# 1. a) Global trend across all sites, same wiggliness [equivalent to GS in Pederson et al. 2019]
mGS <- gam(mass ~ s(year, bs = "tp", k = 10, m = 2) +                       # global trend, m = 2
                  s(year, site, bs = "fs", k = 10, m = 2),                  # site trend, m = 2
            data = df, method = "ML", control = gam.control(nthreads = 4))  # used ML instead of REML for model comparison

# 1. b) Global trend across all sites, different wiggliness [equivalent to GI in Pederson et al. 2019]
mGI <- gam(mass ~ s(year, bs = "tp", k = 10) +                              # global trend
                  s(year, by = site, bs = "tp", k = 10, m = 1) +            # site trend, m = 1
                  s(site, bs = "re"),                                       # random effect for site
            data = df, method = "ML", control = gam.control(nthreads = 4))  

However, for the second set of models, instead of the global trend across all sites, I want to set a regional trend, and I'm not entirely sure I've set the m parameter correctly, especially in 2b) where m = 1 for both the regional and site smooths.

# 2. a) Regional trend for nested sites, same wiggliness
mRS <- gam(mass~ s(year, by = region, bs = "tp", k = 10, m = 1) +              # regional smooth, m = 1 similar to site trend in model 1b)
                 s(region, bs = "re") +                                        # random effect for region similar to 1b)
                 s(year, site, bs = "fs", k = 10, m = 2),                      # site trend, m = 2 similar to site trend in model 1a)
            data = means, method = "ML", control = gam.control(nthreads = 4))  

# 2. b) Regional trend for nested sites, different wiggliness
mRI <- gam(mass ~ s(year, by = region, bs = "tp", k = 10, m = 1) +             # regional smooth, m = 1 similar to site trend in model 1b)
                  s(region, bs = "re") +                                       # random effect for region similar to model 1b)
                  s(year, by = site, bs = "tp", k = 10, m = 1) +               # site trend, m = 1 similar to site trend in model 1b)
                  s(site, bs = "re"),                                          # random effect for site
            data = means, method = "ML", control = gam.control(nthreads = 4))  

Answer 1

The reason for the m = 1 (first derivative penalty) on the factor by smooths is to make the factor by smooths identifiable from what we called the "global" smooth. If you don't have a global smooth then you really don't need m = 1 and are best to use the default second order derivative penalty.

To be clear, when I'm talking about a global smooth I mean the term s(x) in the following model:

y ~ f + s(x) + s(x, by = f, m = 1)

In both of your regional models you don't have the equivalent of s(x) so I would suggest that m = 2 be used otherwise you'll be fitting functions that aren't visually smooth.

Since we wrote the HGAM paper, Simon Wood has added a constrained factor smooth basis that provides a better way to fit a HGAM with a global smooth plus (the equivalent of) by factor smooths, the bs = "sz" basis. This basis is made properly orthogonal to the global smooth, so we don't need to worry about the distinctions of what m to use.

Such a model would be fitted as

y ~ s(x) + s(x, f, bs = "sz")

noting that we don't need the parametric (or the random intercept) for factor f.

So your model 1b could be fitted as

mGI <- gam(mass ~ s(year, bs = "tp", k = 10) +
                  s(year, site, bs = "sz", k = 10),
            data = df, method = "ML", control = gam.control(nthreads = 4))

And your model 2b might run into issues because you have multiple smooths of year, which might not be so identifiable. Instead, you could use the "sz" basis to fit a global smooth and then properly orthogonal "difference" smooths for each site and region:

mRI <- gam(mass ~ s(year) +
                  s(year, region, bs = "sz", k = 10) +
                  s(year, site, bs = "sz", k = 10),
            data = means, method = "ML", control = gam.control(nthreads = 4))

While we described these models in terms of

1. global + subject specific smooth, or
2. subject specific smooths with no common smooth

each model is trying to estimate the smooth effect of year on the response for each level of a factor. The different models decompose the variation in the response into different sets of smooth functions. Option 1. can be efficient if there is common shape to the subject-specific effects of year and if that is part of the research question. But you can use either model to fit your data and even if you don't think there is a global effect, you can use that model in the 2b setting if it helps solve computational issues with the ideal or prefered approach.