Hi everyone,
This is my first time posting to this forum but I am hopeful someone here can help me understand the difference between the functions posterior_epred() and posterior_predict() in the context of mixed effects modelling.
These functions are used in various Bayesian R packages (e.g., rstanarm, brms, marginaleffects, brmsmargins) but it is not clear to me how they differ and, most importantly, when to use one versus the other in a mixed effects modelling context.
I tried to get some insights via Twitter, which were helpful, but I find I am still confused about this. (My mind likes to complicate things and it’s entirely possible that this issue is not as complicated as I make it out to be.)
My confusion started after reading the help file for a close relative of these functions, posterior_linpred():Posterior distribution of the (possibly transformed) linear predictor — posterior_linpred.stanreg • rstanarm (mcstan.org). According to this help file, the purpose of posterior_linpred() is to:
“Extract the posterior draws of the linear predictor, possibly transformed by the inverselink function. This function is occasionally useful, but it should be used sparingly: inference and model checking should generally be carried out using the posterior predictive distribution (i.e., using posterior_predict
).”
Prior to reading the above statement, I thought that posterior_epred() is used to compute the expected value of a response variable on the natural response scale for given predictor values and with all random effects in the model set to 0. Furthermore, I thought that posterior_predict() is used to predict new response values based on known predictor values and specified values of the random effects. But the strong wording about using the posterior_epred() function sparingly stopped me in my tracks.
Now I wonder things like:

Is posterior_epred() a special case of posterior_predict() or is it a totally different animal?

Does posterior_epred() incorporate information not only on WHAT is being estimated but also on HOW uncertainty is computed for the quantity being estimated? In particular, does posterior_epred() estimate a mean response value AND the corresponding uncertainty by completely ignoring random error AND random effects?

Does posterior_predict() incorporate information not only on WHAT is being predicted but also on HOW uncertainty is computed for the quantity being predicted? In particular, does posterior_predict() predict an individual response value AND the corresponding uncertainty by accounting for random error AND random effects (unless we instruct it otherwise)?

If posterior_epred() is a different animal than posterior_predict(), when should we use one versus the other in a mixed effects modelling context?
(Note that I am using terminology like estimation and prediction that I feel familiar with in a Frequentist context  in a Bayesian context, it seems to me that the distinction between the two is lost and everything is referred to as prediction.)
Inspired by Twitter, I went as far as to put together a table which lists several estimation and prediction scenarios I could think of in the context of a relatively simple linear mixed effects model. As a practitioner, I would really love to know whether I should posterior_epred() or posterior_predict() for each of these particular scenarios:
Any answers on the above questions and clarifications on the above table (i.e., what function to use for each task) would be much appreciated.
At this point, I am starting to believe that I shouldn’t even bother with posterior_epred() and should just consider posterior_predict() for every single scenario in my table (based on the help file mentioned above and on some of the responses I received on Twitter). But I would like to understand a bit more why if that is at all possible.
Thank you very much!
Isabella