Is it possible to add standard error or confidence interval to a plot of a predicted vs observed values derived from a multiple regression model? I believe that I have seen such plots as an output in Statistica, but am unsure how to create them in R.
I believe I have a solution (below), but am unsure that I have done this correctly. Basically, I have created a new
dataframe with predictor variable in the range of their possible values. My worry with such an approach is that the prediction is based on the rows of data, and does not really address situations where the variables are randomly selected.
Many thanks for your help.
set.seed(1) n <- 200 x1 <- rnorm(n, mean=10, sd=3) x2 <- rnorm(n, mean=20, sd=5) e <- rnorm(n, mean=10, sd=3) y <- 5 + 2*x1 + 0.5*x2 + e fit <- lm(y ~ x1 + x2) summary(fit) #plot of predicted vs observed pred1 <- predict(fit, se.fit=TRUE) plot(pred1$fit ~ y) abline(0,1, col=8, lwd=2) #new dataframe sequence of each predictor variable in their range df.new <- data.frame(x1=seq(min(x1), max(x1),,100), x2=seq(min(x2), max(x2),,100)) pred2 <- predict(fit, df.new, se.fit=TRUE) #plot of predicted vs observed w/ standard error interval? png("pred_vs_obs.png", width=6, height=6, units="in", res=200) plot(pred1$fit ~ y) abline(0,1, col=8, lwd=2) lines(pred2$fit+1.96*pred2$se.fit ~ pred2$fit, col=2, lty=2, lwd=2) lines(pred2$fit-1.96*pred2$se.fit ~ pred2$fit, col=2, lty=2, lwd=2) dev.off()
The following code elaborates on my hesitation with the method that I used. The relationship between Standard Error (SE) and y is not a precise; i.e. various values of y that are relatively close together, have widely differing SE (black symbols in figure below,
pred1), while the above method predicts a single SE for each predicted y (red symbols,
pred2). Furthermore, using several different combinations of x1 and x2 that always result in the same y-value, I get a single (but different!) SE (green symbol,
pred3). What is going on here? Is there a more correct way of doing this with some sort of permutation method?
#? Do different solutions to a given predicted value always give the same standard error? y.tmp <- rep(40,20) x1.tmp <- seq(0,10, length(y.tmp)) x2.tmp <- (y.tmp - fit$coeff - fit$coeff*x1.tmp) / fit$coeff df3 <- data.frame(x1=x1.tmp, x2=x2.tmp) pred3 <- predict(fit, df3, se.fit=TRUE) YLIM <- range(pred1$se.fit, pred2$se.fit, pred3$se.fit) png("fit.se_vs_fit.png", width=6, height=6, units="in", res=200) plot(pred1$se.fit ~ pred1$fit, ylim=YLIM, lwd=2) points(pred2$se.fit ~ pred2$fit, col=2, lwd=2) points(pred3$se.fit ~ pred3$fit, col=3, lwd=2) legend("topright", legend=c("orig. data", "range of x1 & x2", "various comb. of x1 & x2 \nto acheive y=40"), col=1:3, pch=1, lwd=2, lty=0) dev.off()