Summary - In many fields, count variables often arise in a dependent manner, thus requiring a joint estimation method. Particularly in health care systems, providing reliable estimates of the number of events is crucial to evaluate health care costs. For instance, in the Brazilian health insurance system, actuaries seek to predict the number of doctor appointments and medical exams, as well as other groups of procedures, in order to establish fair premiums. Since a medical exam is performed when prescribed, it is naturally correlated with doctor appointments, making joint analysis of these variables essential. Beyond this intrinsic correlation, unobserved heterogeneity also plays a role. In the previous example, such heterogeneity may arise from portfolio aging, socioeconomic and demographic factors, technological availability, or epidemiological conditions, among others. While some of these factors can be included in regression models as covariates, many relevant variables are difficult to measure, such as the quality of health service providers. Motivated by these challenges, we propose a new class of regression models and study its properties and adequacy. Specifically, we introduce a flexible class of bivariate mixed Poisson (BMP) regression models, which incorporate an exponential-family (EF) distributed component to account for unobserved heterogeneity. This framework addresses overdispersion, a common feature of count data, and provides flexibility in terms of the correlation structure.