I'm just finishing the answer by Marco Cantarini, with $$\newcommand{\li}{\operatorname{li}}\newcommand{\dilog}{\operatorname{Li}_2}\int_1^\infty\li(x)x^{s-1}\,dx=\frac{\log(-1-s)}{s}$$ (for $\Re s<-1$) touched in the past (which calls for the Mellin transform of $\li(x)^2$ via multiplication/convolution). It turns out that the latter has a closed form in terms of the dilogarithm function (which suggests doing integrals with $\li(x)^3$ or $\li(x)^4$ the same way).
As a result, for $a\in\mathbb{C}$ and $0<c<\Re a$ we have $$\int_1^\infty\frac{\li(x)^2}{x^{3+a}}\,dx=\frac1{2\pi i}\int_{c-i\infty}^{c+i\infty}\frac{\log s}{1+s}\frac{\log(a-s)}{1+a-s}\,ds.$$
The integral above is transformed into a principal-value integral along the negative real axis: $$\int_1^\infty\frac{\li(x)^2}{x^{3+a}}\,dx=\text{p.v.}\int_0^\infty\frac{\log(a+x)}{1+a+x}\frac{dx}{x-1}.$$ This is done via the Cauchy integral theorem applied to the boundary of $$\{s\in\mathbb{C}:r<|s|<R\land\Re s<c\land(\Re s>0\lor|\Im s|>r)\}$$ (the part of $|s|<R$ to the left of $\Re s=c$, with a notch around the negative real axis), and taking $R\to\infty$ and $r\to 0$. The latter should be done with care, since $s=-1$ is a singularity on the branch cut of the integrand, but its contributions (at the edges of the cut) cancel each other.
A straightforward (but tedious) way to get a closed form of the last integral (in terms of the dilogarithm) is to write $\text{p.v.}\int_0^\infty=\lim\limits_{R\to\infty}\lim\limits_{r\to 0}\left(\int_0^{1-r}+\int_{1+r}^R\right)$ and evaluate the integrals individually. After simplification using dilogarithm identities, if we put $a=s-1$, we get $$f(s):=\int_1^\infty\frac{\li(x)^2}{x^{2+s}}\,dx=\frac1{1+s}\left[\frac{\pi^2}6+\log(s)^2+2\dilog\left(\frac1s\right)\right]$$ for $\Re s\geqslant 1$. The values of $\dilog(1)$ and $\dilog(1/2)$ yield $f(1)=\pi^2/4$ and $f(2)=\pi^2/9$.
Finally, the value of the given integral is $f(1)-f(2)=5\pi^2/36$ as expected.