使用emmeans和hpd.summary提取后的后验图
我有一个来自参与者的数据集，它提供了与不同程度的奖励(因子pval，级别大小为中小/大)和延迟(因子时间，水平延迟/立即)相关的刺激的评分(从0到100)。数据的子集如下所示：
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)) 接下来，我想从后验分布中抽取两个具体对比的样本(即两个配对比较)。首先，我得到了这些对比的估计使用电子显微镜。原则上，我可以使用函数gather_emmeans_draws (来自tidybayes包)从这些对比的后面抽取样本，没有问题。然而，退一步，mean意味着使用中值作为贝叶斯模型的点估计，而我想使用的是平均值。通过在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)工作时，我似乎只能从后面提取，迫使我使用中间值而不是平均值。
我已经尝试使用summary_emm ()将emmGrid对象转换为emmGrid对象，但这不起作用，并给出了以下错误: nrow(V)中的错误:参数"V“缺失，没有默认情况。
我已经查看了这两条错误消息，但还没有找到解决它们的方法。我还确保更新所有使用的软件包，但这也没有帮助。
因此，我希望：
- --一种将summary_emm对象转换为emmGrid对象(或gather_emmeans_draws接受的任何其他对象)的方法，OR，
- 是另一个函数，它允许我以gather_emmeans_draws的方式从emmeans对象的后面提取。不幸的是，来自brms的函数posterior_samples在这种特殊情况下不起作用，因为我对我的摘要模型输出OR、
- 没有兴趣的对比，这是另一个函数，它允许我以model的方式指定成对的比较，以及一个允许我提取这方面的后绘制的函数。
任何想法都是非常感谢的！
回答 1
Stack Overflow用户
发布于 2022-02-14 22:54:58
关于第一个问题:与大多数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包含了想要的对比的后置图，如果你想要对这些图的点估计，你可以使用平均值，中值，或者任何你想要的东西。
https://stackoverflow.com/questions/71114243
复制相似问题
腾讯云开发者
