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

CLRR · January 6, 2021, 6:53am

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”)
[1] ‘2.14.4’

CLRR · January 7, 2021, 11:04am

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

bbbales2 · January 9, 2021, 11:03pm

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")

CLRR · January 12, 2021, 4:08am

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.

paul.buerkner · January 12, 2021, 7:16am

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.

CLRR · January 12, 2021, 12:46pm

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…

bbbales2 · January 12, 2021, 2:29pm

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.

CLRR · January 13, 2021, 1:22pm

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.

bbbales2 · January 13, 2021, 1:38pm

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)
  )
)

Stephen_Martin · January 13, 2021, 8:03pm

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

Topic		Replies	Views
Approach to use bayes_factor in brms Modeling bayes-factor	2	2038	June 14, 2020
Bayes factor favours reduced model without the predictor of interest: how do to explain? Modeling hierarchical-model , bayes-factor , brms	1	754	March 18, 2022
Find the main effect in a bayesian model Modeling	7	624	December 22, 2022
Bayes factors for a moderating effect brms rstan	11	1182	July 5, 2020
Main effect in a bayesian model Modeling	1	466	September 4, 2022

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

Operating System

Interface Version

Related topics