Mean Response and New-Observation Prediction

Two Different Targets at a Predictor Vector

In multivariate regression, a fitted response at a specified predictor vector estimates an average outcome, not the outcome of a particular future case. Those are different targets. The fitted mean is uncertain because the regression coefficients are estimated; a new observation also has its own random error. Separating these two sources of variation is essential when reporting uncertainty or predicting a future multivariate response.

Use the fixed-design model $Y=XB+U$, where $Y$ is an $n\times q$ response matrix, $X$ is a fixed $n\times p$ design matrix of full column rank, $B$ is a $p\times q$ coefficient matrix, and the rows of $U$ are independent, have mean zero, and share covariance matrix $\Sigma$. Write the least-squares estimator as $\hat B=(X^\top X)^{-1}X^\top Y$. At a specified predictor row vector $x_0$ of length $p$, the mean response is $x_0B$, a row vector of length $q$.

Definition: The mean response at $x_0$ is $x_0B$, and its least-squares estimate is $x_0\hat B$. For a future case at the same predictor vector, write its response as $Y_0=x_0B+\varepsilon_0$, where $\varepsilon_0$ is independent of the training errors, has mean zero, and has covariance $\Sigma$. The prediction error is $Y_0-x_0\hat B$.

The covariance of $x_0\hat B$ describes how the estimated mean would vary over repeated training samples, with the design held fixed. The covariance of the prediction error describes variation in the difference between a future response and the fitted mean. It includes both coefficient-estimation uncertainty and the future case's own error.

Covariance of the Estimated Mean Response

The least-squares estimation error has a direct expression in terms of the training errors:

\hat B-B=(X^\top X)^{-1}X^\top U.

Let $A=(X^\top X)^{-1}X^\top$, and let $a_i$ be its $i$th column. Then $x_0(\hat B-B)=\sum_{i=1}^n (x_0a_i)U_i$, where $U_i$ is the $i$th error row. The rows are independent and each has covariance $\Sigma$. Thus the covariance of this weighted sum is the sum of the squared weights times $\Sigma$. The sum of squared weights is $x_0AA^\top x_0^\top=x_0(X^\top X)^{-1}x_0^\top$.

Theorem: Under the fixed-design model above, the covariance of the estimated mean response at $x_0$ is $$ \operatorname{Cov}(x_0\hat B) = \left[x_0(X^\top X)^{-1}x_0^\top\right]\Sigma. $$ The quantity in square brackets is a nonnegative scalar; it scales every entry of the outcome covariance matrix.

This formula requires only the stated first and second moment assumptions and independence of the training error rows; normality is not needed. To see the derivation in matrix form, $A A^\top=(X^\top X)^{-1}$, so

\begin{aligned} \operatorname{Cov}\{x_0(\hat B-B)\} &=\sum_{i=1}^n (x_0a_i)^2\Sigma\\ &=x_0AA^\top x_0^\top\Sigma\\ &=x_0(X^\top X)^{-1}x_0^\top\Sigma. \end{aligned}

The mean of $x_0\hat B$ is $x_0B$, since $E[\hat B]=B$. Therefore, the covariance just derived is also the covariance of the estimated mean response itself. It is a $q\times q$ matrix: its diagonal entries are the variances for the separate outcomes, and its off-diagonal entries describe their estimation-error covariances.

Worked Example: A Two-Outcome Mean Response

Suppose the predictors are an intercept and a centered quantitative predictor with observed values $-2,-1,0,1,2$. Then $X^\top X$ is diagonal with entries 5 and 10, so $(X^\top X)^{-1}=\operatorname{diag}(1/5,1/10)$. Let $x_0=(1,1)$, let $\hat B=\begin{pmatrix}10&20\\2&-1\end{pmatrix}$, and suppose $\Sigma=\begin{pmatrix}4&1\\1&9\end{pmatrix}$. Calculate the estimated mean and its covariance.

First, multiply the predictor vector by the coefficient estimate:

x_0\hat B =(1,1) \begin{pmatrix}10&20\\2&-1\end{pmatrix} =(12,19).

The scalar factor in the covariance is $x_0(X^\top X)^{-1}x_0^\top=1/5+1/10=3/10$. Therefore,

\operatorname{Cov}(x_0\hat B) =\frac{3}{10} \begin{pmatrix}4&1\\1&9\end{pmatrix} = \begin{pmatrix}1.2&0.3\\0.3&2.7\end{pmatrix}.

The fitted mean for the two outcomes is $(12,19)$. Its estimated sampling variances are 1.2 and 2.7, and its estimation errors have covariance 0.3.

Prediction Error for a Future Observation

The future response is not just the unknown mean. It is the mean plus a new error. Subtracting the fitted mean gives $Y_0-x_0\hat B=\varepsilon_0-x_0(\hat B-B)$. The training-sample coefficient error and the future error are independent under the model assumptions. Their cross-covariance is therefore zero, and changing the sign of one term does not change its covariance. Adding the two covariance matrices gives the prediction-error result.

Theorem: If the future error $\varepsilon_0$ is independent of the training errors and has covariance $\Sigma$, then $$ \operatorname{Cov}(Y_0-x_0\hat B) = \left[1+x_0(X^\top X)^{-1}x_0^\top\right]\Sigma. $$

In detail, write $V=x_0(\hat B-B)$. Since $E[\varepsilon_0]=E[V]=0$, independence gives $\operatorname{Cov}(\varepsilon_0,V)=0$. Expanding the covariance of the difference yields

\begin{aligned} \operatorname{Cov}(\varepsilon_0-V) &=\operatorname{Cov}(\varepsilon_0)+\operatorname{Cov}(V) -\operatorname{Cov}(\varepsilon_0,V)-\operatorname{Cov}(V,\varepsilon_0)\\ &=\Sigma+x_0(X^\top X)^{-1}x_0^\top\Sigma\\ &=\left[1+x_0(X^\top X)^{-1}x_0^\top\right]\Sigma. \end{aligned}

