CODEWITHSUNDEEP
rplotregressionlinear-regressionlm
FKG
FKG

Reputation: 285

Bad plot when plotting 95% confidence interval of linear model prediction

I need a little help with plotting prediction with its confidence interval. Consider the following example

library(Hmisc)
data("mtcars") 

mfit = lm(mpg ~ vs + disp + cyl, data = mtcars)

#disp and cyl at their mean
newcar = data.frame(vs = c(0,1), disp = 230, cyl = 6.188)

pmodel <- predict(mfit, newcar, se.fit=TRUE)

I want to plot the effect of vs (when 0 and 1) when all other variables are held constant (mean / mode).

To do this I run this code below:

plot(1:2, pmodel$fit[1:2], ylim=c(0,1), pch=19, xlim=c(.5,2.5), xlab="X",
     ylab = "Predicted values", xaxt = "n", main = "Figure1")
arrows(1:2, (pmodel$fit[1:2] - 1.96 * pmodel$fit[1:2]),
       1:2, (pmodel$fit[1,1] + 1.96 * pmodel$fit[1:2]),
       length=0.05, angle=90, code=3)
axis(1, at=c(1,2), labels=c("0","1"))

What am I doing wrong here? Thanks!

Upvotes: 1

Views: 395

Answers (2)

Zheyuan Li
Zheyuan Li

Reputation: 73265

Note you have ylim = c(0, 1), which is incorrect. When drawing confidence interval, we have to make sure ylim covers the lower and upper bound of CI.

## lower and upper bound of CI
lower <- with(pmodel, fit - 1.96 * se.fit)
upper <- with(pmodel, fit + 1.96 * se.fit)

## x-location to plot
xx <- 0:1

## set `xlim` and `ylim`
xlim <- range(xx) + c(-0.5, 0.5)  ## extends an addition 0.5 on both sides
ylim <- range(c(lower, upper))

## produce figure
plot(xx, pmodel$fit, pch = 19, xlim = xlim, ylim = ylim, xaxt = "n",
     xlab = "X", ylab = "Predicted values", main = "Figure1")
arrows(xx, lower, xx, upper, length = 0.05, angle = 90, code = 3)
axis(1, at = xx)

enter image description here

Several other comments on your code:

  • it is fit - 1.96 * se.fit not fit - 1.96 * fit;
  • you can plot directly on x-location 0:1, rather than 1:2;
  • it is fit[1:2] not fit[1,1].

Upvotes: 1

Sandipan Dey
Sandipan Dey

Reputation: 23101

In ggplot:

df <- data.frame(x=1:2, pred=pmodel$fit[1:2])
df$lower <- df$pred - 1.96 * pmodel$se.fit[1:2] 
df$upper <- df$pred + 1.96 * pmodel$se.fit[1:2]
ggplot(df, aes(x, pred, group=x, col=as.factor(x))) + geom_point(size=2) + 
  geom_errorbar(aes(ymin=lower, ymax=upper), width=0.1)

enter image description here

Upvotes: 1

Related Questions