In this paper, we address discrete linear systems in the frequency domain, where both frequency and random parameters are considered. Sampling such a system many times is computationally challenging if the size is too large, which is often the case for a finite element discretization of partial differential equations.We propose an adaptive method, steered by dual error indicators, which combines rational Arnoldi model order reduction and sparse grid interpolation with hierarchical Leja nodes. At each Leja node, a reduced order model (ROM) is constructed such that the ROM-error in the frequency domain is balanced with the sparse grid error in the random parameter domain. Both the ROM basis and the sparse grid set are enlarged in a hierarchical manner to achieve a prescribed accuracy in statistical moments. Moreover, parameter sensitivities over the frequency range can be easily extracted from the combined reduced order-surrogate model. In the numerical tests considered in the paper, the method employs sampling sets, which are reduced by at least an order of magnitude, compared to Monte Carlo simulation. Additionally, for an example from vibroacoustics, building a ROM reduces the system size by roughly a factor of 100, while still providing an acceptable accuracy.