The claim we set out to prove is the Rothe-Hagen identity
$$\sum_{k=0}^n
\frac{x}{x+kz} {x+kz\choose k}
\frac{y}{y+(n-k)z} {y+(n-k)z\choose n-k}
= \frac{x+y}{x+y+nz} {x+y+nz\choose n}.$$
We prove it for $x,y,z$ positive integers and since the LHS
and the RHS are in fact polynomials in $x,y,z$ (the
fractional terms cancel with the corresponding binomial
coefficients e.g. $\frac{x}{x+kz} {x+kz\choose k} =
\frac{x}{k!} (x+kz-1)^{\underline{k-1}}$ as long as $x+kz\ne
0$ (consult problem statement)) we then have it for
arbitrary values (we also get polynomials when $k=0$ or
$k=n$.)
Consider the generating function $C(v)$ that satisfies the
functional equation again with $z$ a positive integer
$$C(v) = 1 + v C(v)^z.$$
We ask about again with $x$ a positive integer
$$[v^k] C(v)^x = \frac{1}{k} [v^{k-1}] x C(v)^{x-1} C'(v).$$
This is by the Cauchy Coefficient Formula
$$\frac{x}{k\times 2\pi i}
\int_{|v|=\epsilon} \frac{1}{v^k} C(v)^{x-1} C'(v) \; dv.$$
Here we have assumed that $C(v)$ is analytic in a
neighborhood of the origin so it has a power series
there with radius of convergence $\rho$ and we take
$\epsilon\lt\rho.$
Now we put $C(v) = w$ and we have from the functional
equation $$v = \frac{w-1}{w^z}.$$
To determine the image contour we have from the existence of
the power series by calculating the first two coefficients
from the functional equation that the image of
$|v|=\varepsilon$ is given by $w=1+v+\cdots$ so we get a
circle of radius $\varepsilon$ centered at one plus lower
order terms $(\varepsilon/\rho)^2/(1-\varepsilon/\rho).$
Hence we may choose $\gamma$ such that $|w-1|=\gamma$ is
entirely contained in the image of $|v|=\varepsilon$, which
is contained in an annulus centered at $w=1$ which goes to
$|w-1|=\varepsilon$ as $\varepsilon\ll 1$. (If it's
analytic, then its got the power series, if it's got the
power series then the image is in the annulus.)
This yields the integral
$$\frac{x}{k\times 2\pi i}
\int_{|w-1|=\gamma}
\frac{w^{zk}}{(w-1)^k} w^{x-1} \; dw.$$
The important part is that this image contour includes the
pole at $w=1$, which is the only pole because we assumed
that $x$ and $z$ are positive integers. Continuing,
$$\frac{x}{k\times 2\pi i}
\int_{|w-1|=\gamma}
\frac{1}{(w-1)^k}
\sum_{p=0}^{kz+x-1} {kz+x-1\choose p} (w-1)^p\; dw
\\ = \frac{x}{k} {kz+x-1\choose k-1}
= \frac{x}{x+kz} {x+kz\choose k}.$$
Note that this yields the correct value including for $k=0.$
Now starting from the left of the desired identity we find
$$\sum_{k=0}^n [v^k] C_z(v)^x [v^{n-k}] C_z(v)^y
= [v^n] C_z(v)^x C_z(v)^y = [v^n] C_z(v)^{x+y}.$$
This is the claim. We assumed analyticity of $C(v)$ and
obtained a branch with the given binomial coefficient
expansion. Now for this branch it remains to determine
$\rho.$ We have for the ratio between consecutive terms
$$\frac{x+kz}{x+kz+z} {x+kz+z\choose k+1} {x+kz\choose k}^{-1}
\\ = \frac{k}{k+1} {x+kz+z-1\choose k} {x+kz-1\choose k-1}^{-1}
\\ = \frac{1}{k+1} \prod_{q=0}^{z-1} (x+kz+q)
\prod_{q=1}^{z-1} \frac{1}{x+kz-k+q}
\\ = \frac{k}{k+1} \prod_{q=0}^{z-1} (z+(x+q)/k)
\prod_{q=1}^{z-1} \frac{1}{z-1+(x+q)/k}.$$
In the limit this goes to $z^z/(z-1)^{z-1}$ so that we may
conclude, having determined that if the branch is analytic
in a neighborhood of the origin then its series has radius
of convergence $\rho=(z-1)^{z-1}/z^z \lt 1$ and its
coefficients are $\frac{x}{x+kz} {x+kz\choose k}.$
The same result may be obtained using Lagrange
inversion.
For the LIF computation we put $D(v) =
C(v)-1$ so that we get the functional equation
$$D(v) = v (D(v) + 1)^z.$$
Using the notation from Wikipedia on
LIF
we have $\phi(w) = (w+1)^z$ and $H(v) = (v+1)^x$ and obtain
$$\frac{1}{k} [w^{k-1}] (x (w+1)^{x-1} ((w+1)^z)^k)
= \frac{x}{k} [w^{k-1}] (1+w)^{kz+x-1}
= \frac{x}{k} {kz+x-1\choose k-1}.$$
This matches the first result.
