The following integral $$ \int_0^1 e^{a(x-1)}\ln(-\ln(x)) \, dx$$ looks like the one for which the Laplace method is applicable. The function under the integral looks like this:
($a$ increases from red to black), so the leading contribution to the integral should be at the points 1 and 0. I can split this into two integrals from functions having maximum at 0, however direct application of the Laplace method is problematic since the maxima are infinite. Can it be fixed with some variable change, or do I need something beyond the Laplace method (steepest descents or anything else)?

Let
$$ f(x) = \int_0^1 e^{a(x-1)}\ln(-\ln(x))\,dx. $$
By looking at the graph of the integrand we can see that the largest contribution to the integral comes from a neighborhood of $x=1$. The integral away from the endpoint is exponentially small: if $0 < \epsilon < 1$ then
$$ \left| \int_0^\epsilon e^{a(x-1)}\ln(-\ln(x))\,dx \right| \leq e^{a(\epsilon-1)} \int_0^1 |\ln(-\ln(x))|\,dx. $$
Thus
$$ f(x) = \int_\epsilon^1 e^{a(x-1)}\ln(-\ln(x))\,dx + O\!\left(e^{a(\epsilon - 1)}\right). \tag{1} $$
Just as in the usual Laplace method we now approximate the subdominant factors of the integrand. Near $x=1$ we have
$$ -\ln(x) = \sum_{k=1}^{\infty} \frac{(1-x)^k}{k}, $$
so that
$$ \begin{align} \ln(-\ln(x)) &= \ln\!\left( (1-x) + \sum_{k=2}^{\infty} \frac{(1-x)^k}{k} \right) \\ &= \ln(1-x) + \ln\!\left( 1 + \sum_{k=1}^{\infty} \frac{(1-x)^k}{k+1} \right) \\ &= \ln(1-x) + \sum_{k=1}^{\infty} c_k (1-x)^k \end{align} $$
for some coefficients $\{c_k\}$. In particular, we can write
$$ \ln(-\ln(x)) = \ln(1-x) + \sum_{k=1}^{N} c_k (1-x)^k + R_N(x), $$
where there is a constant $C$ such that $|R_N(x)| \leq C(1-x)^{N+1}$ for all $x \in [\epsilon,1]$.
Substituting this into equation $(1)$ yields
$$ \begin{align} f(x) &= \int_\epsilon^1 e^{a(x-1)}\ln(1-x)\,dx + \sum_{k=1}^{N} c_k \int_{\epsilon}^{1} e^{a(x-1)} (1-x)^k\,dx \\ &\qquad + \int_{\epsilon}^{1} e^{a(x-1)} R_N(x)\,dx + O\!\left(e^{a(\epsilon - 1)}\right). \end{align} \tag{2} $$
We can easily estimate the integral containing the remainder term,
$$ \begin{align} \left|\int_{\epsilon}^{1} e^{a(x-1)} R_N(x)\,dx\right| &\leq \int_{\epsilon}^{1} e^{a(x-1)} |R_N(x)|\,dx \\ &\leq C\int_{\epsilon}^1 e^{a(x-1)}(1-x)^{N+1}\,dx \\ &= C \int_{0}^{1-\epsilon} e^{-ay} y^{N+1}\,dy \qquad\qquad [\text{substitute } y=1-x] \\ &< C \int_0^\infty e^{-ay} y^{N+1}\,dy \\ &= \frac{C \operatorname{\Gamma}(N+2)}{a^{N+2}}, \end{align} $$
so that $(2)$ becomes, after substituting $y = 1-x$,
$$ f(x) = \int_{0}^{1-\epsilon} e^{-ay}\ln(y)\,dy + \sum_{k=1}^{N} c_k \int_{0}^{1-\epsilon} e^{-ay} y^k\,dy + O\!\left(a^{-N-2}\right). \tag{3} $$
Arguing similarly, we can show that attaching the tails to the integrals in the sum $\sum_{k=1}^{N}$ introduces only exponentially small errors, so that
$$ f(x) = \int_{0}^{1-\epsilon} e^{-ay}\ln(y)\,dy + \sum_{k=1}^{N} \frac{c_k \operatorname{\Gamma}(k+1)}{a^{k+1}} + O\!\left(a^{-N-2}\right). \tag{4} $$
Finally we tackle the remaining integral. Starting with the substitution $s = ay$, we have
$$ \begin{align} \int_{0}^{1-\epsilon} e^{-ay}\ln(y)\,dy &= a^{-1} \int_0^{a(1-\epsilon)} e^{-s} \bigl[\ln(s) - \ln(a)\bigr]\,ds \\ &= a^{-1} \int_0^{a(1-\epsilon)} e^{-s}\ln(s)\,ds - a^{-1}\ln(a) \int_0^{a(1-\epsilon)} e^{-s}\,ds. \tag{5} \end{align} $$
Now, for $a$ large enough to make $a(1-\epsilon) > 1$,
$$ 0 < \int_{a(1-\epsilon)}^\infty e^{-s}\ln(s)\,ds < \int_{a(1-\epsilon)}^\infty e^{-s}s\,ds = e^{a(\epsilon-1)}(a(1-\epsilon)+1) $$
so that
$$ \begin{align} \int_0^{a(1-\epsilon)} e^{-s}\ln(s)\,ds &= \int_0^{\infty} e^{-s}\ln(s)\,ds - \int_{a(1-\epsilon)}^\infty e^{-s}\ln(s)\,ds \\ &= -\gamma + O\!\left(a e^{a(\epsilon - 1)}\right), \tag{6} \end{align} $$
where $\gamma$ is the Euler-Mascheroni constant. Similarly,
$$ \int_0^{a(1-\epsilon)} e^{-s}\,ds = 1 + O\!\left(e^{a(\epsilon - 1)}\right). \tag{7} $$
In light of $(6)$ and $(7)$, equation $(5)$ becomes
$$ \int_{0}^{1-\epsilon} e^{-ay}\ln(y)\,dy = -\frac{\ln(a) + \gamma}{a} + O\!\left(ae^{a(\epsilon - 1)}\right). \tag{8} $$
Substituting this into $(4)$, we conclude that an asymptotic series for $f(a)$ as $a \to \infty$ is
$$ f(x) \approx -\frac{\ln(a) + \gamma}{a} + \sum_{k=1}^{N} \frac{c_k \operatorname{\Gamma}(k+1)}{a^{k+1}}. $$
Using a CAS to compute the first few coefficients $\{c_k\}$, we find that
$$ f(x) = -\frac{\ln(a) + \gamma}{a} + \frac{1}{2a^2} + \frac{5}{12a^3} + \frac{3}{4a^4} + \frac{251}{120a^5} + \frac{95}{12a^6} + \cdots. $$