4.4 Prediction

The standard error in the expected value of \(\hat{y}\) at some new set of predictors \(X_n\) is

\[SE(\mu_\hat{y}) = \sqrt{\hat{\sigma}^2 (X_n (X'X)^{-1} X_n')}.\]

The standard error increases the further \(X_n\) is from \(\bar{X}\). If \(X_n = \bar{X}\), the equation reduces to \(SE(\mu_\hat{y}) = \sigma / \sqrt{n}\). If \(n\) is large, or the predictor values are spread out, \(SE(\mu_\hat{y})\) will be relatively small. The \((1 - \alpha)\%\) confidence interval is \(\hat{y} \pm t_{\alpha / 2} SE(\mu_\hat{y})\).

The standard error in the predicted value of \(\hat{y}\) at some \(X_{new}\) is

\[SE(\hat{y}) = SE(\mu_\hat{y})^2 + \sqrt{\hat{\sigma}^2}.\]

Notice the standard error for a predicted value is always greater than the standard error of the expected value. The \((1 - \alpha)\%\) prediction interval is \(\hat{y} \pm t_{\alpha / 2} SE(\hat{y})\).

4.4.0.1 Example

What is the expected value of mpg if the predictor values equal their mean values?

R performs this calucation with the predict() function with parameter interval = confidence.

m <-lm(mpg ~ ., data = d[,1:9])
X_new <- data.frame(Const = 1,
                    cyl = factor(round(mean(as.numeric(as.character(d$cyl))),0), levels = levels(d$cyl)), 
                    disp = mean(d$disp),
                    hp = mean(d$hp),
                    drat = mean(d$drat),
                    wt = mean(d$wt),
                    qsec = mean(d$qsec),
                    vs = factor("S", levels = levels(d$vs)), 
                    am = factor("manual", levels = levels(d$am)))
predict.lm(object = m, 
           newdata = X_new, 
           interval = "confidence")

##   fit lwr upr
## 1  21  17  25

You can verify this by manually calculating \(SE(\mu_\hat{y}) = \sqrt{\hat{\sigma}^2 (X_{new} (X'X)^{-1} X_{new}')}\) using matrix algebra.

X2 <- lapply(data.frame(model.matrix(m)), mean) %>% unlist() %>% t()
X2[2] <- contr.poly(3)[2,1]  # cyl linear
X2[3] <- contr.poly(3)[2,2]  # cyl quadratic
X2[9] <- 1
X2[10] <- 1

y_exp <- sum(m$coefficients * as.numeric(X2))
se_y_exp <- as.numeric(sqrt(rse^2 * 
                              X2 %*% 
                              solve(t(X) %*% X) %*% 
                              t(X2)))

t_crit <- qt(p = .05 / 2, df = n - k - 1, lower.tail = FALSE)
me <- t_crit * se_y_exp
cat("fit: ", round(y_exp, 6),
    ", 95% CI: (", round(y_exp - me, 6), ", ", round(y_exp + me, 6), ")")

## fit:  21 , 95% CI: ( 17 ,  25 )

4.4.0.2 Example

What is the predicted value of mpg if the predictor values equal their mean values?

R performs this calucation with the predict() with parameter interval = prediction.

predict.lm(object = m, 
           newdata = X_new, 
           interval = "prediction")

##   fit lwr upr
## 1  21  15  28

se_y_hat <- sqrt(rse^2 + se_y_exp^2)
me <- t_crit * se_y_hat
cat("fit: ", round(y_exp, 6),
    ", 95% CI: (", round(y_exp - me, 6), ", ", round(y_exp + me, 6), ")")

## fit:  21 , 95% CI: ( 15 ,  28 )