Extract variance-covariance matrix of random effect in brms

Interfaces brms
1

Hi,

Let’s say I have 10 correlated measurements of some quantity in 300 participants, and I want to identify characteristics associated with high values. I could run 10 independent models, or use a multivariate model. But I thought of another trick (not sure if it’s common): put the data in “long” format, and model the quantity of any measurement, with random intercepts and random slopes by type of measurement.

One advantage is that I only use one model, with the random intercept becoming the mean value of each type of measurement. Another advantage is that, with the random slope, I am able to estimate an “average” association of each of the covariates with all 10 measurements, while at the same time being able to look at specific associations. I find this average effect very useful to communicate complex results.

I made some simulations and it seems to work exactly as intended:

# settings
n_patients = 300
n_groups = 10
set.seed(123)  
mu = runif(n=n_groups,min=1,max=20)  # mean value per group
sigma = runif(n=n_groups,min=0.5,max=2) # std dev per group
p_female = .5 # proportion of females
beta_female = rnorm(n=n_groups,mean=2,sd=2) # beta of females per group
rho = randcorr(p = n_groups)
Sigma = diag(sigma) %*% rho %*% diag(sigma)

# simulate
set.seed(123)  
betas = tibble(group=paste0("V",1:n_groups),beta_female=beta_female)
data_long = 
  MASS::mvrnorm(n_patients, mu, Sigma) %>% 
  as_tibble() %>% 
  mutate(female=rbinom(n_patients,1,p_female)) %>% 
  pivot_longer(1:n_groups,names_to = "group") %>% 
  left_join(betas,by="group") %>%
  mutate(value=value+female*beta_female)

# model
modB4 = brm(value ~ -1 + female + (1 + female|group), data=data_long)
summary(modB4)
ranef(modB4)

I guess I could do something similar with a multivariate normal model, but the thing is that my real application is a bit more complex with data imputation and a hurdle-gamma model, and I don’t think it would work as well.

The only thing that I could not find is how to extract the variance-covariance matrix for the random intercept (so basically Sigma or rho in my code). I found a similar question here, but no answer either. Is there a way to do that? I would like to avoid recoding everything in stan.

Thanks in advance!

2

I’m not aware of any simple way to do this, but you can make it by hand with a little work. Just keep in mind brm() parameterizes the upper-level variation parameters in a standard-deviation metric, rather than as variances. You can extract all your \sigma parameters, as well as the correlations, with functions like as_draws_df() (which returns the draws for ALL parameters in the model), or with VarCorr() (which returns the \sigma and correlation parameters). If you use VarCorr(), you’ll want to set the summary argument to FALSE if you want the posterior draws rather than the typical summaries. Whether you use the as_draws_df() or VarCorr() method, you can convert the standard deviations and correlations to variances and covariances with the usual formulas.

3

Here’s an example of me doing something like this for one of the models from Singer and Willett’s (2003) text: 4 Doing Data Analysis with the Multilevel Model for Change | Applied longitudinal data analysis in brms and the tidyverse. You want Section 4.5.2.2, Model D: The controlled effects of COA.

1 Like
4

Thanks for the reply! I will try this out.

1 Like