Frequency-domain finite element models that depend on random inputs pose significant challenges to surrogate modeling and uncertainty propagation. Output functions are typically higher order rational functions and standard polynomial methods do not perform well. One possible solution consists in constructing a rational surrogate to capture the frequency-dependent dynamics, which is subsequently coupled with polynomial approximations to handle the high-dimensional parametric space. The challenge here lies in the fact that the dynamics may change rapidly over the parametric domain. Recently, rational barycentric approximation in higher dimensions has received interest, however, the approach is not free of the curse of dimensionality and extensions to very high-dimensional settings are still under investigation. In this work, we are specifically interested in surrogate / reduced order modeling for uncertainty propagation in the context of Helmholtz or vibroacoustic problems. The response surface of a quantity of interest that depends on frequency, a random input vector and the mesh resolution will be approximated efficiently in the sparse multilinear approximation framework. Our main contribution is to embed rational approximation as part of the sparse multilinear scheme. Using rational kernels (e.g. the Szegö kernel) for the frequency variable, we obtain significantly more accurate response surfaces compared to standard kernels and we can demonstrate the improvements for both single level and multilevel schemes. The contribution will further present an extension to the setting where additional model input parameters are considered.