For $k=0,1,...$, let $a_k$ be the rational number defined by the serie expansion

$\frac{1}{e^x -1} = \frac{1}{x}+ \sum_{k=0}^\infty (-1)^k \frac{a_k}{\Gamma(k+1)} x^k$

The computation of the question simply shows that if there would be a "good value" to attach to $\sum_{n=1}^{+ \infty} n^k$, it should be $a_k$. The real question is why this "good value" is the same as the one obtained via the zeta function: $\zeta (-k)$. The answer is that the analytic continuation of the zeta function is defined by essentially the same computation.

The zeta function is defined by the formula $\zeta(s) = \sum_{n \geq 1} n^{-s}$ for $s \in \mathbb{C}$, $Re(s)>1$. To show that this function admits a meromorphic continuation to the all complex plane, a way to do is to use the trick to write

$n^{-s} = \frac{1}{\Gamma(s)} \int_{0}^{+\infty} e^{-nt} t^{s-1} dt$.

(Remark not useful for what follows: this trick is used very often in quantum field theory, where $t$ is called in this context a "Schwinger parameter" or the "Schwinger proper time").

Using this and summing the geometric serie, one finds for $Re(s)>1$.

$ \zeta(s)= \frac{1}{\Gamma(s)} \int_{0}^{+\infty} \frac{t}{e^t - 1} t^{s-2} dt$.

It is possible to do an integration by part to obtain:

$ - \frac{1}{(s-1) \Gamma(s)} \int_{0}^{+\infty} \frac{d}{dt}(\frac{t}{e^t - 1}) t^{s-1} dt$

The point is that this expression makes sense for $Re(s)>0$ so this formula defines a (and so the) analytic continuation of the zeta function to $Re(s)>0$. Doing successive integrations by part, one obtains the "full" zeta function and it is immediate from this that $\zeta(-k)=a_k$.

This computation defining the analytic continuation of the zeta function is really the same as (precisely the Mellin transform of) the OP's computation. It is one of the standard way to guess what $\sum_{n=1}^{+\infty} n^k$ "should be".

Remark: essentially the same computation appears in other ways to "regularize" $\sum_{n=1}^{+\infty} n^k$, for example the exponential regularization $\sum_{n=1}^{+\infty} n^k e^{-n \epsilon}$