The additional $\Sigma$ is the variation of the future case around its own mean. It is present even if the regression coefficients were known exactly. Thus the prediction-error covariance exceeds the mean-response covariance by $\Sigma$, in the positive-semidefinite sense. This comparison concerns covariance matrices; it does not mean every off-diagonal entry must be numerically larger, since entries of $\Sigma$ may be negative.

Worked Example: An Intercept-Only Mean and Prediction

Suppose there are four independent cases, each with an intercept-only predictor, so $X$ is a column of four ones and $X^\top X=4$. At the same predictor vector $x_0=(1)$, suppose $\Sigma=\begin{pmatrix}4&2\\2&16\end{pmatrix}$. Find the covariance of the estimated mean and the prediction error for one new case.

Here $x_0(X^\top X)^{-1}x_0^\top=1/4$. The estimated mean covariance is

\frac14\Sigma = \begin{pmatrix}1&0.5\\0.5&4\end{pmatrix}.

For a new case, add its independent error covariance $\Sigma$:

\left(1+\frac14\right)\Sigma = \frac54 \begin{pmatrix}4&2\\2&16\end{pmatrix} = \begin{pmatrix}5&2.5\\2.5&20\end{pmatrix}.

The covariance of the future response around the estimated mean is larger because it includes the new case's own variation. In particular, the first outcome's variance is 5 for prediction error, not 1, which is the variance of the estimated mean alone.

Using an Estimated Outcome Covariance

In applications, $\Sigma$ is usually unknown. Under the independent, common-covariance model with full-rank $X$, the residual sum-of-squares-and-cross-products matrix is $E=\hat U^\top\hat U$, where $\hat U=Y-X\hat B$. When $n>p$, the usual unbiased estimate is $S=E/(n-p)$. Substituting $S$ for $\Sigma$ gives estimated covariance matrices for the fitted mean and prediction error:

\widehat{\operatorname{Cov}}(x_0\hat B) = \left[x_0(X^\top X)^{-1}x_0^\top\right]S, \qquad \widehat{\operatorname{Cov}}(Y_0-x_0\hat B) = \left[1+x_0(X^\top X)^{-1}x_0^\top\right]S.

These are covariance estimates, not statements that the future response equals its conditional mean. Under multivariate normal errors, the estimator and prediction error also have multivariate normal distributions, with means $x_0B$ and zero respectively and with the corresponding covariances above when $\Sigma$ is known. If $\Sigma$ is estimated, uncertainty from estimating it must also be considered for procedures that require calibrated prediction regions.

Worked Example: Plugging in a Residual Covariance Estimate

Suppose $n=8$, $p=2$, and the residual SSCP matrix is $E=\begin{pmatrix}24&6\\6&36\end{pmatrix}$. At $x_0=(1,2)$, suppose $(X^\top X)^{-1}=\begin{pmatrix}0.25&-0.05\\-0.05&0.10\end{pmatrix}$, and the fitted mean is $(5,2)$. Calculate the estimated covariance matrices.

There are $n-p=6$ residual degrees of freedom, so

S=\frac{E}{6} = \begin{pmatrix}4&1\\1&6\end{pmatrix}.

The scalar predictor factor is $x_0(X^\top X)^{-1}x_0^\top=0.25-0.10+0.40=0.55$? Check the cross terms: each contributes $-0.05(1)(2)=-0.10$, so the correct calculation is $0.25-0.10-0.10+0.40=0.45$.

Thus the estimated mean covariance is

0.45S = \begin{pmatrix}1.8&0.45\\0.45&2.7\end{pmatrix},

and the estimated prediction-error covariance is

1.45S = \begin{pmatrix}5.8&1.45\\1.45&8.7\end{pmatrix}.

Both calculations use the same fitted mean, $(5,2)$; the distinction is in the uncertainty target. The prediction covariance includes one full copy of the estimated outcome covariance $S$ for the future case's new error.

Common Mistake and Practical Interpretation

A frequent error is to use $x_0(X^\top X)^{-1}x_0^\top\Sigma$ as the covariance for a future response minus its fitted mean. That expression accounts only for estimation of the mean. The future observation has its own error, so the correct prediction-error covariance has the additional $\Sigma$ term. Conversely, adding $\Sigma$ when the target is the mean response overstates uncertainty: the target is an average, not a new individual outcome.

These formulas also rely on the stated model. If future errors are dependent on training errors, the cross-covariance terms need not vanish. If error covariance changes across cases, a common $\Sigma$ is not appropriate. The predictor vector must also define an estimable mean; with a rank-deficient design, the corresponding linear function must be estimable for these fitted-mean statements to apply.

Key takeaway: At a predictor vector $x_0$, the fitted mean has covariance $h_0\Sigma$, where $h_0=x_0(X^\top X)^{-1}x_0^\top$. Prediction error for an independent new observation has covariance $(1+h_0)\Sigma$. The added $\Sigma$ is new-observation variation, not coefficient-estimation uncertainty.

Check Your Understanding

Use the distinction between a mean target and a future-observation target to answer these questions.

What is the covariance of the estimated mean response at $x_0$, and what assumptions are needed for the stated formula?
Why does the prediction-error covariance contain an additional copy of $\Sigma$?
If $x_0(X^\top X)^{-1}x_0^\top=0.2$, what are the multipliers of $\Sigma$ for the estimated mean covariance and prediction-error covariance?
For a two-outcome model, what do the off-diagonal entries of either covariance matrix represent?
Why is it incorrect to use the estimated mean covariance as the prediction-error covariance for a future case?