# How to increase hierarchical GAMM complexity?

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"),
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)),
method = 'fREML',
discrete = TRUE,
family = poisson,
control = list(trace = TRUE))


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.

• Thank you! Quick question, is my random effect smoother "s(fSeason, fCYR, bs = "re")", where both terms are transformed into factors, only used when the fixed effect "s(CYR.std, by = fSeason) + fSeason" isn't needed? I tried to follow the example in "Hierarchical generalized additive models in ecology: an introduction with mgcv", but I had trouble understanding when to use the term like s(taxon, year_f, bs="re"), in the zooplankton example.
– Nate
Oct 29, 2023 at 4:12
• You can add the group means however you want, either though parametric term or via random effects. whether you need f1 * f2 or f1 + f2 (or the ranef equivalents) will depend on how you want to delineate group means. Interactions would be used if you had a smooth for each level of interaction(fSeason, fCYR), but you didn't say you needed anything like that. I really don't understand why you have fCYR when you are also treating the year effect as smooth. (Using both is fine, but they imply different evolutions of time.) Oct 29, 2023 at 11:09
• Aha, thank you! I wasn't sure what I needed. I didn't understand what the syntax meant or that I didn't need both a fixed and random effect for year (I saw "taxon" in both random and fixed form in the paper). There are gaps (100+ days) between our sample events, so a "step-wise" evolution of time makes sense, but would require too many parameter estimates (I'm happy with just one general trend of time per season). Using both forms of year was my advisors suggestion (that the parametric continuous and random factor version did different things and one needed both; reasoning in R comments above)
– Nate
Oct 29, 2023 at 15:58
• Example from paper (pg. 22) -> Model S: zoo_comm_modS <- gam(density_adj ∼ s(taxon, year_f, bs="re") + s(day, taxon, bs="fs", k=10, xt=list(bs="cc")), data=zoo_train, knots=list(day=c(0, 365)), family=Gamma(link="log"), method="REML", drop.unused.levels=FALSE)
– Nate
Oct 29, 2023 at 17:19