Show that $\int_0^1 4 \space\operatorname{li}(x)^3 \space (x-1) \space x^{-3} dx = \zeta(3) $
This solution reduces the expression to an integral that Mathematica knows how to solve. Use a shift in the integrand of Gradshteyn & Rhyzhik 4.351.3 and rewrite the Ei function in terms of li: $$ \frac{1}{2} \Big( \frac{li(x)}{x} \Big)^2 = \int_0^\infty x^t \log{(1+t)}\frac{dt}{2+t} .$$ Use the expansion (B) as given by ComplexYetTrivial. Interchange sums and integral to find $$\int_0^1 \Big( \frac{li(x)}{x} \Big)^2 \Big( \frac{li(x)}{x} \Big) (x-1) dx= 2\int_0^\infty dt \, \frac{\log{(1+t)}}{2+t}\sum_{k=0}^\infty\frac{1}{k+1}\int_0^1 x^t L_k(-\log{x})(x-1)dx $$ Use the well-known evaluation, below, twice $$\int_0^1 x^t L_k(-\log{x}) = t^k/(1+t)^{k+1}$$ to get $$\int_0^1 \Big( \frac{li(x)}{x} \Big)^3(x-1)dx = 2\int_0^\infty dt \, \frac{\log{(1+t)}}{2+t} \sum_{k=0}^\infty\frac{1}{k+1}\Big(\frac{1}{t+1} \big(\frac{t+1}{t+2}\big)^{k+1} - \frac{1}{t} \big(\frac{t}{t+1}\big)^{k+1} \Big)$$ The sums have closed-forms in terms of log. Collect and simplify the integrand to get $$\int_0^1 \Big( \frac{li(x)}{x} \Big)^3(x-1)dx = 2\int_0^\infty dt \, \frac{\log{(1+t)}}{2+t}\Big(\frac{\log{(t+1)}}{t} -\frac{\log{(t+2)}}{t+1}\Big)=\frac{\zeta(3)}{4} $$ where the single integral has been performed by Mathematica.
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.