How to increase hierarchical GAMM complexity?

Question

I'm interested in building a generalized additive mixed model (GAMM) like this as it has nearly the same set-up as mine, but with an additional level of complexity (a 2 level factor called "season", per year). V1-2 are my (time dependent) environmental covariates. I think it makes sense biologically to have a global year smoother with individual site effects that share the same penalty (i.e. s(CYR.std, fSite, bs = 'fs')), but how would I incorporate it with the my random intercept for site? I'm limited by the number of terms I can do as I don't have a lot of data. Therefore, I'm hoping to make this as efficient as possible (Total N=1381).

library(mgcv)

# Is this close to what I need?
m <- gam(count ~ s(V1) + V2 + 
           s(CYR.std, by = fSeason) +
           s(CYR.std, fSite, bs = 'fs'),
         family=poisson, data=df, method = "REML", select = TRUE)

• fSite = factor site (repeated measures design - same 47 sites sampled once per season, every year)
• CYR.std = Continuous year (2008 =1, 2009 =2, etc..)
• fSeason = factor season (2 levels, Wet/Dry)

UPDATE: Does this seem correct? I used the bs = "sz" because I have multiple smooths of year.

# switched to bam() for speed

test <- bam(count ~ s(v1) +
              v2 +
              s(CYR.std) +
             s(CYR.std, fSeason, bs = "sz") +
             s(CYR.std, fSite, bs = "sz"),
            data = toad2, 
            method = 'fREML',
            discrete = TRUE,
            family = poisson,
            select = TRUE,
            control = list(trace = TRUE))

My first try (what I was told these mean):

    mod <- bam(count ~ CYR.std * fSeason + # Easier to explain
               s(v1) + v2 + 
    
               s(fSite, bs = "re") + 
    # No repeated measures per fSite, BUT a random intercept says: each site 
    # can have it's own starting position. If not random, then we are explicitly 
    # saying all sites start with the same abundance. It makes more sense 
    # for abundance to have a random starting point in this case).
                 
               s(fSite, fCYR, bs = "re") + 
    # The station:year interaction captures the correlation (repeated measure) 
    # between the 2 measurements per year, at each site. Captures the broader 
    # level variation among sites and years. 
    
               s(fSite, CYR.std, bs = "re") +
    # Each site can have it's own trend in abundance, over time
                
               s(fSeason, fCYR, bs = "re") + 
    # fSeason within year
               
               offset(log(area_sampled)), 
             data = toad2, 
             method = 'fREML',
             discrete = TRUE,
             family = poisson,
             control = list(trace = TRUE))

Answer 1

The "fs" basis has a fully penalized null space and contains its own "random" intercept and linear terms. As such, you don't need a separate random intercept term when using this basis.

The model you suggest:

m <- gam(count ~ s(V1) + V2 + 
           s(CYR.std, by = fSeason) +
           s(CYR.std, fSite, bs = 'fs'),
         family=poisson, data=df, method = "REML", select = TRUE)

Will include a separate smooth of time per level of fSeason, with potentially different wigglinesses, a smooth of time for each level fsite all sharing the same wiggliness.

You model is incorrectly specified as it doesn't contain the mean for each level of fseason - with factor by smooths you have to include the by factor as a separate parametric effect (or random effect) to model the group means.