Prediction Intervals for Averages of Replicate Outcomes

Why Predict an Average of Future Responses?

A single future response can vary substantially even when its expected value is estimated precisely. In some settings, however, the target is not one new response but the average of several independent responses collected under the same conditions. Examples include averaging repeated assay measurements or repeated outcomes at a specified treatment setting. Averaging reduces the variability of the future observations, but it does not remove uncertainty about the mean estimated from the training data.

The interval must therefore account for two sources of uncertainty: error in the fitted mean at the design point, and random variation among the future responses around that mean. The number of future replicates reduces the second source, not the first.

Model and Target

Use the fixed-design normal linear model. Let \(\mathbf y\) be the \(n\)-vector of observed responses, and let \(\mathbf X\) be a fixed \(n\times p\) design matrix of full column rank, with \(n>p\). Assume \(\mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon\), where \(\boldsymbol\varepsilon\sim N_n(\mathbf 0,\sigma^2\mathbf I_n)\). The unknown regression coefficient is \(\boldsymbol\beta\), and \(\sigma^2>0\) is the common error variance.

At a fixed design point represented by a \(p\)-vector \(\mathbf x_0\), suppose the future responses are \(Y_{01},\ldots,Y_{0m}\), where \(Y_{0j}=\mathbf x_0^\mathsf T\boldsymbol\beta+\varepsilon_{0j}\). Assume the future errors are mutually independent, each has distribution \(N(0,\sigma^2)\), and they are independent of the training errors. The target is their average \(\bar Y_0=m^{-1}\sum_{j=1}^m Y_{0j}\), not the unknown mean by itself.

