Construction of the approximation formula
Let's first consider the symmetric standard situation of an interval $[-r,r]$. Let $F$ be an integral function of $f$ and consider the Taylor expansion of $F(0)$ at $x=\pm r$, so somehow reverse of the usual situation.
$$
F(0)=F(\pm r)\mp f(\pm r)r+\frac12f'(\pm r)r^2\mp\frac16f''(\pm r)r^3+\frac1{24}f'''(0)r^4\mp\frac1{120}f^{(4)}(\pm r)r^5+O(r^6).
$$
In the difference we get
$$
0=\int_{-r}^r f(x)\,dx - [f(r)+f(-r)]r + \frac12[f'(r)-f'(-r)]r^2 - \frac16[f''(r)+f''(-r)]r^3+\frac1{24}[f'''(r)-f'''(-r)]r^4 - \frac1{120}[f^{(4)}(r)+f^{(4)}(-r)]r^5+O(r^6)
$$
Now one can consider $\frac12[f''(r)+f''(-r)](2r)$ as trapezoidal formula for its integral $f'(r)-f'(-r)$,
$$
[f''(r)+f''(-r)]r = [f'(r)-f'(-r)]+\frac12[f'''(r)-f'''(-r)]r^2-\frac16[f^{(4)}(r)+f^{(4)}(-r)]r^3+O(r^4)
$$
The same way,
$$
f'''(r)-f'''(-r)=[f^{(4)}(r)+f^{(4)}(-r)]r+O(r^2)
$$
Now insert backwards
\begin{align}
[f''(r)+f''(-r)]r &= [f'(r)-f'(-r)]+\frac13[f^{(4)}(r)+f^{(4)}(-r)]r^3+O(r^4)
\\
\int_{-r}^r f(x)\,dx &= [f(r)+f(-r)]r - \left(\frac12-\frac16\right)[f'(r)-f'(-r)]r^2 +\left(\frac1{18}-\frac1{30}\right)[f^{(4)}(r)+f^{(4)}(-r)]r^5+O(r^6)
\\
&=[f(r)+f(-r)]r -\frac13[f'(r)-f'(-r)]r^2+\frac1{45}[f^{(4)}(r)+f^{(4)}(-r)]r^5+O(r^6)
\end{align}
Exact derivation of the approximation error
One could now conjecture that a more precise treatment of the remainder terms will give a last term $\frac{2}{45}f^{(4)}(\rho)r^5$ for some $\rho\in (-r,r)$. Another way to confirm this is to take the remainder formula and compute its Taylor expansion.
Set
$$
e(r)=\int_{-r}^r f(x)\,dx - [f(r)+f(-r)]r + \frac13[f'(r)-f'(-r)]r^2.
$$
Then its derivatives are
\begin{align}
e'(r)&=-\frac13[f'(r)-f'(-r)]r+\frac13[f''(r)+f''(-r)]r^2,&e'(0)&=0,
\\
e''(r)&=-\frac13[f'(r)-f'(-r)]+\frac13[f''(r)+f''(-r)]r+\frac13[f'''(r)-f'''(-r)]r^2,&e''(0)&=0,
\\
e'''(r)&=[f'''(r)-f'''(-r)]r+\frac13[f^{(4)}(r)+f^{(4)}(-r)]r^2,&e'''(0)&=0,
\\
e^{(4)}(r)&=[f'''(r)-f'''(-r)]+\frac53[f^{(4)}(r)+f^{(4)}(-r)]r+\frac13[f^{(5)}(r)-f^{(5)}(-r)]r^2,&e^{(4)}(0)&=0,
\\
e^{(5)}(r)&=\frac83[f^{(4)}(r)+f^{(4)}(-r)]+\frac73[f^{(5)}(r)-f^{(5)}(-r)]r
+\frac13[f^{(6)}(r)+f^{(6)}(-r)]r^2,&e^{(5)}(0)&=\frac{16}3f^{(4)}(0),
\end{align}
Employing the extended mean value theorem one finds intermediate points $r>r_1>r_2>r_3>|r_4|$ with
\begin{align}
\frac{e(r)}{r^5}&=\frac{e'(r_1)}{5r_1^4}=\frac{e''(r_2)}{20r_2^3}=\frac{e'''(r_3)}{60r_3^2}
\\
&=\frac{f'''(r_3)-f'''(-r_3)}{60r_3}+\frac{f^{(4)}(r_3)+f^{(4)}(-r_3)}{180}
\\
&=\frac{f^{(4)}(r_3)+6f^{(4)}(r_4)+f^{(4)}(-r_3)}{180}
\end{align}
By the intermediate value theorem of continuous functions there is some $\rho\in(-r,r)$ so that
$$f^{(4)}(r_3)+6f^{(4)}(r_4)+f^{(4)}(-r_3)=8f^{(4)}(ρ)
\implies
e(r)=\frac2{45}f^{(4)}(ρ)r^5
$$
Now note that in shifting and scaling this formula to the subdivision of the interval we have $r=\frac h2$, giving a combined error term of
$$
\frac{h^5}{720}\sum_{k=1}^nf^{(4)}(ρ_k)=\frac{h^5}{720}nf^{(4)}(\eta)=\frac{h^4(b-a)}{720}f^{(4)}(\eta),
$$
the next-to-last again by the intermediate value theorem for
$$
y=\frac1n\sum_{k=1}^nf^{(4)}(ρ_k)\in[\min_{x\in[a,b]}f^{(4)}(x),\max_{x\in[a,b]}f^{(4)}(x)].
$$