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 , \\
\iint \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 , \\
\iint 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) , \\
\iint 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} = \iint 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 discrete eigenbasis $\{ \rvert \lambda_m \rangle \}_{m \in \mathbb{N}}$ (if it has one),
$$ \hat \rho = \sum_m \lambda_m \rvert \lambda_m \rangle \langle \lambda_m \lvert . $$
Given that $\hat \rho^\dagger = \hat \rho, \; \mathrm{Tr}(\hat \rho) = 1,$ and $\langle \lambda_m \lvert \hat \rho \rvert \lambda_m \rangle \ge 0 \; \forall \; \lambda_m$, then $0 \le \lambda_m \le 1 \; \forall \; \lambda_m$. For an orthonormal eigenbasis we have, $\forall \: \lambda_m$,
$$ \hat \rho^2 = \sum_m \lambda_m^2 \rvert \lambda_m \rangle \langle \lambda_m \lvert , \quad 0 \le \lambda_m^2 \le \lambda_m \le 1 , $$
and similarly
$$ \sqrt{\hat \rho} \equiv \sum_m \sqrt{\lambda_m} \rvert \lambda_m \rangle \langle \lambda_m \lvert , \quad 0 \le \lambda_m \le \sqrt{\lambda_m} \le 1 . $$
We choose $\sqrt{\hat \rho}$ to be positive semidefinite because if $\hat \rho = \sqrt{\hat \rho}$, then $\sqrt{\hat \rho}$ can be used to represent a pure state $\left ( \lambda_m = \delta_{mn} , \; n \in \mathbb{N} \right)$. Concluded from this are the properties
$$ \sqrt{\hat \rho}^\dagger = \sqrt{\hat \rho}, \quad \sqrt{\hat \rho} \sqrt{\hat \rho} = \left ( \sqrt{\hat \rho} \right )^2 = \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 (or completely) mixed state $\left ( \lambda_m = 1/N \; \forall \; \lambda_m \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 , \quad 0 \le \varepsilon \le 1 .
\end{align*} $$
The terms $\hat \rho^{1-\varepsilon}$ and $\hat \rho^\varepsilon$ in the last inequality are determined using matrix functions, which also further justifies our choice in $\sqrt{\hat \rho}$.
A way to express $R(q,p)$ should be
$$ \begin{align*}
R(q,p) & = \frac{1}{2 \pi \hbar} \iint \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} \iint \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*} $$