11

How can I calculate and plot a confidence interval for my regression in r? So far I have two numerical vectors of equal length (x,y) and a regression object(lm.out). I have made a scatterplot of y given x and added the regression line to this plot. I am looking for a way to add a 95% prediction confidence band for lm.out to the plot. I've tried using the predict function, but I don't even know where to start with that :/. Here is my code at the moment:

x=c(1,2,3,4,5,6,7,8,9,0)
y=c(13,28,43,35,96,84,101,110,108,13)

lm.out <- lm(y ~ x)

plot(x,y)

regression.data = summary(lm.out) #save regression summary as variable
names(regression.data) #get names so we can index this data
a= regression.data$coefficients["(Intercept)","Estimate"] #grab values
b= regression.data$coefficients["x","Estimate"]
abline(a,b) #add the regression line

Thank you!

Edit: I've taken a look at the proposed duplicate and can't quite get to the bottom of it.

Max Lester
  • 133
  • 1
  • 1
  • 8

2 Answers2

18

You have yo use predict for a new vector of data, here newx. 

x=c(1,2,3,4,5,6,7,8,9,0)

y=c(13,28,43,35,96,84,101,110,108,13)

lm.out <- lm(y ~ x)
newx = seq(min(x),max(x),by = 0.05)
conf_interval <- predict(lm.out, newdata=data.frame(x=newx), interval="confidence",
                         level = 0.95)
plot(x, y, xlab="x", ylab="y", main="Regression")
abline(lm.out, col="lightblue")
lines(newx, conf_interval[,2], col="blue", lty=2)
lines(newx, conf_interval[,3], col="blue", lty=2)

EDIT

as it is mention in the coments by Ben this can be done with matlines as follow:

plot(x, y, xlab="x", ylab="y", main="Regression")
abline(lm.out, col="lightblue")
matlines(newx, conf_interval[,2:3], col = "blue", lty=2)
Alejandro Andrade
  • 2,196
  • 21
  • 40
12

I'm going to add a tip that would have saved me a lot of frustration when trying the method given by @Alejandro Andrade: If your data are in a data frame, then when you build your model with lm(), use the data= argument rather than $ notation. E.g., use 

lm.out <- lm(y ~ x, data = mydata)

rather than 

lm.out <- lm(mydata$y ~ mydata$x)

If you do the latter, then this statement

predict(lm.out, newdata=data.frame(x=newx), interval="confidence", level = 0.95)

seems to either ignore the new values passed using newdata= or there's a silent error. Either way, the output is the predictions from the original data, not the new data.

Also, be sure your x variable gets the same name in the new data frame that it had in the original. That's easier to figure out because you do get an error, but knowing it ahead of time might save you a round of debugging.

Note: Tried to add this as a comment, but don't have enough reputation points.

acullum
  • 121
  • 1
  • 5