Worked Prediction Intervals at Different Leverage Values

Why Compare Prediction Intervals Across Design Points?

A prediction interval answers a question about one future outcome, not just the average outcome at a specified design point. Even when the same model and the same sample are used, prediction intervals can differ across design points: the fitted value may change, and the uncertainty in estimating the mean may change with leverage. The future outcome also has its own random error, so its prediction interval is wider than an interval for the mean response at the same point.

This tutorial works through those differences under the fixed-design normal linear model. The examples compare a point near the center of the observed design with an extrapolation point, examine a model with unequal group sizes, and show the corresponding calculation when the error variance is known.

Definition: Suppose $Y=X\beta+\varepsilon$, where $X$ is a fixed $n\times p$ matrix of rank $p$, and $\varepsilon\sim N_n(0,\sigma^2 I_n)$. At a design vector $x_0\in\mathbb{R}^p$, the fitted mean is $x_0^\top\hat{\beta}$. The leverage of this point is $h_0=x_0^\top(X^\top X)^{-1}x_0$. For a future observation $Y_0=x_0^\top\beta+\varepsilon_0$, assume $\varepsilon_0\sim N(0,\sigma^2)$, independently of the sample errors.

The leverage formula applies to a new design vector, whether or not that vector is one of the rows of $X$. For an observed row, it agrees with the corresponding diagonal entry of the projection matrix (also called the hat matrix). Since $X^\top X$ is positive definite, $h_0\geq 0$. Larger leverage means greater sampling uncertainty in the fitted mean at that point.

The Exact Prediction Interval

Write $\hat{\beta}=(X^\top X)^{-1}X^\top Y$, let $SSE=(Y-X\hat{\beta})^\top(Y-X\hat{\beta})$, and estimate the error variance by $s^2=SSE/(n-p)$. The following result gives the exact interval when $\sigma^2$ is unknown.

Theorem: Under the fixed-design normal model above, with $n>p$, a prediction interval with coverage $1-\alpha$ for one future observation at $x_0$ is $$ x_0^\top\hat{\beta}\ \pm\ t_{n-p,\,1-\alpha/2}\,s\sqrt{1+h_0}. $$ Here $t_{\nu,\,1-\alpha/2}$ is the $1-\alpha/2$ quantile of the Student $t$ distribution with $\nu$ degrees of freedom. Conditional on $X$, the interval covers $Y_0$ with probability $1-\alpha$.

Proof. The prediction error is $$ Y_0-x_0^\top\hat{\beta} =\varepsilon_0-x_0^\top(\hat{\beta}-\beta). $$ Because $\hat{\beta}-\beta=(X^\top X)^{-1}X^\top\varepsilon$, its second term is a linear combination of normal sample errors with variance $\sigma^2 h_0$. It is independent of $\varepsilon_0$, which has variance $\sigma^2$. Thus the prediction error is normal with mean zero and variance $\sigma^2(1+h_0)$.

To establish the required independence from the variance estimate, let $P=X(X^\top X)^{-1}X^\top$ and $M=I_n-P$. The fitted error is determined by $P\varepsilon$, while the residual vector is $M\varepsilon$. Since $PM=0$, these two vectors have zero covariance; because they are jointly normal, they are independent. Also, $\varepsilon_0$ is independent of the sample errors. The prediction error is therefore independent of the residual sum of squares. The matrix $M$ is an idempotent projection of rank $n-p$, so $SSE/\sigma^2=\varepsilon^\top M\varepsilon/\sigma^2$ has a chi-square distribution with $n-p$ degrees of freedom. Consequently, $$ \frac{Y_0-x_0^\top\hat{\beta}}{s\sqrt{1+h_0}} \sim t_{n-p}. $$ Using the two-sided $t$ quantiles gives the stated interval and its conditional coverage. $\square$

The factor $1+h_0$ has two parts. The $h_0$ part is the variance, in units of $\sigma^2$, from estimating the mean at $x_0$; the $1$ is the variance of the new observation's own error. This decomposition explains why the prediction interval is not just an interval for the fitted mean.

Worked Comparisons

Worked Example: A Central Point and an Extrapolation Point

Consider a simple linear regression with an intercept and slope, using $x=-2,-1,0,1,2,3$. The observed responses are $7,6,11,12,14,16$, respectively. Find 95% prediction intervals at $x_0=0.5$, near the center of the design, and at $x_0=4$, beyond its observed range.

The fitted line is $\hat{y}=10+2x$. To verify the fit, subtract $10+2x$ from each response: the residuals are $1,-2,1,0,0,0$. They sum to zero, and their inner product with the $x$-values is $(-2)(1)+(-1)(-2)+(0)(1)=0$. Thus they are orthogonal to both columns of the design matrix, as required for the least-squares fit. Their squared sum is $SSE=1+4+1=6$, so $s^2=6/(6-2)=1.5$.

The design mean is $\bar{x}=0.5$, and $$ S_{xx}=\sum_{i=1}^6(x_i-\bar{x})^2 =\sum_{i=1}^6x_i^2-6\bar{x}^2 =19-1.5=17.5. $$ For an intercept-and-slope model, $h_0=1/n+(x_0-\bar{x})^2/S_{xx}$. At $x_0=0.5$, $h_0=1/6$, and the fitted value is $10+2(0.5)=11$. The prediction variance estimate is $s^2(1+h_0)=1.5(7/6)=1.75$.

