Recently, in the context of https://math.stackexchange.com/a/3392284/198592 , I stumbled over the integral
$$f(n) = \int_0^{\infty } \left(\frac{1-e^{-q}}{q}\right)^n \, dq\tag{1}$$
and wanted to find a closed expession for it.
The integral is convergent $n\gt 1$. Indeed, close to $q=0$ the integrand behaves as $1-\frac{q}{2}+\frac{q^2}{6}+O(q^3)$ so that there is no songularity there. For $q\to \infty$ the integrand becomes $\frac{1}{q^n}$. Hence the integral in convergent under the conditions stated.
Here we confine ourselves to the case of integer powers, i.e. $n=2, 3, ...$.
The first few values can be easily calculated
$$f(n=2..5) = \left\{\log (4),\log \left(\frac{81 \sqrt{3}}{64}\right),\frac{88 \log (2)}{3}-18 \log (3),\frac{5}{24} (-544 \log (2)+162 \log (3)+125 \log (5))\right\}\tag{2}$$
A general formula seemed hard to guess from these cases.
Alternatively, a direct attack on the integral seems to suffer from the singular negative power of $q$ after a binomial expansion of the integrand.
I suggest you try for yourself to find it.
Here is my result
$$f(n) = \frac{1}{(n-1)!}\sum _{k=1}^n (-1)^{k+n} k^{n-1} \log (k) \binom{n}{k}\tag{3}$$
I did not know it in advance. For instance I haven't found it in Gradshteiyn/Ryshik. (EDIT After finishing the OP I found that 3.411.19 is related).
But now we have it it should not be difficult for the reader to prove it.
Notice that
$$\frac{1-e^{-q}}{q} = \int_0^1 e^{-qx}dx$$
So for $n>1$ an integer, we can use Fubini's theorem to rewrite the integral in the following way:
$$\int_0^\infty \left(\int_0^1 e^{-qx}dx\right)^n dq = \int_{[0,1]^n} \int_0^\infty e^{-(x_1+\cdots+x_n)q}dq dx_1\cdots dx_n $$
$$= \int_{[0,1]^n}\frac{1}{x_1+\cdots+x_n}dx_1\cdots dx_n$$
which seems like the kind of integral we should ask a statistician friend about.
We can use $x_1+\cdots+x_n = n\left(1-\left(1-\frac{x_1+\cdots+x_n}{n}\right)\right)$ to get a geometric series in terms of $\left(1-\frac{x_1+\cdots+x_n}{n}\right)$. Using multinomial expansion to expand the geometric series powers and integrating the subsequent product of monomial terms leaves us with the following double sum:
$$\frac{1}{n}\sum_{k=0}^\infty \sum_{m_0+\cdots +m_n=k} {k \choose m_0,\cdots,m_n} \left(-\frac{1}{n}\right)^{k-m_0} \prod_{j=1}^n \frac{1}{m_j+1}$$