Using an example similar to the one from
predict, simulate some independent variable ($x$) data, map them to an observed dependent variable ($y$) using a linear function for the true value (I chose here $3 x - 1$), plus a random normal error. Then fit a linear model to the data.
set.seed(8298) x <- runif(15, -3, 3) x <- sort(x) y <- 3 * x - 1 + rnorm(15) lm1 <- lm(y ~ x)
Let $\hat y = f(x)$ be the linear model.
As explained by the
R documentation on
predict, $\hat y$ and its lower and upper limits can be calculated as:
new <- data.frame(x = seq(-3, 3, 0.5)) pred.w.plim <- predict(lm1, new, interval = "prediction") pred.w.clim <- predict(lm1, new, interval = "confidence")
Then, if I want to represent the uncertainty on the predictions, I can plot $\hat y$ vs $x$, with lines that connect the upper limits of the intervals, and the same for the lower limits:
matplot(new$x, cbind(pred.w.clim, pred.w.plim[,-1]), lty = c(1,2,2,3,3), col = c(1,2,2,3,3), type = "l", lwd = 2, xlab = "x", ylab = "predicted y") points(x, y, pch = 16)
Note that I added the experimental $y$ as points, to show how close the prediction (black line) is to the observed value (black point).
Now, imagine instead a case where $x$ is not a simple real variable, but a complicated vector of several descriptors, and the true relationship between $y$ and $x$ is no longer linear.
You can still make a model that outputs, for each object with a specific $x$, a prediction $\hat y$ with its upper and lower limits (of the 'prediction' type, in analogy with the linear case).
However, now you can no longer plot $\hat y$ vs $x$.
In such cases, we are told it is possible to plot $y$ (observed, experimental) vs $\hat y$:
y_pred <- predict(lm1) plot(y ~ y_pred, pch = 16, xlab = "predicted y", ylab = "experimental y") abline(0, 1)
In this case, to make an example, I will still calculate the prediction limits using the linear model, and with the 'confidence' method to show some variation between them, but you need to imagine that these limits come from some more complex calculation and they are actual prediction intervals.
old <- data.frame("x" = x) pred.old.clim <- predict(lm1, old, interval = "confidence")
And this is where I am stuck, hence my post.
How do I correctly represent the uncertainty or interval on the predictions in such a plot?
Initially I hastily thought that I just needed to put vertical error bars on each point:
ints <- pred.old.clim[,"upr"] - pred.old.clim[,"lwr"] segments(x0 = y_pred, y0 = y - ints / 2, y1 = y + ints / 2)
But this cannot be right, can it, as it uses an experimental value as the center for a distribution that should instead be centered on the predicted value.
The only two alternatives that seemed to make some sense were plotting the bars horizontally from each experimental point:
plot(y ~ y_pred, pch = 16, xlab = "predicted y", ylab = "experimental y") abline(0, 1) segments(y0 = y, x0 = pred.old.clim[,"lwr"], x1 = pred.old.clim[,"upr"])
or vertically, but from the unit line, i.e. from the point where prediction and experiment would agree:
plot(y ~ y_pred, pch = 16, xlab = "predicted y", ylab = "experimental y") abline(0, 1) segments(x0 = y_pred, y0 = pred.old.clim[,"lwr"], y1 = pred.old.clim[,"upr"]) matlines(x = pred.spr.clim[,1], y = pred.spr.clim[,-1], lty = 2, col = 2) matlines(x = pred.spr.plim[,1], y = pred.spr.plim[,-1], lty = 3, col = 3)
What do you think is the correct way of doing this? Any of the above, something else...?