Bivariate Poisson regression is a method for modeling two correlated count response variables. However, standard Poisson models often assume equidispersion, which is frequently violated in real-world data due to overdispersion. To address this issue, the Bivariate Poisson Log-Normal Regression (BPLNR) model is employed, which incorporates random effects to account for variability beyond that captured by the Poisson distribution. This study applies the BPLNR model to analyze the number of leprosy cases in Indonesia in 2021, categorized by the World Health Organization (WHO) into Paucibacillary (PB) and Multibacillary (MB). These two types are known to be correlated and exhibit overdispersion, rendering standard Bivariate Poisson models inadequate. This research contributes by applying BPLNR to leprosy data in Indonesia—an area that has been underexplored in prior studies, which largely employed univariate or standard Poisson approaches and ignored the correlation and overdispersion structure. Data were obtained from the 2021 Indonesian Health Profile and the Central Statistics Agency. Parameter estimation was conducted using Maximum Likelihood Estimation (MLE) with the Newton-Raphson algorithm, and hypothesis testing was performed using the Maximum Likelihood Ratio Test (MLRT). The results confirm that BPLNR effectively models the joint distribution of PB and MB cases while accounting for overdispersion. Key factors influencing both types of leprosy include population density, poverty rate, access to proper sanitation and drinking water, and availability of medical personnel and health facilities. A limitation of this study is the use of aggregate provincial-level data, which may obscure local variation and spatial effects. Future research could integrate spatial modeling techniques or individual-level data to enhance inference.