I am trying to understand the implications of expressing a group-level intercept model in a multivariate syntax with
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,
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 <- brm(y ~ x1 + (1|grp), data=df)
and the following multivariate formula specification (
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