Prediction Intervals For Gamma Responses

What a Prediction Interval Describes

A fitted gamma model gives a conditional mean for a response at specified covariate values. But a future response is not generally equal to that mean: even if the mean were known exactly, individual outcomes would still vary. A prediction interval describes a range for one future response. It is wider in purpose than an interval for the conditional mean, which describes uncertainty about the mean itself.

For gamma response models, the dispersion determines the amount of response variation relative to the mean. Once a model supplies a fitted conditional mean and a dispersion value, gamma quantiles give a direct way to construct a prediction interval. That calculation is conditional on the fitted parameter values; when those values are estimated, the resulting interval does not automatically include all parameter-estimation uncertainty.

Definition: Suppose \(Y\), conditional on covariates \(x\), has a gamma distribution with mean \(\mu(x)>0\) and dispersion \(\phi>0\), so that \(\operatorname{Var}(Y\mid x)=\phi\mu(x)^2\). Its shape is \(1/\phi\) and its scale is \(\phi\mu(x)\). A prediction interval for a future response \(Y_{\mathrm{new}}\) is an interval intended to contain that response at a specified coverage probability.

Here, the gamma scale is the multiplier of a unit-scale gamma variable. In particular, if \(G\) has a gamma distribution with shape \(a\) and scale 1, then \(\phi\mu G\) has shape \(1/\phi\), scale \(\phi\mu\), and mean \(\mu\). Write \(q_p(a)\) for the \(p\) quantile of that unit-scale gamma distribution:

P\{G\le q_p(a)\}=p,\qquad G\sim\operatorname{Gamma}(a,\text{ scale }1).

This notation separates the shape-dependent quantile calculation from the fitted mean. It also makes the role of dispersion visible: changing \(\phi\) changes both the gamma shape and the scale.

Constructing a Plug-In Prediction Interval

Suppose the fitted model gives \(\hat{\mu}(x_0)\) at the covariate values \(x_0\), and let \(\hat{\phi}\) be the dispersion estimate used for prediction. For a central interval with nominal coverage \(1-\alpha\), the plug-in endpoints are:

\left[ \hat{\phi}\hat{\mu}(x_0)\,q_{\alpha/2}(1/\hat{\phi}), \quad \hat{\phi}\hat{\mu}(x_0)\,q_{1-\alpha/2}(1/\hat{\phi}) \right].

The interval is central because it leaves probability \(\alpha/2\) below the lower endpoint and \(\alpha/2\) above the upper endpoint under the fitted gamma distribution. It need not be symmetric around \(\hat{\mu}(x_0)\): gamma distributions are generally right-skewed, especially when their shape is small. The interval’s endpoints are positive, as required for a gamma response.

The formula is a plug-in calculation: it acts as though the fitted mean and dispersion were the true parameter values. If those values were known and the gamma model were correct, the stated coverage would be exact. When they are estimated from data, the interval is useful but its repeated-sampling coverage is not generally exactly \(1-\alpha\).

Worked Example: An Exponential Special Case

At a specified covariate value, suppose \(\hat{\mu}=10\) and \(\hat{\phi}=1\). Find the central 90% plug-in prediction interval.

When \(\phi=1\), the shape is \(1/\phi=1\) and the scale is \(\phi\hat{\mu}=10\). A gamma distribution with shape 1 is exponential. Its unit-scale quantile is \(q_p(1)=-\log(1-p)\), since its distribution function is \(1-e^{-z}\) for \(z\ge 0\).

The lower unit-scale quantile is \(q_{0.05}(1)=-\log(0.95)\approx0.05129\). The upper one is \(q_{0.95}(1)=-\log(0.05)\approx2.99573\). Multiplying each by the fitted scale 10 gives:

[10(0.05129),\,10(2.99573)]\approx[0.513,\,29.957].

Under the fitted exponential distribution, 5% of future responses lie below the lower endpoint and 5% lie above the upper endpoint. The interval is not centered at 10, because the exponential distribution is strongly right-skewed.

Worked Example: Dispersion Changes the Interval

Suppose \(\hat{\mu}=12\) and \(\hat{\phi}=0.5\). Find the central 90% plug-in prediction interval using the unit-scale gamma quantiles for shape 2.

