Cool discovery!
Your observation is probably best understood via the Gamma function. From that Wikipedia page, we can see $k! = \int_0^\infty x^{k} e^{-x} dx$, so your approximation is actually saying $$\int_0^\infty x^{k} e^{-x} dx \approx \sum_{n=1}^\infty n^k e^{-n}$$
When we put it like that, the approximation makes a lot of sense: the summation you found is basically an approximation of that integral using rectangles.
As for carefully proving bounds on exactly how accurate this approximation is as $k$ grows: I see this older MSE question which addresses the same question. The answer on that question goes to very messy lengths to produce a pretty gross error bound involving Bernoulli numbers. So if you want a bound for some specific $k$ then maybe you should look into those results further; but if you're just satisfied with having a better understanding of why your approximation is quite good, then maybe you can spare yourself the hardship of reading through that derivation.
One last comment: The error in this approximation seems quite small, but I'm not sure I'd call it "better than the Stirling formula". After all, the new approximation still involves an infinite sum, and numerically approximating the infinite sum may not be much easier than numerically approximating the original integral. Part of the value of the Stirling formula is that it uses just a few simple functions, so it's easy to plug it into equations or compute values using just a few floating point operations. Basically I'm saying you should just feel good about finding a cool thing, but maybe don't expect it to suddenly replace the role Stirling normally plays when people analyze expressions involving factorials :)
