I think knzhou is right on the money in their answer regarding the core reason for this to be the case. Most of these systems are being probed in the WKB region, which means that, as a rule of thumb, going up by one energy eigenstate means that the state's accessible phase-space area goes up by $\hbar$. The second half of the argument is that, for the simple systems on exhibit, the shape of the accessible phase-space area is mostly independent of the size of the area, or it grows in a pretty predictable fashion by simple scaling on one or both axes. This scaling will also affect $\Delta x$ and $\Delta p$, and it will do so in such a way that their product $\Delta x\,\Delta p$ also increases by $\hbar$.
This insight also points the way at finding systems that don't obey this heuristic, where the shape of the phase-space curve changes with energy, and particularly where the fraction
$$
\frac{\text{volume of phase space at energy }E\leq E_n}{\text{volume of phase-space square box that contains said volume}}
$$
changes as a function of $E_n$.
To get a system like that, I looked for systems whose accessible regions develop 'spikes' as you go up in energy, and the easiest I could think of was systems of the form
$$H=(p-x^3)^2.$$
This needs to be fiddled with a bit to get working - I ended up using $H=p^2-px+x^6$ - and as a function of $p$ and $x$, this Hamiltonian looks like the sort of asymmetric bowl we're trying to create:

This Hamiltonian, in its quantum form $\hat H=\hat p^2-\frac12(\hat p\hat x+\hat x\hat p) + \hat x^6$, is probably amenable to some form of analytical treatment, but the easiest way is to just rough-and-tumble some numerics to get a feel for the structure. In Mathematica, this looks something like
J = 1000; L = 10.; dx = (2 L)/(2 J);
T=-1/2(DiagonalMatrix[Table[1, {2 J}], 1]+DiagonalMatrix[Table[1, {2 J}], -1] -
DiagonalMatrix[Table[2, {2 J + 1}]])/dx^2;
P=-I(DiagonalMatrix[Table[1, {2 J}], 1]-DiagonalMatrix[Table[1, {2 J}], -1])/dx;
X = DiagonalMatrix[Table[j dx, {j, -J, J}]];
and then
SortBy[Eigensystem[2 T - (X^3.P + P.X^3)/2 + X^6]\[Transpose], First]\[Transpose]
The spectrum looks something like this

with eigenfunctions looking like this for the first 20,

and like this for eigenfunctions 280-300,

(which is beginning to look a bit raggedy but still has pretty well-resolved oscillations). This leads us, then, to the uncertainty product $\Delta x\,\Delta p$, calculated as
Table[
Chop[Sqrt[v\[Conjugate].X.X.v] Sqrt[v\[Conjugate].P.P.v]]
, {v, vecs[[1 ;; 300]]}
]
which looks like this:

And this, finally, does not seem to have a linear behaviour. There's more to investigate, for sure, but the overall answer seems to be negative.