# How to add standard error or confidence interval to a plot of a predicted vs observed values?

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.

Example:

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() Edit:

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() • It is strange to see this done with a plot of predicted vs. fit: it makes more sense to see the intervals in a plot of predicted vs. explanatory variables. The reason is that (except in the simplest case of a straight line fit to one explanatory variable) the SE does not depend on the predicted value: it depends on the values of the explanatory variables. Thus the plot isn't even well-defined as a curve; it should be thought of as (at best) a collection of points, each with its own little vertical error band. – whuber Dec 17 '12 at 16:38

1. You want pred2\$fit +1.96*pred2\$se.fit to have 95% confidence bands.
3. I am not sure why you have created an equally spaced set of X's for your confidence limits. Surely you could have run ${\tt lines()}$ on the original values. However, what you did is correct.
• Thank you for your comment. I have added an edit to the question in order to illustrate my hesitation with this approach - SE does not seem to be a constant value for a given y. Any additional clarification would be greatly appreciated. – Marc in the box Dec 17 '12 at 13:30
• I guess that's the point - does it really make sense to display a confidence interval on a predicted y versus observed y plot (i.e. without any explanatory variable, as @whuber mentions). Perhaps there is another way of doing this given the SE associated with the coefficients (summary(fit)\$coeff[,2]) – Marc in the box Dec 18 '12 at 8:43