Multiple outcome logistic regression in R I have a large dataset with multiple binary dependent variables (e.g. y1, y2, y3, y4). yi variables are people's response to their decision making behaviour in different situations. There are many explanatory variables (e.g. x1, x2, ...) which are mainly sociodemographic variables.
I can model each y based on the explanatory variables (e.g. y1~x).
But, all binary dependent variables are correlated and I would like their correlations to be considered in my model. In other words, I would like to model yi ~ xj for all i simultaneously. It can be called a multiple outcome logistic regression model.
I have already tried looking up on the net, glm , and brms packages. None worked for me.
Does anybody know a package that is capable of performing such model?

Maybe nnet::multinom for multinominal logistic regression fits your problem.
If your dichotomous dependent variables are not from a categorical one, perhaps a structural equation model with correlated dichotomous dependent variables is a solution. That would work with lavaan.
