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

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

2 голосов
/ 25 марта 2020

Я хотел бы вычислить контрасты между различными уровнями параметра из моих байесовских моделей смешанных эффектов в 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)

Результат:

enter image description here

 ####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 ИЛИ SOA 200 - SOA 300 <0 с BF = 5,2112. </p>

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

Однако мне не удалось запустить # 1 односторонний (с направлением = "влево" или "вправо") 

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.

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

...