Here's my take on this integral.
$$
I(s, z) = \int_0^\infty t^s (e^t - 1)^{-z}dt
$$
To start off, we'll use the following substitution.
$$u = e^{-t}\\ t = -\ln u \\ dt = -u^{-1} du$$
Which then we'll have
$$
-\int_1^0 (-\ln u)^s (u^{-1} - 1)^{-z} u^{-1}du \\
= (-1)^s \int_0^1 \ln^s u (u^{-1} - 1)^{-z} u^{-1} du \\
= (-1)^s \int_0^1 \ln^s u [u^{-1}(1 - u)]^{-z} u^{-1} du \\
= (-1)^s \int_0^1 \ln^s u (1 - u)^{-z} u^{z-1} du
$$
Now let's take a look at something similar, namely the Beta function, which is defined as the following integral.
$$
\text{B}(x, y) = \int_0^1 t^{x-1} (1-t)^{y-1} dt
$$
What I noticed is that if we take the function's partial derivative with respect to $x$, we'll get
$$
\partial_x \text{B}(x, y) = \partial_x \int_0^1 t^{x-1} (1-t)^{y-1} dt
$$
Using the Leibniz integral rule, we'll get
$$
\partial_x \text{B}(x, y) = \int_0^1 \partial_x t^{x-1} (1-t)^{y-1} dt = \int_0^1 t^{x-1} \ln t (1-t)^{y-1} dt
$$
Similarly, if we take the second order derivative, we'll get
$$
\partial_x^2 \text{B}(x, y) = \int_0^1 t^{x-1} \ln^2 t (1-t)^{y-1} dt
$$
And if we continue this to the $s$-th order derivative, we'll get
$$
\partial_x^s \text{B}(x, y) = \int_0^1 t^{x-1} \ln^s t (1-t)^{y-1} dt
$$
Plug in $x=z$, $y=1-z$ and multiply the whole thing with $(-1)^s$, we'll get back our original integral.
$$
(-1)^s \partial_x^s \text{B}(x, 1-z) \Bigg|_{x = z} = (-1)^s \int_0^1 t^{z-1} \ln^s t (1-t)^{-z} dt = I(s, z)
$$
And for simplicity in the future, I'll use the result from this to turn the Beta function into a series for the sake of ease in the future.
$$
(-1)^s \partial_x^s \sum_{k=0}^\infty \binom{k+z-1}{k} (x+k)^{-1} \Bigg|_{x = z}
$$
The condition for this function (which is also our integral) to converge is for
$$
\begin{cases}
\Re(x) > 0 \\
\Re(y) > 0
\end{cases}
$$
Which means we'll have
$$
\begin{cases}
\Re(z) > 0 \\
\Re(1-z) > 0
\end{cases}
\Rightarrow
\begin{cases}
\Re(z) > 0 \\
\Re(z) < 1
\end{cases}
\Rightarrow 0 < \Re(z) < 1
$$
For the $s$-th order derivative, I noticed that
$$
\partial_x (x+k)^{-1} = -1!(x+k)^{-2} \\
\partial_x^2 (x+k)^{-1} = 2!(x+k)^{-3} \\
\partial_x^3 (x+k)^{-1} = -3!(x+k)^{-4} \\
\partial_x^4 (x+k)^{-1} = 4!(x+k)^{-5}
$$
And if we continue this to the $s$-th order, we'd get
$$
\partial_x^s (x+k)^{-1} = (-1)^s s! (x+k)^{-(s + 1)}
$$
But this only works for natural $s$. However, there's a way around that. By using the Gamma function, we can extend the formula to the real numbers.
$$
\partial_x^s (x+k)^{-1} = (-1)^s \Gamma(s+1) (x+k)^{-(s + 1)} \\
$$
Plugging this back into our original integral, we'll get
$$
(-1)^s \sum_{k=0}^\infty \binom{k+z-1}{k} (-1)^s \Gamma(s+1) (x+k)^{-(s+1)} \Bigg|_{x = z} \\
= (-1)^{2s} \Gamma(s+1) \sum_{k=0}^\infty \binom{k+z-1}{k} (z+k)^{-(s+1)}
$$
There's one problem though, this whole thing diverges for all negative integer $s$. I tried some other ways to evaluate for that case but nothing worked, so we'll just ignore that case in the mean time. I may or may not update this answer.
Now, if you just stopped there, you'd have a good answer for real $s$. But since you asked for complex $s$, my idea is simply to just use the same answer, but instead for real $s$ only, you also use it for complex $s$
Which means our final result for real $s$ is
$$
I(s, z) = (-1)^{2s} \Gamma(s+1) \sum_{k=0}^\infty \binom{k+z-1}{k} (z+k)^{-(s + 1)} \\
\text{for $s \in \mathbb{R}\backslash\mathbb{Z}^-$, $z \in \mathbb{C}$ and $0 < \Re(z) < 1$}
$$
And our final result for complex $s$ is
$$
I(s, z) = (-1)^{2s} \Gamma(s+1) \sum_{k=0}^\infty \binom{k+z-1}{k} (z+k)^{-(s + 1)} \\
\text{for $s \in \mathbb{C}$, $\Re(s) \not\in \mathbb{Z}^-$, $z \in \mathbb{C}$ and $0 < \Re(z) < 1$}
$$
$\int_0^\infty \frac{t^{s-1}}{(e^t-1)^z}dt = \int_1^\infty \frac{t^{s-1}}{(e^t-1)^z}dt
$ where the RHS extends meromorphically to $\Re(s-z) >-K,\Re(z) >0$
– reuns Jan 22 '23 at 22:56