ppc_error_scatter_avg() should be able to plot residuals as function of predicted y, not y

[kruschke]

(Notation below: $y$ is data value, $y_{pred}$ is predicted value, with $y_{pred}$ computed as $stat(y_{rep})$ from posterior draws.)

The usual residual analysis plots $y - y_{pred}$ on the vertical axis with $y_{pred}$ on the horizontal axis.

But ppc_error_scatter_avg() plots $y - y_{pred}$ on the vertical axis with $y$ on the horizontal axis. This is confusing and difficult to interpret (for me, but I'm not alone: https://stats.stackexchange.com/a/146002).

I would be great if ppc_error_scatter_avg() plotted $y_{pred}$ on the horizontal axis, either by default or with that as an option. (Or maybe there's already an easy way to do that; sorry if I missed it.)

Thanks for considering!

P.S. My comment here assumes that $y_{pred}$ is computed as $stat(y_{rep})$ from posterior draws. But the related thread regarding residuals #349 seems to suggest that residuals are computed as $stat(y - y_{rep})$, with $stat(y_{rep})$ not separately, explicitly computed. Hmmm...?

[jgabry]

Thanks for opening the issue.

I would be great if ppc_error_scatter_avg() plotted $y_{pred}$ on the horizontal axis, either by default or with that as an option. (Or maybe there's already an easy way to do that; sorry if I missed it.)

That's a good point. However, we can't just plot y_pred since there's an entire distribution of y_pred for each single y value. We'd need to plot a summary of y_pred, right? In which case we'd be plotting stat(y_pred) vs stat(y - y_pred).

seems to suggest that residuals are computed as $stat(y - y_{rep})$, with $stat(y_{rep})$ not separately, explicitly computed

Yeah that's right. We're defining a predictive error to be y - y_rep and then computing a summary of the errors. This lets us compute any summary of the full distribution of predictive errors. This seems like the right thing to do to me.

@avehtari what do you think about both of these questions?

[kruschke]

@jgabry

That's a good point. However, we can't just plot y_pred since there's an entire distribution of y_pred for each single y value. We'd need to plot a summary of y_pred, right? In which case we'd be plotting stat(y_pred) vs stat(y - y_pred).

Right. My comment uses notation $y_{pred} = stat(y_{rep})$, so the plot would have the residual $y - stat(y_{rep})$ or $stat(y - y_{rep})$ on the vertical axis, and the predicted value $stat(y_{rep})$ on the horizontal axis.

I would prefer the default $stat()$ to be the median, not the mean, but that's a separate issue...

Thanks!

[avehtari]

I think Andrew's comments on this are relevant. See a recent blog post https://statmodeling.stat.columbia.edu/2025/05/11/plotting-truth-vs-predicted-value/ and section 11.3 of Regression and Other Stories: “A confusing choice: plot residuals vs. predicted values, or residuals vs. observed values?”

[jgabry]

Right. My comment uses notation y p r e d = s t a t ( y r e p ) , so the plot would have the residual y − s t a t ( y r e p ) or s t a t ( y − y r e p ) on the vertical axis, and the predicted value s t a t ( y r e p ) on the horizontal axis.

Sorry, I was working too quickly and misread your notation.

I would prefer the default s t a t ( ) to be the median, not the mean, but that's a separate issue...

Yeah I know what you mean (no pun intended). When initially developing the package many years ago the mean seemed like what users would expect the default to be for functions that computed a statistic, which is why it's always been the default for functions like ppc_stat(). However, in hindsight this wasn't the best decision, so we actually now have a warning about this if a user calls ppc_stat() without changing the default. Not optimal, but we were hesitant to change the defaults that users had come to expect. I'd be more inclined to change the defaults if we a bayesplot v2.0, which is certainly a possibility.

[avehtari]

Mean in ppc_stat() is completely different thing from mean in ppc_error_scatter_avg(), and the reason which makes mean in ppc_stat() a bad choice, does not make it make in ppc_error_scatter_avg().

[jgabry]

I think Andrew's comments on this are relevant. See a recent blog post https://statmodeling.stat.columbia.edu/2025/05/11/plotting-truth-vs-predicted-value/ and section 11.3 of Regression and Other Stories: “A confusing choice: plot residuals vs. predicted values, or residuals vs. observed values?”

@avehtari I think that chapter in RAOS argues for doing what @kruschke suggests, right? It says

Figure11.7 shows theresiduals from this model, plotted in two different ways: (a)residuals versus
fitted values, and (b) residuals versus observed values. Figure 11.7a looks reasonable: the residuals
are centered around zero for allfitted values. But Figure 11.7b looks troubling.
It turns out that the firstplot is what we should belooking at,and the second plot is misleading.

In those figures is the residual computed as y - stat(y_rep) or stat(y - y_rep)? That is, it is y minus a point prediction or it a point summary of the distribution of errors stat(y - y_rep)? I guess if using the mean/median it doesn't matter, but with the ppc_error_* functions we're calculating stat(y - y_rep) to summarize the distribution of predictive errors with any stat function. What do you think we should be plotting on each axis?

[jgabry]

Mean in ppc_stat() is completely different thing from mean in ppc_error_scatter_avg(), and the reason which makes mean in ppc_stat() a bad choice, does not make it make in ppc_error_scatter_avg().

Good point

[avehtari]

@avehtari I think that chapter in RAOS argues for doing what @kruschke suggests, right?

Yes.

In those figures is the residual computed as y - stat(y_rep) or stat(y - y_rep)?

As y - stat(y_rep), and y_pred is the same as stat(y_rep). y_rep presents the predictive distribution, y_pred is a point prediction and computed with mean in that example, and residual is y - y_pred.

As the mean of the posterior predictive distribution is the most common point prediction, I think mean is still the best default. It is good to add an option to allow computation of the point prediction with some other function. TJ did also make plots with sd, but that is not sensible point prediction and the language of error or residual is then confusing.

[jgabry]

As y - stat(y_rep), and y_pred is the same as stat(y_rep). y_rep presents the predictive distribution, y_pred is a point prediction and computed with mean in that example, and residual is y - y_pred.

Ok yeah thanks, just double checking.

As the mean of the posterior predictive distribution is the most common point prediction, I think mean is still the best default. It is good to add an option to allow computation of the point prediction with some other function. TJ did also make plots with sd, but that is not sensible point prediction and the language of error or residual is then confusing.

That sounds reasonable. I've added a comment on #349 reflecting where we ended up with this discussion.

Mean in ppc_stat() is completely different thing from mean in ppc_error_scatter_avg(), and the reason which makes mean in ppc_stat() a bad choice, does not make it make in ppc_error_scatter_avg().

Yeah this makes complete sense. I've been trying to get too many things done today and moving too quickly and not stopping to think enough apparently!