I'm wondering how to interpret the bs=fs and by= terms and how two versions of a similarly coded model might differ (conceptually) in a hierarchical generalized additive model.

For example, could s(CYR.std) + s(CYR.std, fSite, bs = "fs") be used if "survey sites are located in the same habitat type (the same sandy beach) and in a similar region of the world, therefore they probably share a global (annual) trend with some minor site specific deviations from that shared trend"? Whereas s(CYR.std, by = fSite) + fSite is more like "survey sites are in different habitats and different parts of the world, so we should model this as having no common global trend in time (s(CYR.std)), because sites are so far away and probably don't have the same general effect on the response (?)

So, is the by= term similar to a treatment level effect (i.e. each is level has an independent effect from the other) and the fs term is more akin to a deviation from an overall trend? If this seems reasonable, what similarly stated real-world scenario might be analogous to both terms being present; like this s(CYR.std) + fSeason + s(CYR.std, fSite, bs = "fs", by = fSeason)?

# Species counts
df <- as.data.frame(rnbinom(1000, mu = 0.6971, size = 1))

df$year <- rep(2011:2020, each=100)
df$CYR.std <- df$year - min(df$year)
df$fCYR <- as.factor(df$year)

df$site <- seq(1, 50, 1)
df$fSite <- as.factor(df$site)

df$season <- rep(c("DRY", "WET"), each=50)
df$fSeason <- as.factor(df$season)

# Salinity (continuous covariate)
df$sal <- sample(0.5:40, 1000, replace = TRUE)

# Area surveyed
df$area <- 3

names(df)[1] <- "count"

What are the conceptual differences between these two models? Are they improperly coded?

m <- bam(count ~ s(CYR.std) + # overall trend
           fSeason + # seasonal means fixed effect
           s(CYR.std, fSite, bs = "fs", by = fSeason) + 
           # Deviations form common annual and seasonal trend for each site..?
           s(fCYR, bs = "re"), # year-to-year effects (average abundance higher for some years)
         offset(log(area)), # 3m^2 area
         method = "fREML",
         discrete = TRUE,
         family = nb(link = "log"),
         data = df)

m2 <- bam(count ~ s(CYR.std) + 
           fSeason +
           s(CYR.std, fSite, bs = "fs") + 
           s(CYR.std, by = fSeason) + 
           s(fCYR, bs = "re"), 
         discrete = TRUE,
         family = nb(link = "log"),
         data = df)
  • $\begingroup$ Have you read the documentation for mgcv::factor.smooth? $\endgroup$
    – PBulls
    Commented Jan 22 at 5:48
  • $\begingroup$ @PBulls the documentation was difficult for me to understand, but according to my reading of this answer (stats.stackexchange.com/questions/482509/…), when the two terms are separate, it means one thinks the effect of the random slopes (the bs=fs factor) doesn't differ between the levels of the by=fac term, I think. As in, sites can differ slightly in their annual trend, but annual trends by site don't differ by season, only the global trend differs by season (?) $\endgroup$
    – Nate
    Commented Jan 22 at 13:08
  • $\begingroup$ Model m could be something like figure 10 in Petersen et al. 2019 (Hierarchical generalized additive models in ecology), but with 2 versions of the random slopes panel? One for each season? $\endgroup$
    – Nate
    Commented Jan 22 at 14:51

1 Answer 1


The discussion below expands upon that in this previous Q&A.

The main difference between bs = "fs" and factor by smooths is in the number of smoothing parameters used for the set of smooths generated by the respective invocations in a formula passed to mgcv::gam() or mgcv::bam(). There are some more mechanistic differences relating to whether factor level means are included in the basis or not.

The "fs" basis

This basis, which we can add to a GAM formula via

s(x, f, bs = "fs", k = 10)

