Is the variation in the residual standard deviation (on sample) accounted for when one builds a prediction interval (PI)?

Question

This question is somehow related to Is the residual, e, an estimator of the error, $\epsilon$?

I also found some information here: Confidence interval of RMSE

Let's say, I got a model that explains sampled X by mean(X).

Using suggestions from Mr. @whuber I reproduced calculations of the PI from predict.lm in R.

## lm PI

x_dat <- data.frame(x = rnorm(100, 0, 1))

lm_model <- lm(data = x_dat, x ~ 1)

summary(lm_model)

lm_preds <- predict.lm(lm_model
                       , x_dat
                       , interval = "prediction"
                       , level = 0.95
                       )

beta.hat <- lm_model$coefficients

# var-covariance matrix

V <- vcov(lm_model)

# mean prediction

Xp <- model.matrix(~ 1, x_dat)
pred <- as.numeric(Xp %*% beta.hat)[1]

# mean prediction variance:

pred_var <- unname(rowSums((Xp %*% V) * Xp))[1]

# confidence

t_stat <- qt((1 - 0.95)/2, df = lm_model$df.residual)

# residual MSE

res_mse <- sum(lm_model$residuals ^ 2) / lm_model$df.residual

# PI

PI <- pred + c(t_stat, -t_stat) * sqrt(pred_var + res_mse)

print(PI)

print(lm_preds[1, ])

> print(PI)
[1] -2.043175  2.572508
> 
> print(lm_preds[1, ])
       fit        lwr        upr 
 0.2646666 -2.0431747  2.5725079 

I only have 2 questions.

Is it right that we assume that the true model error variance is that of sampled residual variance in order to make an unbiased estimate?

Given that we don't know what the true error variance is, does it mean we make a biased estimate of PI, in particular, by not adjusting for an additional dispersion of the residual variance? If so, can we assume that variance of the residual variance follows chi-square distribution in order to get an upper quantile of the value that can be supplied inside for an exact calculation?

predict.lm(...,
        pred.var = 
)

Answer 1

Referencing the initial post, for Q1:  yes, the error variance is estimated (without bias if the model is correct) by an appropriate adjustment to the sample residuals. The bias is corrected with division by a value smaller than $df=n-1$, based on the parameters in the model. (E.g., this is a question students ask a lot: ¿why do we have to divide by $n-2$ instead of $n-1$ for the variance/standard deviation of the residuals for a bivariate model?)  Thus, we have a “reasonable” estimate for the variance

for Q2:  First, if there were bias, then the entire interval would be shifted by the bias...so adjustment via chi-square bounds probably is not appropriate.  Second, to account for the amount of uncertainty, the prediction interval (PI) is slightly elongated.  Again, this references back to the smaller df used to obtain the critical value for this margin of error.

I hope this answers the query...happy to provide more detail if it seems appropriate.

Addendum #1
The standard error for a predicted value from a regression (using bivariate example here) is $$\hat\sigma \sqrt{\frac{1}{n} + \frac{(x_0-\bar{x})^2}{(n-1)s_x^2}}$$ where $\hat\sigma^2$ is the estimated residual variance, $n$ is the sample size, $x_0$ is the value at which you want the prediction, and $s_x^2$ is the variance of the $x$ values.  But to clarify, this is the standard error for the predicted value $\hat{y}$ at $x_0$.  As mentioned in Kleinbaum et al. (2014), this can be used to calculate the confidence interval for the parameter (which in this case is the average predicted value at $x_0$).  To account for the fact that we often have to use an estimate for $\sigma$ ($\hat\sigma$ in this case), to obtain the margin of error, we multiple this by a slightly larger multiplication factor to accommodate that extra element of uncertainty.  (In this case, we use a smaller $df=n-2$ which in turn makes $t_{c.v.}$ slightly larger.)

When we move to the next level of estimation, and ask for an interval in which we might observe values of $y$ at $x_0$, then we need to account for the variance of the parameter estimate and the variance of the distribution of errors.  However, by doing so, we are no longer talking about a confidence interval, as we are talking about a random variable and not a parameter; we call this a prediction interval.  Essentially, the variance used in calculating the “margin of error” would be the sum of the variance associated with the parameter estimate and the variance associated with the residuals.  This is captured by adding a one under the radical: $$\hat\sigma \sqrt{1 + \frac{1}{n} + \frac{(x_0-\bar{x})^2}{(n-1)s_x^2}}$$ Conceptually, the 1 “adds” the estimate for $\sigma$ to the standard error.

Answer 2

I think the truth is that I thought incorrectly.

The following simulation shows that under normal error assumption (normal X actually) the proportion of examples falling out of a 95% PI on a population is not different from 5%.

If I make something stupid, please correct / post an answer.

start_time <- Sys.time()

popul <- 1000000
sam_size <- 25
sim_num <- 10000

big_x <- rnorm(popul, 0, 1)

outofpi <- numeric()

for(i in 1:sim_num)
{

     sam_ind <- sample(popul, sam_size, replace = F)

     x_dat <- data.frame(x = big_x[sam_ind])

     lm_model <- lm(data = x_dat, x ~ 1)

     lm_preds <- predict.lm(lm_model
                            , x_dat
                            , interval = "prediction"
                            , level = 0.95
                            )[1, ]

     newd <- big_x[-sam_ind]

     outofpi[i] <- length(newd[newd < lm_preds[2] | newd > lm_preds[3]]) / length(newd)

}

Sys.time() - start_time

hist(outofpi, breaks = 'fd')

t.test(x = outofpi
       , alternative = "greater"
       , mu = 0.05
       , conf.level = 0.95
       )


    One Sample t-test

data:  outofpi
t = 0.54362, df = 9999, p-value = 0.2934
alternative hypothesis: true mean is greater than 0.05
95 percent confidence interval:
 0.04961823        Inf
sample estimates:
 mean of x 
0.05018844