Plotting confidence intervals from Python statsmodels

Question

I'm trying to plot a statsmodels GLM regression fit together with the original data and the uncertainty of the fit.

The regression is a Gamma regression with log-link (for nonnegative data).

I'm using the following code with statsmodels 0.13.2:

import statsmodels.api as sm

rg = sm.GLM.from_formula("y ~ x", family=sm.families.Gamma(link=sm.families.links.log()), data=df).fit(scale="X2")
fig  = sm.graphics.plot_fit(rg, "x")

This shows the observed and fitted values, but I don't get vertical lines that show the uncertainty. Instead I get an error:

AttributeError                            Traceback (most recent call last)
<ipython-input-69-923641bb2221> in <module>
----> 1 fig  = sm.graphics.plot_fit(rg, "x")

~/.Venvs/ml/lib/python3.8/site-packages/statsmodels/graphics/regressionplots.py in plot_fit(results, exog_idx, y_true, ax, vlines, **kwargs)
    156     ax.plot(x1, results.fittedvalues[x1_argsort], 'D', color='r',
    157             label='fitted', **kwargs)
--> 158     if vlines is True:
    159         _, iv_l, iv_u = wls_prediction_std(results)
    160         ax.vlines(x1, iv_l[x1_argsort], iv_u[x1_argsort], linewidth=1,

~/.Venvs/ml/lib/python3.8/site-packages/statsmodels/sandbox/regression/predstd.py in wls_prediction_std(res, exog, weights, alpha)
     94     #predvar = res3.mse_resid + np.diag(np.dot(X2,np.dot(covb,X2.T)))
     95     # predication variance only
---> 96     predvar = res.mse_resid/weights + (exog * np.dot(covb, exog.T).T).sum(1)
     97     predstd = np.sqrt(predvar)
     98     tppf = stats.t.isf(alpha/2., res.df_resid)

AttributeError: 'GLMResults' object has no attribute 'mse_resid'

Calling sm.graphics.plot_fit(rg, "x", vlines=False) makes the error go away, so it is clearly associated with the uncertainty plot.

• Is there a way to make plot_fit work with this kind of GLM model?

• Is there another way to plot the uncertainty that allows for more flexibility in the visualization (e.g. lines or a "uncertainty band" grey area instead of the vertical lines from plot_fit?