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