How to plot a 95% confidence interval for lm(y~x1+x2)

Question

I followed the post by Max but doesn't seem to understand how to apply the solution given by Alejandro to a different regression equation with form in lm(y~x1+x2).

Below are data for you to reproduce and what I did.

y=c(139.31449, 105.17776, 105.38411,  99.27608,  92.29064,  91.55114,  84.44251,  78.40453,  74.66656,  73.33242,  72.42429,  77.08666)

x1=c(0.04, 0.00, 0.00, 0.00, 0.00, 0.00, 0.04, 0.00, 0.00, 0.00, 0.00, 0.00)

x2=c(0.00, 0.08, 0.10, 0.12, 0.15, 0.20, 0.00, 0.08, 0.10, 0.12, 0.15, 0.20)

lm1 <- lm(y ~ x1+x2)

newx = seq(min(x1+x2),max(x1+x2),by = 0.05)

conf_interval <- predict(lm1, newdata=data.frame(x=newx), interval="confidence",
                         level = 0.95)

plot(x1+x2, y, xlab="x", ylab="y")

abline(lm1, col="lightblue")

matlines(newx, conf_interval[,2:3], col = "blue", lty=2)

I'm not sure if I'm doing the right thing for newx. I had issue with adding regression line or I shouldn't use abline? R not running conf_interval properly as well.

I tried looking for similar topic to mine but struggle to find one. Can someone help please? Thanks.

Answer 1

Your model y~x1+x2 is not a simple linear regression(SLR), so it's confidence interval(CI) cannot be visualized like SLR.

There are several way to plot CI of this model.

First, using predict3d::ggPredict(), for a fixed x2,

ggPredict(lm1, digits = 1, se = TRUE)

Second, by usling plotly::plot_ly and some more to just plot 3-Dimensional confidence plane(?).

xgrid <- seq(0,0.04 , length.out = 30)
ygrid <- seq(0, 0.15, length.out = 30)
newdat <- expand.grid(xgrid, ygrid)
colnames(newdat) <- c("x1", "x2")
predicted <- predict(lm1, newdat, se = TRUE)

ymin <- predicted$fit - 1.96 * predicted$se.fit
ymax <- predicted$fit + 1.96 * predicted$se.fit
fitt <- predicted$fit
z <- matrix(fitt, length(xgrid))
ci.low <- matrix(ymin, length(xgrid))
ci.up <- matrix(ymax, length(xgrid))
library(plotly)
plot_ly(x = xgrid, y = ygrid) %>%
  add_surface(z = z, colorscale = list(c(0,1), c("red", "blue"))) %>%
  add_surface(z = ci.low, opacity = 0.5, showscale = FALSE, colorscale = list(c(0,1),c("grey","grey"))) %>% 
  add_surface(z = ci.up, opacity = 0.5, showscale = FALSE, colorscale = list(c(0,1),c("grey","grey")))

note that x and y are x1 and x2 of your data and z is predicted y.