The fitted shape is \(1/0.5=2\), and the fitted scale is \(0.5(12)=6\). For shape 2, the unit-scale gamma distribution has cumulative distribution function \(F(z)=1-e^{-z}(1+z)\) for \(z\ge0\). Solving \(F(z)=0.05\) and \(F(z)=0.95\) gives \(q_{0.05}(2)\approx0.35536\) and \(q_{0.95}(2)\approx4.7439\).

Rescale these quantiles by 6:

[6(0.35536),\,6(4.7439)]\approx[2.132,\,28.463].

The fitted mean is still 12, but the dispersion differs from the previous example, so the shape and the interval differ too. The relative standard deviation here is \(\sqrt{\phi}=\sqrt{0.5}\), rather than 1.

Why the Quantile Formula Works

The gamma parameterization makes the formula a rescaling of a unit-scale distribution. If \(G\) has shape \(1/\phi\) and scale 1, then \(Y=\phi\mu G\) has mean \(\mu\) and dispersion \(\phi\). Since the multiplier is positive, it preserves the order of values and therefore rescales each quantile by \(\phi\mu\).

Result: If \(\mu>0\) and \(\phi>0\) are known and \(Y_{\mathrm{new}}\mid x_0\) follows the gamma model with mean \(\mu\) and dispersion \(\phi\), then its central \(1-\alpha\) interval is \([\phi\mu q_{\alpha/2}(1/\phi),\,\phi\mu q_{1-\alpha/2}(1/\phi)]\), for \(0<\alpha<1\). Its coverage is exactly \(1-\alpha\).

Standardize.
Write \(Y_{\mathrm{new}}=\phi\mu G\), where \(G\) has shape \(1/\phi\) and scale 1.

Transform the quantiles.
Because \(\phi\mu>0\), the event that \(Y_{\mathrm{new}}\) lies between the proposed endpoints is exactly the event that \(G\) lies between \(q_{\alpha/2}(1/\phi)\) and \(q_{1-\alpha/2}(1/\phi)\).

Calculate coverage.
The gamma distribution is continuous, so its quantiles satisfy \(P\{G\le q_p(1/\phi)\}=p\). Subtracting the lower-tail probability from the upper-tail probability gives \(1-\alpha/2-\alpha/2=1-\alpha\). Thus the interval has the stated coverage. \(\square\)

This result is exact only when the model parameters are known. Substituting \(\hat{\mu}(x_0)\) and \(\hat{\phi}\) gives the practical plug-in interval, not a proof of exact coverage after fitting. That distinction matters most when the training sample is small or the prediction point is far from the covariate patterns represented in the data.

Future-Response Variation and Mean Uncertainty

Even if the true conditional mean \(\mu(x_0)\) were known, a new response would have variance \(\phi\mu(x_0)^2\). Estimating the model adds another source of uncertainty: the fitted mean varies from sample to sample. An interval for the mean targets the latter uncertainty; a prediction interval must address the variation of a future response as well.

For example, under a log link, \(\log\mu(x)=x^\mathsf{T}\beta\). If the estimated coefficient vector has estimated covariance matrix \(\widehat{\operatorname{Cov}}(\hat{\beta})\), then the estimated variance of the fitted linear predictor at \(x_0\) is \(s_\eta^2=x_0^\mathsf{T}\widehat{\operatorname{Cov}}(\hat{\beta})x_0\). The delta method gives the approximate standard error of the fitted mean as \(\hat{\mu}(x_0)s_\eta\), since the derivative of the exponential function at \(\log\hat{\mu}(x_0)\) is \(\hat{\mu}(x_0)\). Such a standard error quantifies uncertainty in the fitted mean; it does not replace the gamma response distribution in the prediction interval.

Result: If the training data and the future response are independent conditional on fixed model parameters, then the prediction error \(Y_{\mathrm{new}}-\hat{\mu}(x_0)\) has variance \(\operatorname{Var}(Y_{\mathrm{new}}\mid x_0)+\operatorname{Var}(\hat{\mu}(x_0))\). This separates future-response variation from fitted-mean variation.

To verify the variance statement, let \(Y_{\mathrm{new}}\) and \(\hat{\mu}(x_0)\) be independent conditional on the fixed parameters. For any independent random variables \(A\) and \(B\), expanding the variance gives \(\operatorname{Var}(A-B)=\operatorname{Var}(A)+\operatorname{Var}(B)-2\operatorname{Cov}(A,B)\). Independence makes the covariance zero. Applying this identity yields the result. In practice, replacing the unknown response variance by \(\hat{\phi}\hat{\mu}(x_0)^2\) and adding an estimated variance for the fitted mean gives a variance approximation, not a gamma distribution for the prediction error. \(\square\)

