# How can I keep the order of main effects in an interaction term, when one of the main effect is dropped?

I would like to check the Bayes Factor (BF) of a main effect (e.g. `zBase`), comparing the following two model: (1) the model which has that main effect (`count ~ zBase + Trt + zBase:Trt`) and (2) the model which does not contain the main effect (`count ~ Trt + zBase:Trt`).

``````library(brms)

## model 1
fit <- brm(
count ~
zBase +
Trt +
zBase:Trt +
(1 | patient),
data = epilepsy,
family = poisson(),
prior = c(
prior(normal(0, 1),    class = b, coef = zBase),
prior(normal(0, 1),    class = b, coef = Trt1),
prior(normal(0, 1),    class = b, coef = zBase:Trt1),
prior(normal(0, 1),    class = sd)
)#,
#backend = "cmdstanr"
)

## model 2
fit.no.M.ZBase <- brm(
count ~
#zBase +
Trt +
zBase:Trt +
(1 | patient),
data = epilepsy,
family = poisson(),
prior = c(
#prior(normal(0, 1),    class = b, coef = zBase),
prior(normal(0, 1),    class = b, coef = Trt1),
prior(normal(0, 1),    class = b, coef = zBase:Trt1),
prior(normal(0, 1),    class = sd)
)#,
#backend = "cmdstanr"
)
``````

Specifying the prior of each parameter, I am able to fit the first model but not the second. Although I overtly set the priors, the programme cannot find the prior specification of the interaction:

``````Error: The following priors do not correspond to any model parameter:
b_zBase:Trt1 ~ normal(0, 1)
Function 'get_prior' might be helpful to you.
``````

After I ran `get_prior`, it turns out that the order of the two main effects (`zBase` and `Trt`) in the interaction term of the second model flips from `zBase:Trt1` to `Trt1:zBase`, automatically and unexpectedly.

So, dropping one of the main effects unexpectedly alter the order of these main effects in an interaction term, as shown below, and interrupts my prior specification. Then, how can I keep the order of main effects in an interaction term, when one of the main effect is dropped?

``````## for model 2, see the interaction term of the formula and returned value
get_prior(
count ~
#zBase +
Trt +
zBase:Trt +
(1 | patient),
data = epilepsy,
family = poisson()
)

>                 prior     class       coef   group resp dpar nlpar bound       source
(flat)         b                                               default
(flat)         b Trt0:zBase                               (vectorized)
(flat)         b       Trt1                               (vectorized)
(flat)         b Trt1:zBase                               (vectorized)
student_t(3, 1.4, 2.5) Intercept                                               default
student_t(3, 0, 2.5)        sd                                               default
student_t(3, 0, 2.5)        sd            patient                       (vectorized)
student_t(3, 0, 2.5)        sd  Intercept patient                       (vectorized)
``````

Note but, the order of the two main effects in the returned value of `get_prior` is `zBase:Trt1`, same as I overtly specified in the formula, when the all main effects are specified in the formula.

``````## for model 1
get_prior(
count ~
zBase +
Trt +
zBase:Trt +
(1 | patient),
data = epilepsy,
family = poisson()
)

>                 prior     class       coef   group resp dpar nlpar bound       source
(flat)         b                                               default
(flat)         b       Trt1                               (vectorized)
(flat)         b      zBase                               (vectorized)
(flat)         b zBase:Trt1                               (vectorized)
student_t(3, 1.4, 2.5) Intercept                                               default
student_t(3, 0, 2.5)        sd                                               default
student_t(3, 0, 2.5)        sd            patient                       (vectorized)
student_t(3, 0, 2.5)        sd  Intercept patient                       (vectorized)
``````

## Operating System

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

## Interface Version

packageVersion(“brms”)
 ‘2.14.4’

Does anybody tell me whether the behaviour is by design or not?

Sorry for the slow response, I’m not super familiar with how this stuff works, but the second thing by itself seems weird enough to make an issue about. brms issues are here: https://github.com/paul-buerkner/brms/issues/new

As a guess, maybe try specifying the interaction in `coef` as a string? Like:

`prior(normal(0, 1), class = b, coef = "zBase:Trt1")`

Thank you for your reply. I’m not sure, too, whether this ‘problem’ is issue-worthy… The reason is that, actually, the computation is done when I manually flip the order of the main effects in the interaction term, though that manipulation itself is tedious and might be error-prone.

Unfortunately, string notation for the interaction term in `coef` does not work, as shown below:

``````library(brms)

## string notation works when all the main effects are specified
fit <- brm(
count ~
zBase +
Trt +
zBase:Trt +
(1 | patient),
data = epilepsy,
family = poisson(),
prior = c(
prior(normal(0, 1),    class = b, coef = "zBase"),
prior(normal(0, 1),    class = b, coef = "Trt1"),
prior(normal(0, 1),    class = b, coef = "zBase:Trt1"),
prior(normal(0, 1),    class = "sd")
)
)

## string notation does NOT work when some main effects are dropped
fit.no.M.ZBase <- brm(
count ~
#zBase +
Trt +
zBase:Trt +
(1 | patient),
data = epilepsy,
family = poisson(),
prior = c(
#prior(normal(0, 1),    class = b, coef = zBase),
prior(normal(0, 1),    class = b, coef = "Trt1"),
prior(normal(0, 1),    class = b, coef = "zBase:Trt1"),
prior(normal(0, 1),    class = sd)
)
)
## Error: The following priors do not correspond to any model parameter:
## b_zBase:Trt1 ~ normal(0, 1)
## Function 'get_prior' might be helpful to you.
``````

This is standard R formula behavior that I have no control over from the brms side. It is not bug but intended. I forgot the exact logic behind this unfortunately.

Thank you for your reply. Now I understand the behaviour is by design. Ummm, but this is counter-intuitive for me… If you remember the logic of R, please let me know… I want to learn why…

Oooo check this: r - Variable order in interaction terms - Stack Overflow

I think that means R arranges order of the interaction according to which bit of the interaction it encounters first.

So for `zBase + Trt + zBase:Trt`, `zBase` is the first thing it encounters so `zBase` goes first in the pair.

Try either putting `Trt` always first or the interaction itself always first and see if that gives you consistent ordering.

The example of `lm()` in the link you shared did not work for `brm()`. Both `Trt + zBase:Trt` and `zBase:Trt + Trt` in the formula consistently return the same error: `Function 'get_prior' might be helpful to you`.

So only the last resort works: to change the order of the main effect terms in the interaction (at least in the prior setting section) manually so that the remained main effect always comes first in the interaction term.

Also change the prior order too: (edit: missed word)

``````make_stancode(
count ~ Trt +
zBase +
zBase:Trt +
(1 | patient),
data = epilepsy,
family = poisson(),
prior = c(
prior(normal(0, 1),    class = b, coef = "zBase"),
prior(normal(0, 1),    class = b, coef = "Trt1"),
prior(normal(0, 1),    class = b, coef = "Trt1:zBase"),
prior(normal(0, 1),    class = "sd")
)
)

## string notation does NOT work when some main effects are dropped
make_stancode(
count ~ Trt +
zBase:Trt +
(1 | patient),
data = epilepsy,
family = poisson(),
prior = c(
#prior(normal(0, 1),    class = b, coef = zBase),
prior(normal(0, 1),    class = b, coef = "Trt1"),
prior(normal(0, 1),    class = b, coef = "Trt1:zBase"),
prior(normal(0, 1),    class = sd)
)
)
``````
1 Like

You could, of course, just make the interaction term(s) yourself. The interaction is literally just the variables multiplied together.