Combinatorial proof that $\sum_{k=0}^n \binom{n}{k} \frac{(-1)^k}{(k+1)^2} = \frac{H_{n+1}}{n+1}$

Question

This recent question contains two proofs that $$\sum_{k=0}^n \binom{n}{k} \frac{(-1)^k}{(k+1)^2} = \frac{H_{n+1}}{n+1}.$$ One, by Antonio Vargas, uses a double integral. The other, by me, uses the absorption identity twice. I'd like to see a combinatorial proof, but I haven't managed to come up with one yet.

Some thoughts so far:

1. Since the right-hand side is clearly not an integer for $n \geq 1$, we may need to interpret the identity as a probability. I would accept a combinatorial-based probability argument.
2. Alternatively, we could turn both sides into integers by multiplying by $(n+1)!^2$. The new identity would be $$ \sum_{k=0}^n \binom{n}{k} \frac{(-1)^k (n+1)!^2}{(k+1)^2} = (n+1)! \, n! \, H_{n+1}.$$ I would accept a combinatorial proof of this version of the identity.
3. We need a combinatorial interpretation of the harmonic numbers. There are two in this paper by Benjamin, Preston, and Quinn. The first is that $n! \, H_n$ is the number of permutations on $n+1$ elements that contain exactly two cycles. That implies that $\dfrac{H_n}{n+1}$ (very close to the right-hand side of the identity we want to prove) is the probability that a random permutation on $n+1$ elements contains exactly two cycles. The second interpretation in the paper is that a random permutation on $n$ elements has, on average, $H_n$ cycles.
4. We likely need inclusion/exclusion to interpret the left-hand side.

Answer 1

I found a probabilistic proof with a combinatorial flavor. I'll post it here so this question is officially answered, as well as for anyone who may be interested in seeing it.

Select points $(x_1, y_1), (x_2, y_2), \ldots, (x_{n+1}, y_{n+1})$ independently from the bivariate uniform distribution on the unit square.

Each side is the probability that, for all $i \in \{2, 3, \ldots, n+1 \}$, either $x_1 < x_i$ or $y_1 < y_i$.

Left side: Let $A_i$ denote the event that either $x_1 < x_i$ or $y_1 < y_i$. We want $\cap_{i=2}^{n+1} A_i$, which can be calculated via inclusion-exclusion. To use inclusion-exclusion, we need the probability that $x_1 > x_i$ and $y_1 > y_i$ for any collection of $k$ of the points in $\{(x_2, y_2), \ldots, (x_{n+1}, y_{n+1})\}$. This is the probability that $x_1$ is the largest of $k+1$ randomly-chosen points and $y_1$ is the largest of $k+1$ randomly-chosen points. Since the $x_i$'s and $y_i$'s are independent, this is $\dfrac{1}{(k+1)^2}$. Therefore, one way of calculating the probability we're after is $$1 - \sum_{k=1}^n \binom{n}{k} \frac{(-1)^{k-1}}{(k+1)^2} = \sum_{k=0}^n \binom{n}{k} \frac{(-1)^k}{(k+1)^2}.$$

Right side: Condition on the position of $x_1$ in the ordering of the $x_i$'s from smallest to largest. The probability that $x_1$ is the smallest is $\dfrac{1}{n+1}$, the probability that $x_1$ is the second-smallest is $\dfrac{1}{n+1}$, and so forth, regardless of the position that $x_1$ takes in the ordering of the $x_i$'s. Then, given that $x_1$ is the $k+1$ smallest of the $x_i$'s, we must have $y_1 < y_j$ for each $j$ such that $x_j < x_1$. Since there are $k$ of these $y_j$, the probability that $y_1$ is smaller than all $k$ of them is $\dfrac{1}{k+1}$. Summing up, we get that the probability we're after is also $$\frac{1}{n+1} \sum_{k=0}^n \frac{1}{k+1} = \frac{H_{n+1}}{n+1}.$$