I think the correct Wigner function for the eigenstates is half yours, so take it to be 
\begin{equation}
W_n(x,p) = \frac{(-1)^n}{\hbar \pi} e^{- z/2} L_n(z ),
\end{equation}
where $z=4 H/\hbar\omega  $. 

You *know* that, since the resolution of the identity must be 
$$
\sum_n  W_n= \frac{1}{2\pi \hbar}=1/h,
$$
which indeed holds trivially by the [generating function of the Laguerre polynomials](https://en.wikipedia.org/wiki/Laguerre_polynomials#Recursive_definition,_closed_form,_and_generating_function), 
$$
\sum_n t^n L_n(z)= \frac{e^{-tz/(1-t)}} {1-t}  ~~.
$$

Your sum then collapses to 
$$
 \sum_n \frac{\bar{n}_\mathrm{th}^n}{(1+\bar{n}_\mathrm{th})^n} W_n(x,p)= 
\frac{e^{-z/2}}{\pi \hbar}  \sum_n \left (\frac{- \bar{n}_\mathrm{th}}{1+ \bar{n}_\mathrm{th}} \right )^n  L_n  (z) = \frac{ (1+\bar{n}_\mathrm{th})}{\pi \hbar  (1+2\bar{n}_\mathrm{th})} ~ e^{-z / 2 (1+2 \bar{n}_\mathrm{th}) } ,
$$
or something of the sort.