Obtain predicted probabilities from rstanarm in ordinal regression

Question

How can I generate the posterior probability distribution for each outcome for each predictor in an ordinal regression?

e.g. what I am looking for is this:

library(rstanarm)
fit_f <- MASS::polr(tobgp ~ agegp, data = esoph)
predict(fit_f,newdata=data.frame(agegp=factor(levels(esoph$agegp))),type = "probs")

Now with rstanarm I do:

fit <- stan_polr(tobgp ~ agegp, data = esoph, method = "logit",
                    prior = R2(0.2, "mean"), init_r = 0.1, seed = 12345)

But how do I obtain the distribution for the individual outcomes/predictors? I do get a distribution of probabilities using epred, but I don't understand for which outcome/predictor?

posterior_epred(fit, newdata=data.frame(agegp=factor(levels(esoph$agegp))))

Answer 1

The easiest way to do this in rstanarm is to use the posterior_predict function to obtain posterior predictions and then calculate the proportion of predictions that fall in each outcome category by observation. In code,

PPD <- posterior_predict(fit) # uses esoph
probs <- t(apply(PPD, MARGIN = 2, FUN = table) / nrow(PPD))

The matrix called probs has rows equal to the number of observations (in esoph) and columns equal to the number of categories in tobgp and each of its rows sums to 1.

head(probs)
  0-9g/day   10-19   20-29     30+
1  0.26400 0.26250 0.22875 0.24475
2  0.25650 0.26750 0.23050 0.24550
3  0.25175 0.27975 0.22450 0.24400
4  0.25575 0.26000 0.24025 0.24400
5  0.26350 0.26625 0.23575 0.23450
6  0.28275 0.26025 0.21500 0.24200