# What are the formulae used in R by predict.lm when interval= a) 'none', b) 'prediction', and c) 'confidence'?

The references provided in the R documentation for predict.lm, taken together, actually leave open a number of possibilities for the formulae for confidence and prediction intervals (including the possible use of tolerance intervals as prediction intervals, though they are not the same thing). My unfamiliarity with the innards of R prevents me from extracting the formulae from the code for predict. Can anyone please provide the formulae actually used?

• +1 for trying to read the source code before asking. Although I think you ought to know how the none case works... – shadowtalker May 27 '15 at 12:21

Prediction "none" is clearly no prediction. To understand the other options, the relevant part of predict.lm is this:

tfrac <- qt((1 - level)/2, df)
hwid <- tfrac * switch(interval, confidence = sqrt(ip),
prediction = sqrt(ip + pred.var))


Where ip is the square of the fitted standard error and pred.var is the variance of the residuals (under the default options).

Consider the following example:

mdl <- lm(hp~disp,mtcars)
predict(mdl,newdata=list(disp=300),se.fit=TRUE,interval="confidence")
$fit fit lwr upr 1 177.0003 159.2944 194.7062$se.fit
[1] 8.669712

$df [1] 30$residual.scale
[1] 42.6459


One can recreate the confidence interval (approximately due to rounding of parameters for brevity) with:

177.003 + qt(0.025,30)*c(1,-1)*8.669712
[1] 159.2971 194.7089


Similarly the prediction interval can be computed with:

177.003 + qt(0.025,30)*c(1,-1)*sqrt(8.669712^2+var(mdl$residuals)) [1] 89.51433 264.49167 predict(mdl,newdata=list(disp=300),interval="prediction") fit lwr upr 1 177.0003 88.12423 265.8764  • Thanks, though I'm disappointed (especially for the confidence intervals) that these are single-point intervals. Draper & Smith and Miller use F instead of t to provide a confidence band for the whole line. Is there an R package with a function that does that? – icc May 28 '15 at 9:48 • Yes, at least for confidence intervals: investr – icc Jun 1 '15 at 7:16 • The (slight) discrepancy between your prediction interval and the one computed by predict.lm arises as you're computing the prediction variance as var(mdl$residuals). If you use sum(mdl$residuals ^ 2) / mdl$df.residual) instead, you should get the exact same answer as R. – Kay Brodersen Jul 31 '15 at 11:38

I can show you some steps how the predict() function calculates for 3 different intervals. I will be taking mtcars as an example and predicting intervals (none, confidence and prediction) at wt = 3 (3000lb)

library(ggplot2)
library(UsingR)
data(mtcars)
y<-mtcars$mpg x<-mtcars$wt

fit<-lm(y~x)
none<-predict(fit,newdata = data.frame(x=3),interval = "none")
con<-predict(fit,newdata = data.frame(x=3),interval = "confidence")
pred<-predict(fit,newdata = data.frame(x=3),interval = "prediction")

g<-ggplot(data=mtcars,aes(x=wt,y=mpg))
g<-g+geom_point(size=2, color = "black", alpha = 0.4)
g<-g+geom_smooth(method = "lm")
g<-g+ggtitle("Prediction vs Confidence using predict()")+theme(plot.title = element_text(hjust = 0.5))
g<-g+labs(x ="wt", y = "mpg")
g<-g+geom_point(aes(x=3,y=pred[3]),colour = "red") #95% prediction level
g<-g+geom_point(aes(x=3,y=con[3]),colour = "blue") #95% confidence level
g

#From the predict() function
none
con
pred

ssx<-round(sum((x-mean(x))^2),2)

#From the manual calculations using formula
con[1] + c(-1,1) * qt(.975,df = fit$df) * sigma(fit) * sqrt(1/n + (3- mean(x))^2/ssx) #confidence interval pred[1] + c(-1,1) * qt(.975,df = fit$df) * sigma(fit) * sqrt(1 + 1/n + (3-mean(x))^2/ssx) #prediction interval


In calculations, which are the last 2 steps in the code above, the only difference is in the standard error calculations (for none we dont output any intervals so understood). The small variations which I am getting may be the case due to rounding issues. So the formula used are pretty much from the theory which are:

standard error used for Confidence interval:

standard error used for Prediction interval: