I am currently learning how to use GAMs and try to model a species response (resp) to different environmental parameters (x1, x2, x3) that were sampled at different locations (random effect). My understanding is that I should allow each predictor to have its own smoother per species level (by=species) but I struggle to understand what this would do in the random effect term and how to best specify it. The example model is of the form:

gam(resp ~ s(x1, by=species) + s(x2, by=species) + 
           s(x3, by=species) + s(location, bs="re")

and I try to understand what is happening if I include either of these three terms:

(1) s(location, bs="re")
(2) s(location, by="species", bs="re")
(3) s(location, species, bs="re")

1 Answer 1


Your base model is incorrectly specified; factor by smooths must have the by factor included as a parametric categorical term in the model, hence you need:

gam(resp ~ species +
      s(x1, by = species) +
      s(x2, by = species) + 
      s(x3, by = species) +
      s(location, bs = "re")

This allows for the mean of the response in each level of species to be included in the model as a specific term, while the smooths by species, which are centred, are used to model smooth variation about these species specific mean responses.

That said...

s(location, bs = "re")

s(location, bs = "re") would indicate an i.i.d. Guassian random intercept term in the model. In other words we get an intercept for each level of the factor location (and it must be a factor to get this interpretation), and those intercepts are shrunk towards 0 as part of the penalized likelihood used in fitting. This is mathematically the same as including (1 | location) in a mixed effects model using lme4 notation for random effects.

s(location, by = species, bs = "re")

(Note that it is by = species, not "species".) This is the same as the above, except the random effect is nested in the levels of the factor species. What this means practically is that the amount of shrinkage over all the locations can vary between species. This is because each smooth generated by by = species gets its own smoothness parameter, $\lambda$, and it is this that controls the degree of shrinkage (penalisation) that each "smooth" is subject to. We're just using the equivalence of smooths and random effects to fit the equivalence of i.i.d. Gaussian random effects in these GAMs.

s(location, species, bs = "re")

Assuming that both location and species are factors, then you would get a random intercept for the groups formed by the interaction of location and species. Consider the following two factors, each with two levels:

levs <- c("A", "B")
df <- data.frame(f1 = factor(sample(levs, 50, replace = TRUE)),
                 f2 = factor(sample(levs, 50, replace = TRUE)))
with(df, table(f1, f2))

The last bit shows the two-way table

f1   A  B
  A 17 11
  B 11 11

So we have the four combinations:

  1. "A:A",
  2. "A:B",
  3. "B:A", and
  4. "B:B"

and we can think about the combination of these two factors creating 4 groups. The term s(f1, f2, bs = "re") would then yield a random intercept for these four groups.

Technically, it results in the addition of ~ f1:f2 - 1 columns to the model matrix:

r$> head(model.matrix(~ f1:f2 - 1, data = df))                                  
  f1A:f2A f1B:f2A f1A:f2B f1B:f2B
1       0       0       1       0
2       0       0       1       0
3       0       0       0       1
4       0       1       0       0
5       0       0       1       0
6       0       0       0       1

which we can see as binary indicators for the groups defined by the interaction of the two factors. These columns are associated with a diagonal penalty matrix (which is what achieves the shrinkage): this is a ridge penalty on the coefficients associated with these columns of the model matrix, and which yields i.i.d. Gaussian random effects.

Continuous variables in s(...., bs = "re")

Adding a continuous variable to a random effect spline will get you a "random slope" term. These always go with at least one factor, so we might have

~ s(temperature, group, bs = "re")

with temperature a continuous variable and group a factor (explicitly coded as a factor in R!). Thus we will get terms ~ temperature:group - 1 added to the model matrix of the model, which we will see as a random effect of temperature within each level of the factor group. Again, these columns are associated with a diagonal(ridge) penalty, leading to i.i.d. Gaussian random effects of temperature.

N <- 10
df <- data.frame(group = factor(sample(c("A", "B"), N, replace = TRUE)),
                 temperature = runif(N, 2, 25))
model.matrix( ~ temperature:group - 1, data = df)

which results in

   temperature:groupA temperature:groupB
1            20.24348            0.00000
2            15.38841            0.00000
3             0.00000           11.35646
4            18.29931            0.00000
5             0.00000           11.15989
6             0.00000            4.20421
7            23.05284            0.00000
8             0.00000           13.64325
9             0.00000           16.17145
10            0.00000           18.00376

and hence we get a separate coefficient for each column, and hence a "random slope" for each group, with the diagonal (ridge) penalty shrinking the coefficients for these columns towards zero, with the amount of shrinkage determined during fitting.

Note that with the formulation in {mgcv}, we cannot achieve the same as (1 + temperature | group), which would imply correlated random effects, where the random intercepts for group and correlated with the random slope of temperature for each group. In {mgcv}, the best we can do is

~ s(group, bs = "re") + s(temperature, group, bs = "re")

which yields uncorrelated random effects, which in {lme4} syntax would be

~ (1 | group) + (0 + temperature | group)

and possibly some other shorthand (IIRC?: (temperature || group)).

  • 1
    $\begingroup$ That is a fantastic answer! Thank you very much. May I ask one more follow-up question: If no a priori knowledge exists about the random structure in a GAM (whether it should just have a random intercept or a random slope), is it safe to select the random structure by AICc score? $\endgroup$ Commented Nov 24, 2021 at 12:32
  • $\begingroup$ @Globoquadrina You shouldn't use AICc; AIC for GAMs is problematic and the AIC computed for GAMs in {mgcv} via function AIC() includes some correction for smoothness selection and is formulated in a way that is appropriate for smooths. Unless the AICc implementation you are using has those features (and I'm not aware of any) you might get mislead by it. I also come from the school of thought that doesn't like this kind of model selection. If I had reason to include a term in the model it should stay in the model & I report the estimated effects accordingly, even if ~0 $\endgroup$ Commented Nov 24, 2021 at 12:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.