Evaluating $I(z,s)=\int_0^1\int_0^1\left(1-\frac{(1-x)(1-y)}{(1-(1-z)x)(1-(1-z)y)}\right)^{s-2}\,\mathrm dx\mathrm dy$ I came across the following double integral in a statistics problem:
For $z>0$
$$
I(z,s)=\int_0^1\int_0^1\left(1-\frac{(1-x)(1-y)}{(1-(1-z)x)(1-(1-z)y)}\right)^{s-2}\,\mathrm dx\mathrm dy.
$$
All I was interested in was the values of $s$ for which $I$ converges $(s>0)$. That said, the integral looked interesting enough to try and evaluate.  I was successful in evaluating the special case $z=1$ (see below) but could not solve for the general solution. Any ideas on how to solve it?

Case for $z=1$:
We have
$$
I(1,s)=\int_0^1\int_0^1\left(1-(1-x)(1-y)\right)^{s-2}\,\mathrm dx\mathrm dy.
$$
Substituting $(u,v)=(1-(1-x)(1-y),x)$, the resulting integral in $v$ is easily evaluated in terms of the natural logarithm; hence,
$$
I(1,s)=-\int_0^1 u^{s-2}\log(1-u)\,\mathrm du.
$$
Observe that $-u^{s-2}\log(1-u)\sim u^{s-1}+\mathcal O(u^s)$ as $u\nearrow 0$, which tells us the integral should converge for $s>0$. Using the series expansion for the logarithm we then integrate termwise to find
$$
I(1,s)=\sum_{k=1}^\infty\frac{1}{k(k+s-1)}=\frac{H_{s-1}}{s-1},
$$
where $H_{s-1}=\psi(s)+\gamma$ is the generalized harmonic number, $\psi(z)=\partial_z\log\Gamma(z)$ is the digamma function and $\gamma=0.577\dots$ is the Euler–Mascheroni constant.

 A: Assuming that $s$ is an integer $\geq 2$, each term in the binomial expansion of the integrand is separable into a product of integrals over $x$ and $y$. Defining $N=s-2$ and $z'=1-z$ for notational simplicity,
