$q,t$-deformed matrix models give rise to representations of the deformed Virasoro algebra and more generally of the quantum toroidal $\mathfrak{gl}_1$ algebra. These representations are described in terms of finite difference equations that induce recursion relations for correlation functions. Under suitable assumptions, these recursions admit unique solutions expressible through "superintegrability" formulas, i.e. explicit closed formulas for averages of Macdonald polynomials. In this paper, we discuss examples arising from localization of 3d $\mathcal{N}=2$ theories, which include $q,t$-deformation of well known classical ensembles: Gaussian, Laguerre and Jacobi. We explain how relations in the quantum toroidal algebra can be used to give a new and universal proof of the known superintegrability formulas, as well as to derive new formulas for models that have not been previously studied in the literature. Finally, we make some remarks regarding the relation between superintegrability and orthogonal polynomials.