We study the existence and azimuthal modulational stability of vortices in the two-dimensional (2D) cubic-quintic nonlinear Schr¨odinger (CQNLS) equation. We use a variational approximation (VA) based on an asymptotically derived ansatz, seeding the result as an initial condition into a numerical optimization routine. Previously known existence bounds for the vortices are recovered by means of this approach. We study the azimuthal modulational stability of the vortices by freezing the radial direction of the Lagrangian functional of the CQNLS, in order to derive a quasi-1D azimuthal equation of motion. A stability analysis is then done in the Fourier space of the azimuthal modes, and the results are analyzed using both the asymptotic variational ansatz and numerically-exact vortices. For unstable vortices, predictions are given for the growth rates of the most unstable azimuthal mode. We also give predictions for the critical value of the frequency, above which all vortices are azimuthally stable. Our predictions are compared to earlier works and corroborated by full 2D simulations. We then briefly study the collisional dynamics between stable vortices of different topological charges.