Worked Example: Comparing Mean and Response Uncertainty

Suppose a log-link model gives \(\hat{\mu}(x_0)=10\), \(\hat{\phi}=0.25\), and a standard error \(s_\eta=0.10\) for the fitted log mean. Compare a 95% approximate interval for the mean with the scale of uncertainty in a future response.

For the fitted mean, the log-scale interval is \(\log(10)\pm1.96(0.10)\). Exponentiating its endpoints gives \(10e^{-0.196}\approx8.22\) and \(10e^{0.196}\approx12.16\). This is an approximate interval for the conditional mean, not for one future observation.

The delta-method standard error of the fitted mean is approximately \(10(0.10)=1\). The fitted gamma response variance is \(\hat{\phi}\hat{\mu}^2=0.25(10^2)=25\), so its standard deviation is 5. Under the independence conditions above, the approximate variance of \(Y_{\mathrm{new}}-\hat{\mu}(x_0)\) is \(25+1=26\), with standard deviation \(\sqrt{26}\approx5.10\). The difference in scale shows why an interval for the mean cannot serve as a prediction interval. The variance calculation alone does not specify a prediction interval; the gamma shape and quantiles are needed for the plug-in interval.

Practical Construction and a Common Pitfall

A practical plug-in calculation follows a short sequence. First, specify the covariate values for the future observation and compute the fitted mean using the model’s link and coefficients. Check that this mean is positive. Next, select the dispersion estimate and calculate the shape \(1/\hat{\phi}\) and scale \(\hat{\phi}\hat{\mu}(x_0)\). Finally, obtain the two gamma quantiles for the desired tail probabilities and multiply them by the scale. If uncertainty in the fitted parameters is important, it requires an additional method rather than a change of interpretation of the plug-in interval.

Fit and evaluate the mean.
At the prediction covariates \(x_0\), calculate \(\hat{\mu}(x_0)\) and verify that it is positive.

Set the gamma parameters.
Use shape \(1/\hat{\phi}\) and scale \(\hat{\phi}\hat{\mu}(x_0)\).

Choose and calculate endpoints.
For central coverage \(1-\alpha\), use the gamma quantiles at \(\alpha/2\) and \(1-\alpha/2\), each multiplied by the fitted scale.

State what is included.
Describe the result as a plug-in prediction interval, and assess separately whether uncertainty in the estimated mean or dispersion needs to be incorporated.

A common mistake is to compute an interval for the mean and call it a prediction interval. For example, a log-scale interval based on the standard error of \(x_0^\mathsf{T}\hat{\beta}\) measures uncertainty in the fitted mean. It does not account for the gamma variation of a new response. The correction is to use the gamma quantiles for future-response variation, while treating uncertainty in estimated parameters as a distinct issue.

When parameter uncertainty is material, one option is a parametric simulation procedure: repeatedly generate model parameter estimates or coefficient vectors according to an uncertainty approximation, calculate the resulting conditional means and dispersions, and then generate a gamma response for each draw. Quantiles of the simulated future responses can reflect both sources of variation. The procedure must match the intended inferential target, and its coverage still depends on the model and uncertainty approximation being adequate.

Key takeaway: The gamma quantile formula describes variation in a future response conditional on the fitted mean and dispersion. An interval for the conditional mean describes a different uncertainty. A plug-in prediction interval should not be presented as though it automatically includes parameter-estimation uncertainty.

Check Your Understanding

Use the gamma mean-dispersion parameterization to answer the following questions.

For a gamma response with mean \(\mu\) and dispersion \(\phi\), what are its shape, scale, and coefficient of variation?
Why are the unit-scale gamma quantiles multiplied by \(\phi\mu\) to obtain response-scale endpoints?
In the exponential special case with mean 10, why is a central prediction interval not symmetric around 10?
What uncertainty does an interval based on the standard error of a fitted log mean describe, and what important source of variation does it omit for a future response?
Under what conditions does the stated central gamma prediction interval have exact coverage, and why does substituting estimated parameters change that claim?