Summary: We introduce a new method for constructing reduced-order models that accurately approximate the frequency domain behavior of large-scale, parametrized, linear dynamical systems. Using multivariate rational interpolation, we compute reduced-order models that match transfer function measurements of the large-scale system. The main tools are a new barycentric formula and Loewner matrices formed directly from measurements. This data-driven approach introduces a new degree of freedom for parametrized model reduction – the ability to choose separate reduced orders for each parameter. More precisely, each reduced order is determined by computing the rank of appropriate Loewner matrices. Moreover, we show how to control the pointwise approximation error through a new formula involving the barycentric form and the smallest singular value of a Loewner matrix. Finally, we also give new state-space realizations for multivariate rational functions and conclude with several numerical examples showcasing the effectiveness of the new method.