First, let's establish that they are a good choice. The SciMLBenchmarks are probably the most comprehensive that there are as of right now for modern methods. This uses the vast number of methods available in DifferentialEquations.jl. Let's check the non-stiff ODEs:
Just to make sure, lets look at some benchmarks against MATLAB, Python, and R to see that these implementations are extremely optimized, outperforming those in MATLAB/Python/R by 50x+. So for fast implementations, across the board Runge-Kutta methods seem to do well on non-stiff ODEs and it's very rare when considering cost/error for multistep methods to really be competitive.
Why?
What compensates that 7-fold increase in the computational cost?
There are two major factors. One of them is stability. The explicit Adams-Bashforth methods are horrifically unstable, losing stability as they grow in order, and thus they are not very competitive in practice.
For this reason, optimized software for non-stiff ODEs which use multistep methods generally do not use the Adams-Bashforth methods, but instead they use the Adams-Bashforth-Moulton coefficients.
But wait, ABM is implicit, doesn't that mean it's for stiff equations? No, that's a very simplistic view. In fact, there's explicit integrators for stiff equations too, and so "implicit = stiff" is too simplistic of a view. Instead, implicit methods can be used for non-stiff equations by using less stable methods to solve the nonlinear system. I.e. you cheat a bit and don't use Newton (which is more stable) and instead use fixed-point iteration (with a predictor-corrector) form. However, once you've moved into this form, an update requires solving a fixed point iteration and thus "multistep methods take one f call to update" looks true from the equations but is not true in practice. Thus the calculus becomes a lot more fuzzy. Even then, the Adams-Bashforth-Moulton is still not very stable.
By 5th order you're at the stability like Euler's method... so while methods like CVODE with Adams coefficients, VCABM, ode113, etc. can technically all go to 13th order and take giant steps with few function iterations, in practice that doesn't happen often because of the stability bound.
For reference, look at the stability region for the Dormand-Prince 8/7 method:
This one isn't actually used in practice (contrary to common belief, the Dormand-Prince 8/7 method is not the method used in dop853, instead it's a different tableau that is only well-documented in the software itself, but I digress), but its stability region is still a lot larger. This means that if things start to become a little stability bound it can be taking larger steps for cost.
But what are all of those other numbers? That takes us to the second major factor: the leading coefficients of the truncation error. In contrast to what most mathematicians would hope, the purpose of numerical ODE solving is not to make $\Delta t \rightarrow 0$. In fact, it's the opposite: you want to make $\Delta t$ as large as possible (subject to hitting the error tolerance that you wished). Thus thinking about asymptotic $\mathcal{O(\Delta t^n)}$ is far too simplistic. People in the Runge-Kutta world noticed that if you optimize the leading coefficient of the truncation error, then the biggest error term is still pretty small even when $\Delta t$ is not zero, and this can allow for taking larger steps. Indeed, the Dormand-Prince method (the tableau behind ode45) is one of these "leading truncation error optimizing" methods. (In reality, it made an assumption and optimized under that constraint, while Tsitouras relaxed this constraint in 2011 and created a tableau that has 20% lower of a coefficient, which if you go back and look at the benchmarks you'll see that's Tsit5 which on average performs a little better than DP5, but that's details). But anyways, the multistep methods have pretty massive leading coefficient values compared to this which means that to hit error tolerance they can require a smaller $\Delta t$ and thus more work in the end.
Those two factors together are some of the big players. There's some other details that can matter too, but I think you can contribute most of the behavior to this.
Now, why are BDF methods like ode15s and CVODE okay for stiff equations? Won't the same leading truncation error argument apply? Yes, it does, but those methods are optimized on stiff equations because they take small steps. That's a whole different story for another time.
