Compound symetery structure (and other) for G matrix of random effects using brms

Hi everyone,

Apologies if this is too trivial. I am new to brms and have been unable to find a straightforward solution to my requirement.

I wish to fit a mixed effects model with multiple random slopes across my groups. Something as shown in the code below. I can’t share my original data hence simulating a sample with similar structure.

library(dplyr)
library(brms)

set.seed(123)

# 7 compositions, replicated thrice across 15 groups. Giving 315 observations
design <- tribble(
  ~x1, ~x2, ~x3,
  1, 0, 0,
  0, 1, 0,
  0, 0, 1,
  0.5, 0.5, 0,
  0.5, 0, 0.5,
  0, 0.5, 0.5,
  1/3, 1/3, 1/3
) %>% slice(rep(1:7, each =3)) %>% 
  mutate(ID = paste0("P", 1:21)) %>% 
  slice(rep(1:21, times = 15)) %>% 
  mutate(group = rep(paste0("G", 1:15), each = 21))

pop_slopes <- c(7, 8, 9)
group_slopes <- tibble(group = paste0("G", 1:15),
                      x1_slopes <- rnorm(15, 0, 1),
                      x2_slopes <- rnorm(15, 0, 2),
                      x3_slopes <- rnorm(15, 0, 3))

sim_data <- design %>% 
  left_join(group_slops, by = "group") %>% 
  mutate(response = x1 * pop_slopes[1] + x1_slopes + 
                    x2 * pop_slopes[2] + x2_slopes +
                    x3 * pop_slopes[3] + x3_slopes +
                    rnorm(n(), 0 , 1))

The current model is have is as follows

# brms model
mod <- brm(formula = bf(response ~ 0 + x1 + x2 + x3 + (0 + x1 + x2 + x3 | group),
                        sigma ~ 0 + group),
           data = sim_data)

I am interested in the covariance matrix (G) of the random effects. To the best of my understanding, by default brms fits an unstructured matrix resulting in a separate variance estimate for each term and a separate covariance between each pair of terms. However, I wish to test out other structures such as the compound symmetry or heterogeneous compound symmetry. I gather that this is made possible via the cosy() and other variants but wasn’t able to find examples that used these structure for models with random slopes.

For reference, I wish to recreate the following glmmTMB model in brms.

glmm_mod <- glmmTMB(response ~ 0 + x1 + x2 + x3 + cs(0 + x1 + x2 + x3 | group),
                    dispformula = ~ 0 + group,
                    data = sim_data) 

I would be very grateful if anyone could provide me an example of fitting such a model with brms or point me towards any relevant reference.

Thanks a lot for any guidance you can provide.

Other relevant information

  • Operating System: Windows
  • brms Version: 2.22.00
1 Like

I have some examples based on Singer & Willett’s (2003) text here: 7 Examining the Multilevel Model’s Error Covariance Structure | Applied longitudinal data analysis in brms and the tidyverse

Also, it looks like this is your first post. Welcome to the Stan forums, @zeke!

Hi @Solomon! Thank you very much for the warm welcome! :)

Thanks also for sharing the link to your book chapter. It was very helpful in understanding how to fit repeated measures models their nuances in brms.

However, I wasn’t able to fully connect the dots to move from the examples in the chapter to my use-case. If I understand correctly, the examples in the chapter were for the variance-covariance matrix of residuals R. However, I am interested in the variance-covariance matrix G of random effects. Usually, the same coding structure applies for both matrices, at least in frequentist packages such as nlme or glmmTMB, but I am unable to replicate these structures in brms.

The main problem I am facing in using the auto-correlation functions such as cosy() and unstr() for my data is due to the multiple random slopes in my models. I am unable to understand where should I specify the multiple random slopes within these functions as they only allow specification for a grouping variable in gr (which would be the group in my example) and a single time variable (I assume) in time.

The closest I could get was to obtain a heterogeneous compound symmetry structure for my random effects variance-covariance matrix using the following code

fit7.2 <- brm(data = sim_data, 
      family = gaussian,
      bf(response ~ 0 + x1 + x2 + x3 + (0 + x1 + x2 + x3 || group) + cosy(gr = group),
         sigma ~ 0 + factor(group)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      seed = 7)

However, I am not sure if this is 100% correct or how to fix the terms x1, x2, and x3 to have a common variance so I can get the compound symmetry structure.

Thanks again for taking the time to answer my query.