calculation of Stefan's constant In the calculation of Stefan's constant one has the integral
$$J=\int_0^\infty \frac{x^{3}}{\exp\left(x\right)-1} \, dx$$
which according to Wikipedia is equal to $\frac{\pi^4}{15}$. 
In this page of Wikipedia there is a (long) method of calculation using  the Taylor expansion of $f(k) = \int_0^\infty \frac{\sin\left(kx\right)}{\exp\left(x\right)-1} \, dx$ obtained with contour integration.
In this article is written that $J=\Gamma(4)\,\mathrm{Li}_4(1) = 6\,\mathrm{Li}_4(1) = 6 \zeta(4)$. How to obtain these equalities?
It is also written that there is number way to obtain the result. Someone know another way?

Let us compute something more general the integrals defined by
$$ I_{n,p}=\int_0^\infty \frac{x^n}{(1-\mathrm{e}^{-x})^{p+1}}\mathrm dx$$
that appear in statistical physics when studying bosonic systems.
The simplest method is to develop the denominator into a series
$$ \frac{1}{(1-\mathrm e^{-x})^{p+1}}=
\sum_{k=0}^\infty \binom{k+p}p\left(\mathrm e^{-x}\right)^k $$
and use the integral 
$$ \int_0^\infty x^n\mathrm{e}^{-kx}\mathrm dx= \frac{n!}{k^{n+1}}.$$
Combining these, one gets
$$I_{n,p}=\frac{n!}{p!} \sum_{k=0}^\infty\frac{(k+p)(k+p-1)\cdots(k+1)}{k^{n+1}}\tag{1}$$
in which we recognize the definition of Riemann's zeta function
$\zeta(k)=\sum_{n=1}^\infty n^{-k}$. The numerator in (1) can be rewritten using Stirling numbers of the first kind thanks to
$$ (k+p)\cdots(k+1)=\sum_{m=0}^{p}\genfrac{[}{]}{0pt}{}{p+1}{m+1}k^m$$
and we obtain finally a closed form
$$ I_{n,p}=\frac{n!}{p!}\sum_{m=0}^p\genfrac{[}{]}{0pt}{}{p+1}{m+1}\zeta(n+1-m).$$
In the case $n=3$, $p=0$ the result is
$$I_{3,0}=3!\,\times\zeta(4)=6\times\frac{\pi^4}{90}=\frac{\pi^4}{15}.$$
