Most sources will take Tao's "identity" as a definition: $dx^1 \wedge \cdots \wedge dx^n$ measures volume in $\mathbf R^n$ in the a way that corresponds to our intuition and previous notions of integration in Euclidean space. It's not clear to me that he's eliminating any "change of variable".
To integrate an $n$-form $\omega$ (let's assume compactly supported so I don't have to worry about convergence) on an $n$-manifold $M$ that is neither (a) embedded in Euclidean space nor (b) parametrized by a single chart, we can do the following. If we first assume that $(U, \varphi)$ is a chart on $M$ containing the support [the points at which the fibre of $\omega$ is non-zero] of $\omega$, then we define
$$
\int_M \omega = \int_{\varphi(U)} (\varphi^{-1})^*\omega
$$
where $(\varphi^{-1})^*\omega$ is the pullback. The integral on the right is computed via Tao's definition.
In general, we take a finite collection of charts $(U_i, \varphi_i)$ covering the support of $\omega$ and a smooth partition of unity $\{\psi_i\}$ of $M$ subordinate to $\{U_i\}$. Then let
$$
\int_M \omega = \sum_i \int_M \psi_i\omega.
$$
This makes sense, because the support of $\psi_i\omega$ is contained in $U_i$. It is not too hard to show that this definition of $\int_M \omega$ is independent of the many choices we've made. For a better discussion, look in any introduction to manifolds, e.g. Lee's Introduction to Smooth Manifolds.
Note. Here the Jacobian is a coordinate-dependent expression hidden in the pullback. I may write out how this goes later if there is interest and time. Also, I would avoid using this to actually compute $\int_M \omega$. There are other ways of computing (which are less theoretically tidy) which are more efficient.
Why I think you formula contains an error $\newcommand{\d}{\mathop{}\!d}$
Let $n = 2$ and $\omega = \d x \wedge \d y$. Then according to your formula
$$I \omega = x \left(\int_0^1 s \d s \right) \d y + y \left( \int_0^1 s \d s\right) \d x = \frac12 (x \d y + y \d x) $$
so $\d (I \omega) = \frac12(\d x \wedge \d y + \d y \wedge \d x) = 0$ and since $\d\omega = 0$ the equation you want to prove does not hold. However if you add a factor of $(-1)^{j-1}$ to the inner sum in the definition of $I \omega$ then $I \omega$ would be $\frac12(x \d y - y \d x)$ in this example and everything works out. The formula for $I$ then is
$$
(I\omega)
= \sum_{i_1\lt\dots\lt i_k}
\sum_{j=1}^kx_{i_j}
\left( \int_0^1s^{k-1}c_{i_1,\dots,i_k}(sx) \d s\right)
(-1)^{j-1} \d x_{i_1}\wedge\dots \wedge \d x_{i_{j-1}}\wedge \d x_{i_{j+1}}\wedge\dots\wedge \d x_{i_k}
$$
Proving $\d(I \omega) + I(\d \omega) = \omega$
Since $\d$ and $I$ are bot additive, it is sufficient to verify this for $\omega = f(x) \d x_{i_1} \wedge \dots \wedge \d x_{i_k}$. By further reordering of the coordinates we can assume $\omega = f \d x_1 \wedge \dots \wedge d x_k$.
Let's denote $\d x_1 \wedge \dots \wedge d x_k$ by $\d x_I$ and $\d x_1 \wedge \dots \wedge \d x_{j-1} \wedge \d x_{j+1} \wedge \dots \wedge d x_k$ by $\d x_{I - j}$.
We have
$$ I \omega = \sum_{j=1}^k x_j \left( \int_0^1 s^{k-1} f(s x) \d s \right) (-1)^{j-1} \d x_{I - j}
$$
and
$$
\begin{align}
\d (I \omega) &= \sum_{i = 1}^n \sum_{j=1}^k \frac{\partial}{\partial x_i} \left[x_j \int_0^1 s^{k-1} f(s x) \d s \right] \d x_i \wedge (-1)^{j-1} d x_{I - j} \\
&= \sum_{i=1}^k \frac{\partial}{\partial x_i} \left[x_i \int_0^1 s^{k-1} f(s x) \d s \right] \d x_I \\
&\phantom{=}+ \sum_{i = k + 1}^n \sum_{j=1}^k x_j \left( \int_0^1 s^k \frac{\partial f}{\partial x_i}(s x) \d s \right) (-1)^{j-1} (-1)^{k-1} d x_{I - j} \wedge d x_i \tag{1}
\end{align}
$$
We can simplify the first sum:
$$
\begin{align}
\sum_{i=1}^k \frac{\partial}{\partial x_i} \left[x_i \int_0^1 s^{k-1} f(s x) \d s \right] &=
\int_0^1 k s^{k-1} f(s x) \d s + \int_0^1 s^k \sum_{i=1}^k x_i \frac{\partial f}{\partial x_i} (s x)
\end{align}
$$
Now
$$
\d \omega = \sum_{i = k + 1}^n \frac{\partial f}{\partial x_i} \d x_i \wedge \d x_I
= (-1)^k \sum_{i = k + 1}^n \frac{\partial f}{\partial x_i} \d x_I \wedge \d x_i
$$
and
$$
\begin{align*}
I(\d \omega) = (-1)^k \sum_{i = k+1}^n &\left[ \sum_{j=1}^k x_j \left( \int_0^1 s^k \frac{\partial f}{\partial x_i}(s x) \right) (-1)^{j-1} \d x_{I - j} \wedge \d x_i \right.\\
&\phantom{[}\left. + x_i \left( \int_0^1 s^k \frac{\partial f}{\partial x_i}(s x) \right) (-1)^k \d x_I\right] \tag{2}
\end{align*}
$$
The first sum in $(2)$ is exactly the negative of the second sum in $(1)$ - $(-1)^k$ vs $(-1)^{k-1}$ - so they cancel when adding $I(\d \omega)$ and $\d (I \omega)$ What remains is
$$
\left( \int_0^1 k s^{k-1} f(s x) \d s + \int_0^1 s^k \sum_{i=1}^n x_i \frac{\partial f}{\partial x_i} (s x) \right) \d x_I
$$
Since the sum on the right hand side is the total derivative of $f(s x)$ with respect to $s$, partial integration shows that this is equal to $f(x) \d x_I = \omega$.
Best Answer
Differential forms cannot be integrated over all subsets of a manifold. You cannot integrate a 2-form over a curve, or a 1-form over a surface. More precisely: integration in measure space is defined over a $\sigma$-algebra. So if integration is defined on two sets $A$ and $B$, the integral is also defined on $A\cap B$. But the same is not true for integration of differential forms. Given a two form $\omega$ in $\mathbb{R}^3$, and take $S$ and $T$ be two smooth surfaces that intersects transversally. $\int_S\omega$ and $\int_T\omega$ make sense, but not $\int_{S\cap T}\omega$.
Generally a differential form, by definition, when taken in a local coordinate representation, is given by functions that are at least continuous. So in fact you don't need the machinery of Lebesgue integration: all the functions you are going to deal with can be treated with Riemann integration.