2
$\begingroup$

In this comment about smooth-factor interactions, Gavin Simpson wrote:

s(x1, by = cat) doesn't create a "fs" class smooth

I understand that

y ~ s(x, by = cat) + cat             (a)
//cat is an unordered factor

creates a separate smooth for each level of the factor as shown in summary() and plot() with a confidence interval around each smooth (and that this model doesn't say anything about whether the smooths for those levels are different).

On the other hand, in

y ~ s(x, cat, bs = 'fs')             (b)

the summary shows one smooth, but the plot() shows two smooths and with no confidence intervals.

So my questions are:

  1. How should the fs smooth in (b) summary() be interpreted? Does the p-value say anything about whether the smooths for cat levels are different?

  2. Why does plot() in (b) shows two curves without confidence intervals?

  3. What is the exact difference between (a) and (b)? When should I use one instead of the other? And how does their interpretation differ?

Many thanks

$\endgroup$
2
$\begingroup$

Q3

The by smooth variant has a separate smoothness parameter for each level of cat. In one sense you can think of these smooths as creating entirely separate smooth functions (they're aren't quite, they share the same basis for example), hence you have nlevels(cat) entries in the summary() output.

fs smooths share a single smoothness parameter, and in that sense are a bit closer to the smooth equivalent of random slopes; all the smooths are being pulled towards 0 (null function), with the smooths for some levels of cat pulled more towards 0 than others.

Q1

You should just treat the fs smooth summary as some measure of overall importance of the set of smooth functions. If the fs smooth is the only one in the model where x is a covariate then you could say interpret this as a test of the amount of variation in the response that can be explained by smooth functions of x, against a null of no effect. It's like an omnibus test in the sense that it is over the set of smooths.

Put another way, it is a bit like a test for a variance term of random slopes in a mixed effects model; this would be on a single variance parameter and you get a single test over the entire "effect" even though it represents $m$ slopes.

Q2

Typically, one would use fs where there where nlevels(cat) was large and you aren't particularly interested in the specific levels, just like you might work with random effects in a mixed effects model.

As there could feasibly be a separate line for the many levels of the factor but we think of this as a single "smooth" you get a single plot with each smooth superimposed to get some idea of the variation within and between individual functions. It would get very messy it if also showed confidence bands for each smooth.

In the by case, because you conditioned on the observed set of levels (sensu fixed effects in a mixed model) it makes sense to view these are being of individual interest, hence a separate plot per level of cat.

Also note that this is just the default plot to get a quick look at the fitted functions.

$\endgroup$
5
  • $\begingroup$ Many thanks, @Gavin! For Q1, you wrote It's [fs] like an omnibus test in the sense that it is over the set of smooths, I guess this means: p-value of fs is not about whether there is an actual interaction between x and cat (i.e., it doesn't say anything about whether the smooths for the factor levels are different in shape) (please correct me - I'm not a statistician). I tried exactly model (b) (note: no separate s(x)) on a dummy data where y depends on x but not on cat, and cat doesn't modify effect of x: the result was fs having a very low p-value. $\endgroup$
    – user291533
    Aug 13 '20 at 14:50
  • $\begingroup$ Here is the code for model (2) that I've tried: x = rnorm(1000, 1, .2); y = rnorm(1000, 1, .1) * x^3 * rnorm(1000, 1, .1) + x^-1; cat = rep(factor(c(1,2)), 500); plot(x, y, col = cat); m = mgcv::gam(y ~ s(x, cat, bs = 'fs'), method = 'REML'); summary(m) $\endgroup$
    – user291533
    Aug 13 '20 at 14:51
  • $\begingroup$ It seems if I want to compare x effect between cat levels, the model should be s(x) + s(x, by = cat) + cat where cat is an ordered factor. (?) $\endgroup$
    – user291533
    Aug 13 '20 at 14:52
  • $\begingroup$ Yes, that's right; it's not a test of different shape, rather a test of a null hypothesis where all the smooths are flat and all the random intercepts are zero. Re your model (2) example, that's expected because x has an effect so the test says the data provide evidence against a zero effect & zero random intercept everywhere null hypothesis. You can do the ordered by but that only compares to a reference level (in your case with two levels that's not an issue) and it can depend on how well you can estimate the reference smooth; flip the order of levels & estimates can change $\endgroup$ Aug 13 '20 at 19:46
  • $\begingroup$ You can use an ordinary factor by variable smooth for this too fromthebottomoftheheap.net/2017/10/10/difference-splines-i but it requires post processing. You could do this with the fs smooth also but you have to handle the random intercept bit of this basis $\endgroup$ Aug 13 '20 at 19:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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