Yes there is an analog of the Laplace-Runge-Lenz vector, even more so! The N-dimensional harmonic oscillator is one of a handful superintegrable systems where you have a maximal number (2 N - 1) of independent constants of the motion, leading to closed trajectories in classical mechanics, and spectrum degeneracies in QM. It's just that higher dimensional oscillators separate and are trivial to solve, so for a long time people ignored their symmetries. OK, the actual symmetry, even for the classical oscillator, is U(N), so geometrical constructions (Saenz) necessitate complex vectors, but let's not dwell on such.
But, recall, in the case of the LRL vector, most of the conserved charges are dependent, so, in the Coulomb potential case, there are only 5
independent charges, even though the symmetry is SO(4), with 6 generators on top of the Hamiltonian, the 3 angular momenta, L, and the 3 components of the LRL vector, A. Still, that vector obeys two constraints with L and E, so the independent charges are 5: E, L, and only one of the 3 components of A. Thus, the trajectory in the 6D phase space is a line, the intersection of 5 hyper-surfaces. Phew! one more and it would have intersected it on a point in phases space and frozen evolution!
We can thus use the entire A, confident that two of its components depend on everybody else.
OK, now, for specificity, and as per your question, let us look at the harmonic oscillator in 3D, so 6D phase space, setting m =1= ω , i.e. absorbing those in our units:
$$
H=(p_x^2 + x^2 + p_y^2 + y^2 + p_z^2 + z^2 )/2.
$$
The symmetry group of it, amazingly, is SU(3), with 8 charges. Two of them are obvious: Beyond E, they are, let us say, $E_x= (p_x^2 + x^2 )/2$ and
$E_y= (p_y^2 + y^2 )/2$, the Cartan subalgebra charges! Of course, $E_z$
is dependent on them, $=E-E_x-E_y$.
We could spend lots of time picking and choosing the rest, but you may convince yourself by taking their PBs with the Hamiltonian above that a rotation around the z-axis,
$$
L_z= x p_y-y p_x,
$$
and one around the x-axis,
$$
L_x= y p_z-z p_y,
$$
are conserved; it might take a bit more work to show that the remaining SU(3) generators are, in fact dependent on these 5 listed, but that's what the pioneers did:
Finally, an indulgence, Curtright & Zachos 2003; online: The charges are organized to trivially demonstrate that this, like all maximally superintegrable systems, is most elegantly expressible in terms of Nambu Brackets, instead of PBs, after a natural reductio ad dimidium -- but it could be argued this is icing on the cake. On the other hand, that is what you are essentially asking for. Separable oscillators are manifestly solved by inspection, and do not need high-powered methods. On the contrary, they may serve to illustrate and streamline such. This they do.
Edit on QM: Quantum mechanically, the SU(3) structure is more recognizable. Using the creation operator triplets $a^\dagger_i$ in the Jordan-Schwinger map, the above 5 invariant charges amount to saturated Gell-Mann matrix bilinears: $a^\dagger \!\!\cdot a; ~
a^\dagger \lambda_3 a; ~a^\dagger \lambda_8 a; ~a^\dagger \lambda_2 a; ~ a^\dagger \lambda_7 a$.