Is it possible to estimate covariance between fixed and random effects in brms? If so, how would I do that?

When we talk about covariances between random effects, we typically mean the covariance between values of `effect_1`

and `effect_2`

across the levels of a common grouping factor. Fixed effects are not estimated across the levels of a grouping factor, and so it isn’t meaningful to wonder how they covary with random terms.

At the same time, something we also might be interested in is the posterior covariance between two parameters across multiple samples from the joint posterior (e.g. two fixed effect terms, or a fixed effect and one single level of a random effect). To evaluate this, you can just take the covariance of the posterior iterations for the two parameters of interest.

Edit: Of course you can have a factor as a fixed predictor in your model, and you can look a the posterior covariance between those terms and n-1 levels of a random effect that shares the same grouping factor. However, this covariance will not and cannot be a model parameter, because the fixed terms by definition have independent priors that are entirely defined by their margins. So looking at this covariance, if it were of interest for some reason, would be a matter of looking at the empirical covariance in the model output, rather than introducing and estimating a covariance parameter in the model (doing so would by definition require converting the fixed terms to random).

For some more context:

Here is my current model

m4_cor ← brm(bf(logIVI ~ 0 + chicksl: year + chickagel:year + (1 | site) + (1|yearsite),

sigma ~ year),

data = dw,

prior = c(set_prior(“normal(0, 0.5)”, class = “b”, coef = “chickagel”),

set_prior(“normal(-0.1, 0.5)”, class = “b”, coef = “chicksl”),

set_prior(“lkj(2)”, class = “cor”),

set_prior(“cauchy(0,1)”, class = “sd”)),

warmup = 1000,

iter = 5000,

chains = 3,

cores = 2,

control = list(adapt_delta = 0.99, max_treedepth = 15),

backend = “cmdstanr”,

threads = threading(2))

I am interested in knowing if there is a correlation between year specific intercept and year specific slope, and if there is correlation between year specific intercept and year specific residual variance… Is that possible with this model structure?

If you treat one or both of the effects as fixed, then you can look at the empirical covariance in the estimates across years. That is, at each posterior iteration, compute the covariance, and then examine the posterior distribution of that quantity. When both of the effects are random, you can directly model the covariance as a parameter in a multivariate random effect.