在序数回归中从 rstanarm 获得预测概率

编程问答 2022-04-22

问题描述

如何为有序回归中的每个预测变量生成每个结果的后验概率分布？

例如 我要找的是这个：

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

现在使用 rstanarm：

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

但是我如何获得单个结果/预测变量的分布？ 我确实使用 epred 获得了概率分布，但我不明白哪个结果/预测变量？

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

解决方法

在 rstanarm 中最简单的方法是使用 posterior_predict 函数获得后验预测，然后通过观察计算落在每个结果类别中的预测比例。在代码中，

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

名为 probs 的矩阵的行等于观察数（在 esoph 中），列等于 tobgp 中的类别数，并且其每一行的总和为 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
bayesianordinalrrrstanarm