In addition to Ira Gessel's answer, let me explain how such Pfaff–Cauchy–Hurwitz type identities are connected to polynomial interpolation techniques and to Raney lemma.

Denote $\mathbf k=\{1,\dotsc,k\}$.

We prove a more general (not polynomial in $f$, but multilinear in $f_1,\dotsc,f_k$) identity
$$
k(f_1\dotsm f_k)^{(k-1)}=\sum_{\emptyset\ne A\subset\mathbf k}\left(\prod_{i\in A} f_i\right)^{(\lvert A\rvert-1)}\left(\prod_{i\notin A} f_i\right)^{(n-\lvert A\rvert)}.
$$
It suffices to prove this identity when all $f_i$'s are power functions $f_i(t)=t^{x_i}$ (by some standard abstract nonsense). In this case we get a polynomial identity in $x_i's$ (where we denote $(x)_n=x(x-1)\dots(x-n+1)$, and also denote $s(A)=\sum_{i\in A}x_i$):
$$
k(s(\mathbf k))_{k-1}=\sum_{A\sqcup B=\mathbf k,A\ne \emptyset}(s(A))_{\lvert A\rvert-1}\,(s(B))_{\lvert B\rvert}.
$$
Both parts are polynomials of degree at most $k-1$, and it suffices to check that they coincide on the following combinatorial simplex $\Delta$: $x_i$ take non-negative integral values such that $\sum x_i\leqslant k-1$.

Note that for any $(x_1,\dotsc,x_k)\in \Delta$ and any partition $\mathbf k=A\sqcup B$, $A\ne \emptyset$, we have $s(A)+s(B)=s(\mathbf k)\leqslant k-1$, thus either $s(A)<\lvert A\rvert-1$ or $s(B)<\lvert B\rvert$ or $s(A)=\lvert A\rvert-1$, $s(B)=\lvert B\rvert$, $s(\mathbf k)=k-1$. Thus if $s(\mathbf k)<k-1$ both parts vanish, and if $s(\mathbf k)=k-1$, we have to prove that
$$
k!=\sum_{A:s(A)=\lvert A\rvert-1} (\lvert A\rvert-1)! (k-\lvert A\rvert)!
$$
I know this identity from Andrei Zelevinsky (2008). His manuscript A peculiar algebraic identity contains a clever proof both of combinatorial and falling factorials identities. Andrei asked not to circulate this text, but I think that this request is out of date….

Consider a permutation $\pi=(c_1,\dotsc,c_k)$ of $1,2,\dotsc,k$. Choose minimal $m$ for which $s(\{c_1,\dotsc,c_m\})<m$. Then by minimality $s(\{c_1,\dotsc,c_m\})=m-1$. For how many permutations did we have $\{c_1,\dotsc,c_m\}=A$, where $A$ is a fixed subset of $\mathbf k$ with $s(A)=\lvert A\rvert-1$? The answer is $(\lvert A\rvert-1)!\lvert B\rvert!$, since for any fixed cyclic order $(c_1,\dotsc,c_m)$ of $A$ exactly one out of $m$ permutations $(c_{i+1},c_{i+2},\dotsc,c_i)$ satisfies $s(c_{i+1},\dotsc,c_{i+t})\geqslant t$ for all $t=1,2,\dotsc,m-1$. This is a variant of Raney's lemma. Thus the identity.

We may avoid combinatorial arguments with Raney's lemma and remain entirely in algebra. The following trick is suggested by Vlad Volkov: we instead check the points for which $(x_1,\dotsc,x_{k-1},x_k+1)\in \Delta$. In this case $\sum x_i\in \{-1,0,1,\dotsc,k-2\}$ and LHS equals 0 unless $x_1=\dotsb=x_{k-1}=0$, $x_k=-1$. In this last case LHS equals $(-1)^{k-1}k!$. Check the same for RHS. Always either $s(A)<\lvert A\rvert-1$ or $s(B)<\lvert B\rvert=k-\lvert A\rvert$, thus if the corresponding summand in RHS is non-zero, one of the sets $A$, $B$ has a sum of $x$'s equal to $-1$, i.e., consists of $k$ and several indices $i$ such that $x_i=0$. After some calculations everything reduces to Vandermonde–Chu identity.

Now a bit about multidimensional interpolation. Consider the space $\Pi(n,k)$ of polynomials $f(x_1,\dotsc,x_k)$ with degree at most $n$. Its dimension equals the number of monomials $x_1^{c_1}\dotsm x_k^{c_k}$, where $c_i$ are non-negative integers such that $\sum c_i\leqslant n$, we write this as $(c_1,\dotsc,c_k)\in \Delta_n^k$. Thus if we manage to find $\lvert\Delta_n^k\rvert$ linearly independent linear forms on $\Pi(n,k)$, then the polynomial $f\in \Pi(n,k)$ is uniquely determined by the values of these forms. I claim that the values at points of $\Delta_n^k$ serve as such linear forms. For proving that they are linearly independent we may simply construct a polynomial from $\Pi(n,k)$ which vanishes at all points of $\Delta_n^k$ but the given point $(t_1,\dotsc,t_k)$. Here it comes:
$$\prod_{i=1}^k \prod_{j=0}^{t_i-1}(x_i-j)\cdot \prod_{\ell=t_1+\dotsb+t_k+1}^n (x_1+\dotsb+x_k-\ell).$$

9more comments