The quantum analogue of Liouville's theorem
It's subtle. The theorem is not there: quantum flows are compressible (Moyal, 1949).
I'll follow Ch. 0.12 of our book, Concise Treatise of Quantum Mechanics in Phase Space, 2014.
The analog of the Liouville density of classical mechanics is the Wigner function in phase space quantum mechanics. Its evolution equation (generalizing Liouville's) is $$ {\partial f \over \partial t} =\{\!\!\{H ,f\}\!\!\}~, $$ where the double (Moyal) brackets indicate a celebrated quantum modification of Poisson brackets by terms of $O(\hbar^2)$, and serve to prove Ehrenfest's theorem for the evolution of expectation values.
For any phase-space function $k(x,p)$ with no explicit time-dependence, $$ \begin{eqnarray} \frac{d\langle k \rangle }{dt} &=& \int\! dx dp~ {\partial f \over \partial t} k \nonumber \\ &=& {1\over i\hbar} \int\! dx dp~ ( H\star f- f\star H)\star k \nonumber\\ &=& \int\! dx dp~ f \{\!\!\{k,H\}\!\!\} = \left \langle \{\!\!\{ k,H\}\!\!\} \right \rangle , \end{eqnarray} $$ where the star product and its manipulations are detailed in said text.
Moyal stressed (discovered?) that his eponymous quantum evolution equation above contrasts to Liouville's theorem (collisionless Boltzmann equation) for classical phase-space densities, $$ {d f_{cl}\over dt}= {\partial f_{cl} \over \partial t} + \dot{x} ~\partial_x f_{cl} + \dot{p}~ \partial_p f_{cl} =0~. $$
Specifically, unlike its classical counterpart, in general, $f$ does not flow like an incompressible fluid in phase space, thus depriving physical phase-space trajectories of meaning, in this context. (Only the harmonic oscillator evolution is trajectoral, exceptionally.)
For an arbitrary region $\Omega$ about some representative point in phase space, the efflux fails to vanish, $$ \begin{eqnarray} {d \over dt}\! \int_{\Omega}\! \! dx dp ~f&=& \int_{\Omega}\!\! dx dp \left ({\partial f \over \partial t}+ \partial_x (\dot{x} f) + \partial_p (\dot{p} f ) \right)\\ &= & \int_{\Omega}\! \!\! dx dp~ (\{\!\!\{ H,f\}\!\!\}-\{H,f\})\neq 0 ~.\nonumber \end{eqnarray} $$
That is, the phase-space region does not conserve in time the number of points swarming about the representative point: points diffuse away, in general, at a rate of O($\hbar^2$), without maintaining the density of the quantum quasi-probability fluid; and, conversely, they are not prevented from coming together, in contrast to deterministic (incompressible flow) behavior.
Still, for infinite $\Omega$ encompassing the entire phase space, both surface terms above vanish to yield a time-invariant normalization for the WF.
The $O(\hbar^2)$ higher momentum derivatives of the WF present in the MB (but absent in the PB---higher space derivatives probing nonlinearity in the potential) modify the Liouville flow into characteristic quantum configurations. So, negative probability regions moving to the left amount to probability flows to the right, etc... Wigner flows are a recondite field, cf. Steuernagel et al, 2013.
For a Hamiltonian $H=p^2/(2m)+V(x)$, the above evolution equation amounts to an Eulerian probability transport continuity equation,
$$\frac{\partial f}{\partial t} +\partial_x J_x + \partial_p J_p=0~,$$
where, for $\mathrm{sinc}(z)\equiv \sin z/~ z$ , the phase-space flux is
$$J_x=pf/m~ ,\\
J_p= -f \mathrm{sinc} \left( {\hbar \over 2} \overleftarrow {\partial _p} \overrightarrow {\partial _x} \right)~~
\partial_x V(x).
$$
Note added. For a recent discussion/proof of the zeros, singularities,and negative probability density features, hence the ineluctable violations of Liouville's theorem in anharmonic quantum systems see Kakofengitis et al, 2017.
The quantum mechanical analogue of the Liouville theorem is given in terms of a density matrix $\rho$ (see https://en.wikipedia.org/wiki/Density_matrix) and states
$$\frac{\partial\rho}{\partial t}=\frac{i}{\hbar}\left[\rho,H\right]$$
This immediately gives us Ehrenfest's theorem, which states that for any observable $A$, the expectation value $\langle A\rangle=\text{tr}(A\rho)$ obeys the equation
$$\frac{\mathrm{d}}{\mathrm{d}t}\langle A\rangle=\frac{i}{\hbar}\langle\left[A,H\right]\rangle$$
Which, in short, says that expectation values obey the classical equations of motion.