How to represent the interval or uncertainty on regression predictions in an 'experimental vs predicted' plot?

Question

Using an example similar to the one from R 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...?

Answer 1

You are close to the usual approach in your second last plot. I would say that the convention is to place the observed (i.e. experimental) values on the X-axis and the fitted values on the Y-axis, assuming that are simply fitting a regression to experimental data. I think the term 'fitted' here is more clear because there is no implication that the values were predicted before the experiment was done. Then you can use vertical error bars to show the uncertainty in your fits. Note that it's also possible to have uncertainty in your observed values, in which case you can use both vertical and horizontal error bars (and you'd want to account for the uncertainty in X-values in the regression as well).