I have 6 binary response variables, and I wish to explore what variables can predict these binary outcomes.

glm(cbind(Y1,Y2)~X1+X2, data = dt) doesn't do the work. It gives me an error: eval(family$initialize) : y values must be 0 <= y <= 1, although all variables have been coded into 1 & 0.

Also, both glm(Y1~X1+X2, data = dt) & glm(Y2~X1+X2, data = dt) work. This means data are in the right class

This answer suggested to use glmer(), but glmer() gives me an error saying that random effects were not specified, and I don't need random effects

This answer didn't provide a way to do multivariate logistic regression

Both Bayesian & frequentist methods are welcomed

  • $\begingroup$ Welcome to Cross Validated! Does it have to be logistic regression, or would multivariate probit regression work? If not, why must you work with logistic regression? $\endgroup$
    – Dave
    Nov 1, 2022 at 14:25
  • $\begingroup$ @Dave thanks for the reply, Dave! Based on my limited knowledge, I think these two can achieve the same goal, so yes it would work! $\endgroup$
    – some
    Nov 1, 2022 at 14:27
  • $\begingroup$ Then there are R packages available for multivariate probit. I don’t think I’ve used this one, but it looks like it does the trick. $\endgroup$
    – Dave
    Nov 1, 2022 at 14:30
  • $\begingroup$ @Dave thanks for the help Dave. I'm currently trying to make the package work. I was wondering what I should do after finding a significant predictor for the responses. Do I run 6 univariate probit regressions on the significant predictor only? $\endgroup$
    – some
    Nov 1, 2022 at 14:40
  • 1
    $\begingroup$ I think this is an R usage/coding problem; outcome must be something like cbind(Y1,Y2), not with a plus sign. See the help page for cbind(). Also see the help page for family(), method 3 for binomial regression : outcome can be "a two-column integer matrix: the first column gives the number of successes and the second the number of failures." Coding per se is off-topic here. If there's a statistical issue beyond coding, please post a new question. $\endgroup$
    – EdM
    Nov 1, 2022 at 15:00

1 Answer 1


For a true multivariate logistic regression with multiple binary outcomes, any of which can have a "success" for any one observation, you can use generalized estimating equations to estimate and take into account the correlations among observations. You effectively get a set of univariate models, but the correlations are taken into account for inference. The R geepack package provides tools for that.

That is typically done for multiple binary observations over time for the same set of individuals, but the basic principle applies for any set of potentially correlated outcomes. This paper illustrates use of geepack with clustered binary outcomes.

The suggestion you saw to use glmer() might have envisioned using the outcome type as a predictor in a logistic model (interacted with the predictors of interest) so that you fit the log-odds of a "success" in part as a function of the outcome type. In that case you would want to include the individuals or cases within which you think there are outcome correlations as random effects. There is no direct handling of multivariate binary outcomes in glmer(), so far as I know.

If instead you have 6 potential mutually exclusive outcomes, only one of which can be a "success" at any one observation, you need to use multinomial or perhaps ordinal regression.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.