Definition: The prediction error for the average of \(m\) future replicates at \(\mathbf x_0\) is \(\bar Y_0-\mathbf x_0^\mathsf T\hat{\boldsymbol\beta}\), where \(\hat{\boldsymbol\beta}=(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf X^\mathsf T\mathbf y\) is the least-squares estimator. Define the leverage for the fitted mean at this point by \(h_0=\mathbf x_0^\mathsf T(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf x_0\).

The fitted mean \(\mathbf x_0^\mathsf T\hat{\boldsymbol\beta}\) has variance \(\sigma^2h_0\). Meanwhile, averaging \(m\) independent future errors gives variance \(\sigma^2/m\). Because the future errors are independent of the training data, these contributions add in the prediction error:

\operatorname{Var}\!\left(\bar Y_0-\mathbf x_0^\mathsf T\hat{\boldsymbol\beta}\right) =\sigma^2\left(h_0+\frac{1}{m}\right).

This variance already shows the role of replication. As \(m\) increases, the future-average noise shrinks, while the uncertainty in the fitted mean remains.

Exact Prediction Interval

Write \(\nu=n-p\) for the residual degrees of freedom, and define the residual variance estimate and residual standard error by \(s^2=\|\mathbf y-\mathbf X\hat{\boldsymbol\beta}\|^2/\nu\) and \(s=\sqrt{s^2}\). The exact interval uses a Student \(t\) quantile, because \(\sigma\) is estimated from the same training data.

Theorem: Under the fixed-design normal linear model above, let \(0<\alpha<1\), and let \(t_{\nu,1-\alpha/2}\) denote the \(1-\alpha/2\) quantile of a Student \(t\) distribution with \(\nu=n-p\) degrees of freedom. Then

\bar Y_0\ \in\ \mathbf x_0^\mathsf T\hat{\boldsymbol\beta} \ \pm\ t_{\nu,1-\alpha/2}\,s\sqrt{h_0+\frac{1}{m}}

has coverage probability exactly \(1-\alpha\) for the future average \(\bar Y_0\).

Proof. Let \(\mathbf P=\mathbf X(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf X^\mathsf T\) be the projection onto the column space of \(\mathbf X\), and let \(\mathbf M=\mathbf I_n-\mathbf P\). Then \(\hat{\boldsymbol\beta}-\boldsymbol\beta=(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf X^\mathsf T\boldsymbol\varepsilon\). Thus the fitted-mean error at \(\mathbf x_0\) is normal with mean zero and variance \(\sigma^2h_0\). The average of the future errors is independent of this error, normal with mean zero, and has variance \(\sigma^2/m\). Consequently,

\frac{\bar Y_0-\mathbf x_0^\mathsf T\hat{\boldsymbol\beta}} {\sigma\sqrt{h_0+1/m}}\sim N(0,1).

The residual sum of squares is \(\boldsymbol\varepsilon^\mathsf T\mathbf M\boldsymbol\varepsilon\). Since \(\mathbf M\) is a symmetric idempotent matrix of rank \(\nu\), normality gives \(\nu s^2/\sigma^2\sim\chi^2_\nu\). Also, \(\mathbf P\mathbf M=\mathbf 0\). Orthogonal projections of a multivariate normal vector are independent when their covariance is zero, so the fitted-mean error is independent of \(s^2\). The future errors are independent of both. Therefore the numerator, standardized by its true standard deviation, is independent of \(s^2\), and

\frac{\bar Y_0-\mathbf x_0^\mathsf T\hat{\boldsymbol\beta}} {s\sqrt{h_0+1/m}} = \frac{\displaystyle \frac{\bar Y_0-\mathbf x_0^\mathsf T\hat{\boldsymbol\beta}} {\sigma\sqrt{h_0+1/m}}} {\displaystyle \sqrt{\frac{\nu s^2/\sigma^2}{\nu}}} \sim t_\nu.

The probability that a \(t_\nu\) random variable lies between \(-t_{\nu,1-\alpha/2}\) and \(t_{\nu,1-\alpha/2}\) is \(1-\alpha\). Rearranging that event gives the stated interval. \(\square\)

Worked Calculations

Worked Example: Replicates in an Intercept-Only Model

Suppose \(n=10\) independent training observations follow an intercept-only normal model. The fitted mean is \(\bar y=50\), and the residual standard error is \(s=4\). Predict the average of \(m=4\) future observations with 95% coverage.

Here \(p=1\), so \(\nu=10-1=9\). In an intercept-only model, \(h_0=1/n=0.1\), and the 95% critical value is \(t_{9,0.975}=2.262\). The standard error for predicting the future average is

s\sqrt{h_0+\frac{1}{m}} =4\sqrt{0.1+\frac{1}{4}} =4\sqrt{0.35} \approx 2.366.

The margin of error is \(2.262(2.366)\approx5.352\). Therefore the interval is \(50\pm5.352\), or approximately \((44.648,55.352)\). The \(1/4\) term reflects the noise in averaging four future observations; the \(0.1\) term remains because the population mean was estimated from only ten training observations.

Worked Example: A Replicate Average at a Regression Design Point

A simple linear regression is fitted using \(n=8\) observations at predictor values \(0,1,\ldots,7\). Suppose the fitted coefficients are intercept \(10\) and slope \(2\), and the residual standard error is \(s=2.5\). At predictor value \(x_0=3.5\), predict the average of \(m=9\) future responses.

The design point is at the mean predictor value, \(3.5\). For this design, \(\sum_i(x_i-\bar x)^2=42\), so the leverage for the fitted mean is \(h_0=1/8+(3.5-3.5)^2/42=1/8\). The fitted mean is \(10+2(3.5)=17\). There are \(\nu=8-2=6\) residual degrees of freedom, and \(t_{6,0.975}=2.447\). The prediction standard error is

2.5\sqrt{\frac18+\frac19} =2.5\sqrt{\frac{17}{72}} \approx 2.5(0.4859) \approx1.215.

The margin is \(2.447(1.215)\approx2.973\), giving the 95% prediction interval \(17\pm2.973\), or approximately \((14.027,19.973)\). Even though nine future responses are averaged, the \(1/8\) contribution from estimating the mean is not divided by nine.

Worked Example: Increasing the Number of Future Replicates

Suppose a fitted regression has \(n=16\), \(p=2\), residual standard error \(s=3\), and leverage \(h_0=1/16\) at the design point of interest. The fitted mean there is \(20\). Compare 95% intervals for the average of one future response and the average of 25 future responses.

Both intervals use \(\nu=16-2=14\) and \(t_{14,0.975}=2.145\). For one future response, the standard error is \(3\sqrt{1/16+1} = 3\sqrt{1.0625}\approx3.092\), so the margin is \(2.145(3.092)\approx6.633\). The interval is approximately \((13.367,26.633)\).

For 25 future responses, the standard error is \(3\sqrt{1/16+1/25}=3\sqrt{0.1025}\approx0.960\). The margin is \(2.145(0.960)\approx2.060\), giving approximately \((17.940,22.060)\). The interval narrows substantially, but it does not collapse to the fitted value \(20\): its width still reflects uncertainty in estimating the mean.

What Replication Can and Cannot Reduce

The interval half-width is \(t_{\nu,1-\alpha/2}s\sqrt{h_0+1/m}\). For fixed training data, design point, and confidence level, it decreases as \(m\) increases, because \(1/m\) decreases. But it cannot decrease below the half-width associated with the fitted-mean uncertainty alone.

Result: Holding the training data, \(h_0\), and confidence level fixed, the prediction half-width decreases with \(m\) and converges to \(t_{\nu,1-\alpha/2}s\sqrt{h_0}\) as \(m\to\infty\). Indeed, \(1/m\) decreases to zero, and the square-root function is continuous on the positive numbers. Thus \(\sqrt{h_0+1/m}\to\sqrt{h_0}\).

This limiting statement concerns averaging more future replicates at the same design point; it does not say that the fitted mean becomes known. The uncertainty in that fitted mean comes from the finite training sample. The pointwise mean-response interval has standard-error factor \(s\sqrt{h_0}\), whereas the prediction interval for an average of \(m\) future responses has factor \(s\sqrt{h_0+1/m}\). The prediction interval is therefore wider for every finite \(m\), and approaches the mean-response interval's width as \(m\) grows.

Common Mistakes

Dividing the leverage by \(m\). Replication reduces the variance of future noise to \(\sigma^2/m\); it does not average away the error in the fitted mean. The correct factor is \(h_0+1/m\), not \((h_0+1)/m\).
Using the single-response noise term for an average. For one future response, \(m=1\), so the factor is \(h_0+1\). For an average of \(m\) independent future responses, the factor is \(h_0+1/m\).
Using a mean-response interval for a future average. The mean-response interval includes fitted-mean uncertainty but omits future noise. For any finite \(m\), a prediction interval must include the additional \(1/m\) term.
Assuming exact coverage without the model assumptions. The stated t calibration depends on normal errors, a correctly specified linear mean, common variance, and independence of training and future errors. Without these assumptions, the formula need not have exact finite-sample coverage.

Key takeaway: To predict an average of \(m\) independent future responses at a fixed design point, center the interval at the fitted mean and use the standard-error factor \(s\sqrt{h_0+1/m}\), with a \(t_{n-p}\) critical value. Replication reduces future-response noise, while uncertainty in the fitted mean remains.

Check Your Understanding

Use the model assumptions and formulas above to answer these questions.

What are the two components of the variance of the prediction error for an average of \(m\) future responses?
In the exact interval, which part of the standard-error factor changes with the number of replicates, and which part does not?
If \(m=1\), what does the prediction standard-error factor become?
Why does the interval width approach a positive fitted-mean uncertainty limit as \(m\) grows, when \(h_0>0\)?
Name two assumptions needed for the stated exact t coverage.