$$
I(1-z',N+2) = \int_{0}^1 dx \int_{0}^1 dy \sum_{k=0}^N \binom{N}{k} \left( -\frac{(1-x)(1-y)}{(1-z'x)(1-z'y)}\right)^k \\
= \sum_{k=0}^N \binom{N}{k} (-1)^k \left\{
\left(
    \int_{0}^1 \left( \frac{1-x}{1-z'x}\right)^k \, dx
\right)^2
\right\}
$$
As a check, when $z=1$ and thus $z'=0$, the inner integral is $1/(1+k)$, and so
$$
I(1,N+2) = \sum_{k=0}^N \binom{N}{k} \frac{(-1)^k}{(1+k)^2}
$$
which sums (as verified by Mathematica) to the same result derived by you. For general $z'$, define
$$
\xi(z) =\int_0^1 \left( \frac{1-x}{1-(1-z)x} \right)^k \, dx
$$
Unfortunately I cannot evaluate this integral, but it can be expanded to linear order around $z=1$ to give
$$
\xi(z) = \frac{1}{1+k} + \frac{k}{(k+1)(k+2)} (1-z) + \mathcal{O}((1-z)^2)
$$
yielding (for integral $s \geq 2$ at least) (edit and with some corrections to converting the expansion of $\xi(z)$ into that of $\xi^2(z)$)
$$
I(z,s) = \frac{\psi(s)+\gamma}{s-1} + 2 \frac{s(\psi(s+1)+\gamma)-2s+1}{s(s-1)} (1-z) + \mathcal{O}((1-z)^2)
$$
ADDITION
As pointed out by the OP, from Mathematica or integral tables (c.f. https://dlmf.nist.gov/15.6#E1) we have
$$
\xi(z) = \int_0^1 \left(\frac{1-x}{1-(1-z)x}\right)^k \, dx = \frac{1}{k+1} F(1,k;k+2;1-z)
$$
(Note that the referenced table represents the identity using the Olver form of the hypergeometric function $\mathbf{F}(a,b;c;z)=F(a,b;c;z)/\Gamma(c)$, and that $F(a,b; c; z)$ is symmetric in $a$ and $b$). This gives the result, for integral $s \geq 2$,
$$
\left. I(z,s)\right|_{s \in \mathbb{N}, s\geq 2} = \sum_{k=0}^{s-2} \binom{s-2}{k} \frac{(-1)^k}{(1+k)^2} F(1,k;k+2,1-z)^2
$$
At least formally, without worrying about convergence, this could be generalized to non-integer $s$ (and $s=0,1$) by extending the series to infinity (with generalized binomial coefficients), but this is probably less useful. The Taylor series of the squared hypergeometric function is
$$
F(1,k;k+2;1-z)^2 = \left( \sum_{n=0}^\infty \frac{(1)_n (k)_n}{(k+2)_n} \frac{(1-z)^n}{n!} \right)^2 \\
= \left( \sum_{n=0}^\infty \frac{k(k+1)}{(k+n)(k+n+1)} (1-z)^n \right)^2 \\
= \sum_{n=0}^\infty \left\{ \sum_{m=0}^n \frac{k(k+1)}{(k+m)(k+m+1)} \frac{k(k+1)}{(k+n-m)(k+n-m+1)} \right\} (1-z)^n
$$
This inner sum is somewhat messy, but can be broken up with partial fractions to yield
$$
\sum_{n=0}^\infty \frac{2k^2(k+1)^2}{2k+n+1} \left(  \frac{H_{k+n}-H_{k-1}}{2k+n} - \frac{H_{k+n+1}-H_k}{2k+n+2}  \right) (1-z)^n
$$
Substituting into the expression for $I(s,z)$ and exchanging the order of summation, the Taylor series $I(z,s) = \sum_{n=0}^\infty C_n (1-z)^n$ has coefficients
$$
C_n = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{2k^2 (-1)^k}{2k+n+1} \left(  \frac{H_{k+n}-H_{k-1}}{2k+n} - \frac{H_{k+n+1}-H_k}{2k+n+2}  \right)
$$
This reduces to the special cases above for $n=0,1$, e.g.
$$
C_0 = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{2k^2 (-1)^k}{2k+1} \left(  \frac{1}{2k^2} - \frac{1}{2(k+1)^2}  \right) = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{k^2 (-1)^k}{(k+1)^2} \\
C_1 = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{2k^2 (-1)^k}{2k+2} \left(  \frac{1}{k(k+1)} - \frac{1}{(k+1)(k+2)}  \right) = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{2k (-1)^k}{(k+1)^2(k+2)}
$$
and so on.
A: This is meant to supplement the solution derived by @jwimberley. The solution derived was (after some algebraic manipulations):
$$
I(z,s)=\sum_{k=0}^{s-2}\frac{(2-s)_k}{k!}\frac{F^2(1,k;k+2;1-z)}{(k+1)^2},
$$
which held for $z>0$ and $s\in\Bbb N$ with $s\geq 2$. I would like to show that this solution in fact holds for all $s\in\Bbb C: \Re s>0$. First note that the solution above is equivalent to
$$
I(z,s)=1+\sum_{k=1}^\infty\frac{(2-s)_k}{(k+1)^2k!}F^2(1,k;k+2;1-z),
$$
since the first term is always one and $s\in\Bbb N\land s\geq 2\implies (2-s)_k=0$ for $k=s-1,s,s+1,\dots$
To show convergence for other values of $s$ we will use the integral representation of the hypergeometric function (DLMF $15.6.1$) to write
$$
F\left({1,k \atop k+2};1-z\right)=\int_0^1\frac{1}{1-(1-z) t} \frac{t^{k-1}(1-t)}{\operatorname{B}(k,2)}\,\mathrm dt,\quad k\in\Bbb N
$$
which can be interpreted as an expected value w.r.t. the $\operatorname{Beta}(k,2)$ distribution, i.e. $\mathsf E[(1-(1-z) X)^{-1}]$ with $X\sim\operatorname{Beta}(k,2)$. By Jensen's inequality we then have
$$
F^2\left({1,k \atop k+2};1-z\right)\leq \mathsf E[(1-(1-z) X)^{-2}]=\int_0^1\frac{1}{(1-(1-z) t)^2} \frac{t^{k-1}(1-t)}{\operatorname{B}(k,2)}\,\mathrm dt.
$$
Now note for $z>0$ and $t\in[0,1]$ that
$$
\min\{1,z^{-2}\}\leq\frac{1}{(1-(1-z) t)^2}\leq\max\{1,z^{-2}\}.
$$
Consequently,
$$
F^2\left({1,k \atop k+2};1-z\right)\leq\max\{1,z^{-2}\}
$$
and
$$
\tag{1}
I(z,s)\leq1+\max\{1,z^{-2}\}\left(\sum_{k=0}^\infty\frac{(2-s)_k}{(k+1)^2k!}-1\right).
$$
By the properties of the Pochhammer symbol we now write $k+1=(2)_k/(1)_k$. Substituting this result into our inequality for $I$ gives
$$
I(z,s)\leq1+\max\{1,z^{-2}\}\left(\sum_{k=0}^\infty\frac{(1)_k(1)_k(2-s)_k}{(2)_k(2)_kk!}-1\right),
$$
which according to DLMF $\S 16.2(\mathrm{iii})$ is absolutely convergent if $\Re s>0$. Furthermore, if we evaluate the series we can obtain an explicit expression for the upper bound on $I$. Upon inspection of $(1)$ we see that
$$
I(z,s)\leq1+\max\{1,z^{-2}\}\left(I(1,s)-1\right).
$$

Hence, for $z>0$ and $s\in\Bbb C:\Re s>0$:
$$
I(z,s)=\sum_{k=0}^\infty\frac{(2-s)_k}{k!}\frac{F^2(1,k;k+2;1-z)}{(k+1)^2}
$$
and
$$
I(z,s)\leq1+\max\{1,z^{-2}\}\left(\frac{H_{s-1}}{s-1}-1\right).
$$

A: Maybe helpful, I don't really know, sometimes a different perspective sparks something. If you center the integration with the changes of variable
$x = (1 + u)/2$, $y = (1 + v)/2$ the integral becomes
\begin{equation}
z^{s - 2}\int_{-1}^{1}\int_{-1}^{1}
\left( \frac{u  v z - 2uv + uz + vz + z + 2}{(uz -u + z +1)(vz - v + z + 1)}\right)^{s-2}\,du\;dv.
\end{equation}
For a given $z$ the denominator vanishes when either
$u$ or $v$ equal $(1 + z)/(1 - z)$.
Interesting stuff!
Added fun fact
If computer algebra systems are to be trusted we have
\begin{equation}
I(z, 3) = 
{\left(z^{2} + z \log z - 3 \, z + 2\right)} {\left(z - \log z - 1\right)} z
{\left(z - 1\right)}^{-4}.
\end{equation}
I haven't checked if that matches with what we've seen so far, but for $z$ near $1$ this looks like $\frac{3}{4} + \frac{1}{6}(z - 1) - \frac{1}{9}(z - 1)^2 + \cdots$.
As you work your way higher there seems to be a pattern of the form
\begin{equation}
I(z, n) = (1 - z)^{-2n + 2}z \left( A_n^{(2n-1)}(z) + B_n^{(2n-3)}(z)\log^2 z + C_n^{(2n-3)}(z)\log z \right)
\end{equation}
where $A_m^{(k)}(z)$, $B_m^{(k)}(z)$, and $C_m^{(k)}(z)$ are polynomials of degree $k$ in $z$. Would be fun to see if a three-term recurrence relation for them is hiding somewhere.
