A group-level intercept model using a multivariate brms formula?

Interfaces brms
1

I am trying to understand the implications of expressing a group-level intercept model in a multivariate syntax with brms.

To illustrate, I’ll begin by simulating some data with hiearchical properties, where there are 20 observations in each of 26 groups grp, the group-level intercept is distributed Normal(0, 2), the beta coefficient of interest associated with variable x1 is fixed at 1, with some Normal(0,1) error, epsilon:

library(tidyverse)
grpEffects <- data_frame(
  grp=LETTERS[1:26], 
  grpEffect=rnorm(26, sd=2)
)
df <- data_frame(
  x1=rnorm(20*26), 
  grp=rep(LETTERS[1:26], each=20)
)
df <- df %>%
  left_join(grpEffects, by="grp") %>%
  bind_cols(epsilon=rnorm(20*26)) %>%
  mutate(y=1*x1 + grpEffect + epsilon)

Let’s say my goal is to estimate the posterior effect associated with x1. What is the difference between fitting this data using a group-level intercept (m1)

m1 <- brm(y ~ x1 + (1|grp), data=df)

and the following multivariate formula specification (m2)?

m2 <- brm(bf(y ~ x1) + 
          bf(x1 ~ 0 + (1|grp)) + 
          set_rescor(FALSE), data=df)

I understand that different default priors may be used in each case. But beyond that, are these models describing the same data-generating process? Superficially it seems so. The estimated parameters are as expected in both cases. It seems that algebraically these models should be identical. Why might one be preferred over the other?

m2 seems to sample much better, with the number of effective samples at the maximum for most parameters (i.e. all iterations are effective), whereas effective samples are much lower in m1. Perhaps this performance difference can be attributed to priors?

Perhaps more of a modelling question than a brms one, but asking here for interface reasons.

  • Operating System: Windows 10
  • brms Version: 2.4.0
2

These two models are fundamentally different! m1 has only one response variable with a varying intercept on y over variable grp. m2 has two response variables and a varying intercept on x1 (not on y). You can’t interprete those formulas literally as if you were plugging in 0 + (1|grp) into x1. That’s a syntax used in brms’ non-linear models but not in the multivariate ones.

3

Okay–understood. So the reason why these two models give virtually identical parameter estimates is likely because there is no covariance between the y and x1 responses that make up the multivariate normal being modelled in m2. The special case of near equivalence is a consequence of how the fake data were made. Is that a fair assessment?

4

The reason for the similar estimates is pure luck I would say ;-)