In this question, user Yrogirg explains that it is possible to derive a chain rule for $GF(2)[\mathbf{x}]$. That is, there exists a finite difference chain rule for $GF(2)^n\to GF(2)$ functions, where the finite difference is given by the notation:
$$(D_u g)(\mathbf{x}) = g(\mathbf{x}+\mathbf{u})+g(\mathbf{x})$$
This rule then takes the form:
$$(D_u f\circ g) (\mathbf{x}) = (D_u f)(g(\mathbf{x}))\cdot (D_u g)(\mathbf{x})$$
At the end of their answer, they mention a problem that I would also like to know the answer to, is there a similar chain rule for functions $GF(2)^n\to GF(2)^m$?
So, I don't have time to work out the details, but I strongly believe that if a similar rule holds for functions from $GF(2)^n$ to $GF(2^m)$ it should be possible to verify under the idea of Hasse derivatives.
https://en.wikipedia.org/wiki/Hasse_derivative
Let $k[X]$ be a polynomial ring over a field $k.$ The r-th Hasse derivative of $X^n$ is
$${\displaystyle D^{(r)}X^{n}={\binom {n}{r}}X^{n-r},}$$ if $n \geq r$ and zero otherwise.
There is a nice discussion below:
https://math.fontein.de/2009/08/12/the-hasse-derivative/
Notes:
If you are only interested in first derivatives things are simpler.
You do have vectors, not scalars. In that case, the univariate representation under which the above quantities are defined may not be directly usable unless $m$ divides $n,$ i.e., we have a subfield relationship between the domain and range.
Also every function over finite fields is uniquely representable as a polynomial by Lagrange interpolation. So perhaps polynomials can play a role in any case.