This is not a complete answer, but just a description of two ideas that might help with the evaluation of the integral
$$ I \equiv 4 \int \limits_0^1 \left(\frac{\operatorname{li}(x)}{x}\right)^3 (x-1) \, \mathrm{d} x \, . $$
They are based on methods that can be applied to find the easier integral
$$ J \equiv \int \limits_0^1 \left(\frac{\operatorname{li}(x)}{x}\right)^2 \, \mathrm{d} x \, . $$
The first approach relies on integration by parts and the series
$$ x-1 = \sum \limits_{k=1}^\infty \frac{1}{k!} \ln^k (x) \, , \, x > 0 \, .$$
In order to evaluate $J$ we can use the antiderivative $x \mapsto 1-\frac{1}{x}$ of $x \mapsto \frac{1}{x^2}$ to avoid problems with the singularity of $\operatorname{li}(x)$ at $x = 1$ . We get
\begin{align}
J &= 2 \int \limits_0^1 \frac{\operatorname{li}(x)}{x} \frac{1-x}{\ln(x)} \, \mathrm{d} x = - 2 \sum \limits_{k=1}^\infty \frac{1}{k!} \int \limits_0^1 \frac{\operatorname{li}(x)}{x} \ln^{k-1} (x) \, \mathrm{d} x\\
&= 2 \sum \limits_{k=1}^\infty \frac{1}{k! k} \int \limits_0^1 \ln^{k-1} (x) \, \mathrm{d} x = 2 \sum \limits_{k=1}^\infty \frac{1}{k! k} (-1)^{k-1} (k-1)! \\
&= 2 \sum \limits_{k=1}^\infty \frac{(-1)^{k-1}}{k^2} = 2 \eta (2) = \zeta(2) = \frac{\pi^2}{6} \, .
\end{align}
Similarly, we can use the antiderivative $x \mapsto \frac{(x-1)^2}{2 x^2}$ of $x \mapsto \frac{x-1}{x^3}$ to find
\begin{align}
I &= - \frac{3}{2} \int \limits_0^1 \left(\frac{\operatorname{li}(x)}{x}\right)^2 \frac{(x-1)^2}{\ln(x)} \, \mathrm{d} x \\
&= \frac{3}{2} \sum \limits_{k=0}^\infty \frac{1}{k!} \int \limits_0^1 \operatorname{li}^2 (x) \frac{1-x}{x} \frac{\ln^{k-1} (x)}{x} \, \mathrm{d} x \, . \tag{A}
\end{align}
We can now integrate by parts once more to obtain at least one term that reduces to a multiple of $\zeta(3)$ as in the simpler case. However, I have not managed to solve the remaining integrals yet. We could of course use the series again to express the remaining $1-x$ in terms of logarithm powers, but that does not seem to solve the problem.
The second suggestion employs the Fourier-Laguerre series
$$ \operatorname{li} (x) = - x \sum_{n=0}^\infty \frac{\mathrm{L}_n (-\ln(x))}{n+1} \, , \, x \in (0,1) \, , \tag{B}$$
of the logarithmic integral. It can be proved by deriving a recurrence relation for the coefficients from that of the Laguerre polynomials.
Using the substitution $x = \mathrm{e}^{-t}$ and the orthogonality relation of the Laguerre polynomials we immediately obtain
$$ J = \sum \limits_{p=0}^\infty \sum \limits_{q=0}^\infty \frac{1}{(p+1)(q+1)} \int \limits_0^\infty \mathrm{L}_p (t) \mathrm{L}_q (t) \mathrm{e}^{-t} \, \mathrm{d} t = \sum \limits_{p=0}^\infty \frac{1}{(p+1)^2} = \zeta(2) = \frac{\pi^2}{6} \, .$$
Similarly, we have
$$ I = 4\sum \limits_{p=0}^\infty \sum \limits_{q=0}^\infty \sum \limits_{r=0}^\infty \frac{1}{(p+1)(q+1)(r+1)} \int \limits_0^\infty \mathrm{L}_p (t) \mathrm{L}_q (t) \mathrm{L}_r (t) (1- \mathrm{e}^{-t}) \mathrm{e}^{-t} \, \mathrm{d} t \, .$$
General formulas for integrals involving three Laguerre polynomials appear to be known (see this paper or this one for a generalisation). I do not know whether they are nice enough to reduce the triple series to a representation of $\zeta(3)$ though.
Remark: After doing some numerical calculations I now suspect that the triple series diverges. This is probably due to the fact the original series $(\mathrm{B})$ only converges in $L^2$, so it cannot be used here. For the simpler integral everything works out though.
It is of course possible to combine the two methods by applying the Laguerre series $(\mathrm{B})$ in equation $(\mathrm{A})$. I do not know if these ideas can be used to get the final result, but maybe they can help someone else to find a way.
