brms模型两种预测均值方法结果不同，如何判断正确性？

我尝试绘制brms模型的密度曲线并添加一条展示预测均值的线。我用两种不同方法计算均值，得到的结果略有差异。我是否忽略了什么？该如何判断哪种方法正确？

library(tidyverse)
library(brms)
library(tidybayes)

set.seed(42)
df <- data_frame(site = rep(LETTERS[1:4], each = 25),
                 sex = rep(c("f", "m"), 50),
                 length = rnorm(100, 20, 2),
                 mass = rnorm(100, 50, 5))

m <- brm(mass ~ length + site + sex, 
         data = df, 
         prior = c(prior(normal(0, 10), class = Intercept),
                   prior(normal(0, 1), class = b),
                   prior(normal(0, 1), class = sigma)))

# 获取后验分布的预测值
pred <- data_frame(site = rep(LETTERS[1:4], each = 2),
                   sex = rep(c("f", "m"), 4),
                   length = 20) %>% 
  add_predicted_draws(m)

# 方法1：从预测值中计算均值
predM <- group_by(pred, site, sex) %>% summarise(mP = mean(.prediction))

# 方法2：从后验汇总结果计算均值
postSum <- posterior_summary(m)
calcM <- data_frame(site = rep(LETTERS[1:4], each = 2),
                    sex = rep(c("f", "m"), 4),
                    length = 20) %>% 
  mutate(mC = postSum[1] + # 截距项
              postSum[2] * length + # length的斜率
              case_when(site == "B" ~ postSum[3], # 不同site的截距调整
                        site == "C" ~ postSum[4],
                        site == "D" ~ postSum[5],
                        TRUE ~ 0) +
              if_else(sex == "m", postSum[6], 0)) # 不同sex的截距调整


ggplot(data = pred, aes(x = .prediction, col = sex)) + geom_density() +
  geom_vline(data = predM, aes(xintercept = mP, col = sex), linetype = "dashed") +
  geom_vline(data = calcM, aes(xintercept = mC, col = sex), linetype = "solid") +
  facet_wrap(~site)

差异原因 #

两种方法的核心区别在于是否包含模型随机误差项（sigma），以及计算逻辑的不同：

• 方法1：add_predicted_draws生成的是**后验预测分布（PPD）**样本，模拟的是真实观测的结果，每个样本都加入了模型的随机噪声。虽然噪声的理论均值为0，但由于后验样本数量有限，实际计算出的PPD均值会和无噪声的理论均值存在微小抽样偏差。
• 方法2：直接用各参数的后验均值代入模型公式，计算的是**线性预测器（μ）**的均值——这是模型拟合的「潜在平均水平」，完全不包含随机误差。

哪种方法正确？ #

取决于你的展示目标：

• 若要展示新观测数据的平均预期位置，选方法1，它更贴近实际观测的分布中心。
• 若要展示模型拟合的理论均值（无噪声的潜在水平），选方法2，它直接反映模型参数的汇总结果。

验证方法 #

如果想让两种方法结果完全一致，可以把方法1中的add_predicted_draws替换为add_fitted_draws——后者生成的是不含随机噪声的线性预测器样本，此时计算出的均值会和方法2的结果几乎完全重合：

# 替换为add_fitted_draws获取线性预测器样本
pred <- data_frame(site = rep(LETTERS[1:4], each = 2),
                   sex = rep(c("f", "m"), 4),
                   length = 20) %>% 
  add_fitted_draws(m)

# 重新计算均值
predM <- group_by(pred, site, sex) %>% summarise(mP = mean(.value))

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