Prediction Intervals in Multiple Regression

Why Multiple Regression Needs a Prediction Interval

A fitted multiple-regression model gives a point prediction for a response at a specified combination of predictor values. But a new individual outcome will usually differ from that fitted value for two reasons: the regression coefficients have been estimated from a sample, and the new outcome has its own random error. A prediction interval must account for both sources of uncertainty.

The basic idea is the same as for a prediction interval in simple regression, but the predictor values are now collected into a design vector. That vector allows the uncertainty in all the estimated coefficients—including their covariances—to contribute to the interval.

Model and Notation

Suppose the fitted model has \(p\) coefficients, including an intercept, and is based on \(n\) observations. For a new case, write its design vector as \(\mathbf{x}_0\), with the intercept entry first. For example, with two quantitative predictors, \(\mathbf{x}_0=(1,x_{01},x_{02})^\mathsf{T}\). For a categorical predictor, its entry or entries follow the coding used to fit the model.

Let \(\hat{\boldsymbol{\beta}}\) be the vector of estimated coefficients and let \(s\) be the residual standard error. The point prediction for the new case is \(\hat{y}_0=\mathbf{x}_0^\mathsf{T}\hat{\boldsymbol{\beta}}\). Let \(\widehat{\operatorname{Cov}}(\hat{\boldsymbol{\beta}})\) denote the estimated covariance matrix of the coefficient estimates. The conventional prediction-interval standard error is

The first term, \(s^2\), represents the residual variation of an individual response around its conditional mean. The quadratic term represents uncertainty in the estimated conditional mean at \(\mathbf{x}_0\). The off-diagonal entries of the coefficient covariance matrix matter: they capture how coefficient estimates vary together.

In ordinary least squares, when the design matrix is \(X\), the usual estimated coefficient covariance matrix is \(s^2(X^\mathsf{T}X)^{-1}\). Thus the same standard error can be written

These are two versions of the same calculation. If the coefficient covariance matrix is supplied directly, use it as the estimated covariance matrix in the first formula; do not multiply its quadratic form by \(s^2\) again.

Why the Two Sources of Variance Are Added

The variance formula follows by comparing the future response with its point prediction. Under the linear model, write the new response as \(Y_0=\mathbf{x}_0^\mathsf{T}\boldsymbol{\beta}+\varepsilon_0\), where \(\boldsymbol{\beta}\) is the true coefficient vector and \(\varepsilon_0\) is the new case’s error. The prediction error is

For this calculation, treat the new predictor values as fixed, or condition on them. Assume the new error is independent of the errors used to fit the model, has mean zero, and has variance \(\sigma^2\). The fitted coefficient vector is unbiased under the correctly specified linear mean model, and its covariance is \(\operatorname{Cov}(\hat{\boldsymbol{\beta}})\). Because the two terms in the prediction error are independent and have mean zero, the variance of their difference is the sum of their variances:

Replacing the unknown variance quantities with their estimates gives the prediction standard error above. For the usual \(t\)-based interval, use the residual degrees of freedom, \(df=n-p\), and the appropriate \(t\) critical value. The conventional procedure relies on the linear model being suitable for the data, independent errors with constant variance, and the usual normal-error condition for \(t\)-based calibration. Its detailed repeated-sampling coverage properties are a separate issue.

Constructing the Interval

For a two-sided \(100(1-\alpha)\%\) prediction interval, use

Here \(t_{1-\alpha/2,\,df}\) is the positive \(t\) critical value with cumulative probability \(1-\alpha/2\) and \(df\) degrees of freedom. The point prediction and both interval endpoints are expressed in the response’s units. The interpretation concerns one new response at the specified predictor values, not the average response of all cases with those values.

Worked Example With a Categorical Predictor

Categorical predictors use the same formula; the design vector must simply reflect the model’s coding. Suppose a hypothetical model predicts exam score from study hours and a tutoring indicator, coded 1 for tutoring and 0 otherwise. The fitted coefficients are \(50\) points for the intercept, \(2.4\) points per study hour, and \(3.0\) points for tutoring. Consider a new student who studied 5 hours and received tutoring, so \(\mathbf{x}_0=(1,5,1)^\mathsf{T}\). Let \(s=4\) points and \(df=36\). Suppose

First calculate the point prediction:

The coefficient-variance contribution is the sum of the diagonal terms \(1+25(0.04)+0.64=2.64\) and the cross terms \(2(1)(5)(-0.10)+2(1)(1)(-0.20)+2(5)(1)(0.01)=-1.30\). Thus the quadratic form is \(2.64-1.30=1.34\), and

For a 90% interval, \(t_{0.95,36}\approx1.688\). The margin is approximately \(1.688(4.164)=7.029\) points, so the interval is

This interval concerns one future student with the specified study-hours and tutoring profile. The indicator value is entered as 1 because the student received tutoring; a student without tutoring would have a different design vector and prediction.

Worked Example: Using Residual Standard Error and \((X^\mathsf{T}X)^{-1}\)

A third hypothetical model has an intercept and two quantitative predictors, with \(n=30\) observations and \(p=3\) fitted coefficients. Thus \(df=n-p=27\). At a new profile \(\mathbf{x}_0=(1,2,1)^\mathsf{T}\), the estimated coefficients are \((8,1.5,-0.5)^\mathsf{T}\), the residual standard error is \(s=3\) response units, and

The point prediction is \(8+(2)(1.5)+(1)(-0.5)=10.5\) response units. For the quadratic form, multiplying the matrix by \(\mathbf{x}_0\) gives \((0.10,0.02,0.03)^\mathsf{T}\). Taking the dot product with \(\mathbf{x}_0\) gives \(0.10+2(0.02)+0.03=0.17\). Therefore,

For a 90% interval, \(t_{0.95,27}\approx1.7033\). The margin is approximately \(1.7033(3.245)=5.527\) response units. The resulting interval is

This calculation uses \(s\sqrt{1+\mathbf{x}_0^\mathsf{T}(X^\mathsf{T}X)^{-1}\mathbf{x}_0}\). Equivalently, the estimated coefficient covariance matrix is \(s^2(X^\mathsf{T}X)^{-1}\), so the direct covariance-matrix formula gives the same standard error.

Common Mistakes and a Practical Check

• Leaving out the new response’s residual variation. Using only \(\sqrt{\mathbf{x}_0^\mathsf{T}\widehat{\operatorname{Cov}}(\hat{\boldsymbol{\beta}})\mathbf{x}_0}\) measures uncertainty in the fitted mean, not prediction uncertainty for an individual outcome. Include \(s^2\) under the square root.
• Counting the intercept incorrectly. The design vector includes an intercept entry, and \(p\) counts the intercept among the fitted coefficients. Hence the usual residual degrees of freedom are \(n-p\).
• Scaling a supplied covariance matrix twice. If software or a problem gives \(\widehat{\operatorname{Cov}}(\hat{\boldsymbol{\beta}})\), use its quadratic form directly. Multiplying that term by \(s^2\) again is incorrect.
• Ignoring covariance terms. The coefficient covariance matrix is not generally diagonal. Omitting its off-diagonal entries can change the standard error.
• Describing the interval as covering a mean. The interval is for one new response at the stated profile. The conditional mean is a different target and does not include the new response’s residual variation.

Before reporting an interval, check that the design vector matches the model’s predictor coding, the covariance input is on the correct scale, the \(s^2\) term is present, and the \(t\) critical value uses the residual degrees of freedom. Then state the response units and specify the new case’s predictor profile.