There appears to be a split between the case where $f$ and $h$ are real and when they can be complex. In particular, $S(g)$ and $I(g)$ will agree where they are well-defined for the real case, but there are counterexamples in the complex case.
I will use this answer on the moments of the absolute value of $h(x)$ to answer my own question for the real case: when $S(g)$ has a nonzero radius of convergence $R$, then for $|g| < R$ and $f(x)$ and $h(x)$ both real, then $S(g)$ and $I(g)$ match.
First we need a quick lemma. Note that by Cauchy-Schwarz $$b(n+1)^2 \equiv \left(\int_{-\infty}^{+\infty} |h(x)|^{n+1}e^{-f(x)}dx\right)^2 \leq \int_{-\infty}^{+\infty} |h(x)|^{n}e^{-f(x)}dx\cdot\int_{-\infty}^{+\infty} |h(x)|^{n+2}e^{-f(x)}dx$$
so $$b(n+1)^2 \leq b(n) b(n+2) \leq \max(b(n), b(n+2))^2,$$
so $$b(n+1) \leq \max(b(n), b(n+2)) \tag{1}\label{eq1}$$
Now consider $S(g)$. When $|g|<R$, $S(g)$ absolutely converges by the ratio test, so in this range,
$$\sum_{n=0}^\infty \frac{|g|^n}{n!} \left|\int_{-\infty}^{\infty}dx h(x)^n e^{-f(x)}\right| < \infty $$
It follows then that summing only the even terms
$$\sum_{n=0}^\infty \frac{|g|^{2n}}{2n!} \left|\int_{-\infty}^{\infty}dx h(x)^{2n} e^{-f(x)}\right| = \sum_{n=0}^\infty \frac{|g|^{2n}}{2n!} \int_{-\infty}^{\infty}dx |h(x)|^{2n} e^{-f(x)} < \infty \tag{2}\label{eq2}$$
is finite, where in the middle I used the reality of $h$ and $f$ to pull the absolute value sign into the integral.
Note that if we defined odd terms with absolute values inside the integral, the lemma $\eqref{eq1}$ yields the first inequality below, the second following from non-negativity of the integrals:
$$\sum_{n=0}^\infty \frac{|g|^{2n+1}}{(2n+1)!} \int_{-\infty}^{\infty}dx |h(x)|^{2n+1} e^{-f(x)} \leq \sum_{n=0}^\infty \frac{|g|^{2n+1}}{(2n+1)!} \max\left(\int_{-\infty}^{\infty}dx |h(x)|^{2n} e^{-f(x)}, \int_{-\infty}^{\infty}dx |h(x)|^{2n+2} e^{-f(x)}\right) \leq \sum_{n=0}^\infty \frac{|g|^{2n+1}}{(2n+1)!} \left(\int_{-\infty}^{\infty}dx |h(x)|^{2n} e^{-f(x)} + \int_{-\infty}^{\infty}dx |h(x)|^{2n+2} e^{-f(x)}\right)$$
The ratio test applied to the last sum above gives a finite answer by virtue of the limiting ratio being the same as the limiting ratio for the even only series in $\eqref{eq2}$. Thus,
$$\sum_{n=0}^\infty \frac{|g|^{2n+1}}{(2n+1)!} \int_{-\infty}^{\infty}dx |h(x)|^{2n+1} e^{-f(x)} < \infty \tag{3}\label{eq3}$$
Putting this all together by adding $\eqref{eq2}$ and $\eqref{eq3}$, for $|g|<R$, we have
$$\sum_{n=0}^\infty \frac{|g|^{n}}{n!} \int_{-\infty}^{\infty}dx |h(x)|^{n} e^{-f(x)} < \infty $$
However, this is a sufficient condition for Fubini's theorem to allow the swapping of the sum and integral sign for the original problem, so we have at last
$$S(g) \equiv \sum_{n=0}^\infty \frac{g^{n}}{n!} \int_{-\infty}^{\infty}dx h(x)^{n} e^{-f(x)} = \int_{-\infty}^{\infty}dx \sum_{n=0}^\infty \frac{g^{n}}{n!} h(x)^{n} e^{-f(x)} = \int_{-\infty}^{\infty}dx e^{-f(x)+g h(x)} \equiv I(g)$$
Thus, when $S(g)$ converges, it is equal to $I(g)$. More carefully stated, when $S(g)$ has a radius of convergence $R$, then $S(g)=I(g)$ when $|g| < R$.
The above assumed $f(x)$ and $h(x)$ real in e.g. $\eqref{eq2}$. The complex case allows for co-opting a famous example from probability as a clear example of $S(g) \neq I(g)$.
Consider the complex functions
$$f(x) = -x^{1/4}+\ln(\sin(x^{1/4})) \text{ for } x\geq 0$$
and
$$h(x) = ix$$
I will take $f(x)$ to be $-\infty$ for $x\leq 0$ for simplicity in the following.
Then the case I am considering is
$$I(g) = \int_{0}^{\infty} e^{-x^{1/4}}\sin(x^{1/4})e^{igx}$$
which is the Fourier transform of the continuous $L^{1}$ function $e^{-x^{1/4}}\sin(x^{1/4}) \theta(x)$, where $\theta(x)$ is Heaviside's step function. This is a unique, invertible Fourier transform. In particular, this Fourier transform cannot merely be the identically $0$ function in $g$, as that would be the Fourier transform of $0$.
However, infamously,
$$\int_{0}^{\infty} e^{-x^{1/4}}\sin(x^{1/4})x^n =0 $$
for all non-negative integers $n$, and so
$$S(g) = 0$$
$S(g)$ is identically zero and hence converges everywhere, but it nevertheless does not agree with $I(g)$ which cannot be zero for all values of $g$ by the invertibility of Fourier transforms. Thus, this gives an example where $S(g)$ converges but to the wrong answer.
