How to calculate prediction interval in GLM (Gamma) / TweedieRegression in Python? I have checked much source from webs about conducting the prediciton interval, especially in GLM function. One of the approaches is about Prediction Intervals for Machine Learning https://machinelearningmastery.com/prediction-intervals-for-machine-learning/ from Jason Brownlee. However,his method targets to the linear regression, and it might not be appropriate to the GLM (Gamma) to some degrees. Another approach I found is to use bootstrapping method to conduct the prediction interval. However, the computation was so time-consuming, and my computer's memory was killed when running the function from the article https://saattrupdan.github.io/2020-03-01-bootstrap-prediction/. I am confused how to conduct the prediction interval in an appropriate way in GLM (Gamma most probably) in Python instead in R. I have found a related package in R, but I do not want to use R to conduct the interval. Another related information I found from the web is Gamma GLM - Derive prediction intervals for new x_i: Gamma GLM - Derive prediction intervals for new x_i.
 A: Its a bit involved, but it should be doable.
As that post says, in order to get a prediction interval you have to integrate over the uncertainty in the coefficients.  That is hard to do analytically, but we can instead simulate it.  Here is some gamma regression data
N = 100
x = np.random.normal(size = N)

true_beta = np.array([0.3])
eta = 0.8 + x*true_beta
mu = np.exp(eta)
shape = 10

#parameterize gamma in terms of shaope and scale
y = gamma(a=shape, scale=mu/shape).rvs()

Now,  I will fit the gamma regression to this data

X = sm.tools.add_constant(x)

gamma_model = sm.GLM(y, X, family=sm.families.Gamma(link = sm.families.links.log()))
gamma_results = gamma_model.fit()

gamma_results.summary()

          Generalized Linear Model Regression Results           
Dep. Variable:  ,y               ,  No. Observations:  ,   100  
Model:          ,GLM             ,  Df Residuals:      ,    98  
Model Family:   ,Gamma           ,  Df Model:          ,     1  
Link Function:  ,log             ,  Scale:             ,0.075594
Method:         ,IRLS            ,  Log-Likelihood:    , -96.426
Date:           ,Mon, 30 Nov 2020,  Deviance:          ,  7.7252
Time:           ,22:45:07        ,  Pearson chi2:      ,  7.41  
No. Iterations: ,7               ,                     ,        
Covariance Type:,nonrobust       ,                     ,        
     ,   coef   , std err ,    z    ,P>|z| ,  [0.025 ,  0.975] 
const,    0.8172,    0.028,   29.264, 0.000,    0.762,    0.872
x1   ,    0.2392,    0.029,    8.333, 0.000,    0.183,    0.296



So long as I have enough data, we can make a normal approximation to the sampling distribution of the coefficients.
The mean and covariance can be obtained from the model summary.
beta_samp_mean = gamma_results.params
beta_samp_cov = gamma_results.cov_params()
dispersion = gamma_results.scale

Now, it is just a matter of sampling fake data using these estimates and taking quantiles.
X_pred = np.linspace(-2, 2)
X_pred = sm.tools.add_constant(X_pred)

num_samps = 100_000
possible_coefficients = np.random.multivariate_normal(mean = beta_samp_mean, cov = beta_samp_cov, size = num_samps)
linear_predictions = [X_pred@b for b in possible_coefficients]


y_hyp = gamma(a=1/dispersion, scale = np.exp(linear_predictions)*dispersion).rvs()

# Here is the prediction interval
l, u = np.quantile(y_hyp, q=[0.025, 0.975], axis = 0)

Its easy to then plot the prediction interval
yhat = gamma_results.predict(X_pred)
fig, ax = plt.subplots(dpi = 120)
plt.plot(X_pred[:,1], yhat, color = 'red', label = 'Estimated')
plt.plot(X_pred[:, 1], np.exp(0.8 + X_pred[:, 1]*true_beta), label = 'Truth')
plt.fill_between(X_pred[:, 1], l, u, color = 'red', alpha = 0.1, label = 'Prediction Interval')

for i in range(10):
    y_tilde = gamma(a=shape, scale=np.exp(0.8 + X_pred[:, 1]*true_beta)/shape).rvs()
    plt.scatter(X_pred[:, 1], y_tilde, s = 1, color = 'k')
plt.scatter(X_pred[:, 1], y_tilde, s = 1, color = 'k', label = 'New Data')


plt.legend()


Math of what is going on
Our data $y$ are distributed according to
$$ y\vert X \sim \mbox{Gamma}(\phi, \mu(x)/\phi) $$
At least I think that is the correct parameterization of the Gamma, I can never get it right.  In any case, assuming we use a log link for the model, this means
$$ \mu(x) = \exp(X\beta)$$
The thing is, we never know $\beta$, we only get $\hat{\beta}$ because we have to estimate the parameters of the model.  The parameters are thus a random variable (because different data can yield different parameters).  Theory says that with enough data, we can consider
$$ \hat{\beta} \sim \mbox{Normal}(\beta, \Sigma) $$
and some more theory says that plugging in our estimate for $\beta$ and $\Sigma$  should be good enough. Let $\tilde{y}\vert X$ be data I might see for observations with covariates $X$. If I could, I would really compute
$$ \tilde{y} \vert X \sim \int p(y\vert X,\beta)p (\beta) \, d \beta $$
and then take quantiles of this distribution.  But this integral is really hard, so instead we just approximate it by simulating from $p(\beta)$ (the normal distribution) and passing whatever we simulated to $p(y\vert X, \beta)$ (in this case, the gamma distribution).
Now, I realize I've been quite fast and loose here, so if any readers want to put a little more rigour into my explanation, please let me know in a comment and I will clean it up.  I think this should be good enough to give OP an idea of how this works.