for numeric x and factor f, will result in a smooth function for level of factor f using k basis functions. This smooth is fully penalized, meaning that it has no penalty null space (no functions for which the penalty has no effect). This means that we don't need to remove the constant function from the span of the basis (and hence why I say there will be 10 basis functions in the above example, each smooth getting it's own 10 coefficients (weights) for those basis functions), so we include the group (levels of f) means. The linear function is also penalized (in normal smooths with the default penalty — or a higher order degree penalty — the linear function would be unpenalized as it's second derivative, curvature, is 0 everywhere). These extra penalties mean we have the usual wiggliness penalty, the penalty on the group means, and the one on the linear functions. The latter two penalties are in the form of identity matrices, so the group means and linear functions act like random intercepts and slopes respectively.

So, we have three smoothing parameters to select while fitting with this specific definition s(x, f, bs = "fs", k = 10). The upshot of this is that all the smooths created in this basis share the same smoothing parameter and hence have roughly the same wiggliness. In this sense then, the smooths are sharing information about the wiggliness, in the same way that random intercepts and slops share information as they can be thought of as being draws from a Gaussian distribution centred on 0 with some standard deviation (which is inversely proportional to the smoothing parameter for the group means or the group slopes [linear functions]).

The key thing here is that the family of smooths created by this basis share wiggliness information and hence we think of them more like random smooths in the terminology of random intercepts and random slopes from mixed models.

The factor by smooth

The equivalent factor by smooth to the "fs" one I showed above is

s(x, by = f, k = 10, m = 2)

These are only equivalent in the sense that both produce a smooth for each level of f. The key differences are

  1. this smooth definition doesn't include the group means. The sum-to-zero constraints are applied to each smooth so they are centred about 0 like most smooths in {mgcv}, and we need include something else in the model for the group means,
  2. each smooth has a penalty null space of dimension 1 (for the default m = 2) which means the linear function is in the span of the basis but it is unpenalized, it's can't be shrunk back to a zero function), and
  3. the family of smooths created by this definition all have separate smoothing parameters, meaning that one of the smooths could be really wiggly and one could a really smooth, and others anywhere in between.

In practice then, we need to include this factor by smooth in the model with an associated parametric factor term

f + s(x, by = f, k = 10, m = 2)

where the f term outside s() will model the group means. A s(f, bs = "re") could also be used here in place of the parametric f term.

Given these definitions then, the distinctions when considering their use in a model are: do you expect the family of smooths to have similar wiggliness or do you want to allow for potentially quite different wigglinesses across the family?

More philosophically, one might view the "fs" basis as being the random effect equivalent random intercepts and random slopes but applied to smooth functions, and factor by smooths as being the equivalent of fixed effects but applied to smooths. This philosophy gets confusing quickly though, as, when fitting using REML or ML smoothness selection all penalized splines have a random effect mirror image.

Because the factor by smooths are not fully penalized, there are some more subtle differences when we come to fit models discussed in Pedersen et al (2019) that we coined Hierarchical GAMs (HGAMs). There we discussed models that might have a average or global smooth effect plus group specific effects.

In our terminology then,

s(x, f, bs = "fs", k = 10)

is model S, while this

f + s(x, by = f, k = 10, m = 2)

is model I for Shared and Individual smoothness respectively. And then for the average or Global versions, which have a common smooth plus group-specific smooths we have model GS:

s(x) + s(x, f, bs = "fs", k = 10)

and model GI

f + s(x) + s(x, by = f, k = 10, m = 1)

where for the GI model we've had to shift the penalty to the integrated first derivative of the smooth to make the effects of the two types of smooth identifiable. The m=1 here also means that the factor by smooth also has 0 penalty null space once the constant term is removed by the sum to zero constraint. This is because a linear function has a constant first derivative that is not 0 (unless it has been shrunk to a $f(x_i) = 0 \; \forall \; i$ (i.e. a zero function). Importantly, the group means are still not included in this basis so we still need to model them explicitly.

