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, $$ \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.