Correlation between random effects

Modeling
,
1

Hi guys, i´m new kind of new in the multilevel modeling in brms world, and actually im litte confuse.

The point is to fit this model in brms

y_i = \beta_1 x_{1i}+ \beta_2 x_{2i}+\epsilon_i

in my understanding this is the way

brm(y~-1+x1+x2,data = Data)

The above is the fixed effects model which is consistent with other frequentist libraries and I have no problem with it.

the confusion comes with the random effects model for which I search for values of \beta_j
present variation with respect to another variable, let’s say “Study”, so with a little search in the brms documentation I found that I could model the random effects so that the slopes vary across “ Study ” by

brm(y~-1+x1+x2+(x1+x2|Study),data=Data)

I don’t know exactly if this is 100% correct although in the same way I find consistent results with respect to other libraries.

But the real point of my post is that I do not yet find, if it exists, a way to model the correlation between the random effects in a way other than by the prior distributions and cosy or cor_fixed but I did not clearly understand how they are implemented

i,e say \theta_1 and \theta_2 the true random effects, the intention is to model Cor(\theta_1,\theta_2) =\frac{1}{2}

Any help or some kind of guidance with brms will be a great help

THANK YOU

1 Like
2

Yes, that would usually be a way to write such a model.

I admit that this is a pretty unusual request, so I would recheck you are sure you actually have good reasons to constrain the correlation in this way. But it turns out it is almost possible. brms recognizes a special type of prior called constant that will pin a model parameter to a specific value. The problem is that the correlation matrix is treated as one parameter with a single prior for the whole matrix, so you either need to fix the whole correlation matrix or keep it all variable. Presuming you want to fix the whole matrix, we need to remove the Intercept from the varying effects (to let us have a pure 2x2 correlation matrix with just one element of interest). We can then run:

## Simple simulated data
data=data.frame(y = rnorm(10), x1 = rnorm(10), x2 = rnorm(10), Study = sample(LETTERS[1:3], 10, replace = TRUE))

brm(y~-1+x1+x2+(0 + x1+x2|Study), data = data, 
  prior = prior(
    "constant(cholesky_decompose([[1,0.5],[0.5,1]]))", 
    class = "cor", group = "Study"))

You can learn what class and group values are available via get_prior. Note that we write the correlation matrix using [] and we also call cholesky_decompose (as brms works with the Cholesky decomposition of the matrix internally).

Best of luck with your model!