Regarding the first candidate I proposed, $\lvert \langle q \vert \hat \rho \rvert p \rangle \rvert^2$ has the properties
$$ \int \lvert \langle q \vert \hat \rho \rvert p \rangle \rvert^2 dp = \langle q \lvert \hat \rho^2 \rvert q \rangle , \quad \int \lvert \langle q \vert \hat \rho \rvert p \rangle \rvert^2 dp = \langle p \lvert \hat \rho^2 \rvert p \rangle , \\
\int \int \lvert \langle q \vert \hat \rho \rvert p \rangle \rvert^2 dq dp = \mathrm{Tr}(\hat \rho^2) \le 1 , $$
which may invalidate the conjecture that accompanied it.
An allure to finding a classical entropy analogue with the Wigner distribution comes from the properties
$$ \int W(q,p) dp = \langle q \lvert \hat \rho \rvert q \rangle , \quad \int W(q,p) dq = \langle p \lvert \hat \rho \rvert p \rangle , \\
\int \int W(q,p) dq dp = \mathrm{Tr}(\hat \rho) = 1 .$$
I think a related question to the one first asked is if there is a joint density distribution $R(q,p)$ with marginal densities $Q(q)$ and $P(p)$ such that
$$ \int R(q,p) dp = \langle q \lvert \hat \rho \rvert q \rangle = Q(q) , \quad \int R(q,p) dq = \langle p \lvert \hat \rho \rvert p \rangle = P(p) , \\
\int \int R(q,p) dq dp = \mathrm{Tr}(\hat \rho) = 1 , $$
and as such with $f(x) = - x \ln x$ also satisfies the property of entropy subadditivity,
$$ S_{qp} \le S_q + S_p , $$
where
$$ S_q = \int f(Q(q)) dq , \quad S_p = \int f(P(p)) dp , \\
S_{qp} = \int \int f(R(q,p)) dq dp = S_{pq} . $$
Under these conditions, I propose another candidate
$$ R(q,p) = \lvert \langle q \vert \sqrt{\hat \rho} \rvert p \rangle \rvert^2 = \lvert \langle p \vert \sqrt{\hat \rho} \rvert q \rangle \rvert^2 \ge 0 . $$
As to what is meant by $\sqrt{\hat \rho}$, consider representing $\hat \rho$ diagonally in its eigenbasis $\{ \rvert \lambda_m \rangle \}_{m \in \mathbb{N}}$,
$$ \hat \rho = \sum_m \lambda_m \rvert \lambda_m \rangle \langle \lambda_m \lvert . $$
Given that $ \hat \rho^\dagger = \hat \rho, \enspace \langle \lambda_m \lvert \hat \rho \rvert \lambda_m \rangle \ge 0 \enspace \forall \enspace \lambda_m$, and $\mathrm{Tr}(\hat \rho) = 1$, then $0 \le \lambda_m \le 1$. For at least an orthonormal eigenbasis, what is then meant is
$$ \sqrt{\hat \rho} \equiv \sum_m \sqrt{\lambda_m} \rvert \lambda_m \rangle \langle \lambda_m \lvert . $$
Concluded from this are the properties
$$ \sqrt{\hat \rho}^\dagger = \sqrt{\hat \rho}, \quad \left ( \sqrt{\hat \rho} \right )^2 = \sqrt{\hat \rho} \sqrt{\hat \rho} = \hat \rho , \\
\mathrm{Tr}(\hat\rho^2) \le \mathrm{Tr}(\hat\rho) = 1 \le \mathrm{Tr}(\sqrt{\hat\rho}). $$
In finite basis of dimension $N$, the inequality above for a totally mixed density matrix $\left ( \rho_{mn} = \delta_{mn}/N \right )$ reads
$$ 1/N \le 1 \le \sqrt N . $$
An example application of the operator $\sqrt{\hat \rho}$ is with fidelity. This clarification was made to compare with other works.
Some comparisons to the first candidate I proposed. From the Cauchy-Schwarz inequality we have
$$ \begin{align*}
R(q,p) & \le \langle q \lvert \sqrt{\hat \rho} \rvert q \rangle \langle p \lvert \sqrt{\hat \rho} \rvert p \rangle , \\
{} & \le \langle q \lvert \sqrt{\hat \rho} \sqrt{\hat \rho} \rvert q \rangle \langle p \vert p \rangle = Q(q) \delta(0) = \infty , \\
{} & \le \langle q \vert q \rangle \langle p \lvert \sqrt{\hat \rho} \sqrt{\hat \rho} \rvert p \rangle = \delta(0) P(p) = \infty , \\
{} & \le \langle q \lvert \hat \rho^{1-\varepsilon} \rvert q \rangle \langle p \lvert \hat \rho^\varepsilon \rvert p \rangle , \qquad 0 \le \varepsilon \le 1 .
\end{align*} $$
A way to express $R(q,p)$ should be
$$ \begin{align*}
R(q,p) & = \frac{1}{2 \pi \hbar} \int \int \langle q \vert \sqrt{\hat \rho} \rvert q_1 \rangle \langle q_2 \vert \sqrt{\hat \rho} \rvert q \rangle \exp \left (\frac{i p (q_1-q_2)}{\hbar} \right ) dq_1 dq_2 & \\
{} & = \frac{1}{2 \pi \hbar} \int \int \langle p_1 \vert \sqrt{\hat \rho} \rvert p \rangle \langle p \vert \sqrt{\hat \rho} \rvert p_2 \rangle \exp \left (\frac{i (p_1 - p_2) q}{\hbar} \right ) dp_1 dp_2 . &
\end{align*} $$
