`posterior_epred` following `stan_polr` is not an expectation

Summary:

posterior_epred(stan_polr(...)) yields a $S \times N$ matrix of numbers between $0$ and $1$ that are not the expectation of anything

Description:

This was an unanticipated but anticipatable consequence of introducing posterior_epred as an alias for posterior_linpred(*, transform = TRUE). In the case of an ordinal model estimated by stan_polr, calling posterior_linpred(*, transform = TRUE) did not make sense (with more than two categories), but at least it was somewhat clear that it was applying a transformation (such as plogis, pnorm, etc.) to the linear predictor. With posterior_epred, it still does not make sense (with more than two categories) and is not at all clear that the result is not the expectation of the outcome, which does not really exist unless you reconceptualize the outcome as a one-hot vector of size $J$ in which case its expectation is a simplex.

We should probably throw an error in this case (unless there are only two categories) or else return a $S \times N \times J$ array

Reproducible Steps:

library(rstanarm)
example(stan_polr) # creates fit
str(posterior_epred(fit))

RStanARM Version:

2.21.1 but applies to several versions before that

R Version:

Doesn't matter

Operating System:

Doesn't matter

Start by running the reported example(stan_polr) reproduction and inspecting how posterior_epred delegates to posterior_linpred for stan_polr. Trace the behavior for binary and multi-category outcomes, then determine whether the intended completed behavior is an error or an S × N × J result and add coverage for that decision.