Has anyone seen this combinatorial identity involving the Bernoulli and Stirling numbers? Does anyone know a nice (combinatorial?) proof and/or reference for the following identity?

$$\left( \frac{\alpha}{1 - e^{-\alpha}} \right)^{n+1} \equiv \sum_{j=0}^n \frac{(n-j)!}{n!} |s(n+1, n+1-j)| \alpha^j \bmod \alpha^{n+1}.$$

Here
$$\frac{\alpha}{1 - e^{-\alpha}} = \sum_{i=0}^{\infty} (-1)^i \frac{B_i \alpha^i}{i!}$$
is one of the generating functions for the Bernoulli numbers, and $|s(n+1, n+1-j)|$ is an unsigned Stirling number of the first kind.
Motivation (feel free to ignore): this identity comes from two different computations of the Todd class of $\mathbb{CP}^n$. One uses the Euler sequence. The other involves computing the holomorphic Euler characteristic $\chi(\mathcal{O}(k))$ of the line bundles $\mathcal{O}(k)$ using that the higher cohomology of $\mathcal{O}(k)$ vanishes for $k$ large enough and that for $k \ge 0$, $H^0(\mathcal{O}(k))$ is the dimension of the space of homogeneous polynomials of degree $k$ in $n+1$ variables, which is ${k+n \choose n}$, then working out what the Todd class must be using Hirzebruch-Riemann-Roch. This is a bit indirect to say the least, and I have no idea how to convert it into combinatorics. 
 A: The coefficients $B_j^{(r)}$ defined by$$\sum_{j = 0}^\infty B_j^{(r)} {{x^j}\over{j!}} = \left({x\over{e^x - 1}}\right)^r$$are usually called higher order Bernoulli numbers, so your identity is a formula for $B_j^{(r)}$ for $j < r$.
Let $c(n, k) = |s(n, k)| = (-1)^{n - k}s(n, k)$. This is a fairly standard notation, used, for example, in Stanley's "Enumerative Combinatorics". Then$${{(-\log(1 - x))^k}\over{k!}} = \sum_{n = k}^\infty c(n, k) {{x^n}\over{n!}}.$$Differentiating this equation gives$${{(-\log(1 - x))^k}\over{(1 - x)k!}} = \sum_{n = k}^\infty c(n + 1, k + 1) {{x^n}\over{n!}}.\tag*{$(1)$}$$For any formal Laurent series $f = f(\alpha)$ we define the residue of $f$, denoted $\text{res}\,f$, to be the coefficient of $\alpha^{-1}$ in $f$. So the coefficient of $\alpha^k$ in $f$ is $\text{res}\,f/\alpha^{j + 1}$.
We will apply Jacobi's change of variables formula for residues, which is a form of the Lagrange inversion formula. See e.g. Gessel's survey of Lagrange inversion at https://arxiv.org/abs/1609.05988, Theorem 4.1.1.
Suppose that $f(\alpha)$ is a formal Laurent series in $\alpha$ and that $g(\alpha) = g_1 \alpha + g_2\alpha^2 + \ldots$ is a formal power series in $\alpha$ with $g_1 \neq 0$. Then Jacobi's formula says that$$\text{res}\,f(\alpha) = \text{res}\,f(g(\alpha))g'(\alpha).$$We apply Jacobi's formula with$$f(\alpha) = \left({\alpha\over{1 - e^{-\alpha}}}\right)^{n + 1} \alpha^{-j - 1}$$and$$g(\alpha) = -\log(1 - \alpha).$$Then the coefficient of $\alpha^j$ in$$\left({\alpha\over{1 - e^{-\alpha}}}\right)^{n + 1}$$is\begin{align*}
\text{res}\,f(\alpha) & = \text{res}\,f(g(\alpha))g'(\alpha) \\ & = \text{res}\, {{(-\log(1 - \alpha)/\alpha)^{n + 1}}\over{(-\log(1 - \alpha))^{j + 1}(1 - \alpha)}} \\ & = \text{res}\,{{(-\log(1 - \alpha))^{n - j}}\over{\alpha^{n + 1}(1 - \alpha)}}.
\end{align*}This is the coefficient of $\alpha^n$ in$${{(-\log( 1- \alpha))^{n - j}}\over{1 - \alpha}}$$which by $(1)$ is$${{(n - j)!}\over{n!}} c(n + 1, n - j + 1).$$
