Как вычислить контрасты между уровнями параметра в байесовских моделях со смешанными эффектами и получить байесовские факторы в R?

Question

Я хотел бы вычислить контрасты между разными уровнями параметра из моих байесовских моделей смешанных эффектов в R и создать байесовские факторы. Мой результат (Jud) является двоичным (1= Да / Синхронно, 0= Нет / Несинхронно), а параметр SOAsF - это фактор с 6 уровнями (0, 100, 200, 300, 400, 500).

Следуя различным руководствам / функциям [#1],[# 2] и [#3], вот мой код с 3 разными способами:

 library(emmeans)
    library(brms)
    library(modelbased)

    brm_acc_1<-brm(Jud ~ SOAsF +(1|pxID),data =dat_long, family=bernoulli("logit"), prior = set_prior('normal(0,10)'), iter = 2000, chains=4,  save_all_pars = TRUE)
    summary(brm_acc_1)
    brms::conditional_effects(brm_acc_1)

    ####1           
    groups <- emmeans(brm_acc_1, ~ SOAsF)
    group_diff <- pairs(groups)
    (groups_all <- rbind(groups, group_diff))        
    bayesfactor_parameters(groups_all, prior = brm_acc_1, direction = "two-sided", effects = c("fixed", "random", "all"))

    ####2   
    ppc <- pp_check(brm_acc_1, type = "stat_grouped", group = "SOAsF")
    #contrast 200 - 300
    contrast_300_200 <- ppc$data$value[ppc$data$group == "200"] - ppc$data$value[ppc$data$group == "300"]
    quantile(contrast_300_200*100, probs = c(.5, .025, .975))

    ####3
    h_1 <- hypothesis(brm_acc_1, "SOAsF200 < SOAsF300")
    print(h1, digits = 4)
    h2 <- hypothesis(brm_acc_1, "SOAsF200 > SOAsF300")
    print(h2, digits = 4)

Итог:

 ####1
# Bayes Factor (Savage-Dickey density ratio)

    Parameter    |       BF
    -----------------------
    0, .         | 8.08e-03
    100, .       |     0.61
    200, .       | 7.29e+03
    300, .       |    67.77
    400, .       |    21.81
    500, .       |     0.28
    ., 0 - 100   |     2.75
    ., 0 - 200   | 1.90e+05
    ., 0 - 300   |   410.42
    ., 0 - 400   |   570.11
    ., 0 - 500   |     1.03
    ., 100 - 200 |      0.5
    ., 100 - 300 |     0.05
    ., 100 - 400 |     0.02
    ., 100 - 500 | 7.13e-03
    ., 200 - 300 |     0.01
    ., 200 - 400 |     0.01
    ., 200 - 500 |     1.11
    ., 300 - 400 | 7.21e-03
    ., 300 - 500 |      0.1
    ., 400 - 500 |     0.04

    * Evidence Against The Null: [0]

  ####2
   50%      2.5%     97.5% 
 1.988631 -3.707585  7.680694 

  ####3

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper
1 (SOAsF200)-(SOAsF... < 0   0.0881    0.0903  -0.0612   0.2372
  Evid.Ratio Post.Prob Star
1     0.1919     0.161     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper
1 (SOAsF200)-(SOAsF... > 0   0.0881    0.0903  -0.0612   0.2372
  Evid.Ratio Post.Prob Star
1     5.2112     0.839     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Итак, возьмем пример контрастов 200 и 300. Являются ли звуки, представленные в SOA 200, одинаково оцениваемыми как синхронизированные (да), по сравнению со звуками, представленными в SOA 300?

Способ №1, кажется, обеспечивает свидетельство нулевой гипотезы SOA 200 - SOA 300 = 0 с BF = 0,01; так они кажутся одинаково оправданными как синхронными?

Способ №2, кажется, мало доказывает нулевую гипотезу SOA 200 = SOA 300 с доказательством 1,988631% 95% ДИ [-3,707585, 7,680694].

Способ № 3, кажется, дает доказательства альтернативной гипотезы SOA 200 > SOA 300 OR SOA 200 - SOA 300 < 0 с BF = 5,2112.

Нахожу ли я различия, потому что №1 двусторонний, а №3 односторонний?

Однако мне не удалось запустить №1 односторонний (с direction = "left" или "right")

bayesfactor_parameters(groups_all, prior = brm_acc_1, direction = ">",  effects = c("fixed", "random", "all") )
Computation of Bayes factors: sampling priors, please wait...
Error in `$<-.data.frame`(`*tmp*`, "ind", value = 8L) : 
  replacement has 1 row, data has 0

Или № 3 двусторонний (гипотеза (brm_acc_1, " SOAsF200 - SOAsF300 = 0 "))

Hypothesis Tests for class b:
               Hypothesis Estimate Est.Error CI.Lower CI.Upper
1 (SOAsF200-SOAsF300) = 0   0.0881    0.0903  -0.0914   0.2682
  Evid.Ratio Post.Prob Star
1         NA        NA     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Я застрял, любая помощь будет принята с благодарностью. Спасибо.

Answer 1

У меня есть частичный ответ на вопрос, почему у вас есть разница между №1 и №3. В № 1 вычисляется коэффициент Байеса, что означает, что функция дает вам свидетельство H1, разделенное на свидетельство H0. В № 3 вы получаете доказательства вашей гипотезы >0 (H1), но это не байесовский фактор. Если вы хотите факторизовать, вам следует вычислить доказательства для гипотезы =0 (H0) и сделать коэффициент H1/H0. (для БФ10).

Ваше здоровье.