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 Mö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

rational function identity

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

https://math.stackexchange.com/questions/3753212/is-there-any-geometrical-intuition-for-the-factorials-in-taylor-expansions/3753600#3753600