生成多预测因子同步变化序列 实现贝叶斯后验预测

你原有代码里的预测输入集仅包含-1SD、均值、+1SD三个离散点位，所以预测结果只有三个独立点，绘制的拟合线自然不平滑。要实现三个预测因子同步在-1SD到+1SD区间连续变动，不需要做全因子网格组合（全因子组合适用于变量独立变动的场景），只需要先生成统一的偏移进度序列，再按每个变量自身的均值、标准差换算为实际取值即可，你常用的seq()或modelr::seq_range()都可以用来生成这个梯度序列。

核心实现逻辑 #

• 先设定梯度点密度，一般取50-200个点即可满足绘图平滑需求，点数过大会徒增预测计算耗时
• 生成0到1区间的等距序列作为统一偏移比例：比例0对应所有变量取-1SD值，比例0.5对应所有变量取均值，比例1对应所有变量取+1SD值
• 对每个变量按公式 取值 = (均值 - 1*SD) + 偏移比例 * 2*SD 换算得到对应位置的实际值，保证三个变量完全同步偏移
• 将生成的等长向量合并为tibble格式的预测输入集，传入epred_draws()即可得到连续梯度的后验预测结果

修改后的可复现代码 #

library(brms)
library(tidybayes)
library(ggplot2) 
library(ggthemes) 
library(dplyr)

# 生成示例数据集 加随机种子方便结果复现
set.seed(123)
data <- tibble(
  outcome = rnorm(100, 2, 2), 
  var_1 = rnorm(100, 5, 2), 
  var_2 = rnorm(100, 8, 2), 
  var_3 = rnorm(100, 10, 2)
) 

# 拟合模型 关闭采样过程的冗余日志输出
m1 <- brms::brm(outcome ~ var_1 + var_2 + var_3, data, refresh = 0)

# 提前计算每个变量的均值、标准差锚点
v1_mean <- mean(data$var_1)
v1_sd <- sd(data$var_1)
v2_mean <- mean(data$var_2)
v2_sd <- sd(data$var_2)
v3_mean <- mean(data$var_3)
v3_sd <- sd(data$var_3)

# 生成连续梯度预测数据集
n_point <- 100 # 梯度点数量 可按需调整
shift_seq <- seq(0, 1, length.out = n_point) # 统一偏移比例序列
new_data <- tibble(
  var_1 = (v1_mean - v1_sd) + shift_seq * 2 * v1_sd,
  var_2 = (v2_mean - v2_sd) + shift_seq * 2 * v2_sd,
  var_3 = (v3_mean - v3_sd) + shift_seq * 2 * v3_sd
)

# 生成后验预期值预测
pred_1 <- m1 %>% 
  tidybayes::epred_draws(new_data) 

# 绘制平滑预测图
plot_1 <- ggplot(pred_1, aes(x = var_1, y = .epred)) +
  stat_lineribbon() + 
  scale_fill_brewer(palette = "Reds") +
  labs(x = "三个变量同步偏移幅度", y = "结果变量预期值",
       fill = "可信区间") +
  ggthemes::theme_pander() +
  theme(legend.position = "bottom") + 
  scale_x_continuous(
    breaks = c(v1_mean - v1_sd, v1_mean, v1_mean + v1_sd), 
    labels =  c("-1 SD", "均值", "+1 SD")
  )

print(plot_1)

场景拓展提示 #

如果你后续需要做变量独立变动的预测（比如固定另外两个变量在均值，仅观测单个变量变动的边际效应），可以用tidyr::expand_grid()给每个变量传入单独的梯度序列生成全组合数据集，不要用base版本的expand.grid()，它生成的普通data.frame容易丢失变量属性，和tidyverse系列函数兼容性差。

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