After looking at Antonio's comment, I thought a little bit more about how to find the expansion of my integral. I don't know if it's in poor taste to answer one's own question, but it was interesting enough, to me anyway, that I thought it was worth posting.
First, as Antonio points out it is sufficient to only consider $f(x)=1$.
To expand $I_1(\epsilon)$ let $x=\epsilon e^{-t}$ and define $L:= \log(1/\epsilon)$, after the dust clears the integral becomes
$$
I_1(\epsilon) = \int_{-L}^\infty \frac{ L+t }{\sqrt{1+e^{t}}} dt.
$$
The singular behavior comes from the contribution to this integral near $t=-\infty$. We can subtract off the limiting divergent behavior to write
\begin{align*}
I_1(\epsilon) &= \frac{1}{2}L^2 + \int_{-L}^0 \frac{ L+t }{\sqrt{1+e^{t}}} - (L+t) dt
+ \int_0^\infty \frac{ L+t }{\sqrt{1+e^{t}}} dt \\
&= \frac{1}{2}L^2 + L \int_{-\infty}^{\infty} \left( \frac{1}{\sqrt{1+e^{t}}} -
\chi_{\mathbb{R}_-}(t) \right) dt
+ \int_{-\infty}^\infty \left( \frac{t}{\sqrt{1+e^{t}}} -
t \chi_{\mathbb{R}_-}(t) \right) dt \\
&\qquad- \int_{-\infty}^{-L} \frac{ L+t }{\sqrt{1+e^{t}}} - (L+t) dt.
\end{align*}
Here $\chi_{\mathbb{R}_-}(t)$ is the indicator function of the negative half line:
$$
\chi_{\mathbb{R}_-}(t) = \begin{cases} 1 & x<0 \\ 0 & x>0. \end{cases}
$$
Each of the integrals is convergent. The first two are $\epsilon$ independent so that the formula really is an asymptotic series, and they can be computed numerically or using Mathematica (I simplified further using dilogarithm identities which Mathematica didn't use). The last integral is vanishingly small as can be verified by expanding the denominator in powers of $e^t$ and integrating term by term. The result is
\begin{align*}
I_1(\epsilon) &= \frac{1}{2} \log(1/\epsilon)^2 + 2\log(2) \log(1/\epsilon) + \frac{\pi^2}{3} + 2 \log^2 (2) + \mathcal{O}({\epsilon \log \epsilon}) \\
&= \frac{1}{2} \log(4/\epsilon)^2 + \frac{\pi^2}{3} + \mathcal{O}({\epsilon \log \epsilon}).
\end{align*}
For general $f(x)$ using Antonio's comment the expansion is given by
$$
I_f(\epsilon) = f(0) \left[\frac{1}{2} \log(4/\epsilon)^2 + \frac{\pi^2}{3} \right] +
\int_0^1 \log(1/x) \frac{ f(x) - f(0) }{x} dx + o(1).
$$