根据brm

Question

我有一个

brm

模型，我想为一个参数（长度）的固定值生成预测值。我认为我需要使用

posterior_predict

函数（或者可能是

posterior_epred

函数）来执行此操作，但是这些函数都会生成没有标签的大型矩阵，因此不可能知道实际生成的内容。

library(tidyverse)
library(brms)
set.seed(42)
df <- data.frame(mass   = rnorm(100, mean = 50, sd = 5),
                 length = rnorm(100, mean = 20, sd = 2),
                 sex    = sample(c("f", "m"), 100, replace = TRUE),
                 site   = sample(LETTERS[1:3], 100, replace = TRUE))

m <- brm(mass ~ length * sex + site, data = df, prior = prior(normal(0, 10)))

newdata <- data.frame(site   = rep(LETTERS[1:3], each = 30),
                      sex    = rep(c("f", "m"), 45),
                      length = 20)
pred1 <- bind_cols(newdata, predict(m, newdata)) # results in a usable dataframe
pred2 <- posterior_predict(m, newdata)           # results in a large matrix with no labels
pred3 <- posterior_epred(m, newdata)             # results in a large matrix with no labels

ggplot(data = pred1, aes(x = Estimate, y = site, col = sex)) + geom_violin()

Answer 1

posterior_predict

输出说明

posterior_predict()

和

posterior_epred()

都会生成 4000x90 矩阵。为什么是 4000 行？它有 4000 行，因为

brm()

拟合的模型使用 MCMC 算法对参数的联合后验分布进行采样。默认情况下，它使用四个马尔可夫链来执行此操作，每个马尔可夫链返回 1000 个样本，总共 4000 个。为什么是 90 列？因为本例中的

newdata

有 90 行。

pred1

和

pred2

的列对应于

newdata

的行，顺序相同。

顺便说一句，

posterior_predict()

和

posterior_epred()

之间的差异类似于频率分析中预测区间和均值周围置信区间之间的差异。

posterior_predict()

为您提供后验预测分布的绘图，该分布包含残差，因此具有更高的方差。相反，

posterior_epred()

给出了后验期望值的分布。

无需任何额外软件包的解决方案

要获得

newdata

中每一行的后验预测分布的中位数和 95% 等尾可信区间，您可以这样做：

posteriorpredictive_CI <- bind_cols(
  newdata,
  t(apply(pred2, 2L, quantile, probs = c(.5, .025, .975)))
)

使用 tidybayes 包的附加解决方案

您可能会发现 tidybayes 包使这变得更容易。函数

add_predicted_draws()

（用于后验预测分布）和

add_epred_draws()

（用于后验期望值）输出长格式的“整洁”数据帧。它们有 360,000 行，因为每行

newdata

都重复 4000 次，每次后验抽取一次。然后，从

ggdist

导入的函数 median_qi()（它是 tidybayes 的依赖项）可用于以等尾可信区间进行汇总。这将为您提供与之前的代码类似的结果，其中包含一些额外的 ID 列，用于标识间隔的宽度以及摘要的每一行来自

newdata

的哪一行。

library(tidybayes)

posteriorpredictive_CI <- newdata |>
  add_predicted_draws(m) |>
  median_qi()

间隔的宽度默认为 0.95，但可以使用

.width

的

median_qi()

参数进行更改。