# 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.

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



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)

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.

• Thank you for the reply, but I didn't get the part of num_samps = 100_000 and possible_coefficients = np.random.multivariate_normal(mean = beta_samp_mean, cov = beta_samp_cov, size = num_samps) Are these codes corresponding to the hat of Beta ~ Normal(Beta, sigma) by using the simulation method like resampling or bootstrapping? linear_predictions = [X_pred@b for b in possible_coefficients] Does this code mean the sampling 100,000 X_pred? tbc Dec 1 '20 at 9:09
• y_hyp = gamma(a=1/dispersion, scale = np.exp(linear_predictions)*dispersion).rvs() This code I'm confused and get lost. Thanks. Dec 1 '20 at 9:09
• If I would like to apply the way in my real dataset, I think I need to make some changes in y_hyp, but I didn't get this part(Is the y_hyp called the 'simulated distribution' that I need to do the quantiles?) , which is different from yhat. Dec 1 '20 at 9:26
• Pardon my replies too many times. I tried your methods in my real dataset, and it seems working well. However, I would like to make sure that I really undrerstand the rationale behind the problem. Based on my comprehension so far, the y_hyp, is the way to simulate the integral,right? So that's why when I modify you codes: serveral lines of codes where I pasted to the above comments, the method works? Appreciate. Dec 1 '20 at 10:18
• Ok, there is a lot to unpack here. I'm not going to explain the code, it should be fairly self explanatory to someone familiar with python. In essence, we draw from a gamma distribution which is parameterized to have np.exp(linear_predictor) as its mean. Since linear_predictor is a random quantity (it depends on the possible_coefficients, the sampling dist of the coefficients), y_hyp is an approximation to the integral. Dec 1 '20 at 13:59