Exact Prediction Ellipsoids for Multiple Future Responses

Why Predict Several Responses Jointly?

A collection of future outcomes is often more useful as a unit than as separate predictions. For example, a researcher may want to predict responses at several design points and assess whether the entire observed vector is plausible. Separate prediction intervals do not, by themselves, provide a stated joint coverage probability. A joint prediction ellipsoid addresses that question by accounting for the covariance among prediction errors and calibrating the region for the vector as a whole.

The setting here is the normal linear model with a fixed design matrix. The key distinction from a region for mean responses is that future observations contain both uncertainty from estimating the mean and new response noise. The ellipsoid must include both sources.

Model and Prediction-Error Covariance

Definition: Let $Y\in\mathbb{R}^n$ satisfy $Y=X\beta+\varepsilon$, where $X$ is a fixed $n\times p$ matrix of full column rank, $\beta\in\mathbb{R}^p$, and $\varepsilon\sim N_n(0,\sigma^2 I_n)$, with $\sigma^2>0$. Let $Y_0\in\mathbb{R}^q$ be the vector of future responses at specified design points, so $Y_0=X_0\beta+\varepsilon_0$, where $X_0$ is a fixed $q\times p$ matrix. Assume $\varepsilon_0\sim N_q(0,\sigma^2 I_q)$ and that $\varepsilon_0$ is independent of $\varepsilon$. Write $\hat\beta=(X'X)^{-1}X'Y$, $s^2=\|Y-X\hat\beta\|^2/(n-p)$, and $A=I_q+X_0(X'X)^{-1}X_0'$.

The natural joint predictor is $X_0\hat\beta$. Its prediction error is $E_0=Y_0-X_0\hat\beta$. Since $\hat\beta-\beta=(X'X)^{-1}X'\varepsilon$,

E_0=\varepsilon_0-X_0(\hat\beta-\beta), \qquad \operatorname{Cov}(E_0)=\sigma^2\left[I_q+X_0(X'X)^{-1}X_0'\right] =\sigma^2 A.

The identity matrix in $A$ represents the new response noise. The other term represents uncertainty from estimating $\beta$. In general, $A$ is not diagonal: prediction errors at different design points share uncertainty from the same fitted coefficients. Also, $A$ is positive definite because for every nonzero $u\in\mathbb{R}^q$,

u'Au=u'u+\left[X_0'u\right]'(X'X)^{-1}\left[X_0'u\right]>0.

Thus $A^{-1}$ exists, even if some rows of $X_0$ are identical.

Worked Example: Covariance at Two Design Points

Suppose a simple linear regression is fitted at $x=-2,-1,0,1,2$, with a column of ones and a column of $x$-values in $X$. Find the matrix $A$ for future responses at $x=0$ and $x=1$.

The design is centered, so $X'X$ is diagonal:

X'X= \begin{pmatrix} 5&0\\ 0&10 \end{pmatrix}, \qquad (X'X)^{-1}= \begin{pmatrix} 1/5&0\\ 0&1/10 \end{pmatrix}.

The future design matrix and the contribution from estimating the regression mean are

X_0= \begin{pmatrix} 1&0\\ 1&1 \end{pmatrix}, \qquad X_0(X'X)^{-1}X_0'= \begin{pmatrix} 1/5&1/5\\ 1/5&3/10 \end{pmatrix}.

Therefore

A= \begin{pmatrix} 6/5&1/5\\ 1/5&13/10 \end{pmatrix}.

The prediction errors have variances $6\sigma^2/5$ and $13\sigma^2/10$, and covariance $\sigma^2/5$. The nonzero covariance is why the joint region cannot generally be constructed by treating the two errors as unrelated.

The Exact Prediction Ellipsoid

Theorem: Under the model and independence assumptions above, let $\nu=n-p>0$. For $0<\alpha<1$, the region $$ \mathcal{P}_{1-\alpha} = \left\{y_0\in\mathbb{R}^q: (y_0-X_0\hat\beta)'A^{-1}(y_0-X_0\hat\beta) \le q s^2 F_{q,\nu}(1-\alpha) \right\} $$ has exact conditional coverage $P(Y_0\in\mathcal{P}_{1-\alpha}\mid X,X_0)=1-\alpha$. Here $F_{q,\nu}(1-\alpha)$ denotes the $1-\alpha$ quantile of an $F$ distribution with $q$ and $\nu$ degrees of freedom.

The quantile is scaled by $q s^2$, not just $s^2$. This factor follows from the standard $F$-statistic convention: an $F_{q,\nu}$ variable is the ratio of two independent chi-squared variables divided by their degrees of freedom.

Why the Calibration Is Exact

Standardize the prediction error.
Because $E_0\sim N_q(0,\sigma^2 A)$ and $A$ is positive definite, $Z=A^{-1/2}E_0/\sigma\sim N_q(0,I_q)$. Consequently, $E_0'A^{-1}E_0/\sigma^2=Z'Z\sim\chi^2_q$.

Identify the residual variance distribution.
Let $H=X(X'X)^{-1}X'$ and $M=I_n-H$. Then the residual vector is $MY$, and its squared length is $\varepsilon'M\varepsilon$, because $MX=0$. The matrix $M$ is a symmetric idempotent projection of rank $\nu=n-p$. Orthogonal diagonalization of $M$, together with the spherical normal distribution of $\varepsilon/\sigma$, gives $\nu s^2/\sigma^2=\varepsilon'M\varepsilon/\sigma^2\sim\chi^2_\nu$.

Establish independence.
The fitted-coefficient error depends on $H\varepsilon$, whereas the residual depends on $M\varepsilon$. These are jointly normal and have zero covariance because $HM=0$, so they are independent. The future error $\varepsilon_0$ is independent of the training errors by assumption. Hence $E_0$ and $s^2$ are independent.

Form the $F$ ratio.
By the preceding steps, $E_0'A^{-1}E_0/\sigma^2$ and $\nu s^2/\sigma^2$ are independent chi-squared variables with $q$ and $\nu$ degrees of freedom. Therefore $$ \frac{E_0'A^{-1}E_0}{q s^2} = \frac{(E_0'A^{-1}E_0/\sigma^2)/q}{(\nu s^2/\sigma^2)/\nu} \sim F_{q,\nu}. $$ The probability that this ratio is at most $F_{q,\nu}(1-\alpha)$ is exactly $1-\alpha$, which is precisely the event that $Y_0$ belongs to the stated ellipsoid. $\square$

This argument also explains why the result is finite-sample exact: it uses normal distributions, orthogonal projections, and an exact $F$ ratio, not a large-sample approximation.

Worked Example: Checking a Future Vector Against the Ellipsoid

Take six regression design points $x=-3,-2,-1,1,2,3$, with an intercept and slope. Then $n=6$, $p=2$, and $\nu=4$. For future responses at $x=0$ and $x=1$, suppose the fitted line is $\hat y=10+2x$ and $s^2=4$. Find the ellipsoid threshold at 95% coverage and check the vector $y_0=(20,12)'$.

Here $X'X$ has diagonal entries $6$ and $28$, since the design values sum to zero and their squared values sum to $28$. Therefore

A= \begin{pmatrix} 7/6&1/6\\ 1/6&101/84 \end{pmatrix}, \qquad \det(A)=\frac{11}{8}.

For $F_{2,4}$, its upper-tail probability at $x\geq0$ is $(1+x/2)^{-2}$. Setting this equal to $0.05$ gives $F_{2,4}(0.95)=2(\sqrt{20}-1)\approx6.944$. Thus the squared-radius threshold is

q s^2 F_{q,\nu}(0.95) =2(4)(6.944)\approx55.55.

The predicted vector is $X_0\hat\beta=(10,12)'$, so the prediction error for the proposed $y_0$ is $e=(10,0)'$. The first diagonal entry of $A^{-1}$ is $(8/11)(101/84)=202/231$, giving

e'A^{-1}e=100\left(\frac{202}{231}\right)\approx87.45>55.55.

Therefore $y_0=(20,12)'$ lies outside the 95% prediction ellipsoid. The conclusion uses the full covariance matrix and the joint threshold, not a separate comparison of each coordinate with an unadjusted interval.

Worked Example: Replicated Future Responses at One Point

Consider an intercept-only model with $n=5$ training observations and $q=3$ future observations at the same design point. Then $X$ is a column of five ones and $X_0$ is a column of three ones. Since $X'X=5$,

A=I_3+\frac{1}{5}\mathbf{1}\mathbf{1}' = \begin{pmatrix} 6/5&1/5&1/5\\ 1/5&6/5&1/5\\ 1/5&1/5&6/5 \end{pmatrix}.

Each future prediction error has variance $6\sigma^2/5$, and every pair has covariance $\sigma^2/5$. The future response errors themselves are independent, but all three prediction errors use the same estimated intercept. That shared estimation uncertainty creates their positive covariance. The exact joint ellipsoid uses this $A$ and the calibration $3s^2F_{3,4}(1-\alpha)$.

Geometry and a Common Pitfall

The ellipsoid is centered at $X_0\hat\beta$. Its shape is determined by $A$: directions with larger prediction-error variance have greater spread, and off-diagonal entries encode dependence among coordinates. The factor $s^2F_{q,\nu}(1-\alpha)$ sets its overall size. If $q=1$, the region reduces to the familiar two-sided exact prediction interval, since $F_{1,\nu}(1-\alpha)=t_{\nu}(1-\alpha/2)^2$.

Common mistake: using only the covariance of the fitted mean. The matrix $X_0(X'X)^{-1}X_0'$ describes uncertainty in estimating the mean-response vector. Using it alone omits $I_q$, which represents the new response noise. The resulting region is too small for predicting actual future observations. A second error is to use separate scalar $t$ cutoffs and call the resulting collection a joint region with the same nominal coverage. Those intervals address coordinatewise coverage; the theorem’s $F$ calibration addresses coverage of the whole vector.

Key takeaway: For multiple future normal responses, use the full matrix $A=I_q+X_0(X'X)^{-1}X_0'$ and the exact threshold $q s^2F_{q,n-p}(1-\alpha)$. The identity term accounts for future noise; the remaining term accounts for fitted-mean uncertainty.

Check Your Understanding

Answer each question using the model and assumptions stated above.

Which part of $A$ accounts for new response noise, and which part accounts for estimating the regression coefficients?
Why can prediction errors at distinct future design points be correlated even when the future response errors are independent?
State the exact distribution of $(Y_0-X_0\hat\beta)'A^{-1}(Y_0-X_0\hat\beta)/(q s^2)$.
What would be omitted if $X_0(X'X)^{-1}X_0'$ were used as the prediction-error covariance factor?
In the exact calibration, what are the numerator and denominator degrees of freedom of the $F$ distribution?