Perhaps this will be quick for some of you but I would like to represent the standard error of the regression using plot().
So, if I have a data like this:
x1 <- 1:500
b0 <- 17
b1 <- 0.5
sigma <- 7
er <- rnorm(x1,0,sigma)
y <- b0 + b1*x1 + er
model1 <- lm(y~x1)
plot(x1,y)
abline(model1,col="red",lwd=5)
How could I represent the standard error, in lines, for that regression?
Thanks in advance!
You may want to use ggplot() instead of plot(). You have a slightly more modern layout and much more possibilities.
Add the following code after your example code. For better visibility I have changed the original data a bit (See line: y <- b0 + b1*x1 + er*5.)
# ------ test case ggplot ----------------
library (ggplot2)
# --- change some data fields ------------
y <- b0 + b1*x1 + er*5
df <- data.frame(y,x1)
ggplot(data = df, aes(x1,y)) +
geom_point() +
geom_smooth(method="lm", color="red", fill = "blue")
You can do it manually calculating the regression without the error and from that put the general lines of the sd() of all the residual errors:
x1 <- 1:500
b0 <- 17
b1 <- 0.5
sigma <- 7
er <- rnorm(x1,0,sigma)
y <- b0 + b1*x1 + er
GeneralStandardDev<-sd(model1$residuals)
UpperLine<- model1$coefficients[1]+model1$coefficients[2]*x1 + GeneralStandardDev
LowerLine<- model1$coefficients[1]+model1$coefficients[2]*x1 - GeneralStandardDev
model1 <- lm(y~x1)
plot(x1,y)
abline(model1,col="red",lwd=5)
lines(x1, UpperLine, col = "blue")
lines(x1, LowerLine, col = "blue")
predict(model1, se.fit = TRUE) if you want to calculate the standard error of predictions.
lines(x1[order(x1)], predict(model1, interval = "confidence")[, "lwr"][order(x1)]); lines(x1[order(x1)], predict(model1, interval = "confidence")[, "upr"][order(x1)])