Lang's Jacobian identity: slicker, elementary proof? In Jeffrey Lang, A Jacobian identity in positive characteristic, J. Commut. Algebra, Volume 7, Number 3 (2015), pp. 393--409, the following result is proven:

Theorem 1. Let $p$ be a prime. Let $\mathbf{k}$ be a commutative $\mathbb{F}_p$-algebra. Let $n$ be a nonnegative integer. Let $R$ be the polynomial ring $\mathbf{k}\left[x_1, x_2, \ldots, x_n\right]$. Let $f_1, f_2, \ldots, f_n$ be $n$ polynomials in $R$. Let $\nabla$ be the differential operator $\prod\limits_{i=1}^n \left(\dfrac{\partial}{\partial x_i}\right)^{p-1}$ on $R$. Let $M \in R^{n\times n}$ be the Jacobian matrix of $f_1, f_2, \ldots, f_n$; this is the $n\times n$-matrix over $R$ whose $\left(i,j\right)$-th entry is $\dfrac{\partial f_i}{\partial x_j}$. Then,
  \begin{align}
& \sum\limits_{\left(i_1, i_2, \ldots, i_n\right) \in \left\{0,1,\ldots,p-1\right\}^n} f_1^{i_1} f_2^{i_2} \cdots f_n^{i_n} \nabla\left(f_1^{p-1-i_1} f_2^{p-1-i_2} \cdots f_n^{p-1-i_n} \right) \\
& = \left(-1\right)^n \left(\det M\right)^{p-1} .
\end{align}

(I have taken the liberty to correct the typo in the paper where the sum ranged over $\left\{1,2,\ldots,p-1\right\}^n$ instead of $\left\{0,1,\ldots,p-1\right\}^n$. Note that later in the paper, in Proposition 1.5, there is a "$0 \leq j \leq n-1$" that should be a "$1 \leq j \leq n-1$"; this suggests blaming a miscommunication between author and editors about which $0$ to replace by a $1$.)
Theorem 1 generalizes Glynn's determinant formula (see Section 6 of Hendrik Lenstra, The unit theorem for finite-dimensional algebras, arXiv:1703.07273v1); indeed, it is easy to see that we can obtain the latter formula from Theorem 1 by setting $f_j = \sum\limits_{i=1}^n a_{ij} x_i$.

Question. Is there an elementary (e.g., combinatorial, inductive, or Hopf-algebraic) proof of Theorem 1?

The proof in Lang's paper relies on a different paper, which in turn relies on some algebraic geometry. I am not sure what is actually used in the proof, but the whole construct seems rather indirect.
One approach that looks promising is to consider the $\mathbf{k}$-linear map $\widetilde{\nabla} : R \otimes R \to R$ (all tensors are over $\mathbf{k}$) that sends each pure tensor $a \otimes b$ to $a \nabla\left(b\right)$. Then, the left hand side of Theorem 1 is
\begin{equation}
\widetilde{\nabla}\left( \prod_{i=1}^n \left(1 \otimes f_i - f_i \otimes 1\right)^{p-1} \right)
\end{equation}
(where the product is in $R \otimes R$). This is because of the classical fact that $\left(X-Y\right)^{p-1} = \sum\limits_{i=0}^{p-1} X^i Y^{p-1-i}$ in $\mathbb{F}_p\left[X,Y\right]$. But the question is whether $\widetilde{\nabla}$ has any good properties with respect to the Hopf algebra $R$.
Also, the operator $\nabla$ can be described rather explicitly on monomials: For any $n$ nonnegative integers $a_1, a_2, \ldots, a_n$, we have
\begin{align}
& \nabla \left(x_1^{a_1} x_2^{a_2} \cdots x_n^{a_n}\right) \\
& = \begin{cases}
\left(-1\right)^n x_1^{a_1-\left(p-1\right)} x_2^{a_2-\left(p-1\right)} \cdots x_n^{a_n-\left(p-1\right)}, & \text{ if } a_i \equiv p-1 \mod p \text{ for all } i ; \\
0, & \text{ otherwise }.
\end{cases}
\end{align}
 A: Awesome question! I haven't looked at Lang's paper yet, so I can't comment on whether this will be a different approach, but it is elementary. I will make use of Glynn's determinant formula at some point later on, so in order to keep this self contained I will start by giving a combinatorial proof of it.
