How good is this approximation of $n$ factorial?

For any $\alpha > 0$ we have that $\sum_{k\geq 1}e^{-\alpha k}=\frac{1}{e^{\alpha}-1}$, hence by differentiating both sides $n$ times with respect to $\alpha$ we get that $$ (-1)^n\sum_{k\geq 1} k^n e^{-\alpha k} = \frac{d^n}{d\alpha^n}\,\frac{1}{e^\alpha-1}=\frac{d^n}{d\alpha^n}\left[\frac{1}{\alpha}+\sum_{m\geq 1}\frac{B_m}{m!}\alpha^{m-1}\right] \tag{1}$$ hence: $$ (-1)^n \sum_{k\geq 1} k^n e^{-\alpha k} = \frac{n!(-1)^n}{\alpha^{n+1}}+\sum_{m\geq n+1}\frac{B_m}{m!}\cdot\frac{(m-1)!}{(m-n-1)!}\alpha^{m-n-1}\tag{2} $$ and by evaluating at $\alpha=1$ it comes the magic: $$ \sum_{k\geq 1}\frac{k^n}{e^k} = \color{red}{n!}+(-1)^n \sum_{m\geq n+1}\frac{B_m}{m(m-n-1)!}\tag{3} $$ The behaviour of Bernoulli numbers (I am going to talk about Bernoulli numbers with even index, since every Bernoulli number with odd index is zero, with the exception of $B_1$) is quite erratic: till $B_{12}$ they are all less than one in absolute value, then their absolute value starts growing pretty fast: $|B_{2 n}| \sim 4 \sqrt{\pi n} \left(\frac{n}{ \pi e} \right)^{2n}$ for large values of $n$. Since the remainder series in $(3)$ converges pretty fast and the first Bernoulli numbers are essentially negligible, for small values of $n$ (namely $n\leq 16$) $S_n$ and $n!$ are very close, as conjectured.

By the Euler-Maclaurin formula,

$$S_n\approx n!-\frac1{2e}$$

You can get better approximations by taking more terms.

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} \left.\sum_{k = 1}^{\infty}{k^{n} \over \expo{k}}\right\vert_{\ n\ \geq\ 0} & = -\delta_{n0} + \sum_{k = 0}^{\infty}{k^{n} \over \expo{k}} = -\delta_{n0} + \int_{0}^{\infty}{k^{n} \over \expo{k}}\,\dd k + \left.{1 \over 2}{k^{n} \over \expo{k}}\right\vert_{\ k\ =\ 0} -2 \int_{0}^{\infty}{\Im\pars{\bracks{\ic x}^{n}/\expo{\ic x}} \over \expo{2\pi x} - 1}\,\dd x \\[5mm] & = -\,{1 \over 2}\,\delta_{n0} + n! - 2\left\{\begin{array}{lcl} \ds{-\pars{-1}^{n/2}\int_{0}^{\infty}{x^{n}\sin\pars{x} \over \expo{2\pi x} - 1}\,\dd x} & \mbox{if} & \ds{n}\ \mbox{is}\ even \\[3mm] \ds{\phantom{-\,}\pars{-1}^{\pars{n - 1}/2}\int_{0}^{\infty}{x^{n}\cos\pars{x} \over \expo{2\pi x} - 1}\,\dd x} & \mbox{if} & \ds{n}\ \mbox{is}\ odd \end{array}\right. \end{align}

---

$$ \left.\vphantom{\large A}n!\,\right\vert_{\ n\ \geq\ 0} = \sum_{k = 1}^{\infty}{k^{n} \over \expo{k}} + {1 \over 2}\,\delta_{n0} + 2\left\{\begin{array}{lcl} \ds{-\pars{-1}^{n/2}\int_{0}^{\infty}{x^{n}\sin\pars{x} \over \expo{2\pi x} - 1}\,\dd x} & \mbox{if} & \ds{n}\ \mbox{is}\ even \\[3mm] \ds{\phantom{-\,}\pars{-1}^{\pars{n - 1}/2}\int_{0}^{\infty}{x^{n}\cos\pars{x} \over \expo{2\pi x} - 1}\,\dd x} & \mbox{if} & \ds{n}\ \mbox{is}\ odd \end{array}\right. $$

Relative Error plot:

Absolute Error plot: It's worsening as $\ds{n > 16}$. Such behaviour is fully discussed in @Jack D'Aurizio answer.