On-shell Poisson brackets and time derivative In classical statistical mechanics, the information about a given system is given by a distribution of probability over phase space $\rho(p,q,t)$.
Let $H(p,q, t)$ be the hamiltonian of the system and $S^t(p,q)$ the associated hamiltonian flow $\forall t>0$.
According to Liouville's equation, on the solutions of Hamilton's equation$^1$ (i.e. the hamiltonian flow) we have.
$$\frac{\partial\rho}{\partial t}=-\{H,\rho\}\tag{1}$$
Now, let $d\Gamma=dpdq$ and
$$Tr()=\int d\Gamma$$
where the integral is extendend over all phase space.
The averages value of an observable $F(p,q)$ at time t is given by
$$\langle F\rangle_t=Tr(\rho(t)F)$$
because the time dependence is encoded in the distribution function. I wrote $\rho(t)$ because $(p,q)$ are dummy integration variables. The time derivative of the average value is then given by
$$\frac{d\langle F\rangle_t}{dt}=Tr\left(\frac{\partial\rho}{\partial t}F \right) \tag{$\star$}$$
Now here is the problem. According to the book I'm reading$^2$, this is equivalent to the quantum Schroedinger picture (density matrix evolves and observables do not). To find the equivalent of Heisenberg picture, he uses $(1)$ in $(\star)$. Although the conclusions are correct, it looks somehow wrong. $(1)$ is evaluated over the solutions of Hamilton's equations, while $(p,q)$ are just integrated variables representing the coordinates of phase space. So, why does this work?
My guess is that what underlies that is a substition in the integral given by the hamiltonian flow, whose jacobian is one (it is a canonical transformation)
$$(\star)=\int d\Gamma\partial_t\rho(p,q,t)F(p,q)=\int\underbrace{J_{S^t}}_{=1}d\Gamma\partial_t\rho(S^t(p,q), t)F(S^t(p,q))$$
and now $(1)$ can be used.
Is this wrong? Or maybe $(1)$ can also be used off-shell?

$^1$This is the crucial point in my question.
$^2$ Principles of Statistical Mechanics, Amnon Katz (1967).
 A: According to Liouville's theorem
$$\frac{d\rho}{dt}(S^t(p,q),t)=0\iff\frac{\partial\rho}{\partial t}(S^t(p,q),t)+\{H,\rho\}=0\tag{A}$$
where the Poisson Brackets are evaluated on the Hamiltonian flow as well. It actually turns out$^1$ that a quantity is a constant of motion iff this holds in all of phase space, meaning in this case the equation holds off-shell.
$(A)$ is equivalent to
$$\rho(S^t(p,q), t)=\rho(p,q,0)\tag{B}$$
Nevertheless, the point was to move the time dependence from the distribution function to the observable, as it happens in quantum mechanics for Heisenberg picture, so at some point the substition in the integral was necessary. On the other hand, the time derivative of the expected value in $(\star)$ (see OP) was actually unneeded, for one only needs to act on $\langle F\rangle_t$ as it follows
$$\langle F\rangle_t=\int d\Gamma\rho(p,q,t)F(p,q)=^2\int d\Gamma|J_{S^t}|\rho(S^t(p,q), t)F(S^t(p,q))\underbrace{=}_{(B)\land J_{S^t}=1}\int d\Gamma\rho(p,q,0)F((S^t(p,q)).$$

$^1$ Valter Moretti, Meccanica Analitica; Remark 12.31.
$^2$ The integration is over all phase space, so I used the same notation and multiplied by Jacobian determinant.
