AbstractPatient‐reported outcomes (PROs) are currently being increasingly used as primary outcome measures in observational and experimental studies since they inform clinicians and researchers about the health‐status of patients and generate data to facilitate improved care. PROs usually appear as discrete and bounded with U, J, or inverse J shapes, and hence, exponential family members offer inadequate distributional fits. The beta‐binomial distribution has been proposed in the literature to fit PROs. However, the fact that the beta‐binomial distribution does not belong to the exponential family limits its applicability in the regression model context, and classical estimation approaches are not straightforward. Moreover, PROs are usually measured in a longitudinal framework in which individuals are followed up for a certain period. Hence, each individual obtains several scores of the PRO over time, which leads to the repeated measures and defines the correlation structure in the data. In this work, we have developed and proposed an estimation procedure for the analysis of correlated discrete and bounded outcomes, particularly PROs, by a beta‐binomial mixed‐effects model. Additionally, we have implemented the methodology in thePROregpackage in R. Because there are similar approaches in the literature to address the same issue, this work also incorporates a comparison study between our proposal and alternative methodologies commonly implemented in R and shows the superior performance of our estimation procedure. This paper was motivated by the analysis of the health‐status of patients with chronic obstructive pulmonary disease, where the main objective is the assessment of risk factors that may affect the evolution of the disease. The application of the proposed approach in the study leads to clinically relevant results.