Bivariate Generalized Linear Mixed Meta-Analysis model with {brms}

Modeling
1

Hi, I want to fit a Bivariate Generalized Linear Mixed Meta-Analysis model using a dichotomous outcome with brms.

The model is described in Appendix B in this article. In brief, it’s a binomial model with the logit link, in which the average control and treatment groups log odds are treated as fixed effects. Moreover, there are two study-specific parameters (control and treatment groups log odds, random-effects), and these are assumed to follow a bivariate normal distribution with covariance matrix “E”. This matrix includes the between-study variances for the baseline and treatment log odds + the correlation between the baseline and treatment risks in the logit scale. Unfortunately, the authors used SAS and didn’t provide the code.

After some search, I found the same model along with R code using {lme4} in another article. It’s the one depicted in “Model 6: the 'Van Houwelingen bivariate” here. The code they provided is:

lme4::glmer(cbind(event,n-event)~factor(treat)+(control+treat-1|study), data=thedata1, family=binomial(link="logit"))

The brms translation seems very straightforward. My code (please let me know if it’s wrong):

brm(event | trials(n) ~ treat + (control + treat - 1| study),   data = thedata1, family = binomial(link = "logit"))

In addition to fitting this model in brms, I want to generate a figure provided in the first article I mentioned above.

fig
fig750×1378 133 KB

The authors provide multiple formulas in their Appendix B on how to estimate marginal and conditional effects from this model. I don’t know how to translate it to {brms}/{emmeans}/{marginaleffects} synthax.

In summary, I would like to know if my brms translation is correct and how to estimate multiple odds ratios and risk differences while conditioning on multiple baseline risks, so I can plot it (as shown in the figure above).

The full code for data + model fitting can be found below. Thanks!

# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")

# Install/load packages
pacman::p_load(tidyverse,
               lme4,
               brms)

pacman::p_install_gh("stan-dev/cmdstanr")


# Data
thedata1 = 
  dplyr::tribble(
    ~study, ~treat, ~n, ~event, ~control,
    1, 0, 377, 113, 1,
    1, 1, 377, 128, 0,
    2, 0, 40, 4, 1,
    2, 1, 41, 6, 0,
    3, 0, 100, 20, 1,
    3, 1, 101, 22, 0,
    4, 0, 1010, 201, 1,
    4, 1, 1001, 241, 0
  )

# ----

m1 = 
  lme4::glmer(
    cbind(event, n - event) ~ factor(treat) + (control + treat - 1|study),
    data=thedata1,
    family=binomial(link="logit"))

summary(m1)


b1 = 
  brms::brm(
    formula = event | trials(n) ~ treat + (control + treat - 1| study),
    
    data = thedata1,
    family = binomial(link = "logit"),
    cores = parallel::detectCores(),
    backend = "cmdstanr",
    chains = 4,
    warmup = 2000, 
    iter = 4000, 
    seed = 123)

b1