Compute an integral containing a product of powers of logarithms.

Let $q \ge 1$ and $p \ge 0$ be integers. Consider a following integral: \begin{equation} {\mathcal I}^{(q,p)} := \int\limits_0^1 [\log(1-\eta)]^q [\log(\eta)]^p \frac{1}{\eta} d \eta \end{equation} Clearly that integral is proportional to the Nielsen generalized polylogarithm at unity. Now by using Euler's beta function integral it is easy to see that : \begin{equation} {\mathcal I}^{(q,p)} = \frac{\partial^p}{\partial \theta_1^p} \frac{\partial^q}{\partial \theta_2^q} \left. \left[ \frac{(\theta_1-1)! \theta_2!}{(\theta_1+\theta_2)!} \right] \right|_{\theta_1=\theta_2=0} \end{equation} We have computed the derivative with respect to $\theta_2$ using the Faa di Bruno's formula then we have set $\theta_2=0$ then we differentiated the result $p$ times using Mathematica and finally set $\theta_1=0$. As a result we discovered the following relations: \begin{eqnarray} &&1!\cdot{\mathcal I}^{(q,0)} = - \Psi^{(q)}(1) \\ &&2! \cdot{\mathcal I}^{(q,1)} = - \Psi^{(q+1)}(1) + \sum\limits_{j=1}^{q-1} \binom{q}{j} \Psi^{(j)}(1) \Psi^{(q-j)}(1) \\ &&3! \cdot {\mathcal I}^{(q,2)} = -2 \Psi^{(q+2)}(1) + 3 \cdot 1_{q\ge 2}\cdot\sum\limits_{j=1}^{q-1} \binom{q}{j}\cdot\left[ \Psi^{(j+1)}(1) \Psi^{(q-j)}(1)+\Psi^{(j)}(1) \Psi^{(q+1-j)}(1) \right]+ \\ &&-2 \cdot 1_{q \ge 3}\cdot\sum\limits_{1 \le j < j_1 \le q-1} \binom{q}{j,j_1-j,q-j_1} \Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(q-j_1)}(1)\\ &&4! \cdot {\mathcal I}^{(q,3)} = -6 \Psi^{(q+3)}(1)+\\ &&12 \cdot\sum\limits_{j=1}^{q-1} \binom{q}{j} \left[ \Psi^{(j)}(1) \Psi^{(q-j+2)}(1) + \frac{3}{2} \Psi^{(j+1)}(1) \Psi^{(q-j+1)}(1)+\Psi^{(j+2)}(1) \Psi^{(q-j+0)}(1)\right]+\\ &&-12 \cdot \sum\limits_{1 \le j < j_1 \le q-1} \binom{q}{j,j_1-j,q-j_1} \left[\Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(q-j_1+1)}(1)+\Psi^{(j)}(1) \Psi^{(j_1-j+1)}(1) \Psi^{(q-j_1)}(1)+\Psi^{(j+1)}(1) \Psi^{(j_1-j)}(1) \Psi^{(q-j_1)}(1)\right]+\\ &&6\cdot \sum\limits_{1\le j < j_1 < j_2 \le q-1} \binom{q}{j,j_1-j,j_2-j_1,q-j_2} \Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(j_2-j_1)}(1) \Psi^{(q-j_2)}(1) \end{eqnarray} where $\Psi^{(j)}(1)$ is the polygamma function at unity. Now the question would be how do we find the result for $p \ge 3$? The multitude of terms that appear in the Faa di Bruno formula is difficult to deal with. Is there a more elegant way of arriving at the result?

We have $$ \log^q(1-\eta) = q!\sum_{n\geq q}(-1)^q{\,n\, \brack q}\frac{\eta^n}{n!}\tag{1} $$ hence $$ \mathcal{I}^{(q,p)}=p!q!\sum_{n\geq q}\frac{(-1)^{p+q}}{n!\,n^{p+1}}{\,n\,\brack q}\tag{2} $$ and the problem boils down to computing some Euler sums, once the Stirling numbers of the first kind are converted into combinations of generalized harmonic numbers. In this context values of $p$ or $q$ greater than $3$ lead to intractable problems by hand: that is a good moment for invoking the help of a CAS, without shame.