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.


1 Answer 1


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()


          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')


enter image description here

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.

  • $\begingroup$ 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 $\endgroup$
    – Anaconda
    Dec 1, 2020 at 9:09
  • $\begingroup$ y_hyp = gamma(a=1/dispersion, scale = np.exp(linear_predictions)*dispersion).rvs() This code I'm confused and get lost. Thanks. $\endgroup$
    – Anaconda
    Dec 1, 2020 at 9:09
  • $\begingroup$ 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. $\endgroup$
    – Anaconda
    Dec 1, 2020 at 9:26
  • $\begingroup$ 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. $\endgroup$
    – Anaconda
    Dec 1, 2020 at 10:18
  • 2
    $\begingroup$ 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. $\endgroup$ Dec 1, 2020 at 13:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.