在使用 emmeans 和 hpd.summary 后从后验中提取绘图
Extracting draws from posterior after using emmeans and hpd.summary
我有一个来自参与者的数据集，它提供了与不同幅度的奖励（因子 pval，水平 small/medium/large）和延迟（因子时间，级别 delayed/immediate)。数据的一个子集如下所示：
structure(list(ppnr = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("7", "8"), class = "factor"),
time = structure(c(2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L), .Label = c("del", "imm"), class = "factor"), pval = structure(c(1L,
1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L), .Label = c("pval_L",
"pval_M", "pval_S"), class = "factor"), rating = c(66, 55,
81, 11, 30, 11, 18, 28, 85, 61, 6, 76), stimJPG = structure(c(1L,
2L, 3L, 5L, 4L, 6L, 5L, 1L, 3L, 2L, 6L, 4L), .Label = c("pStim01.jpg",
"pStim02.jpg", "pStim03.jpg", "pStim04.jpg", "pStim05.jpg",
"pStim06.jpg"), class = "factor")), row.names = 283:294, class = "data.frame")
ppnr time pval rating stimJPG
283 7 imm pval_L 66 pStim01.jpg
284 7 del pval_L 55 pStim02.jpg
285 7 imm pval_M 81 pStim03.jpg
286 7 del pval_M 11 pStim05.jpg
287 7 imm pval_S 30 pStim04.jpg
288 7 del pval_S 11 pStim06.jpg
289 8 imm pval_L 18 pStim05.jpg
290 8 del pval_L 28 pStim01.jpg
291 8 imm pval_M 85 pStim03.jpg
292 8 del pval_M 61 pStim02.jpg
293 8 imm pval_S 6 pStim06.jpg
294 8 del pval_S 76 pStim04.jpg
为了调查评分是否受到与线索相关的奖励的时间和幅度的影响，我 运行 brms 中的以下模型：
n_chains <- 4
n_threads <- 2
options(contrasts = c("contr.sum", "contr.poly"))
model <- brm(rating ~ time*pval + (1 + time + pval | ppnr) + (1 + time * pval | stimJPG), data = data, backend = "cmdstanr", chains = n_chains, cores = n_chains, threads = threading(n_threads), iter = 4000, warmup = 2000, control = list(adapt_delta = 0.9999, max_treedepth = 15))
接下来，我想从后验分布中抽取样本用于两个特定的对比（即两个成对比较）。首先，我使用 emmeans 获得了这些对比的估计值。我原则上可以使用函数 gather_emmeans_draws（来自 tidybayes 包）从这些对比的后部抽取样本而不会出现问题。然而，退后一步，emmeans 使用中位数作为贝叶斯模型的点估计，而我想使用均值。通过在 emmeans 对象上使用 hpd.summary 可以获得平均值。但是，这会将 emmeans 创建的 emmGrid 对象转换为 summary_emm 对象。不幸的是，gather_emmeans_draws() 不接受 summery_emm 对象，但似乎只接受 emmGrid 对象（或一般的 S4 对象）。参见：
emm_withmedian <- emmeans(model, pairwise ~ pval * time)$contrasts
emm_withmean <- hpd.summary(emm_withmedian, point.est = mean)
#This results in all pairwise comparisons, but I am only interested in 2 of these:
contrast estimate lower.HPD upper.HPD
pval_L del - pval_M del 14.79 3.87 26.17
pval_L del - pval_S del 26.67 11.69 42.55
pval_L del - pval_L imm -5.85 -17.98 7.67
pval_L del - pval_M imm 9.51 -3.17 22.61
pval_L del - pval_S imm 17.70 4.23 31.75
pval_M del - pval_S del 11.89 -1.45 26.84
pval_M del - pval_L imm -20.64 -33.83 -7.43
pval_M del - pval_M imm -5.28 -18.19 6.96
pval_M del - pval_S imm 2.91 -9.40 16.33
pval_S del - pval_L imm -32.53 -47.46 -18.05
pval_S del - pval_M imm -17.16 -29.95 -3.68
pval_S del - pval_S imm -8.98 -22.43 5.10
pval_L imm - pval_M imm 15.36 4.28 26.58
pval_L imm - pval_S imm 23.55 9.58 39.50
pval_M imm - pval_S imm 8.19 -4.94 22.43
Point estimate displayed: mean
HPD interval probability: 0.95
#I then want to draw from the posterior, but that's where it goes wrong:
posteriorsamples <- gather_emmeans_draws(emm_withmean)
Error in as_tibble(object@grid) :
trying to get slot "grid" from an object (class "summary_emm") that is not an S4 object
#Just for comparison's sake, if I would do the following, it would be no problem, because it uses an emmGrid object as input:
posteriorsamples2 <- gather_emmeans_draws(emm_withmedian)
因此，如果我直接从 emmGrid 对象 (emm_withmedian) 工作，我似乎只能从后验绘制，迫使我使用中位数而不是平均值。
我已经尝试使用 as.emmGrid() 将 summary_emm 对象转换为 emmGrid 对象，但这不起作用，并出现以下错误：nrow 错误(V) : 缺少参数“V”，没有默认值。
我已经查看了这两条错误消息，但还没有找到解决它们的方法。我还确保更新了所有使用的包，但这也没有帮助。
因此，我正在寻找：
- 一种将 summary_emm 对象转换为 emmGrid 对象（或 gather_emmeans_draws 接受的任何其他对象）的方法，
或者，
- 另一个函数，允许我以 gather_emmeans_draws 的方式从 emmeans 对象的后部绘制。不幸的是，来自 brms 的函数 posterior_samples 在这种特定情况下不起作用，因为我在摘要模型输出中没有感兴趣的对比，
或者，
- 另一个允许我以 emmeans 的方式指定成对比较的函数，以及一个允许我提取后验图的函数。
非常感谢任何想法！
关于第一个问题：与大多数 summary 方法一样，返回的对象只是一个摘要，它不包含将其转换回对象的信息总结。但是，原始 emmGrid 对象确实包含所有需要的内容。
另一个障碍是试图从你不想要的对比中工作，而不是得到你想要的对比。通常最好分两个单独的步骤进行均值和对比。做起来很简单：
EMM <- emmeans(model, ~ pval * time)
CON <- contrast(EMM, list(c1 = c(<desired coefs>), c2 = c(<desired coefs>)))
draws <- coda::as.mcmc(CON)
后者在内部进行相同的计算，使 linfct 槽成为单位矩阵，post.beta 槽等于估计值。
无论哪种方式，您看到的摘要显示中位数这一事实并不重要 - draws 包含所需对比的后验图，如果您想对这些图进行点估计，您可以使用均值、中位数或您想要的任何其他值。
我有一个来自参与者的数据集，它提供了与不同幅度的奖励（因子 pval，水平 small/medium/large）和延迟（因子时间，级别 delayed/immediate)。数据的一个子集如下所示：
structure(list(ppnr = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("7", "8"), class = "factor"),
time = structure(c(2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L), .Label = c("del", "imm"), class = "factor"), pval = structure(c(1L,
1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L), .Label = c("pval_L",
"pval_M", "pval_S"), class = "factor"), rating = c(66, 55,
81, 11, 30, 11, 18, 28, 85, 61, 6, 76), stimJPG = structure(c(1L,
2L, 3L, 5L, 4L, 6L, 5L, 1L, 3L, 2L, 6L, 4L), .Label = c("pStim01.jpg",
"pStim02.jpg", "pStim03.jpg", "pStim04.jpg", "pStim05.jpg",
"pStim06.jpg"), class = "factor")), row.names = 283:294, class = "data.frame")
ppnr time pval rating stimJPG
283 7 imm pval_L 66 pStim01.jpg
284 7 del pval_L 55 pStim02.jpg
285 7 imm pval_M 81 pStim03.jpg
286 7 del pval_M 11 pStim05.jpg
287 7 imm pval_S 30 pStim04.jpg
288 7 del pval_S 11 pStim06.jpg
289 8 imm pval_L 18 pStim05.jpg
290 8 del pval_L 28 pStim01.jpg
291 8 imm pval_M 85 pStim03.jpg
292 8 del pval_M 61 pStim02.jpg
293 8 imm pval_S 6 pStim06.jpg
294 8 del pval_S 76 pStim04.jpg
为了调查评分是否受到与线索相关的奖励的时间和幅度的影响，我 运行 brms 中的以下模型：
n_chains <- 4
n_threads <- 2
options(contrasts = c("contr.sum", "contr.poly"))
model <- brm(rating ~ time*pval + (1 + time + pval | ppnr) + (1 + time * pval | stimJPG), data = data, backend = "cmdstanr", chains = n_chains, cores = n_chains, threads = threading(n_threads), iter = 4000, warmup = 2000, control = list(adapt_delta = 0.9999, max_treedepth = 15))
接下来，我想从后验分布中抽取样本用于两个特定的对比（即两个成对比较）。首先，我使用 emmeans 获得了这些对比的估计值。我原则上可以使用函数 gather_emmeans_draws（来自 tidybayes 包）从这些对比的后部抽取样本而不会出现问题。然而，退后一步，emmeans 使用中位数作为贝叶斯模型的点估计，而我想使用均值。通过在 emmeans 对象上使用 hpd.summary 可以获得平均值。但是，这会将 emmeans 创建的 emmGrid 对象转换为 summary_emm 对象。不幸的是，gather_emmeans_draws() 不接受 summery_emm 对象，但似乎只接受 emmGrid 对象（或一般的 S4 对象）。参见：
emm_withmedian <- emmeans(model, pairwise ~ pval * time)$contrasts
emm_withmean <- hpd.summary(emm_withmedian, point.est = mean)
#This results in all pairwise comparisons, but I am only interested in 2 of these:
contrast estimate lower.HPD upper.HPD
pval_L del - pval_M del 14.79 3.87 26.17
pval_L del - pval_S del 26.67 11.69 42.55
pval_L del - pval_L imm -5.85 -17.98 7.67
pval_L del - pval_M imm 9.51 -3.17 22.61
pval_L del - pval_S imm 17.70 4.23 31.75
pval_M del - pval_S del 11.89 -1.45 26.84
pval_M del - pval_L imm -20.64 -33.83 -7.43
pval_M del - pval_M imm -5.28 -18.19 6.96
pval_M del - pval_S imm 2.91 -9.40 16.33
pval_S del - pval_L imm -32.53 -47.46 -18.05
pval_S del - pval_M imm -17.16 -29.95 -3.68
pval_S del - pval_S imm -8.98 -22.43 5.10
pval_L imm - pval_M imm 15.36 4.28 26.58
pval_L imm - pval_S imm 23.55 9.58 39.50
pval_M imm - pval_S imm 8.19 -4.94 22.43
Point estimate displayed: mean
HPD interval probability: 0.95
#I then want to draw from the posterior, but that's where it goes wrong:
posteriorsamples <- gather_emmeans_draws(emm_withmean)
Error in as_tibble(object@grid) :
trying to get slot "grid" from an object (class "summary_emm") that is not an S4 object
#Just for comparison's sake, if I would do the following, it would be no problem, because it uses an emmGrid object as input:
posteriorsamples2 <- gather_emmeans_draws(emm_withmedian)
因此，如果我直接从 emmGrid 对象 (emm_withmedian) 工作，我似乎只能从后验绘制，迫使我使用中位数而不是平均值。
我已经尝试使用 as.emmGrid() 将 summary_emm 对象转换为 emmGrid 对象，但这不起作用，并出现以下错误：nrow 错误(V) : 缺少参数“V”，没有默认值。
我已经查看了这两条错误消息，但还没有找到解决它们的方法。我还确保更新了所有使用的包，但这也没有帮助。
因此，我正在寻找：
- 一种将 summary_emm 对象转换为 emmGrid 对象（或 gather_emmeans_draws 接受的任何其他对象）的方法， 或者，
- 另一个函数，允许我以 gather_emmeans_draws 的方式从 emmeans 对象的后部绘制。不幸的是，来自 brms 的函数 posterior_samples 在这种特定情况下不起作用，因为我在摘要模型输出中没有感兴趣的对比， 或者，
- 另一个允许我以 emmeans 的方式指定成对比较的函数，以及一个允许我提取后验图的函数。
非常感谢任何想法！
关于第一个问题：与大多数 summary 方法一样，返回的对象只是一个摘要，它不包含将其转换回对象的信息总结。但是，原始 emmGrid 对象确实包含所有需要的内容。
另一个障碍是试图从你不想要的对比中工作，而不是得到你想要的对比。通常最好分两个单独的步骤进行均值和对比。做起来很简单：
EMM <- emmeans(model, ~ pval * time)
CON <- contrast(EMM, list(c1 = c(<desired coefs>), c2 = c(<desired coefs>)))
draws <- coda::as.mcmc(CON)
后者在内部进行相同的计算，使 linfct 槽成为单位矩阵，post.beta 槽等于估计值。
无论哪种方式，您看到的摘要显示中位数这一事实并不重要 - draws 包含所需对比的后验图，如果您想对这些图进行点估计，您可以使用均值、中位数或您想要的任何其他值。