I am not sure if I understand your motivation. Is the goal to get formula for rate of change from given initial density matrix to some specified set of density matrices?
Formally, I do not know, but here is an obvious problem with such an idea. The golden rule gives the number of transitions, per unit time, from a single initial state to some state in the energy band considered.
The formula uses a sum over all Hamiltonian eigenstates in the energy band. Superpositions of the eigenstates have to be ignored in the rule's derivation in order to get the well-known formula.
If instead of Hamiltonian eigenvectors one uses density matrices as allowed states, one has much greater set of allowed states. How would you restrict this set of final states so that one is able to derive formula similar to the golden rule? If you allow only those matrices that have all entries zero and single 1 on a diagonal, then we're back to the original description and using density matrices brings only notational complications - how would one define matrix element $V_{ik}$ if one has to use density matrix instead of integral over $\psi$ functions?