Worked Prediction Ellipsoids for Replicated Outcomes

Why Predict Replicated Outcomes Jointly?

Suppose a fitted model is used to predict several future responses, all measured at the same design point. Their future errors are independent, but their predictions share the same estimated regression coefficients. Consequently, the prediction errors are correlated. A joint prediction ellipsoid accounts for both sources of variation and gives a region for the entire vector of future responses, not just separate ranges for each response.

The construction below uses the fixed-design normal linear model. The exact cutoff depends on estimating the error variance from the original sample, so it is a finite-sample \(F\) cutoff rather than a normal quantile.

Model and Exact Prediction Region

Definition: Let \(Y=X\beta+\varepsilon\), where \(Y\) is an \(n\)-vector, \(X\) is a fixed \(n\times p\) matrix of rank \(p\), and \(\varepsilon\sim N(0,\sigma^2 I_n)\). Assume \(n>p\). For \(m\) future responses at fixed design rows collected in an \(m\times p\) matrix \(X_*\), assume \(Y_*=X_*\beta+\varepsilon_*\), where \(\varepsilon_*\sim N(0,\sigma^2 I_m)\) and is independent of \(\varepsilon\). Write \(\hat\beta=(X'X)^{-1}X'Y\), \(s^2=(Y-X\hat\beta)'(Y-X\hat\beta)/(n-p)\), and \(A=I_m+X_*(X'X)^{-1}X_*'\).

The \(I_m\) term in \(A\) represents the independent future errors. The other term represents uncertainty from estimating the regression coefficients. The matrix \(A\) is positive definite: for any nonzero \(u\in\mathbb R^m\), \(u'Au=u'u+(X_*'u)'(X'X)^{-1}(X_*'u)>0\).

Theorem: Under the assumptions above, for \(0<\alpha<1\), an exact conditional \(1-\alpha\) prediction region for \(Y_*\), given the fixed design matrices \(X\) and \(X_*\), is

$$ \left\{y_*: (y_*-X_*\hat\beta)'A^{-1}(y_*-X_*\hat\beta) \leq m s^2 F_{m,n-p}(1-\alpha) \right\}, $$ where \(F_{m,n-p}(1-\alpha)\) is the \(1-\alpha\) quantile of the \(F\) distribution with \(m\) and \(n-p\) degrees of freedom.

The condition \(n>p\) is essential: it ensures positive residual degrees of freedom and makes \(s^2\) and the stated \(F\) cutoff well-defined.

Why the Cutoff Is Exact

The prediction error is \(Y_*-X_*\hat\beta=\varepsilon_*-X_*(X'X)^{-1}X'\varepsilon\). Its mean is zero, and independence of the original and future errors gives

\operatorname{Cov}(Y_*-X_*\hat\beta) =\sigma^2\left[I_m+X_*(X'X)^{-1}X_*'\right] =\sigma^2 A.

Because the errors are jointly normal, this linear transformation is multivariate normal. Since \(A\) is positive definite, whitening by \(A^{-1/2}\) gives \(A^{-1/2}(Y_*-X_*\hat\beta)/\sigma\sim N(0,I_m)\). Its squared Euclidean length therefore has a \(\chi^2_m\) distribution.

To establish independence from the estimated variance, let \(P_X=X(X'X)^{-1}X'\). The residual vector is \((I_n-P_X)\varepsilon\). The covariance between the part of the prediction error involving \(\varepsilon\) and this residual vector is zero, since \(X'(I_n-P_X)=0\). The future error is independent of both. These quantities are jointly normal, so zero covariance implies independence. Also, \(I_n-P_X\) is an idempotent projection of rank \(n-p\); hence \((n-p)s^2/\sigma^2\sim\chi^2_{n-p}\).

Thus the squared whitened prediction error and the scaled residual sum of squares are independent chi-square variables. Their ratio, each divided by its degrees of freedom, has an \(F_{m,n-p}\) distribution:

\frac{(Y_*-X_*\hat\beta)'A^{-1}(Y_*-X_*\hat\beta)} {m s^2} \sim F_{m,n-p}.

Taking the \(1-\alpha\) quantile of this pivot gives the stated region, with conditional coverage exactly \(1-\alpha\). \(\square\)

Replicates at One Design Point

When all \(m\) future responses share the same design row \(x_0'\), \(X_*\) consists of \(m\) copies of \(x_0'\). Define the leverage \(h_0=x_0'(X'X)^{-1}x_0\). Then \(A=I_m+h_0\mathbf 1\mathbf 1'\), where \(\mathbf 1\) is the \(m\)-vector of ones. The off-diagonal entries \(h_0\) describe the shared coefficient-estimation contribution to covariance; the diagonal entries \(1+h_0\) include both that contribution and the individual future error.

The geometry can be seen by splitting any prediction-error vector \(d\) into its average component and its deviations from that average. Write \(\bar d=m^{-1}\mathbf 1'd\) and \(d_{\perp}=d-\bar d\mathbf 1\), so \(\mathbf 1'd_{\perp}=0\). The matrix \(A\) has eigenvalue \(1+mh_0\) in the direction of \(\mathbf 1\), and eigenvalue \(1\) in every direction perpendicular to \(\mathbf 1\). Therefore,

d'A^{-1}d =d_{\perp}'d_{\perp}+\frac{m\bar d^{\,2}}{1+mh_0}.

This decomposition follows by applying \(A\) to each component: \(A\mathbf 1=(1+mh_0)\mathbf 1\), while \(Ad_{\perp}=d_{\perp}\). The common shift of all replicates is discounted by the larger uncertainty in their shared fitted mean; differences among replicates are not affected by coefficient estimation at that common design point.

Worked Example: Two Replicates at One Design Point

Fit a straight-line model with intercept and slope using \(x=-1,0,1,2\) and responses \(y=9,9,11,15\). Predict two future responses at \(x_0=1\), using a 90% prediction region.

The design matrix has columns \(1\) and \(x\), so \(X'X=\begin{pmatrix}4&2\\2&6\end{pmatrix}\) and \((X'X)^{-1}=\frac1{20}\begin{pmatrix}6&-2\\-2&4\end{pmatrix}\). The fitted line is \(\hat y=10+2x\): its fitted values at the observed \(x\)'s are \(8,10,12,14\), with residuals \(1,-1,-1,1\). Their sum is zero and their product with the \(x\)-values sums to \(-1+0-1+2=0\), confirming the least-squares residual orthogonality. The residual sum of squares is \(4\), so \(n-p=2\) and \(s^2=4/2=2\).

At \(x_0=1\), the fitted mean is \(12\) and the leverage is \(h_0=(1,1)(X'X)^{-1}(1,1)'=0.3\). For two replicates,

A=I_2+0.3\mathbf 1\mathbf 1' =\begin{pmatrix}1.3&0.3\\0.3&1.3\end{pmatrix}.

The 90% cutoff uses \(F_{2,2}(0.90)=9\). Indeed, the \(F_{2,2}\) distribution has CDF \(G(q)=q/(1+q)\), so \(G(9)=9/10\). The ellipsoid is centered at \((12,12)'\), and its right-hand cutoff is \(m s^2 F=2\cdot2\cdot9=36\). The eigenvalues of \(A\) are \(1.6\) in the common direction and \(1\) in the contrast direction. Thus the semiaxis lengths are \(\sqrt{36(1.6)}=\sqrt{57.6}\approx7.59\) and \(\sqrt{36}=6\), respectively.

Worked Example: Two Future Responses at Different Design Points

Use the same fitted line and residual variance as above, but predict one response at \(x=-1\) and one at \(x=2\), again with 90% coverage. The fitted center is \((8,14)'\). The two design rows are \((1,-1)\) and \((1,2)\). Using the displayed inverse of \(X'X\), their leverage values are both \(0.7\), and their cross-product through \((X'X)^{-1}\) is \(-0.2\). Therefore,

A=\begin{pmatrix}1.7&-0.2\\-0.2&1.7\end{pmatrix}, \qquad \det(A)=1.7^2-0.2^2=2.85=\frac{57}{20}, \qquad A^{-1}=\frac1{57}\begin{pmatrix}34&4\\4&34\end{pmatrix}.

For the proposed future vector \(y_*=(9,12)'\), the deviation from the center is \(d=(1,-2)'\). Its ellipsoid distance is

d'A^{-1}d =(1,-2)\frac1{57}\begin{pmatrix}34&4\\4&34\end{pmatrix} \begin{pmatrix}1\\-2\end{pmatrix} =\frac{154}{57}\approx2.702.

The cutoff remains \(36\), so this proposed vector lies inside the 90% joint prediction ellipsoid. The calculation uses the covariance between the two prediction errors; treating the two responses as if their prediction errors were independent would omit the off-diagonal entries of \(A\).

Worked Example: Inside Both Coordinate Ranges, Outside the Ellipsoid

Return to the two replicates at \(x=1\). Consider the proposed future vector \((17,7)'\), whose deviation from the center \((12,12)'\) is \(d=(5,-5)'\). This is entirely a contrast direction: its components sum to zero. In that direction \(A\) has eigenvalue \(1\), so \(d'A^{-1}d=d'd=25+25=50\). Since \(50>36\), the vector is outside the 90% prediction ellipsoid.

Yet both coordinates are within the coordinate ranges obtained by projecting the ellipsoid onto either axis. Each diagonal entry of \(A\) is \(1.3\), so the projected half-width is \(\sqrt{36(1.3)}=\sqrt{46.8}\approx6.84\). The coordinate deviations \(5\) and \(-5\) each fall within that half-width. This confirms that lying within every projected coordinate range does not guarantee lying inside the joint region.

Common Mistake and Practical Interpretation

A common mistake is to create \(m\) separate prediction intervals and assume that their combined rectangle is the joint prediction region. Separate intervals address one coordinate at a time; they do not use the ellipsoid’s covariance structure or provide the stated joint coverage. Conversely, checking that every coordinate falls within the ellipsoid’s projected ranges is not sufficient to establish that the whole vector is inside it, as the last example demonstrates. To check joint membership, calculate the single quadratic form \(d'A^{-1}d\) and compare it with \(m s^2F_{m,n-p}(1-\alpha)\).

The ellipsoid is useful when the scientific target is a batch of future replicate outcomes—for example, several measurements at a fixed treatment setting—and the joint behavior of that batch matters. Replication reduces uncertainty about the average of the future outcomes, but it does not remove the individual future errors. The covariance matrix and the exact cutoff keep those two sources of uncertainty distinct.

Key takeaway: For a vector of future responses, combine future-error variation and coefficient-estimation variation in \(A\), then use the \(F_{m,n-p}\) cutoff when \(\sigma^2\) is estimated. The resulting quadratic inequality describes joint, not coordinate-by-coordinate, coverage.

Check Your Understanding

Use the model assumptions and formulas in this tutorial to answer the following.

Why is \(n>p\) needed for the exact prediction region when the error variance is estimated?
For two replicates at a design point with leverage \(h_0\), what are the eigenvalues of \(A\) in the common and contrast directions?
What does the identity term in \(A=I_m+X_*(X'X)^{-1}X_*'\) represent?
For the vector \((17,7)'\) in the final worked example, why does being within both projected coordinate ranges not imply membership in the joint ellipsoid?
Which degrees of freedom appear in the finite-sample \(F\) cutoff, and what determines each one?