Prediction and confidence intervals - large number of predictions

Question

I'm using a regression model to predict one quantity, $y$, given another, $x$. I'm trying to estimate the error in future predictions of $y$, but I'm wondering in which scenarios I can fairly use the prediction interval (PI) vs. the confidence interval (CI). Let's start with this from Reddy (2011):

A regression equation can be used to predict future values of $y$ provided the $x$ value is within the domain of the original data from which the model was identified. One differentiates between the two types of predictions (...):

(a) mean response or standard error of regression where one would like to predict the mean value of $y$ for a large number of repeated $x_0$ values. The mean value is directly deduced from the regression equation while the variance is:

$$ \sigma^2(\hat{y}_0) = MSE\left[\frac{1}{n} + \frac{(x_0 -\bar{x})^2}{s_{xx}}\right] $$

(b) individual or specific response or standard error of prediction where one would like to predict the specific value of y for a specific value $x_0$. This error is larger than the error in the mean response by an amount equal to the RMSE. Thus,

$$ \sigma^2(\hat{y}_0) = MSE\left[1 + \frac{1}{n} + \frac{(x_0 -\bar{x})^2}{s_{xx}}\right] $$

Finally, the 95% CL for the individual response at level $x_0$ is:

$$ y_0 = \hat{y}_0 \pm t_{0.05/2} \cdot\sigma(\hat{y}_0) $$

Given the above, here's what I interpret:

1. To express the error in a single predicted value (one $y$ value given one $x_0$) , use the PI.
2. If a large number of predictions are made for a single $x_0$, you can use the CI.

Questions: Suppose I am summing up a large number of $y$ predictions, for a wide range of $x$ values that are never quite the same:

$$ y_{tot} = \sum \hat{y}(x_i) $$

and suppose I estimate the error in $y_{tot}$ by assuming:

$$ e_{y_{tot}} = \sum e_\hat{y} $$

where $e_\hat{y}$ is the interval $t_{0.05/2} \cdot\sigma(\hat{y}_0)$. Given that I'm making a large number of $y$ predictions,

1. Would it be fair to consider this a "mean" response and use the CI rather than PI, even if it is not a repeated prediction of the same value of $x_0$? The PI seems like an overly conservative estimate for this application.

2. Is my assumption of simply summing errors in $y$ (not in quadrature) an ok one? We typically do this in my field to be conservative, as summing in quadrature assumes independence in the errors, correct? Does this depend on the model's residuals and if they're normally distributed?

Reference:

Reddy, T. A. (2011). Applied data analysis and modeling for energy engineers and scientists. Springer Science & Business Media.

Answer 1

The model is based on a vector of responses $y$ associated with explanatory variables arranged into a parallel "model matrix" $X$ with one row per observation. Assuming, as is usual, that the responses differ from "true" values $X\beta$ by independent, identically distributed errors with zero mean and (unknown) variance $\sigma^2,$ the Ordinary Least Squares estimate of the parameters is

$$\hat\beta = (X^\prime X)^{-1}X^\prime y\tag{*}$$

(at least when $X$ is of full rank).

Similarly, assume the future responses $z$ are governed by the same model and will correspond to a matrix $Z,$ again with one row for each set of explanatory values of the future responses.

To complete the specification of the problem, let $\gamma$ represent a linear combination of the future responses, which will be the (random) number $\gamma^\prime z.$

According to the OLS estimate $(*),$ the individual predictions (at the rows of $Z$) are the random variables

$$\hat z = Z\hat\beta + e$$

where $e$ is the vector of (unobservable) errors. The expectation of the desired linear combination of these predictions is

$$E[\gamma^\prime\hat z] = E[\gamma^\prime (Z\hat\beta + e)] = Z\gamma^\prime E[\hat\beta] + \gamma^\prime E[e] = \gamma^\prime Z\beta + \gamma^\prime 0 = \gamma^\prime Z\beta,$$

the $\gamma$ linear combination of the "true" responses $Z\beta.$ Its variance is

