<p>Multivariate binary data, often encountered in various fields including epidemiology, biology, and finance, pose unique challenges due to their discrete nature and complex dependencies. Copulas, which provide a flexible framework for capturing and modelling dependencies between random variables, offer a promising approach for analyzing such data. In this work, we propose a copula-based joint regression approach to investigate the relationship between a vector of binary response variables and a set of covariates, with aim to provide a unified and flexible modelling of the predictive probability of success for multivariate binary outcomes. To gain modelling flexibility, the proposed approach assumes that i) the marginal distributions of the binary response variables depend on the covariates through copula models, and ii) the outcomes dependence is also captured through a copula model. Several joint regression models fall within the proposed framework, including multivariate latent probit and logistic models. We provide closed-form expression for the estimator of the predictive probability of multivariate binary outcomes, and derive its asymptotic properties (weak convergence and i.i.d. representation). We validate the performance of the proposed methodology by conducting simulation studies and analyzing the Western Collaborative Group Study (WCGS) cohort. To evaluate the risk of coronary heart disease and corneal arcus (two binary outcomes) and a set of covariates, including age, BMI, systolic and diastolic blood pressure and smoking habits.</p>