5
$\begingroup$

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.

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

enter image description here

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[1] - fit$coeff[2]*x1.tmp) / fit$coeff[3]

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

enter image description here

$\endgroup$
  • 5
    $\begingroup$ 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. $\endgroup$ – whuber Dec 17 '12 at 16:38
2
$\begingroup$

That looks fairly good. A few points:

  1. You want pred2\$fit +1.96*pred2\$se.fit to have 95% confidence bands.
  2. I see that you have drawn confidence lines for the expected mean, not predicted observations. Is that what you wanted?
  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.
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Marc in the box Dec 17 '12 at 13:30
  • $\begingroup$ The standard error of the predicted mean varies with X. The further you are from the mean of the X's, the larger the SE. Your second graph makes a lot of sense. $\endgroup$ – Placidia Dec 17 '12 at 16:11
  • 1
    $\begingroup$ I'm not sure that I understand all of your points. What do you mean by "different combinations of X and Y" etc? The SE depends on X1 and X2, as well as the Y values of the training sample. I guess if the same Y occurs at different values of X1 and X2, you would get different SE's. The SE is not a function of Y alone. $\endgroup$ – Placidia Dec 17 '12 at 16:14
  • $\begingroup$ 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]) $\endgroup$ – Marc in the box Dec 18 '12 at 8:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.