In operation, rotating bladed disks (blisks) are often subject to high levels of dynamic loading, resulting in large amplitudes of forced vibrations especially at resonance. Moreover, variations in structural properties of individual sectors, referred to as mistuning, can lead to strain energy localization and can amplify forced responses. To prevent damages caused by high cycle fatigue, various frictional damping sources are introduced to dissipate vibration energy. Due to the nonlinear behavior of frictional contacts, conventional methods to study the dynamics of the blisk–damper systems are based often on numerical time integration, which is time-consuming and can be computationally prohibitive due to the large sizes of commercial blisk models. Existing techniques for model reduction either rely heavily on cyclic symmetry of the blisk–damper system or are based on component mode synthesis (CMS). However, in the presence of mistuning, cyclic symmetry no longer exists. Also, mistuning is random and best studied statistically. Repetitive CMS condensation for a large amount of random mistuning patterns can lead to a computationally formidable task. This paper presents a reduced-order modeling (ROM) technique to efficiently capture the nonlinear dynamic responses of blisk–damper systems with both small perturbations in blade material properties (small mistuning) and significant changes in the blisk geometries (large mistuning). The ROMs are formed by projecting the blisk–damper systems onto a novel mode basis that mimics the contact behavior. This mode basis contains normal mode shapes of the mistuned blisk–damper systems with either sliding or sticking conditions enforced on the contact surfaces. These mode shapes are computed through the N-PRIME method, a technique recently developed by the authors to efficiently obtain mode shapes for blisks with simultaneous large and small mistuning. The resulting modal nonlinear equations of motion (EOM) are solved by a hybrid frequency/time (HFT) domain method with continuation. In the HFT method, the contact status and friction forces are determined in the time domain by a quasi-two-dimensional contact model at each contact point, whereas the modal EOM are solved in the frequency domain according to a harmonic balance formulation. The forced responses computed by the proposed ROMs are validated for two systems with distinct mistuning patterns. A statistical analysis is performed to study the effectiveness of the frictional dampers under random mistuning patterns.