[proofplan]
We first prove the formula for indicator functions, where it is exactly the definition of the pushforward measure $\mu_X=\mathbb P\circ X^{-1}$. Linearity of the nonnegative integral then gives the formula for nonnegative simple functions. Finally, we approximate an arbitrary nonnegative measurable function from below by nonnegative simple functions and pass to the limit on both sides using the [Monotone Convergence Theorem](/theorems/???).
[/proofplan]
[step:Verify the identity for indicators of measurable sets]
Let $A\in\mathcal E$, and define the indicator function
\begin{align*}
\mathbb 1_A:E&\to\{0,1\}\\
x&\mapsto
\begin{cases}
1,&x\in A,\\
0,&x\notin A.
\end{cases}
\end{align*}
Since $X:(\Omega,\mathcal F)\to(E,\mathcal E)$ is measurable, the preimage
\begin{align*}
X^{-1}(A)=\{\omega\in\Omega:X(\omega)\in A\}
\end{align*}
belongs to $\mathcal F$. For every $\omega\in\Omega$,
\begin{align*}
(\mathbb 1_A\circ X)(\omega)=\mathbb 1_A(X(\omega))=\mathbb 1_{X^{-1}(A)}(\omega).
\end{align*}
Using the definition of the integral of an indicator function and the definition of $\mu_X$, we obtain
\begin{align*}
\int_\Omega \mathbb 1_A(X(\omega))\,d\mathbb P(\omega)
&=\int_\Omega \mathbb 1_{X^{-1}(A)}(\omega)\,d\mathbb P(\omega)\\
&=\mathbb P(X^{-1}(A))\\
&=\mu_X(A)\\
&=\int_E \mathbb 1_A(x)\,d\mu_X(x).
\end{align*}
[/step]
[step:Extend the identity to nonnegative simple functions]
Let $m\in\mathbb N$, let $a_1,\dots,a_m\in[0,\infty)$, let $A_1,\dots,A_m\in\mathcal E$, and define the nonnegative simple function
\begin{align*}
s:E&\to[0,\infty)\\
x&\mapsto \sum_{j=1}^m a_j\mathbb 1_{A_j}(x).
\end{align*}
Then $s$ is $\mathcal E$-measurable, and for every $\omega\in\Omega$,
\begin{align*}
s(X(\omega))=\sum_{j=1}^m a_j\mathbb 1_{A_j}(X(\omega)).
\end{align*}
By finite linearity of the nonnegative integral and by the indicator case already proved,
\begin{align*}
\int_\Omega s(X(\omega))\,d\mathbb P(\omega)
&=\sum_{j=1}^m a_j\int_\Omega \mathbb 1_{A_j}(X(\omega))\,d\mathbb P(\omega)\\
&=\sum_{j=1}^m a_j\int_E \mathbb 1_{A_j}(x)\,d\mu_X(x)\\
&=\int_E s(x)\,d\mu_X(x).
\end{align*}
[guided]
The point of this step is that a simple function is a finite linear combination of indicators, and the previous step already handled each indicator exactly.
Let $m\in\mathbb N$, let $a_1,\dots,a_m\in[0,\infty)$, let $A_1,\dots,A_m\in\mathcal E$, and define
\begin{align*}
s:E&\to[0,\infty)\\
x&\mapsto \sum_{j=1}^m a_j\mathbb 1_{A_j}(x).
\end{align*}
This is a nonnegative $\mathcal E$-measurable simple function because each $A_j$ is $\mathcal E$-measurable and the sum is finite. Composing with $X$ gives, for every $\omega\in\Omega$,
\begin{align*}
s(X(\omega))=\sum_{j=1}^m a_j\mathbb 1_{A_j}(X(\omega)).
\end{align*}
Each function $\omega\mapsto \mathbb 1_{A_j}(X(\omega))$ is the indicator of $X^{-1}(A_j)\in\mathcal F$, so it is $\mathcal F$-measurable.
Now use finite linearity of the nonnegative integral. Since all coefficients $a_j$ are nonnegative, no integrability hypothesis is needed:
\begin{align*}
\int_\Omega s(X(\omega))\,d\mathbb P(\omega)
&=\sum_{j=1}^m a_j\int_\Omega \mathbb 1_{A_j}(X(\omega))\,d\mathbb P(\omega).
\end{align*}
The indicator case applies to each $A_j\in\mathcal E$, so
\begin{align*}
\int_\Omega \mathbb 1_{A_j}(X(\omega))\,d\mathbb P(\omega)
=
\int_E \mathbb 1_{A_j}(x)\,d\mu_X(x)
\end{align*}
for every $j\in\{1,\dots,m\}$. Substituting these identities and using finite linearity again on $(E,\mathcal E,\mu_X)$ gives
\begin{align*}
\int_\Omega s(X(\omega))\,d\mathbb P(\omega)
&=\sum_{j=1}^m a_j\int_E \mathbb 1_{A_j}(x)\,d\mu_X(x)\\
&=\int_E \sum_{j=1}^m a_j\mathbb 1_{A_j}(x)\,d\mu_X(x)\\
&=\int_E s(x)\,d\mu_X(x).
\end{align*}
Thus the pushforward formula holds for every nonnegative simple function.
[/guided]
[/step]
[step:Approximate the nonnegative measurable function from below by simple functions]
Let $g:E\to[0,\infty]$ be $\mathcal E$-measurable. By the [Simple Approximation Theorem](/theorems/???), applied to the nonnegative measurable function $g$ on $(E,\mathcal E)$, there exists a sequence $(s_k)_{k\in\mathbb N}$ of nonnegative $\mathcal E$-measurable simple functions
\begin{align*}
s_k:E&\to[0,\infty)
\end{align*}
such that
\begin{align*}
0\le s_k(x)\le s_{k+1}(x)\quad\text{for every }x\in E\text{ and every }k\in\mathbb N,
\end{align*}
and
\begin{align*}
\lim_{k\to\infty}s_k(x)=g(x)
\end{align*}
for every $x\in E$. Define
\begin{align*}
h_k:\Omega&\to[0,\infty)\\
\omega&\mapsto s_k(X(\omega))
\end{align*}
for each $k\in\mathbb N$, and define
\begin{align*}
h:\Omega&\to[0,\infty]\\
\omega&\mapsto g(X(\omega)).
\end{align*}
Because $X$ is $\mathcal F/\mathcal E$-measurable and each $s_k$ and $g$ is $\mathcal E$-measurable, the functions $h_k$ and $h$ are $\mathcal F$-measurable. Moreover,
\begin{align*}
0\le h_k(\omega)\le h_{k+1}(\omega)
\end{align*}
for every $\omega\in\Omega$ and every $k\in\mathbb N$, and
\begin{align*}
\lim_{k\to\infty}h_k(\omega)=h(\omega)
\end{align*}
for every $\omega\in\Omega$.
[/step]
[step:Pass to the limit with monotone convergence on both measure spaces]
The sequence $(h_k)_{k\in\mathbb N}$ consists of nonnegative $\mathcal F$-measurable functions on $(\Omega,\mathcal F,\mathbb P)$ and increases pointwise to $h$. Therefore the [Monotone Convergence Theorem](/theorems/???) gives
\begin{align*}
\int_\Omega g(X(\omega))\,d\mathbb P(\omega)
=
\int_\Omega h(\omega)\,d\mathbb P(\omega)
=
\lim_{k\to\infty}\int_\Omega h_k(\omega)\,d\mathbb P(\omega).
\end{align*}
For every $k\in\mathbb N$, the simple-function case gives
\begin{align*}
\int_\Omega h_k(\omega)\,d\mathbb P(\omega)
=
\int_\Omega s_k(X(\omega))\,d\mathbb P(\omega)
=
\int_E s_k(x)\,d\mu_X(x).
\end{align*}
The sequence $(s_k)_{k\in\mathbb N}$ consists of nonnegative $\mathcal E$-measurable functions on $(E,\mathcal E,\mu_X)$ and increases pointwise to $g$. Applying the [Monotone Convergence Theorem](/theorems/???) again gives
\begin{align*}
\lim_{k\to\infty}\int_E s_k(x)\,d\mu_X(x)
=
\int_E g(x)\,d\mu_X(x).
\end{align*}
Combining the preceding three displayed identities yields
\begin{align*}
\int_\Omega g(X(\omega))\,d\mathbb P(\omega)
=
\int_E g(x)\,d\mu_X(x).
\end{align*}
This is the desired pushforward formula.
[guided]
We now pass from simple functions to an arbitrary nonnegative measurable function by monotone convergence. The approximation step gave a sequence $(s_k)_{k\in\mathbb N}$ of nonnegative $\mathcal E$-measurable simple functions increasing pointwise to $g$, and it also gave the composed sequence
\begin{align*}
h_k:\Omega&\to[0,\infty)\\
\omega&\mapsto s_k(X(\omega))
\end{align*}
increasing pointwise to
\begin{align*}
h:\Omega&\to[0,\infty]\\
\omega&\mapsto g(X(\omega)).
\end{align*}
We first apply the [Monotone Convergence Theorem](/theorems/???) on the measure space $(\Omega,\mathcal F,\mathbb P)$. Its hypotheses are satisfied: each $h_k$ is nonnegative and $\mathcal F$-measurable, and the sequence satisfies $h_k(\omega)\le h_{k+1}(\omega)$ for every $\omega\in\Omega$ and every $k\in\mathbb N$, with pointwise limit $h$. Hence
\begin{align*}
\int_\Omega g(X(\omega))\,d\mathbb P(\omega)
=
\int_\Omega h(\omega)\,d\mathbb P(\omega)
=
\lim_{k\to\infty}\int_\Omega h_k(\omega)\,d\mathbb P(\omega).
\end{align*}
For each fixed $k\in\mathbb N$, the function $s_k$ is a nonnegative simple function on $E$. The simple-function case therefore applies to $s_k$ and gives
\begin{align*}
\int_\Omega h_k(\omega)\,d\mathbb P(\omega)
=
\int_\Omega s_k(X(\omega))\,d\mathbb P(\omega)
=
\int_E s_k(x)\,d\mu_X(x).
\end{align*}
This identity is the bridge between the two measure spaces: it transfers the integral of the approximation after composition with $X$ to the integral of the approximation against the pushforward measure.
Now apply the [Monotone Convergence Theorem](/theorems/???) on the measure space $(E,\mathcal E,\mu_X)$. Its hypotheses are satisfied because each $s_k$ is nonnegative and $\mathcal E$-measurable, the sequence is pointwise increasing, and its pointwise limit is $g$. Thus
\begin{align*}
\lim_{k\to\infty}\int_E s_k(x)\,d\mu_X(x)
=
\int_E g(x)\,d\mu_X(x).
\end{align*}
Combining the limit identity on $\Omega$, the simple-function identity for each $s_k$, and the limit identity on $E$, we obtain
\begin{align*}
\int_\Omega g(X(\omega))\,d\mathbb P(\omega)
=
\int_E g(x)\,d\mu_X(x).
\end{align*}
This proves the formula for every nonnegative $\mathcal E$-measurable function $g:E\to[0,\infty]$.
[/guided]
[/step]