tl;dr From a thermodynamic perspective, the virial equation of state is a general form into which all smooth equations of states can be cast. From a statistical mechanics perspective, we can derive it from the Mayer cluster expansion, where the $n$th virial coefficient encodes (after some hand-waving) $n$-particle collision events.
An equation of state defines a particular thermodynamic system. Thermodynamics as a theory takes an equation of state as an input and tells you what happens when you add work (be it volumetric, magnetic, etc.) or heat. From the perspective of thermodynamics, any (sufficiently smooth) equation of state whatsoever is pretty much acceptable, in that we can then go on to do the manipulations that will tell us about work and heat along paths in parameter space. From a more mathematical perspective, we can think of the space of variables $(p, T, \rho)$ as a three-dimensional manifold, and the equation of state ($p = \rho R T$ for the ideal gas) selects a particular two-dimensional submanifold as the space of states. If you want to carry this perspective forward, you can use differential geometry to define all of the regular notions of thermodynamics using differential forms (see here).
The equations of state that you mention are all important because their empirical parameters heuristically distill some important effects that are the leading corrections to ideal gas behavior in the relevant regimes. For van der Waals, the two parameters account for finite sized particles and two-particular interactions.
These equations are useful, because it's common to go tabulate all of their parameter values from measurements and then use them to predict gas behavior. But there's an important generalization that's sometimes slightly more unwieldy but has a nice theoretical interpretation: the virial equation of state:
$$
\frac{p}{RT} = A(T) \rho + B(T) \rho^2 + C(T) \rho^3 + \ldots.
$$
This is at its core a Taylor expansion of the pressure as a function of density, which means that if we believe that our equation of state is smooth (analytic), it should always be able to be expressed in this form. $A(T)$ is always 1, because all gases are ideal in the low-density limit. The higher virial coefficients are the corrections as the density increases, and in general they are allowed to have any temperature dependence. If you want, you can also expand $B(T)$ and $C(T)$ in powers of $T$, and then you can fully express an analytic equation of state. The linked Wikipedia page shows how the van der Waals and several other empirical equations of state can be cast into a virial form.
Statistical mechanics is the theory that starts with a microscopic description of a system and tries to abstract away to connect it to the thermodynamic properties. In the context of equations of state of gases, we first find the partition function (which is in general very difficult but straightforward for the ideal gas) and then take derivatives to find the pressure, which will be expressed as a function of the temperature and density as $p = \rho R T$. For real gases, the complication is that the phase space of the entire system has a complicated structure that doesn't factorize nicely into single-particle pieces because particles interact with each other. However, there's a procedure called the Mayer cluster expansion which allows us to compute the partition function in the limit that inter-particle interactions are short-range, so for the most part, molecules fly around ballistically, without interacting---the defining assumption of the gas phase. In statistical field theory (which is what this is quickly becoming), it's dangerous to ascribe too much meaning to formal manipulations and resummings like the cluster expansion. But at a high level, the $n$th term of the expansion represents interactions between a cluster of $n$ particles at once: the third term will represent the contribution from three particles all within range of each others' potential at once, etc. This means that the higher terms only really end up happening when the density is very high and the probability of a several-particle collision is not extremely small. In fact, if we find the pressure from the partition function, what we get is the virial expansion! This leads us to the (hand-wavy) interpretation that the $\rho^n$ term of the virial equation of state represents the effect of $n$-fold collisions between particles.
