Multivariate Logistic Regression with brms

I am a new user of brms and I am exploring the way to conduct multivariate logistic regression with brms. I have six binary response variables and five predictors, one is continuous, one is ordinal, and three others are binary. Based on my understanding I found I could use the bernoulli family.

However, it seems this specification assumed the six response variables to be independent, which is against my initial purpose of using multivariate regression, but I can’t find an alternative ways based on the R documentations of brms. Therefore I was wondering if it’s possible to obtain some help from the peers in this forum.

Helps will be much appreciated!

Here’s the code to reproduce my question:

# Set seed for reproducibility
set.seed(123)

# Number of observations
n <- 1000

# Generate predictors
x1 <- rnorm(n, mean = 50, sd = 10)
x2 <- sample(1:5, n, replace = TRUE)
x3 <- rbinom(n, 1, 0.5)
x4 <- rbinom(n, 1, 0.3)
x5 <- rbinom(n, 1, 0.7)

# Generate response variables
y1 <- rbinom(n, 1, plogis(x1*0.02 + 0.3 + x2*0.1 + x3*0.5 + x4*0.4 + x5*0.6))
y2 <- rbinom(n, 1, plogis(x1*0.01 + 0.2 + x2*0.2 + x3*0.3 + x4*0.6 + x5*0.7))
y3 <- rbinom(n, 1, plogis(x1*0.03 + 0.1 + x2*0.3 + x3*0.4 + x4*0.5 + x5*0.2))
y4 <- rbinom(n, 1, plogis(x1*0.05 + 0.4 + x2*0.4 + x3*0.2 + x4*0.3 + x5*0.1))
y5 <- rbinom(n, 1, plogis(x1*0.02 + 0.6 + x2*0.1 + x3*0.7 + x4*0.8 + x5*0.9))
y6 <- rbinom(n, 1, plogis(x1*0.04 + 0.5 + x2*0.2 + x3*0.1 + x4*0.2 + x5*0.3))

# Combine predictors and response variables into a data frame
data <- data.frame(x1, x2, x3, x4, x5, y1, y2, y3, y4, y5, y6)

# brms model
model = brm(mvbind(y1,y2,y3,y4,y5,y6) ~ x1+x2+x3+x4+x5, data = data, family = bernoulli(link = "logit"))

  • Operating System: Mac OS Monterey 12.4
  • brms Version: 2.19.0

You’re right, brms doesn’t offer residual correlation parameters for multiple binary response variables. You might instead consider fitting a cross-classified multilevel model where responses are nested within dependent variables and within participants. See Paul’s IRT preprint (here) or his paper with Mattie on ordinal models (here).

thanks so much for your reply Solomon! I will read these two papers, but based on my brief reading, I’m not sure this two models would fit my design as my data is cross-sectional and not nested. The dependency among response variables are not within-person but because of the similarities among measured items (but it’s not a factor model either). Also the response variables are not ordinal but categorical. Do you have other suggestions? In any case, thanks for your help and suggestions!

This framework can handle cross-sectional data; it’s just a very different way of thinking about nesting.