0

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))))

Misha
  • 3,114
  • 8
  • 39
  • 60

1 Answers1

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
Ben Goodrich
  • 4,870
  • 1
  • 19
  • 18
  • Sweet—but if I want the entire posterior probability distribution for each cell for plotting- how can this be done? The _epred function- which probability is this? – Misha Jun 25 '21 at 13:30