Predict function for lm object in R

Question

Why are prediction_me and prediction_R not equal? I'm attempting to follow the formula given by Lemma 5 here. Does the predict function use a different formula, have I made a mistake in my computation somewhere, or is it just rounding error? (the two are pretty close)

set.seed(100)
# genrate data
x    <- rnorm(100, 10)
y    <- 3 + x + rnorm(100, 5)
data <- data.frame(x = x, y = y)
# fit model
mod  <- lm(y ~ x, data = data)

# new observation
data2 <- data.frame(x = rnorm(5, 10))

# prediction for new observation
d    <- as.matrix(cbind(1, data[,-2]))
d2   <- as.matrix(cbind(1, data2))
fit  <- d2 %*% mod$coefficients 
t    <- qt(1 - .025, mod$df.residual)
s    <- summary(mod)$sigma
half <- as.vector(t*s*sqrt(1 + d2%*%solve(t(d)%*%d, t(d2))))

prediction_me <- cbind(fit, fit - half, fit + half)

prediction_R <- predict(mod, newdata = data2, interval = 'prediction')


prediction_me
prediction_R

Perhaps you may want to look into [this](https://stackoverflow.com/a/48706431/6574038) how to predict manually. — jay.sf, Feb 21 '18 at 13:07
`half <- as.vector(t*s*sqrt((1 + t(d2[, -1])%*%solve(t(d)%*%d, d2[, -1]))))` returns an error. — nghauran, Feb 21 '18 at 13:31

score 2 · Accepted Answer · answered Feb 21 '18 at 14:19

Your current code is almost fine. Just note that the formula in Lemma 5 is for a single newly observed x. For this reason, half contains not only relevant variances but also covariances, while you only need the former ones. Thus, as.vector should be replaced with diag:

half <- diag(t * s * sqrt(1 + d2 %*% solve(t(d) %*%d , t(d2))))
prediction_me <- cbind(fit, fit - half, fit + half)
prediction_R <- predict(mod, newdata = data2, interval = 'prediction')

range(prediction_me - prediction_R)
# [1] 0 0

Predict function for lm object in R

1 Answers1