3
$\begingroup$

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?

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

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
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – icc May 28 '15 at 9:48
  • $\begingroup$ Yes, at least for confidence intervals: investr $\endgroup$ – icc Jun 1 '15 at 7:16
  • 2
    $\begingroup$ 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. $\endgroup$ – Kay Brodersen Jul 31 '15 at 11:38
0
$\begingroup$

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

enter image description here

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:

enter image description here

standard error used for Prediction interval:

enter image description here

$\endgroup$

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.