$$\operatorname{Var}(\gamma^\prime\hat z) = \operatorname{Var}(\gamma^\prime( Z\hat\beta + e)).\tag{**}$$

Because $e$ is assumed independent of the $y,$ it is independent of $\gamma^\prime Z\hat \beta,$ permitting us to sum these variances separately. Since the components of $e$ are independent,

$$\operatorname{Var}(\gamma^\prime e) = \gamma^\prime \operatorname{Var}(e) \gamma = \sigma^2\gamma^\prime \gamma = \sigma^2\,|\gamma|^2.$$

For instance, when $Z$ is $k$ additional values to be summed, $\gamma^\prime = (1,1,\ldots,1)$ (a vector with $k$ entries) this variance reduces to $k\sigma^2.$

Finally, the variance of the first term in $(**)$ is

$$\operatorname{Var}(\gamma^\prime Z \hat \beta) = \gamma^\prime Z \operatorname{Var}(\hat\beta) Z^\prime \gamma.$$

The variance of the coefficient estimate $\hat\beta$ is obtained, as usual, from $(*)$ as

$$\operatorname{Var}(\hat\beta) = (X^\prime X)^{-1}\sigma^2.$$ This gives a useful expression for the first variance term,

$$\operatorname{Var}(\gamma^\prime Z \hat \beta) = \gamma^\prime Z (X^\prime X)^{-2} Z^\prime \gamma\, \sigma^2.$$

Whence the standard error of the predicted random variable $\hat z$ is

$$\operatorname{se}(\hat z)^2 = (\gamma^\prime Z (X^\prime X)^{-2} Z^\prime \gamma + \gamma^\prime \gamma)\sigma^2.$$

An approximate $1-\alpha$ interval (exact when the errors are Normally distributed) is thereby obtained by setting limits at distances $t_{\alpha/2, \nu}\operatorname{se}(\hat z)$ on either side of the predicted value $\gamma^\prime Z\hat\beta,$ after substituting an estimate $\hat\sigma^2$ of $\sigma^2.$ (The usual OLS estimate is the unbiased one, but this method applies even to other estimates, such as the maximum likelihood estimate (the mean squared residual).) The term $\nu$ is the "degrees of freedom," equal to the number of rows of $X$ minus its rank, and $t$ represents the percentile of that Student t distribution. You can use the standard Normal distribution when $\nu$ is large; typically, $\nu \ge 20$ permits this.

To interpret this interval, contemplate the entire procedure of observing $y,$ fitting the model, constructing the prediction interval from that, observing $z,$ and comparing $\gamma^\prime z$ to that interval. The chance this value lies in the interval (is "covered by" the interval) is approximately $1-\alpha.$ Thus, the risk that the interval fails to cover $\gamma^\prime z$ is only $\alpha$ -- and this risk doesn't depend on actually observing $z.$ Thus, if all you have available are the data $X$ and $y,$ in the sense just described you can have $100(1-\alpha)\%$ "confidence" that $\gamma^\prime z$ lies in your interval.

Here is an example. Observations were made of the model shown at the left at ten equally spaced points $x_i$ between $0$ and $1.$ The "future" observations $z$ were made at the $k=5$ red points. The chosen linear combination was, as in the question, the sum of their values.

In a simulation, 5000 datasets were created according to the OLS model (using Normal errors with $\sigma=1/4$) and, independently, 5000 sets of values of $z$ were also randomly created. In each iteration a $100(1-0.05) = 95\%$ prediction interval for the sum of the $z$ was computed and compared to the observed sum. The plot at the right shows (1) the true sum (horizontal dark line), which was constant throughout; (2) the lower (blue) and upper (red) prediction limits for the first 500 iterations, sorted in decreasing order; (3) the predicted values $\gamma^\prime\hat z$ (cyan); and (4) the observed values of the $\gamma^\prime z,$ displayed as gray dots. The ones that fell beyond their prediction interval are shown darker and larger than the others.