A better way to fit the GI model now (for models where the factor doesn't appear in smooths of different numeric covariates) is the "sz" basis (sz standing for sum to zero constraint, which is how it does what it does):

s(x) + s(x, f, bs = "sz", k = 10)

where now the group means are now in the basis (so me don't need to model them explicitly via a parametric f term), but the linear function is in the basis and hence unpenalized (the group means are unpenalized too IIRC). By construction, this "sz" basis produces models that are identifiable without changing the order of the penalty on the group-level smooths. As with the factor by, the smooths for each level of f have their own smoothing parameter, so wigglinesses can be different across the family of smooths.

Now that we have cleared all that up, here's something that will confuse everything again: the family of smooths created by the factor by smooth and the "sz" basis can actually share a single smoothing parameter if you wish; this is done through the id mechanism:

# Model GS via "fs" basis
s(x) + s(x, f, bs = "fs", k = 10)

# Model GS via factor by with same smoothing parameter
f + s(x) + s(x, by = f, k = 10, m = 1, id = 1)

# Model GS via the "sz" basis
s(x) + s(x, f, bs = "sz", k = 10, id = 1)

At this point then, the difference is really just

  1. are the group means in the basis for the smooths? And if they are, are the group means penalized?
  2. are the linear functions penalized?
  3. what is the order of the penalty on the group-specific smooths?
  4. group-level smooths share a smoothing parameter

For the `"fs" basis we answer

  1. YES (YES),
  2. YES,
  3. 2nd order derivative penalty
  4. YES

For the factor by we answer

  1. NO (NO if added via parametric term f, YES if you added them via s(f, bs = "re"))
  2. YES (with m = 1)
  3. 1st order derivative penalty
  4. NO (but with id=1 YES)

For the "sz" basis we answer

  1. YES (?? I'm not sure)
  2. NO
  3. 2nd order derivative penalty
  4. NO (but with id=1 YES)

Your choice really boils down to how you set up an identifiable model (the m=1 we have to add to the factor by version can make smooths not be visually smooth as the factor by smooth is piece-wise linear, so there is a consequence of using it) and then whether you want the group-specific smooths to share a smoothing parameter (have similar wiggliness) or each to have it's own smoothing parameter (wiggliness is expected to be different across the groups).

That's really the choice being made here; you can coax the different definitions into fitting the same model (ish) with shared wigglinesses, but the factor by and the "sz" basis do allow for different group-level wigglinesses. Note however that all of these are parameterized differently internally so don't expect each of these models to return the same fitted values for covariate x and f:

# Model GS via "fs" basis
s(x) + s(x, f, bs = "fs", k = 10)

# Model GS via factor by with same smoothing parameter
f + s(x) + s(x, by = f, k = 10, m = 1, id = 1)

# Model GS via the "sz" basis
s(x) + s(x, f, bs = "sz", k = 10, id = 1)

How you, the analyst, interprets these disparate model definitions is not a statistical question. You can say that the "fs" basis should be used when all the sites come from the same habitat (where f codes for plots all on the same habitat) but the "sz" basis should be used when the habitats differ (where f codes for plots that are on different habitats). If the effect of x is similarly smooth over x across the different habitats, the "fs" will capture this well, but it won't if the effect of x is wiggly over x for some habitats but not others.

The difference between models with global or average smooths and those without is also something to consider but don't conflate it with the distinctions between the bases / smooth definitions! You can fit model S using the factor by or "fs" basis while Model I can be fitted by the factor by smooth. Whether you go for models with or without global smooths depends on

whether you want to know if there is some common pattern to the smooth effects of x over factor levels f and to estimate it?

It really is as simple as that.

If you want to get philosophical about it (but beware! That way, there be dragons!) we could think of the "fs" basis as being closer in spirit to random effects and the "sz" and factor by models being closer to fixed effects. But really the differences are more mechanistic than that and choosing between them is really down to how much wiggliness is expected to vary across the group-level smooths.

Now, if you aren't interested in the common pattern it can still be efficient to use a GS or GI model if there is common pattern. It is more efficient (in terms of model complexity) to fit one of these GS or GI:

# GS
s(x) + s(x, f, bs = "fs", k = 10)
f + s(x) + s(x, by = f, k = 10, m = 1, id = 1)
s(x) + s(x, f, bs = "sz", k = 10 id = 1)

# GI
f + s(x) + s(x, by = f, k = 10, m = 1)
s(x) + s(x, f, bs = "sz", k = 10)

if there is common pattern than to fit

# S
s(x, f, bs = "fs", k = 10)
f + s(x, by = f, k = 10, m = 1, id = 1)

# I
f + s(x, by = f, k = 10, m = 2)

because in the G versions, we can spend just the 10 coefficients estimating the common pattern via term s(x) and then group-specific smooths probably won't need to be a complex (wiggly) as they are just fine-tuning the common pattern to fit each groups data. As such many of the coefficients can be shrunk back towards 0 and we end up with lower overall EDF for the model. In the S and I models, however, we might need to use most of the 10 coefficients per group just to model the same pattern in each group. This is wasteful an inefficient use of the data.

If you're having a hard time visualising this, consider insect abundances measured throughout the year, which will vary seasonally. Some species come out earlier in the year than others, but they all share roughly the same shape of phenological curve — the peak varies and the shapes are a little different, but they all increase slowly at first, then quickly to a peak, then decline back to zero again. It might be efficient to model this common shape with a GS or GI model. But it the peaks vary too much over the season and the shapes of the phenological curves vary too much across the species, we might not gain that much using the GS or GI models and the S or I models would be easier to specify and use fewer smoothing parameters and have a smaller model matrix (so be more efficient to fit from a computational point of view).

As for your model m and m2, the main difference is in

s(CYR.std, fSite, bs = "fs", by = fSeason)


s(CYR.std, fSite, bs = "fs") + 
s(CYR.std, by = fSeason)

and (apart from the missing m=1 on the second version's factor by smooth) the main conceptual difference is that the first one says tha

all sites from the same season share the same wiggliness (smoothing parameter) in the between-year trend,

while the second form says

there is some between-year common trend to all sites in a particular season, but all sites share roughly the same site-specific trend wiggliness, only differing in their wiggliness through the common trend for each season.

  • 2
    $\begingroup$ Awesome, thank you! $\endgroup$
    – Nate
    Commented Jan 23 at 13:12

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.