How can I reproduce R's prediction band calculation?

Question

Suppose I fit a linear regression model in R like so

> head(mtcars[, c("mpg", "wt")])
                   mpg    wt
Mazda RX4         21.0 2.620
Mazda RX4 Wag     21.0 2.875
Datsun 710        22.8 2.320
Hornet 4 Drive    21.4 3.215
Hornet Sportabout 18.7 3.440
Valiant           18.1 3.460

# Fit linear regression model
> lrmodel <- lm(mpg ~ wt, data = mtcars)

> print(lrmodel)

Call:
lm(formula = mpg ~ wt, data = mtcars)

Coefficients:
(Intercept)           wt  
     37.285       -5.344

I can make predictions on new data (or in this case, the same data) with prediction bands.

# Make predictions with lower and upper prediction bands @ 0.99 level
> preds <- predict(lrmodel, mtcars, interval = "prediction", level = 0.99)
> head(preds)
                       fit      lwr      upr
Mazda RX4         23.28261 14.72715 31.83807
Mazda RX4 Wag     21.91977 13.39748 30.44207
Datsun 710        24.88595 16.26877 33.50313
Hornet 4 Drive    20.10265 11.59662 28.60868
Hornet Sportabout 18.90014 10.38722 27.41307
Valiant           18.79325 10.27904 27.30747

Given my fitted model, how can I reproduce these prediction bands manually?

---

What I've tried

According to this article prediction bands are given by

$$ \hat{y}_h \pm t_{(1-\alpha/2, n-2)} \times \sqrt{MSE \times \left(1+ \frac{1}{n} + \dfrac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)} $$

where

• $x_h$ is the x location of the point of interest
• $\hat{y}_h$ is the predicted value at the point of interest
• $t_{(1-\alpha/2, n-2)}$ is the t-multiplier with n-2 degrees of freedom

Plugging and chugging..

# Lower band calculation
xh <- lrmodel$model$wt[1]          # 2.62
n <- length(lrmodel$residuals)     # 32
tval <- dt(x = 0.99/2, df = n - 2) # 0.3487666
xbar <- mean(lrmodel$model$wt)     # 3.21725
mse <- mean(lrmodel$residuals^2)   # 8.697561
yhat <- preds[1, "fit"]            # 23.28261
yhat - tval * sqrt(mse * (1 + 1/n + (xh - xbar)^2/sum((lrmodel$model$wt - xbar)^2)))
# 22.23202

My result, 22.23202, is quite different than R's result for the same $x$ value, 14.72715. Where did I go wrong?

Answer 1

Thanks to @dipetkov's hints in the comments, I was able to correct my mistakes, and reproduce the calculation as follows:

# Lower band calculation
xh <- lrmodel$model$wt[1]               # 2.62
n <- length(lrmodel$residuals)          # 32
tval <- qt(p = 0.01/2, df = n - 2)      # -2.749996
xbar <- mean(lrmodel$model$wt)          # 3.21725
mse <- sum(lrmodel$residuals^2)/(n - 2) # 9.277398
yhat <- preds[1, "fit"]                 # 23.28261

# upper band (31.83807)
yhat - tval * sqrt(mse * (1 + 1/n + (xh - xbar)^2/sum((lrmodel$model$wt - xbar)^2)))

# lower band (14.72715)
yhat + tval * sqrt(mse * (1 + 1/n + (xh - xbar)^2/sum((lrmodel$model$wt - xbar)^2)))