# What are the standard errors of the predictions from predict.lm in R?

In R, ?predict says:

If the logical se.fit is TRUE, standard errors of the predictions are calculated.

An example:

> predict(lm(mpg ~ wt + cyl, data = mtcars), se.fit=TRUE)$se.fit  0.6011667 0.4976294 0.7252444 0.4602669 0.7752706 0.5178496 0.7267482 1.0000172 0.9793969  0.5108741 0.5108741 0.6544576 0.6819424 0.6718159 1.1525645 1.2633704 1.2125441 0.7270859  0.8820281 0.7988791 0.7380797 0.7442464 0.7773252 0.6623197 0.6616629 0.7700721 0.7322422  0.9282190 0.9023791 0.5342815 0.7267482 0.8176667  How are those standard errors defined? How are they calculated? I looked at the predict.lm code, and it has a lot of branches. An excerpt of the code is something like: ip[, i] <- if (any(iipiv > 0L)) as.matrix(X[, iipiv, drop = FALSE] %*% Rinv[ii, , drop = FALSE])^2 %*% rep.int(res.var, p) # ... se later returned as "se.fit" se <- sqrt(ip)  Wikipedia defines the standard error of the betas (slope of a 1-dimensional predictor), but not the standard error of the predictions. How is se.fit defined, using some standard notation? Related (or the same?), how is it computed? ## 1 Answer As suggested by its name, se.fit returns the standard error of the fit. This is the standard error associated with the estimated mean value of the response variable at given values of the predictor variables included in a linear regression model fitted with the lm() function in R. Here is some R code that verifies that for a simplified model with just wt as a predictor of mpg: # model m = lm(mpg ~ wt, data = mtcars) # model summary summary(m) # residual degrees of freedom summary(m)$df

# 95% confidence interval for the mean value of mpg for ALL cars
# represented by the ones in the mtcars dataset which have
# wt = 3.2
p <- predict(m, newdata = data.frame(wt = 3.2), se.fit=TRUE,
interval="confidence")

p

# reported interval (lwr, upr) SHOULD be constructed as
# fit +/- t(alpha/2,n-2)*se.fit
# where t(alpha/2,n-2) is a critical value from the
# t distribution with n-2 degrees of freedom and
# alpha = 0.05

# compute the critical value
crit <- qt(p= 0.05/2, df = summary(m)\$df, lower.tail = FALSE)
crit

# check that computing the half-length of the interval and
# dividing it by the critical value gives the same result
# as that reported by se.fit
(p$$fit[,"upr"] - p$$fit[,"fit"])/crit


This excellent post illustrates the difference between the *standard error of fit" and the "standard error of prediction" and gives the formulas used to derive the two quantities.

• What "excellent post"? Jan 4 at 23:06