With $n-p=4$ degrees of freedom, $t_{4,0.975}\approx2.776$. The interval's half-width is $2.776\sqrt{1.75}\approx3.673$, giving $[11-3.673,\ 11+3.673]$, or approximately $[7.327,\ 14.673]$.

At $x_0=4$, the fitted value is $18$, while $$ h_0=\frac{1}{6}+\frac{(4-0.5)^2}{17.5} =\frac{1}{6}+0.7=\frac{13}{15}. $$ The estimated prediction variance is $1.5(1+13/15)=2.8$. Its half-width is $2.776\sqrt{2.8}\approx4.646$, giving the interval $[13.354,\ 22.646]$, approximately. The extrapolation interval is wider because its leverage is higher; the interval is also centered at a different fitted value.

Worked Example: Unequal Group Sizes and Group-Specific Leverage

Suppose a model has an intercept and a binary treatment indicator. Three control outcomes are $8,10,12$, and two treatment outcomes are $13,15$. Calculate a 95% prediction interval for one future outcome in each group.

The fitted value for each group is its sample mean: $10$ for control and $14$ for treatment. The control residuals are $-2,0,2$, with squared sum $8$; the treatment residuals are $-1,1$, with squared sum $2$. Thus $SSE=10$, $n=5$, $p=2$, and $s^2=10/3$.

In this model, the fitted value for each group is its group mean. The variance of a group mean is $\sigma^2/n_g$, so the leverage for a point in a group of size $n_g$ is $1/n_g$. It is $1/3$ for control and $1/2$ for treatment. Using $t_{3,0.975}\approx3.182$, the control prediction variance estimate is $(10/3)(1+1/3)=40/9$. Its standard error is $\sqrt{40/9}\approx2.108$, so its half-width is $3.182(2.108)\approx6.708$. The interval is approximately $[10-6.708,\ 10+6.708]=[3.292,\ 16.708]$.

For treatment, the prediction variance estimate is $(10/3)(1+1/2)=5$. Its standard error is $\sqrt{5}\approx2.236$, and its half-width is $3.182(2.236)\approx7.115$. The interval is approximately $[14-7.115,\ 14+7.115]=[6.885,\ 21.115]$. The treatment interval is wider because the treatment group is smaller and therefore has greater leverage. Its center is also higher because the fitted treatment mean is higher.

Worked Example: A Known Error Variance

A validated process provides a known error variance $\sigma^2=4$. A simple linear regression with an intercept uses $x=-2,-1,0,1,2$ and responses $4,4.5,5,5.5,6$, giving the fitted line $\hat{y}=5+0.5x$. Find a 95% prediction interval at $x_0=1$.

Here $\bar{x}=0$ and $S_{xx}=4+1+0+1+4=10$. The leverage is $h_0=1/5+1^2/10=0.3$, and the fitted value is $5.5$. Since $\sigma$ is known, the standardized prediction error is normal rather than Student $t$, with standard deviation $\sqrt{\sigma^2(1+h_0)}=\sqrt{4(1.3)}=\sqrt{5.2}$.

The 95% normal critical value is $z_{0.975}\approx1.96$. The half-width is $1.96\sqrt{5.2}\approx4.469$, so the exact interval is approximately $[5.5-4.469,\ 5.5+4.469]=[1.031,\ 9.969]$. The residuals happen to be zero here, but they do not determine the interval: the calculation uses the stipulated known error variance.

Why Leverage Changes the Width

For a fixed model, sample, and confidence level, the estimated-variance interval has width $$ W(h_0)=2t_{n-p,\,1-\alpha/2}s\sqrt{1+h_0}. $$ All factors other than $h_0$ are fixed and nonnegative, and $\sqrt{1+h_0}$ increases as $h_0$ increases. Therefore, if two design points have leverages $h_A<h_B$, their prediction interval widths satisfy $$ W(h_A)<W(h_B). $$ This comparison concerns width, not the interval's location: different fitted values move the centers. The inequality also assumes the same sample, model, and coverage level. Across different fitted models, both $s$ and the critical value may change, so leverage alone need not determine which interval is wider.

A common mistake is to use $s\sqrt{h_0}$ as the prediction standard error because it describes uncertainty in the estimated mean, then call the result a prediction interval. That omits the independent error in the future observation. For one future outcome the correct estimated standard error is $s\sqrt{1+h_0}$; for an interval about the mean response, the additional $1$ is not present. A second mistake is to interpret a wider interval as evidence that the fitted mean is larger. Width reflects uncertainty, while the center is the fitted value.

Key takeaway: For one future observation, leverage contributes to uncertainty in the fitted mean, and the future outcome contributes an additional unit of error variance. With an estimated error variance, combine both through $s\sqrt{1+h_0}$ and use the appropriate $t$ critical value.

Check Your Understanding

Work through these questions using the definitions and calculations above.

In the fixed-design normal model, what assumptions make the prediction-error statistic independent of the residual sum of squares?
Why does the prediction variance contain $1+h_0$, rather than just $h_0$?
For two design points in the same fitted model, what can be concluded about their prediction-interval widths if one has strictly higher leverage?
In a model with an intercept and a binary group indicator, which group has higher leverage when its sample size is smaller? Why?
When the error variance is known, which critical distribution is used for the exact prediction interval, and why?

Worked Prediction Intervals at Different Leverage Values

What You'll Learn