It is well-known that: $$ \int_{0}^{+\infty}\frac{x}{e^{x}-1}\,dx = \zeta(2) = \sum_{n\geq 1}\frac{1}{n^2} \tag{1}$$
but what is known about $$ I_2 = \int_{0}^{+\infty}\frac{x^2}{e^x-1-x}\,dx \tag{2} $$ or $$ I_n = \int_{0}^{+\infty}\frac{x^n}{e^x-\sum\limits_{k=0}^{n-1}\frac{x^k}{k!}}\,dx \tag{3}$$ ?
They seem pretty hard to tackle through the usual geometric series approach or the residue theorem, but it would be interesting to find some sharp bounds, too. An approximation that came to my mind is, for instance: $$ I_n \approx \int_{0}^{+\infty} n! e^{-\frac{z}{n+1}}\,dz = (n+1)!$$ but it does not seem so tight. A better approximation is: $$ \exp\left(-\frac{x}{n+1}-\frac{x^2}{2(n+1)^2}\right)\leq \frac{\frac{x^n}{n!}}{\sum_{m>n}\frac{x^m}{m!}}\leq\exp\left(-\frac{x}{n+1}-\frac{nx^2} {2(n+1)^2(n+2)}\right)$$ but that just improves the previous approximation by a constant factor.
The previous inequality has a nice side effect. From $\frac{x}{e^x-1}\approx\exp\left(-\frac{x}{2}-\frac{x^2}{24}\right)$ it follows that: $$ \int_{0}^{+\infty}\left(-\frac{1}{2}-\frac{x}{12}\right)\frac{x}{e^x-1}\,dx \approx \int_{0}^{+\infty}\left(-\frac{1}{2}-\frac{x}{12}\right)\exp\left(-\frac{x}{2}-\frac{x^2}{24}\right)\,dx = -1,$$ or:
$$ \frac{\zeta(2)}{2!}+\frac{\zeta(3)}{3!}\approx 1.$$