I have a question about mgcv's formula syntax for interactions with smooths s() (I am actually using this syntax within brms).

Say I have two interacting factors X1 and X2 and I want to estimate how this interaction is moderated by a smoothed term X3 (a "varying coefficients" model). According to the mgcv help, when the by variable is a factor, it should be included as a main effect, but I am not sure how that applies when the interaction term is also a by variable. So, would I write that as

Y ~ X1 * X2 + s(X3, by=X1) + s(X3, by=X2) + s(X3, by=interaction(X1,X2)) 

or as

Y ~ X1 * X2 + s(X3, by=interaction(X1,X2)) 

Now, to get to the heart of my question, let us suppose that X2 is a continuous predictor, not a factor (according to the mgcv help, when the by variable is a continuous predictor, it should not be included as a main effect). So, which would be the correct syntax

Y ~ X1 + s(X3, by=X1) + s(X3, by=X2) + s(X3, by=interaction(X1,X2))


Y - X1 + s(X3, by=interaction(X1,X2))

or are both wrong?


For the two factors, I would use the second form

df <- transform(df, X1X2 = interaction(X1, X2, drop = TRUE))
Y ~ X1 * X2 + s(X3, by = X1X2)

(where df is the object passed to data.) You need to be careful with interaction() creating all possible combinations rather than the combinations in the data.

I might even go so far as doing

df <- transform(df, X1X2 = factor(paste(X1, X2, sep = '_')))
Y ~ X1 * X2 + s(X3, by = X1X2)

The terms produced by the intercept, X1, and X2 will code for the means of each group, while X1X2 should give a factor also coding for each group such that we get a smooth for each combination of X1 and X2. It doesn't really matter how the group means get code, just so long as we have a unique factor for each combination of them and have this coded in the dataset. The smooths and the factor fixed effects in the formula are somewhat independent in how we represent the effects in the model; so long as we get a smooth per unique combination and we represent the group means in some way we are good.

I'm not sure it is strictly possible to do what you want with a continuous X2. This came up recently on Twitter and there Eric Pedersen suggested

Y ~ X1 + te(X3, X2, bs = c('cr', 're'), by = X1)

which would give a tensor product of a cubic regression spline in X3 and a random slope (i.e. a linear effect) of X2.

In the Twitter conversation where Eric suggested that solution, the only other way we could think to do it would be by writing a new smooth.construct function that would simply add an unpenalized, linear term to the model matrix. But no one went so far as to implement that as yet.

  • $\begingroup$ Thank you. I am unsure what the "coefficients" would mean in Y ~ X1 * X2 + s(X3, by = X1X2). Suppose X1 and X2 each have 2 levels and are dummy coded; for the following 4 possible cases of X1 and X2, are the these predictions correctly specified: for $X1=0$, $X2=0$, predictions are $b_0 + f(X3)$; for $X1=1$, $X2=0$, predictions are $b_0 + b_{X1} + f(X3)$; for $X1=0$, $X2=1$, predictions are $b_0 + b_{X2} + f(X3)$; for $X1=X2=1$, predictions are $b_0 + b_{X1} + b_{X2} + b_{X1X2} + f(X3)$? So, in each case, it is the same X3 "curve", just differently offset? $\endgroup$ May 11 '20 at 5:22
  • 1
    $\begingroup$ No, the model now has $f_j(X3)$ where $j$ indexes the levels of X1X2. Each of the $j$ smooths has it's own smoothness parameter, so the smooths can have different wiggliness. Basically the model matrix gets $j$ x k new columns, whose values are all 0 except for the rows related to the $j$th smooth, ie whose combination of X1 and X2 is that of the $j$th group. The parametric terms for X1 and X2 do shift the curve up and down; the $j$ smooths are centred so we need these terms to get the group means. But this model also specifies $j$ separate smooths. $\endgroup$ May 11 '20 at 14:31
  • $\begingroup$ Multiple separate smooths -- excellent, now I understand (I was mentally multiplying instead of concatenating and so incorrectly thinking X1X2 had just two levels)! In the last formula Y ~ X1 + te(X3, X2, bs = c('cr', 're'), by = X1), I don't yet understand how the 're' argument turns X2 into a slope (a penalized linear effect) but I think having these slopes penalized makes sense in my current use case, so this is a good solution. For X3, is there any reason to prefer 'cr' over, for example, 'tp'? $\endgroup$ May 12 '20 at 1:18
  • 1
    $\begingroup$ For the random effect basis bs = 're', the mgcv-convention is to add random intercept for a factor variable or a random (linear) slope for a continuous variable. However, as you are only passing the one continuous variable you simply get a shrunk linear coefficient. When you do a tensor product of this with a smooth function, the effect is to have the shrunk linear coefficient for X2 vary as a smooth function of X3. A cubic regression spline basis is the default in te() so I just used that, because you have to provide something for the first marginal smooth. $\endgroup$ May 12 '20 at 2:02

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.