1

我有一个来自参与者的数据集，它提供了与不同幅度的奖励（因子 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))

接下来，我想从两个特定对比（即两个成对比较）的后验分布中抽取样本。首先，我使用 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 的函数terior_samples 在这种特定情况下不起作用，因为我的摘要模型输出没有感兴趣的对比，或者，
  • 另一个允许我以 emmeans 的方式指定成对比较的函数，以及一个允许我提取其后验图的函数。

任何想法都受到高度赞赏！

4

1 回答 1

0

关于第一个问题：与大多数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包含所需对比的后验，如果您想要对这些平局进行点估计，您可以使用均值、中位数或任何其他想。

于 2022-02-14T22:54:58.067 回答