Remarkable logarithmic integral $\int_0^1 \frac{\log^2 (1-x) \log^2 x \log^3(1+x)}{x}dx$ We have the following result ($\text{Li}_{n}$ being the polylogarithm):
$$\tag{*}\small{ \int_0^1 \log^2 (1-x) \log^2 x \log^3(1+x) \frac{dx}{x} = -168 \text{Li}_5(\frac{1}{2}) \zeta (3)+96 \text{Li}_4(\frac{1}{2}){}^2-\frac{19}{15} \pi ^4 \text{Li}_4(\frac{1}{2})+\\ 12 \pi ^2 \text{Li}_6(\frac{1}{2})+8 \text{Li}_4(\frac{1}{2}) \log ^4(2)-2 \pi ^2 \text{Li}_4(\frac{1}{2}) \log ^2(2)+12 \pi ^2 \text{Li}_5(\frac{1}{2}) \log (2)+\frac{87 \pi ^2 \zeta (3)^2}{16}+\\ \frac{447 \zeta (3) \zeta (5)}{16}+\frac{7}{5} \zeta (3) \log ^5(2)-\frac{7}{12} \pi ^2 \zeta (3) \log ^3(2)-\frac{133}{120} \pi ^4 \zeta (3) \log (2)-\frac{\pi ^8}{9600}+\frac{\log ^8(2)}{6}- \\ \frac{1}{6} \pi ^2 \log ^6(2)-\frac{1}{90} \pi ^4 \log ^4(2)+\frac{19}{360} \pi ^6 \log ^2(2) }$$
This is extremely amazing: almost all other similar integrals are not expressible via ordinary polylogarithm.
The solution is however non-trivial. There are two methods: first is to find enough linear relations between similar integrals, once the rank is high enough, solving the system gives $(*)$; second method is to convert the integral into multiple zeta values, then use known linear relations between them. None of these methods can explain the result's simplicity.

Question: Is there a simpler method to prove (*), or a conceptual explanation of its elegance?

Any thought is welcomed. Thank you very much.

I wrote a Mathematica package, it can calculate the integral in subject and many similar ones. The following command calculates $(*)$:
MZIntegrate[Log[1-x]^2*Log[x]^2*Log[1+x]^3/x, {x,0,1}]

It can also solve some other integrals.
The package can be obtained here. I hope it can benefit those interested in related integral/series.

Remarks on the question:

*

*It's known that $\zeta(\bar{3},1,\bar{3},1)$ is
very reminiscent to the RHS of $(*)$. But both the simplicity of
$\zeta(\bar{3},1,\bar{3},1)$ and its connection to the integral are
elusive to me.

*(Added by Iridescent) This contains nearly all known general formulas
of these log integrals. However it does not help much on solving OP's
problem.

 A: Here are some ideas towards explaining the form of the right hand side. I'm a bit stuck and my main approach hasn't worked out. This may just be rephrasing things in terms of other log-integrals, but hopefully this is a useful way of looking at the problem.
Taking the integral
$$
I = \int_0^1 \log^2(1-x) \log^2(x) \log^3(1+x) \frac{dx}{x}
$$
we can also rewrite this as
$$
I = \int_0^\infty \log^2(1-e^{-x}) \log^2(e^{-x}) \log^3(1+e^{-x}) \; dx
$$
which is suited for interpretation as a Mellin transform. Specifically, the power of $x$, is controlled by the power on $\log(x)$ in the original integral format as
$$
I = \int_0^\infty x^2 \log^2(1-e^{-x})\log^3(1+e^{-x}) \; dx
$$
according to Mathematica we have in general a result for the Mellin transform of the other components
$$
\mathcal{M}[\log^n(1\pm e^{-x})](s) = (-1)^n n! \Gamma(s) S_{s,n}(\mp 1)
$$
invoking the Neilsen Generalisation of the polylogarithm, $S_{s,n}$. This does recreate the series expansion for $\log(1+e^{-x})$ but the series for $\log(1-e^{-x})$ has a $\log(x)$ term, which might be causing a problem.
We could toy with the idea of a formal series via the Ramanujan Master Theorem, using these Mellin transforms
$$
\log^n(1\pm e^{-x}) = \sum_{k=0}^\infty \frac{(-1)^{k+n} n!}{k!} S_{-k,n}(\mp 1)x^k
$$
and then the Cauchy product
$$
\log^a(1 + e^{-x})\log^b(1 - e^{-x}) = \left( \sum_{k=0}^\infty \frac{(-1)^{k+a} a!}{k!} S_{-k,a}(-1)x^k \right)\left( \sum_{k=0}^\infty \frac{(-1)^{k+b} b!}{k!} S_{-k,b}(1)x^k \right)
$$
$$
\log^a(1 + e^{-x})\log^b(1 - e^{-x}) = \sum_{k=0}^\infty \left(\sum_{l=0}^k \frac{(-1)^{a+b+k} a! b!}{l!(k-l)!} S_{-l,a}(-1) S_{l-k,b}(1)\right) x^k
$$
alternatively
$$
\log^a(1 + e^{-x})\log^b(1 - e^{-x}) = \sum_{k=0}^\infty \frac{(-1)^k}{k!} \left(\sum_{l=0}^k (-1)^{a+b} a! b! \binom{k}{l} S_{-l,a}(-1) S_{l-k,b}(1)\right) x^k
$$
plausibly leading to (via RMT)
$$
\mathcal{M}\left[ \log^a(1 + e^{-x})\log^b(1 - e^{-x})\right](s) = \Gamma(s) \sum_{l=0}^{-s} (-1)^{a+b} a! b! \binom{-s}{l} S_{-l,a}(-1) S_{l-k,b}(1)
$$
then we would conceptually have (with some dodgy negative parts) an answer for the integral as a sum over (four?) pairs of generalized Polylogs, specifically in the case that $s=3$.
This motivates an expression in terms of pairs of $S_{n,k}(z)$, we can guess a term and quickly find
$$
-8\cdot3 \cdot 19 S_{2,2}(1)S_{1,3}(-1) = -\frac{19}{15} \pi ^4 \text{Li}_4\left(\frac{1}{2}\right)-\frac{133}{120} \pi ^4 \zeta (3) \log (2)+\frac{19 \pi ^8}{1350}+\frac{19}{360} \pi ^6 \log ^2(2)-\frac{19}{360} \pi ^4 \log ^4(2)
$$
this covers a few of the terms in your expression R.H.S. It is likely that other terms contribute to $\pi^8$ for example. I can't get an explicit value for $S_{2,3}(-1)$ to explore this further, but I would assume this holds a $\mathrm{Li}_5(1/2)$ term among others, and the other factor is $S_{1,2}(1) = \zeta(3)$. Perhaps your linear combinations method can be rephrased in terms of the generalized polylogarithm?
