A triple integral involving $\text{Li}_2$ I am grateful to Cornel Ioan Valean for sharing this nice problem.

Find a closed form in terms of the values of the $\zeta$ function for the following triple integral:
  $$ \iiint_{(0,1)^3}\frac{\text{Li}_2(1-xyz)}{1-xyz}\,dx\,dy\,dz.$$

My approach has been the following one


*

*By exploiting the dilogarithm reflection formula, it is enough to compute the above integral where $\text{Li}_2(1-xyz)$ is replaced by $1,\log(xyz)$ and $\text{Li}_2(xyz)$;

*By Fubini's theorem these integrals boil down to Euler sums of weight $5$, which can be computed in terms of values of the $\zeta$ function through Theorems $2.2$ and $3.1$ of Flajolet-Salvy.


I am looking for an alternative/more elementary approach.
 A: We use standard notation. 
\begin{eqnarray}
\int\limits_{(0,1)^3} \frac{Li_2(1-x y z)}{1-x y z} dx dy dz =\\ \int\limits_{(0,1)^2} \frac{\zeta(3)-Li_3(1- y z )}{y z} dy dz=\\\int\limits_0^1 \frac{-\log(z)\left(Li_3(1-z)-\zeta(3)\right)-1/2\left(Li_2(1-z)^2-\zeta(2)^2\right)}{z}dz=\\
\frac{1}{2} \int\limits_0^1 \log(z)^2 \frac{Li_2(1-z)}{1-z} dz=\\
\frac{1}{2} \int\limits_0^1 \log(z)^2 \frac{\left(\frac{\pi^2}{6} - \log(z)\log(1-z)-Li_2(z) \right)}{1-z} dz=\\
\frac{1}{2}\left(2 \zeta(2)\zeta(3)+\int\limits_0^1 \log(z)^2 \frac{\left( -\log(z) \log(1-z)-Li_2(z)\right)}{1-z} dz\right)=\\
\frac{1}{2}\left(2 \zeta(2)\zeta(3)-\frac{1}{2}\left(-\Psi^{(4)}(1) + 6 \Psi^{(1)}(1) \Psi^{(2)}(1)\right)+\int\limits_0^1 \log(z)^2 \frac{\left(-Li_2(z)\right)}{1-z} dz\right)=\\
\frac{1}{2}\left(2 \zeta(2)\zeta(3)-\frac{1}{2}\left(-\Psi^{(4)}(1) + 6 \Psi^{(1)}(1) \Psi^{(2)}(1)\right)-\int\limits_0^1 \left( \frac{Li_2(z)^2}{z}-\frac{Li_1(z)^2}{z} \log(z)^2\right)dz\right)=\\
\frac{1}{2}\left(2 \zeta(2)\zeta(3)-\frac{1}{2}\left(-\Psi^{(4)}(1) + 6 \Psi^{(1)}(1) \Psi^{(2)}(1)\right)-\left(2{\bf H}^{(2)}_3(1)-2 \zeta(5)\right)\right)=\\
\frac{1}{2}\left(2 \zeta(2)\zeta(3)-\frac{1}{2}\left(-\Psi^{(4)}(1) + 6 \Psi^{(1)}(1) \Psi^{(2)}(1)\right)-\left(6 \zeta(2) \zeta(3)-11 \zeta(5)\right)\right)=\\
\zeta(2) \zeta(3)-1/2 \zeta(5)
\end{eqnarray}
Here ${\bf H}^{(p)}_q(x):=\sum\limits_{m=1}^\infty H_m^{(p)}/m^q x^m$ and $\Psi^{(q)}(\cdot)$ is the polygamma function.
Update: Let $d\ge 1$ and $q\ge 1$. As a matter of fact the following result holds:
\begin{equation}
{\mathcal I}_q^{(d)}:=\int\limits_{(0,1)^d} \frac{L_q(1-\prod\limits_{\xi=1}^d x_\xi)}{1-\prod\limits_{\xi=1}^d x_\xi} \prod\limits_{\xi=1}^d d x_\xi = \frac{1}{(d-1)!} \int\limits_0^1 [\log(1/z)]^{d-1} \frac{Li_q(1-z)}{1-z} dz
\end{equation}
which is readily seen by introducing a unity $1= \int\limits_{\mathbb R} \delta(z- \prod\limits_{\xi=1}^d x_\xi) dz$ into the integrand on the left hand side and then integrating over the $x$-variables. Now, by using the reflection formula for the di-logarithm we readily derive the following result:
\begin{equation}
{\mathcal I}_2^{(d)} = \zeta(d+2)+\zeta(2)\zeta(d)+1/2 \frac{(-1)^d}{(d-1)!} \left(-\Psi^{(d+1)}(1)+\sum\limits_{j=1}^{d-1} \binom{d}{j} \Psi^{(j)}(1) \Psi^{(d-j)}(1) \right)- {\bf H}^{(2)}_d(1)
\end{equation}
for $d\ge 3$.
Unfortunately if $q>2$ it is not quite clear for me how to do the integral since, according to my knowledge, there are no reflection formulae for polylogarithms of order strictly greater than two.
