rstanarm有序回归中获取各预测因子对应各结局的后验概率分布

posterior_epred针对stan_polr返回的结果是三维数组，各维度的对应规则如下：

• 第一维度：后验抽样迭代次数，长度等于模型设置的总有效后验抽样数（默认配置下为4条链各1000次，共4000次）
• 第二维度：和你传入newdata的行顺序完全对应，也就是你设置的预测因子各水平的顺序
• 第三维度：和原始数据中结局变量的因子水平顺序完全对应

1. 基础提取方法 #

首先统一存储预测数据集和后验结果，避免顺序错乱：

# 单独存储要预测的预测因子水平
new_dat <- data.frame(agegp = factor(levels(esoph$agegp)))
# 生成后验预测概率
epred <- posterior_epred(fit, newdata = new_dat)

如果要提取单个预测因子水平+单个结局的后验概率分布，直接按维度取切片即可，示例如下：

# 示例：提取agegp第2个水平（即new_dat第2行）、tobgp第3个结局类别的后验概率分布
single_post_dist <- epred[, 2, 3]

# 可直接查看分布统计特征或可视化
summary(single_post_dist)
hist(single_post_dist, main = "agegp=35-44时tobgp第三类的后验概率分布")

2. 批量规整为数据框 #

如果需要批量筛选不同预测因子、结局的分布，可将三维数组转换为带标签的长格式数据框：

library(dplyr)
library(tidyr)

epred_df <- as.data.frame.table(epred) %>%
  # 重命名各维度列
  setNames(c("iter_id", "pred_row_id", "outcome_level_id", "post_prob")) %>%
  # 匹配预测因子的实际水平
  mutate(agegp_level = new_dat$agegp[as.integer(pred_row_id)],
         # 匹配结局的实际类别
         tobgp_level = levels(esoph$tobgp)[as.integer(outcome_level_id)]) %>%
  # 保留需要的字段
  select(iter_id, agegp_level, tobgp_level, post_prob)

转换完成后可直接按标签筛选所需数据，示例：

# 筛选agegp为"65+"对应所有结局的后验概率分布
age_65_all_outcome <- epred_df %>%
  filter(agegp_level == "65+")

3. 结果验证 #

你之前用MASS::polr得到的概率是点估计，对应后验概率的均值，可通过以下代码得到和MASS输出结构完全一致的结果：

# 对每个预测因子水平、结局类别计算后验均值，和MASS的predict输出完全对应
apply(epred, c(2,3), mean)

内容的提问来源于stack exchange，提问作者Misha