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$

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.

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). $$