r - R事后测试brms

Question

我有以下模型：

prior1 <- c(
  prior(normal(0, 50), class = b),
  prior(exponential(0.1), class = sd),
  prior(exponential(0.1), class = sigma))

BMvpa <- brm (
  RT ~ 1 + GroupC*PrimeC*CongC*EmoC*SexC
  + (1 + CongC *PrimeC*EmoC || ID) 
  + (1 + GroupC*PrimeC*CongC*EmoC*SexC || Target),
  data = df1,
  family = exgaussian(),
  prior = prior1,
  warmup = 2000, iter = 5000,
  chains = 3,
  cores= parallel::detectCores(),
  sample_prior = TRUE
)

这是我的输出的一部分：

                                   Predictors Estimate     SE   Lower   Upper  Rhat   BF01
1                                   Intercept  700.175 13.470 674.330 726.617 1.000   <NA>
2                                       Group   49.027 33.542 -17.261 115.122 1.000   0.51
3                                       Prime  -12.197  2.816 -17.799  -6.655 1.000  0.017
4                                  Congruency  -15.879  2.798 -21.507 -10.435 1.001 <0.001
5                                     Emotion   17.092  6.860   3.740  30.373 1.000   0.32
6                                         Sex    5.362 24.381 -42.347  52.871 1.001  2.031
...
21                    Group x Prime x Emotion   22.339  8.509   5.464  39.136 1.001   0.24
...

每个变量（Group、Sex、Prime、Congruency、Emotion）都是二分法并编码为 -0.5（TD、M、LSF、ICG、Joy）、+0.5（ASD、F、HSF、CG、Anger）。

我想更详细地了解交互组 x Prime x Emotion（如图所示），并想知道 Prime 对每个组和每种情绪的影响的后验分布。

我想了2个策略。

1/首先使用emmeans：

BMGPE_emm <- emmeans(BMvpa, ~ GroupC:PrimeC:EmoC)
BMGPE_fac <- update(BMGPE_emm, levels =list(GroupC= c("TD","ASD"), PrimeC= c("LSF","HSF"), EmoC=c("joy","anger")))
contRT1 <- as.data.frame(contrast(BMGPE_fac, method = "pairwise", by = c("GroupC","EmoC")))

输出：

1 LSF - HSF     TD   joy  8.706978 -0.3446188  18.33661
2 LSF - HSF    ASD   joy 20.487280 10.2622944  30.46115
3 LSF - HSF     TD anger 15.029082  6.2702713  24.67623
4 LSF - HSF    ASD anger  4.412052 -5.6749680  14.60393

我对此不确定，因为我只期望负估计（HSF 素数减少响应时间）。此外，是否有可能在这里计算贝叶斯因子？

2/其次，使用hypothesis函数（我阅读了Matti Vuorre 的博文）。我认为这将是最好的，但我得到的结果更奇怪，我想我可能犯了一个错误（我预计只有负估计）：

> hypothesis(BMvpa,c(qAJ = "PrimeC + 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
  Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1        qAJ     3.77     17.35   -30.51    37.79       3.56      0.78     
---

> hypothesis(BMvpa,c(qAA = "PrimeC + 0.5*GroupC + 0.5*EmoC = 0"))
Hypothesis Tests for class b:
  Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1        qAA    20.86     17.46   -13.54    55.14       1.63      0.62     
---

> hypothesis(BMvpa,c(qTJ = "PrimeC - 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
  Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1        qTJ   -45.26     17.34   -78.93   -10.78       0.15      0.13    *
---

> hypothesis(BMvpa,c(qTA = "PrimeC - 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
  Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1        qTA   -45.26     17.34   -78.93   -10.78       0.15      0.13    *
---

所以我的问题是：我怎么能得到每个子组（和BF）中素数效应的后验分布。

标签： rbayesianposthocemmeans