In this simulation, 94.3% of the prediction intervals covered the values of $\gamma^\prime z.$

Much can be learned by varying the specifics of this simulation. Here, then, is the R code that produced it. It shows how all the formulas can be expressed in a matrix-oriented programming language. All the work is done in just three lines of code, following the comment "The prediction of the sum and its variance."

(I set the number of iterations in the simulation to 5,000 because that takes less than a second to complete and provides a precise estimate of the coverage.)

#
# The model.
#
f <- function(x, b=beta) c(cbind(1, x) %*% b)
generate <- function(x, b=beta) {
  y <- f(x, b)
  y + rnorm(length(y), 0, sigma)
}
#
# True coefficients.
#
beta <- c(2, -1)
ff <-  ~ .
#
# Prediction interval size.
#
alpha <- 0.05
#
# Common SD.
#
sigma <- 1/4
#
# Number of future observations
#
k <- 5
set.seed(17)
#
# X values.
#
n <- 10
X.df <- data.frame(x = seq(0, 1, length.out=n))
X <- model.matrix(ff, X.df)
df <- nrow(X) - ncol(X)
t.alpha <- qt(alpha/2, df)
#
# X values for future observations.
#
zlim <- c(-1/2, 3/2)
Z.df <- data.frame(x = runif(k, zlim[1], zlim[2]))
Z <- model.matrix(ff, Z.df)
#
# The linear combination to use.
#
gamma <- rep(1, nrow(Z.df))
#
# Plot the x values and model.
#
xlim <- range(c(X.df$x, Z.df$x))
x <- seq(xlim[1], xlim[2], length.out=3)
y <- f(x)
par(mfrow=c(1,2))
plot(x, y, type="l", main="Model and Sample Points")
points(X.df$x, f(X.df$x), pch=21, bg=gray(0.5))
points(Z.df$x, f(Z.df$x), pch=21, bg="Red")

sim <- replicate(5e3, {
  #
  # The fit.
  #
  y <- generate(X.df$x)
  fit <- lm.fit(X, y)
  fit.s <- summary(fit)
  #
  # The prediction of the sum and its variance.
  #
  z.hat <- sum(gamma * f(Z.df$x, fit$coefficients))
  se <- c(sqrt(crossprod(solve(crossprod(X), t(Z) %*% gamma)) + nrow(Z))) * sigma
  #
  # The prediction interval.
  #
  PI <- t.alpha * se * c(1, -1) + z.hat
  #
  # The future observations.
  #
  z.sum <- sum(gamma * generate(Z.df$x))
  c(Covers = prod(PI - z.sum) <= 0,
    LPL = PI[1],
    UPL = PI[2],
    Prediction = z.hat,
    Actual = z.sum,
    se = se)
})
#
# Summarize the results.
#
m <- mean(sim["Covers", ])
s <- sqrt(m * (1-m) / length(sim))
c(`1-alpha`=1-alpha, Coverage=m, `SE(Coverage)`=s)

if (ncol(sim) > 500) sim <- sim[, 1:500] # Limit the display
i <- order(sim["UPL", ])# - sim["LPL", ])
sim <- sim[, rev(i)]

plot(c(1, ncol(sim)), range(sim[c("LPL", "UPL", "Prediction"), ]), type="n",
     sub=paste0(round(100*m, 1), "% Coverage"),
     xlab="Iteration", ylab="Value",
     main="Prediction Intervals")
abline(h = sum(gamma * f(Z.df$x)), lwd=2, col="Black")
points(sim["Prediction", ], pch=19, 
       cex=0.5, col=hsv(1/2, alpha=1/8))
points(sim["UPL", ], pch="-", col="Red")
points(sim["LPL", ], pch="-", col="Blue")
points(sim["Actual", ], pch=19, 
       cex=ifelse(sim["Covers",], 0.5, 0.8),
       col=ifelse(sim["Covers",], gray(0, alpha=0.125), gray(0)))

par(mfrow=c(1,1))