Using the alternative "Poisson-like" parameterization of the generalized negative binomial distribution in terms of its mean $\lambda$, it follows that the PMF can be written $$ \frac{\Gamma(r+n)}{\Gamma(r)} \frac{r^r}{(r+\lambda)^{r+n}} \frac{\lambda^n}{n!} \,,$$ and therefore its mean $\lambda$ must equal $$\sum_{n=0}^{\infty} n \frac{\Gamma(r+n)}{\Gamma(r)} \frac{r^r}{(r+\lambda)^{r+n}} \frac{\lambda^n}{n!} = \lambda \,.$$ Pulling out all terms in the above sum which don't depend on $n$ leads to the identity:
$$\Gamma(r) \left(\frac{r+\lambda}{r}\right)^r \lambda = \sum_{n=1}^{\infty} \frac{\Gamma(r+n)}{(n-1)!} \left(\frac{\lambda}{r + \lambda} \right)^n \,. $$
Question: How could one prove the above identity "directly"/"non-probabilistically", presumably using properties of the $\Gamma$ function?
My guess is that because $\frac{\Gamma(r+n)}{\Gamma(r)} \approx r^n$ for $n$ large due to Stirling's approximation. Would someone just plug that into the first expression for the mean, and then try to apply the squeeze theorem? Can that really be good enough when the approximation is only good for large $n$ -- how could one possibly show that the error from smaller $n$ "disappears" as well? Specifically I tried doing this and got several unexpected additional factors.
Let $a:=\lambda/(r+\lambda)$. By monotone convergence, swap integral and sum in the definition of the Gamma Function:
\begin{align*} \sum_{n=1}^\infty \frac{\Gamma(r+n)}{(n-1)!}a^n &= \int_0^{\infty}e^{-u}\sum_{n=1}^\infty \frac{1}{(n-1)!} u^{r+n-1} a^ndu\\ &=\int_0^\infty e^{-u}u^rae^{ua}du\\ &=a(1-a)^{-r-1}\Gamma(r+1)\\ &=\frac{\lambda}{\lambda+r}\left(\frac{r+\lambda}{r}\right)^{r+1}\Gamma(r+1)\\ &=\lambda \left(\frac{r+\lambda}{r}\right)^r\Gamma(r). \end{align*}