The computation of the steady state response of large finite element discretized systems subject to periodic excitations is unfeasible because of excessive run time and memory requirements. One could in principle resort to reduced order models stemming from the high fidelity counterparts, which typically require a solution time orders of magnitude smaller. However, when many simulations are required, as in the case of parametric studies, the overall effort could be still significant and the analysis process could be severely hindered. In this work, we propose a sensitivity approach to assess the influence of model parameters on the nonlinear dynamic response. As opposed to the costly evaluation of reduced order solutions over a range of excitation frequencies and model parameters, the sensitivities of a nominal response allows one to approximate the dynamic response by a simple evaluation of an expansion in the directions spanning the parameter space. Special care must be taken on the closure equation that needs to be appended to the system of equations stemming from the harmonic balance method. We discuss the limitations of the current constant frequency approach and propose an improvement. We demonstrate the merits of the proposed approach on a micro-electro-mechanical system affected by parameterized manufacturing defects. Leveraging from a previous contribution, the nonlinear response and the sensitivities are obtained from a reduced order model which is analytical in the defect parameters. Our procedure is able to deliver accurate probability density functions of quantities of interest (e.g. nonlinear resonance peaks, triple solution bandwidth, etc) against statistical distributions of manufacturing defects at negligible computational cost.