Lemma 1: (Glynn's determinant formula) Suppose $A$ is a matrix with entries in a commutative $\mathbb F_p$-algebra $T$. The coefficient of $x_{1}^{p-1}x_{2}^{p-1}\cdots x_n^{p-1}$ in the expansion of $\prod_{i=1}^{n}\left(\sum_{j=1}^n a_{ij}x_j\right)^{p-1}$ is equal to $(\det A)^{p-1}$.
Proof: By the Macmahon Master Theorem this coefficient is the same as the coefficient of $x_{1}^{p-1}x_{2}^{p-1}\cdots x_n^{p-1}$ in the series $\det (I-AX)^{-1}\in T[[x_1,x_2,\cdots,x_n]]$ where $X$ is the diagonal matrix with $x_i$ on the diagonal. Notice that this coefficient is unchanged if we multiply by $\det(I-AX)^p$ since this doesn't affect monomials where all exponents are $<p$. However this leaves us with $\det(I-AX)^{p-1}$ which is a polynomial of total degree $(p-1)n$ so our coefficient comes from the top homogeneous part, and from here is easily seen to equal $\left(\det A\right)^{p-1}$.

Define the ideal $I=\langle \epsilon_1^p-1,\dots,\epsilon_n^p-1\rangle$ in $R_1=R[\epsilon_1,\dots,\epsilon_n]=\mathbf{k}[x_1,x_2,\dots,x_n,\epsilon_1,\dots,\epsilon_n]$. We can distinguish between two operators on $R_1$:
$$\nabla_{x}=\prod\limits_{i=1}^n \left(\dfrac{\partial}{\partial x_i}\right)^{p-1}, \nabla_{\epsilon}=\prod\limits_{i=1}^n \left(\dfrac{\partial}{\partial \epsilon_i}\right)^{p-1}$$
and notice that since $\frac{\partial}{\partial \epsilon_{i}}$ annihilate $\epsilon_i^p-1$, these are well defined in $R_1/I$ as well.
Lemma 2: We have the following relation in $R_1/I$
$$\begin{align}
& \sum\limits_{\left(i_1, i_2, \ldots, i_n\right) \in \left\{0,1,\ldots,p-1\right\}^n} f_1^{i_1} f_2^{i_2} \cdots f_n^{i_n} \nabla_x\left(f_1^{p-1-i_1} f_2^{p-1-i_2} \cdots f_n^{p-1-i_n} \right) \\
& = \frac{\nabla_{\epsilon}\prod_{i=1}^n(f_i(\epsilon_1 x_1,\dots,\epsilon_n x_n)-f_i(x_1,\dots,x_n))^{p-1}}{x_1^{p-1}x_2^{p-1}\cdots x_n^{p-1}}
\end{align}$$
Proof: We will make use of the identity in $R_1/I$ $$\nabla_x F(x_1,\cdots,x_n)=\frac{\nabla_{\epsilon} F(\epsilon_1x_1,\dots,\epsilon_n x_n)}{x_1^{p-1}x_2^{p-1}\cdots x_n^{p-1}}$$
which is pretty much immediate from the explicit description of $\nabla$ given at the end of your question. By applying this identity to every instance of $\nabla_x$ on the left hand side of our equation and then using the fact that
$$\sum_{i=0}^{p-1} f_j(x)^i f_j(\epsilon x)^{p-1-i}=(f_j(\epsilon x)-f_j(x))^{p-1}$$
we arrive at the right hand side.

Lemma 3: We have
$$\prod_{i=1}^n(f_i(\epsilon_1 x_1,\dots,\epsilon_n x_n)-f_i(x_1,\dots,x_n))^{p-1}=\prod_{i=1}^n \left(\sum_{j=1}^n (\epsilon_j-1)x_j\frac{\partial f_i}{\partial x_j}\right)^{p-1}$$
Proof: We can write 
$$f_i(\epsilon_1 x_1,\dots,\epsilon_n x_n)-f_i(x_1,\dots,x_n)=\sum_{j=1}^n (\epsilon_j-1)x_j\frac{\partial f_i}{\partial x_j}+S_i$$
when expanding in powers of the form $(\epsilon_1-1)^{r_1}(\epsilon_2-1)^{r_2}\cdots (\epsilon_n-1)^{r_n}$ with coefficients in $R$, where $S_i$ consists of "higher order terms", i.e. terms where $r_1+r_2+\cdots+r_n\geq 2$. Since we have $(\epsilon_i-1)^p=0$ we notice that in the expansion of $$\prod_{i=1}^n \left(S_i+\sum_{j=1}^n (\epsilon_j-1)x_j\frac{\partial f_i}{\partial x_j}\right)^{p-1}$$
any term which uses a monomial from any $S_i$ will have degree $r_1+r_2+\cdots +r_n> (p-1)n$ and will therefore vanish, since this implies that some $r_i\geq p$.

We are now ready to complete the proof of Theorem 1 since
$$\nabla_{\epsilon}\prod_{i=1}^n \left(\sum_{j=1}^n (\epsilon_j-1)x_j\frac{\partial f_i}{\partial x_j}\right)^{p-1}$$
is $(-1)^n$ times the coefficient of $\epsilon_1^{p-1}\epsilon_2^{p-1}\cdots \epsilon_n^{p-1}$ in $\prod_{i=1}^n \left(\sum_{j=1}^n (\epsilon_j-1)x_j\frac{\partial f_i}{\partial x_j}\right)^{p-1}$. This last coefficient is of course the same as taking the coefficient of $(\epsilon_1-1)^{p-1}\cdots (\epsilon_n-1)^{p-1}$ in the same expression but expanded in monomials $(\epsilon_1-1)^{r_1}(\epsilon_2-1)^{r_2}\cdots (\epsilon_n-1)^{r_n}$. By Glynn's determinant this last coefficient is $x_1^{p-1}x_2^{p-1}\cdots x_n^{p-1}\left(\det M\right)^{p-1}$ and this completes our proof.
Remark: This proof, and all the lemmas, adapt easily (essentially just change every occurrence of $(p-1)$ to $(p^m-1)$) to give the following generalization of Lang's identity: For any $m\geq 1$, let's denote by $\nabla^{(m)}$ the linear operator that sends monomials of the form $x_1^{a_1}x_2^{a_2}\cdots x_n^{a_n}$ to $x_1^{a_1-p^m+1}x_2^{a_2-p^m+1}\cdots x_n^{a_n-p^m+1}$ when $a_i\cong -1\pmod{p^m}$ for all $1\le i\le n$, and to zero otherwise. The following identity holds:
$$\begin{align}
& \sum\limits_{\left(i_1, i_2, \ldots, i_n\right) \in \left\{0,1,\ldots,p^m-1\right\}^n} f_1^{i_1} f_2^{i_2} \cdots f_n^{i_n} \nabla^{(m)}\left(f_1^{p^m-1-i_1} f_2^{p^m-1-i_2} \cdots f_n^{p^m-1-i_n} \right)= \left(\det M\right)^{p^m-1} .
\end{align}$$
