In one of my other questions, I asked about an approximation to $x = \log_y(y+1)$ of the form $y = f(x)$. One of the answers (by Daniel Schepler) said that $$y = \sum_{n=0}^{\infty}\frac{a_n}{n!}x^{-n}$$ where $a_n$ is the $n$th derivative of $y-(y+1)^x=0$ evaluated at $x = 0, y = 1$.
Following this process using Mathematica, I got $$y = 1+kx^{-1}+\left(\frac{1}{2}k+\frac{1}{2}k^2\right)x^{-2}+\left(\frac{1}{4}k+\frac{5}{8}k^2+\frac{1}{6}k^3\right)x^{-3}+o(x^{-4})$$ where $k = \ln(2)$. Checking a few more evaluations to be sure, I found that each coefficient of $x^{-n}$ was a polynomial in $k$ of degree $n$. I saw that the coefficients of $k^j$ followed a pattern based on $j, n$. For example, the coefficient of $k$ in $x^{-n}$ is $$\frac{1}{2^{n-1}}$$ for $n \ge 1$. As an example, this should mean that the coefficient of $k$ in $x^{-3}$ is $\frac{1}{2^{2}} = \frac{1}{4}$, which is correct.
The coefficient of $k^2$ is $$\frac{(n-1)(n+2)}{2\cdot2^{n}}$$ for $n \ge 2$.
The coefficient of $k^3$ is $$\frac{(n-1)(n-2)(n^2+n+4)}{24\cdot2^{n}}$$ for $n \ge 3$.
Checking up to $j = 5$, which is admittedly a small sample, I found that the coefficient of $k^j$ is of the form $$\frac{P_j(n)*\prod_{r=1}^{j-1}(n-r)}{m_n\cdot 2^n}$$ where $P_j(n)$ is a polynomial of degree $j-1$ with a leading coefficient of $1$ and all coefficients are integers. $m_n$ is a positive integer dependent on $n$. $m_n$ may or may not be equal to $2^{n-2}(n-1)!\cdot n!$, as the values agree for $n = 1$ to $n=5$.
Either expanding on my method or using another, how can I find a closed form for the coefficient of $x^{-n}$ in the approximation of $y$?