7
$\begingroup$

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
 [1] 0.6011667 0.4976294 0.7252444 0.4602669 0.7752706 0.5178496 0.7267482 1.0000172 0.9793969
[10] 0.5108741 0.5108741 0.6544576 0.6819424 0.6718159 1.1525645 1.2633704 1.2125441 0.7270859
[19] 0.8820281 0.7988791 0.7380797 0.7442464 0.7773252 0.6623197 0.6616629 0.7700721 0.7322422
[28] 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?

$\endgroup$
3
$\begingroup$

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[2] 

# 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[2], 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.

$\endgroup$
1
  • 7
    $\begingroup$ What "excellent post"? $\endgroup$ Jan 4 at 23:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.