Use of AccuracyGoal & MachinePrecision in NIntegrate

Question

I just learned about the Inverse Symbolic Calculator and it seems like a very useful tool. For example, to find an analytic solution to

$$ \int_0^{\infty } \frac{1-e^{-x} (x+1)}{\left(e^x-1\right) \left(e^{-x}+e^x\right) x} \, dx, $$

one can solve it numerically to obtain

0.156595806752698829513363962452

and submit this to the Calculator to obtain the analytic expression

$$ -\frac{\gamma }{2}+\frac{\pi }{8}-\frac{1}{4} 3 \log (2)+\frac{\log (\pi )}{2}. $$

(Source)

However, I'm having trouble reproducing this result, i.e., using NIntegrate with the right settings for WorkingPrecision, AccuracyGoal, and PrecisionGoal to obtain the numerical result above.

I eventually found that

NIntegrate[(1 - Exp[-x] (1 + x))/(x (Exp[x] - 1) (Exp[x] + Exp[-x])), 
    {x, 0, Infinity}, WorkingPrecision -> 40, AccuracyGoal -> 25]

yielded

0.1565958067526988295133639624516335266315,

which the Calculator can't find an analytic form for, but if I chop off the last several digits and submit instead

0.1565958067526988295133639624516335

then the Calculator finds the correct analytic expression.

This took a lot of guess-and-check, and I would like to know if there is a more sure-fire way to get a bunch of digits out of NIntegrate to look up analytic expressions like this.

(I've read the docs on PrecisionGoal and AccuracyGoal as well as this MSE link.)