In a Gamma GLM, the statistical model for each observation 𝑖 is assumed to be $Y_i \sim Gamma(shape, scale)$, where $E(Y_i) = \mu_i = f(X_i\beta)$, and $f$ is the link function.

I've used MLE to estimate $\hat{\beta}$ and $\hat{scale}$, and wish to produce a 90% prediction interval on a new point $Y'$ given $X'$.

I can produce the confidence intervals on $E(Y|X') = \mu'$ by using link function $f$ on the normally distributed confidence intervals for $X\hat{\beta}$. Let's say $\hat{\mu'} = 10$ and 90% confidence intervals are [5, 30].

However, we want the intervals from the distribution of $Y'$, not $\mu'$. Intuitively, these intervals should be much wider than the confidence intervals for $\mu'$ I think they should also be wider than the 5th and 95th percentile of a single Gamma distribution with $\mu=\hat{\mu'}$, since the uncertainty around $\hat{\mu'}$ should translate into increased uncertainty around the final distribution, sort of like an vague prior on a bayesian posterior distribution.

What is the correct way to model prediction intervals on the new point $Y'$?

The below schema shows how uncertainty on $\mu'$ translates into many possible gamma distributions and a wide prediction interval for $Y'$

schema of Y'




  • 3
    $\begingroup$ The interval you are asking for is called a prediction-interval. A confidence interval is something else. $\endgroup$
    – whuber
    Aug 21, 2019 at 19:39
  • $\begingroup$ Thanks @whuber, fixed $\endgroup$ Aug 21, 2019 at 19:50

2 Answers 2


The prediction interval for a new observation depends on both the assumed inherent randomness in this case given by the gamma distribution, and the uncertainty coming from the parameters that are estimated and not assumed to be known.

In general there is no analytical or closed form expression for the combination of the two effects. The two main options are to ignore parameter uncertainty and to use simulation methods.

Ignoring parameter uncertainty: If we ignore that the parameters are estimated with some sampling uncertainty, the the distribution of a new observation is just given by the assumed gamma distribution. We can use the mean and scale estimates to compute the relevant prediction intervals using e.g. the distribution methods in scipy.stats. The parameterization might have to be transformed from the GLM parameterization to the scipy.stats parameterization.

Simulation Methods: One possibility is to use full bootstrap on the original sample to simulate new observations. The simpler method is to assume that the asymptotic normal distribution for the parameter estimates are appropriate and simulate the parameters of the gamma distribution from the mean and covariance of the parameter estimates. For each sampled parameter we can sample a new observation and compute a confidence interval from the simulated observations. One problem with this is that GLM only estimates the mean parameters directly, the scale is estimated using deviance or pearson chi2. That is, GLM in statsmodels in other packages does not provide a joint covariance for mean and scale parameter.

Because of these problems, statsmodels currently provides prediction intervals for new observations that take parameter uncertainty into account only for the linear normal case, i.e. OLS.

  • 1
    $\begingroup$ Thanks @Josef. How about sampling the underlying data with replacement, fitting the full Gamma GLM on each new sample and predicting y_i, and then take the distribution of those y_i s ? That way, we get some natural variation in the scale parameter by differing deviances $\endgroup$ Aug 21, 2019 at 20:29
  • 1
    $\begingroup$ Note that predicting y_i usually means to predict the mean or expected value and not an observation. You would need a draw from the predictive distribution for a new observation. Something like that is an option that includes all estimation uncertainty and what I meant with full bootstrap. The disadvantage is that it is costly to compute because the model has to be estimated for each draw of a bootstrap sample. Using the original parameter estimates as starting values will speed up the computations but there is still a large amount of overhead to compute. $\endgroup$
    – Josef
    Aug 22, 2019 at 1:21

This publication answers your question exactly

"A derivation of prediction intervals for gamma regression"



The paper is not publicly available, but if it helps, here's an R code for implementation. Say, you want to model the outcome y which is gamma distributed on two predictors x1 and x2. A prediction and prediction interval is required a new point x0=(x01,x02)=(0.75,-1).

Generate Fake Data

x1   = runif(100,-1,1); x2 = rnorm(100)
Mean = exp(.1+x1+x2/10) 
y    = rgamma(100,5,5/Mean)

Fit a gamma model and extract the model matrix

mod   = glm(y~x1+x2,family=Gamma(link="log"))
alpha = (mod$deviance/mod$df.residual)^-1  # 
X     = model.matrix(mod)
n     = nrow(X); p = ncol(X) 

Assume prediction is required at x1 =0.75 and x2 = -1

x0    = c(1,0.75,-1)
pred  = exp(x0%*%coef(mod)) # Prediction of y at x0 

Find quantiles of G_hat as proposed in eq(3) in the paper

h     = as.vector(t(x0)%*%solve(t(X)%*%X)%*%x0)
a1    = rchisq(10^6,2*alpha) # 10^6 is arbitrary large.
a2    = rt(10^6,n-p)
G     = a1/(a2*sqrt(h/alpha) + 1)
G_25  = quantile(G,0.025); G_975 = quantile(G,0.975)

A 95% prediction interval at x0


[1] 0.5745241 4.8665244

  • $\begingroup$ Paywall, unfortunately $\endgroup$ Sep 10, 2020 at 20:56

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.