This isn't an answer but I have a function which gets you as close to $\frac{1}{4}$ as you desire. Further, I think the method I used could be used to find an even higher value (EDIT: maybe not, see comments at bottom) but I'm not sure (I might come back to this if I can get better results). To this end, I will first give my example and then explain how I found it in the hopes that it might lead to a final answer at some point.
Define the function
$$f(x)=\begin{cases}
\epsilon +\frac{1}{2} & x\leq \frac{1}{2}-\epsilon \\
-\frac{x^3 (1-2 \epsilon )}{\epsilon ^3}+\frac{3 x^2 (\epsilon -1) (2 \epsilon -1)}{2 \epsilon ^3}-\frac{3 x \left(4 \epsilon ^2-4 \epsilon +1\right)}{4 \epsilon ^3}-\frac{-8 \epsilon ^3-6 \epsilon ^2+5
\epsilon -1}{8 \epsilon ^3} & \frac{1}{2}>x>\frac{1}{2}-\epsilon \\
1 & x\geq \frac{1}{2}
\end{cases}$$
for some $0<\epsilon<\frac{1}{2}$. Now, it is rather difficult to show but this function fits the conditions presented above. Further, we have
$$\int_0^1(f(f(x))-f(x))dx=\frac{1}{4}-\frac{3}{4}\epsilon+\frac{1}{2}\epsilon^2 $$
Clearly, we can get within any $\epsilon>0$ of $\frac{1}{4}$ as we desire.
Now, my method to find this function:
$1)$ I started under the assumption that $f(x)$ was a piecewise function
$$f(x)=\begin{cases}
a_1 & x< b \\
a_2 & x>b
\end{cases}$$
(ignoring what happens at $x=b$).
$2)$ I then plugged this function into mathematica and got the result
$$\int_0^1(f(f(x))-f(x))dx=(a_2-a_1)b$$
when $b<a_1$. For other possible cases, I got a negative value and since $f(x)=0$ gives an integral of $0$ I knew that these cases could be discarded.
$3)$ I looked at this and thought about how to maximize it. Obviously, $a_2$ should be set to $1$. This then gave me
$$=(1-a_1)b$$
when $a_1>b$.
$4)$ I then set $a_1=b+\epsilon$ to get
$$=(1-b-\epsilon)b$$
This function is maximized at $b\approx \frac{1}{2}$. Thus, our original function becomes
$$f(x)=\begin{cases}
\frac{1}{2}+\epsilon & x< \frac{1}{2} \\
1 & x>\frac{1}{2}
\end{cases}$$
$5)$ To make our function fit, I used the cubic
$$h(x)=sx^3+tx^2+ux+v$$
and set
$$f(x)=\begin{cases}
\epsilon +\frac{1}{2} & x\leq \frac{1}{2}-\epsilon \\
h(x) & \frac{1}{2}>x>\frac{1}{2}-\epsilon \\
1 & x\geq \frac{1}{2}
\end{cases}$$
All that was left was to solve for $s,t,u,v$ such that $f(x)$ met all the conditions above. I then plugged $f(x)$ back into mathematica to get the answer above.
Commentary: The obvious next step is to try this process with a step function of three levels rather than two. However, this already seems extremely difficult and I fear some clever idea is needed in order to make this work. The hope eventually would be to extend this to a step function of $n$ steps, and then in the limit we would hopefully have a function which maximizes the integral.
EDIT: I tried to extend the method above to a step function with three steps, but I kept running up on the $\frac{1}{4}$ barrier. Don't know if that means it is impossible to get past that, but I'm beginning to think that may be the upper limit. Especially since the linked post has answers which have $\frac{n^2}{4}$ in them as a sort of limiting factor. That is in the limit
$$\lim_{n\to\infty}\frac{n^2}{4}\cdot \frac{1}{g(n)}=1$$
and
$$g(n)<\frac{n^2}{4}$$
(where in that question $g(n)$ was the maximum). This is kind of like the answer that I got.
EDIT 2: After condition $2$ was changed to monotonicly increasing, there exists a function $f(x)$ such that the integral is exactly $\frac{1}{4}$. For
$$f(x)=\begin{cases}
\frac{1}{2} & x< \frac{1}{2} \\
1 & x\geq \frac{1}{2}
\end{cases}$$
we have
$$\int_0^1(f(f(x))-f(x))dx=\frac{1}{4}$$
