For an elementary fact like this, which may have been reinvented a thousand times, it is hard to find the first paper where this appeared. However, let me give some missing context. There is a whole industry in constructive quantum field theory and statistical mechanics about related "smart" interpolation formulas or Taylor formulas with integral remainders. These are used to perform so-called cluster expansions. For the OP's identity, there is no loss of generality in taking $u=(0,0,\ldots,0)$ and $v=(1,1,\ldots,1)$. In this case, via Möbius inversion in the Boolean lattice, the formula comes from the following identity.
Let $L$ be a finite set. Let $f:\mathbb{R}^L\rightarrow \mathbb{R}$, $\mathbf{x}=(x_{\ell})_{\ell\in L}\mapsto f(\mathbf{x})$ be a sufficiently smooth function, and let $\mathbf{1}=(1,\ldots,1)\in\mathbb{R}^L$, then $$ f(\mathbf{1})=\sum_{A\subseteq L}\int_{[0,1]^A}d\mathbf{h} \left[\left(\prod_{\ell\in A}\frac{\partial}{\partial x_{\ell}}\right)f\right](\psi_A(\mathbf{h})) $$ where $\psi_A(\mathbf{h})$ is the element $\mathbf{x}=(x_{\ell})_{\ell\in L}$ of $\mathbb{R}^L$ defined from the element $\mathbf{h}=(h_{\ell})_{\ell\in A}$ in $[0,1]^A$ by the rule: $x_{\ell}=0$ if $\ell\notin A$ and $x_{\ell}=h_{\ell}$ if $\ell\in A$. Of course one needs to 1) apply this to all $L$'s which are subsets of $[p]$, 2) use MoebiusMöbius inversion in the Boolean lattice, and 3) specialize to $L=[p]$, and this gives the OP's identity.
The above formula is the most naive one of its kind used to do a "pair of cubes" cluster expansion. See formula III.1 in the article
A. Abdesselam and V. Rivasseau, "Trees, forests and jungles: a botanical garden for cluster expansions".
It is also explained in words on page 115 of the book
V. Rivasseau, "From Perturbative to Constructive Renormalization".
Now the formula is a particular case of a much more powerful one, namely, Lemma 1 in
A. Abdesselam and V. Rivasseau, "An explicit large versus small field multiscale cluster expansion",
where one sums over "allowed" sequences $(\ell_1,\ldots,\ell_k)$ of arbitrary length of elements of $L$, instead of subsets of $L$. The notion of allowed is based on an arbitrary stopping rule. The above identity corresponds to "allowed"$=$"without repeats", or the stopping rule that one should not tack on an $\ell$ at the end of a sequence where it already appeared. By playing with this kind of choice of stopping rule one can use Lemma 1 of my article with Rivasseau, to prove the Hermite-Genocchi formula, the anisotropic Taylor formula by Hairer in Appendix A of "A theory of regularity structures" and many other things. When $f$ is the exponential of a linear form for instance, one can obtain various algebraic identities as in the MO posts
Identity involving sum over permutations
I forgot to mention, one can use Lemma 1 to derive the Taylor formula from calculus 1. This corresponds to $L$ having one element and defining allowed sequences as the ones of length at most $n$. See