如何在R中计算brms位置尺度模型的位置与尺度预测值？

一、生成位置和尺度部分的预测值 #

brms对位置-尺度模型的多响应部分支持直接通过fitted()或predict()函数指定resp参数，分别计算均值（位置部分，resp="mu"）和标准差（尺度部分，resp="sigma"）的预测值，步骤如下：

1. 构造预测网格数据
   先创建包含所有自变量组合的新数据集，覆盖edu_lvl的所有类别和SSI1Z的全范围，固定wave和S003为典型值（如中位数），后续可通过参数设置边缘化随机效应：

# 构造预测用的网格数据
new_data <- expand.grid(
  edu_lvl = unique(Data$edu_lvl),
  SSI1Z = seq(min(Data$SSI1Z), max(Data$SSI1Z), length.out = 100),
  wave = median(Data$wave),
  S003 = median(Data$S003)
)

2. 计算位置部分（预测均值）
   使用fitted()函数指定resp="mu"，添加re_formula=NA可边缘化随机效应，得到总体层面的边际预测：

# 计算均值的预测值及95%置信区间
pred_mu <- fitted(modelout, 
                  newdata = new_data,
                  resp = "mu",
                  re_formula = NA)
# 合并到新数据集
new_data$mu_est <- pred_mu[, "Estimate"]
new_data$mu_low <- pred_mu[, "Q2.5"]
new_data$mu_high <- pred_mu[, "Q97.5"]

3. 计算尺度部分（预测sigma）
   同理，指定resp="sigma"计算标准差的预测值：

# 计算sigma的预测值及95%置信区间
pred_sigma <- fitted(modelout, 
                     newdata = new_data,
                     resp = "sigma",
                     re_formula = NA)
# 合并到新数据集
new_data$sigma_est <- pred_sigma[, "Estimate"]
new_data$sigma_low <- pred_sigma[, "Q2.5"]
new_data$sigma_high <- pred_sigma[, "Q97.5"]

二、可视化跨层交互效应 #

直接使用ggplot2绘制所需的交互图，灵活性更高：

1. 图1：SSI1Z与预测均值的交互（按教育水平分组）

library(ggplot2)

ggplot(new_data, aes(x = SSI1Z, y = mu_est)) +
  geom_line(aes(color = as.factor(edu_lvl), linetype = as.factor(edu_lvl)), linewidth = 1) +
  geom_ribbon(aes(ymin = mu_low, ymax = mu_high, fill = as.factor(edu_lvl)), alpha = 0.2, color = NA) +
  labs(x = "SSI1Z", 
       y = "political_interestZ预测均值", 
       color = "教育水平", 
       linetype = "教育水平", 
       fill = "教育水平") +
  theme_minimal()

2. 图2：SSI1Z与预测sigma的交互（按教育水平分组）

ggplot(new_data, aes(x = SSI1Z, y = sigma_est)) +
  geom_line(aes(color = as.factor(edu_lvl), linetype = as.factor(edu_lvl)), linewidth = 1) +
  geom_ribbon(aes(ymin = sigma_low, ymax = sigma_high, fill = as.factor(edu_lvl)), alpha = 0.2, color = NA) +
  labs(x = "SSI1Z", 
       y = "political_interestZ预测标准差", 
       color = "教育水平", 
       linetype = "教育水平", 
       fill = "教育水平") +
  theme_minimal()

补充说明 #

• 若需要保留特定随机效应（如针对某个国家或调查波次），移除re_formula=NA即可，但边际预测（边缘化随机效应）通常更适合展示总体的交互模式。
• fitted()返回的是模型拟合的均值/标准差估计，而predict()返回的是因变量的预测值（包含随机误差），此处我们需要的是位置和尺度参数的估计，因此fitted()更合适。

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