Here is a first solution, which uses no complex variable nor any differentiation, just plain real analysis. A second solution follows, which uses Fubini theorem.
For every positive $t$, let $I_t(a,b)$ denote the integral of the same function than in the integral which defines $I(a,b)$, but integrated from $t$ to $+\infty$. Then $ I_t(a,b)=\int_{t}^{+\infty}\frac{\mathrm{d}x}{x\left(\mathrm{e}^{bx}+1\right)}-\int_{t}^{+\infty}\frac{\mathrm{d}x}{x\left(\mathrm{e}^{ax}+1\right)}= \int_{bt}^{at}\frac{\mathrm{d}x}{x\left(\mathrm{e}^{x}+1\right)}, $ where the first decomposition of $I_t(a,b)$ into two parts is legal because we got rid of the problem at $x=0$ and the second formula stems from the changes of variables $x\leftarrow bx$ and $x\leftarrow ax$ in the two parts.
Extracting $1/(2x)$ from the last fraction and integrating, one gets $ I_t(a,b)=\frac12\log\left(\frac{a}b\right)- \int_{bt}^{at}\frac{\mathrm{e}^{x}-1}{2x\left(\mathrm{e}^{x}+1\right)}\mathrm{d}x=\frac12\log\left(\frac{a}b\right)+O(t), $ where the $O(t)$ term when $t\to0$ comes from the fact that the function in the last integral is bounded around $x=0$ and from the fact that this integral is over an interval of length $at-bt=O(t)$. Finally, when $t\to0$, $I_t(a,b)\to I(a,b)$ hence $ I(a,b)=\frac12\log\left(\frac{a}b\right). $
Here is a second solution, based on Fubini theorem. Introduce $ u(c,x)=\frac{\mathrm{e}^{cx}}{(1+\mathrm{e}^{cx})^2}=\frac{\partial}{\partial c}\left(\frac{-1}{x(1+\mathrm{e}^{cx})}\right). $ One sees that $ \frac{\mathrm{e}^{ax}-\mathrm{e}^{bx}}{x\left(\mathrm{e}^{ax}+1\right)\left(\mathrm{e}^{bx}+1\right)}=\frac1{x\left(\mathrm{e}^{bx}+1\right)}-\frac1{x\left(\mathrm{e}^{ax}+1\right)}=\int_b^au(c,x)\mathrm{d}c, $ hence $ I(a,b)=\int_{0}^{\infty}\left(\int_b^au(c,x)\mathrm{d}c\right)\mathrm{d}x. $ Since $u\ge0$, Fubini-Tonelli theorem says one can interchange the order of integration, hence $ I(a,b)=\int_b^aU(c)\mathrm{d}c,\qquad U(c)=\int_{0}^{\infty}u(c,x)\mathrm{d}x. $ Obviously, $ u(c,x)=\frac{\partial}{\partial x}\left(\frac{-1}{c(1+\mathrm{e}^{cx})}\right),\qquad U(c)=\left.\frac{-1}{c(1+\mathrm{e}^{cx})}\right|^{x=+\infty}_{x=0}=\frac1{2c}, $ and the value of $I(a,b)$ follows.