How to compute contrats between levels of a parameter in bayesian mixed-effects models and produce bayes factors in R?

Question

I would like to compute contrats between different levels of a parameter from my bayesian mixed-effects models in R, and produce bayes factors. My outcome (Jud) is binary (1=Yes/In synch, 0=No/Out of synch) and the parameter SOAsF is a factor with 6 levels (0, 100, 200, 300, 400, 500).

Following different tutorials/functions [#1],[#2] and [#3], here is my code with the 3 different manners:

 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)

Outcome:

 ####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.

So taking the example of the contrasts 200 vs 300. Are sounds presented at SOA 200 equally judged as in synch (yes) compared to sounds presented at SOA 300?

Manner #1 seems to provide evidence toward the null hypothesis SOA 200 - SOA 300 = 0 with BF = 0.01; so they seem equally jusged as in synch?

Manner #2 seems to provide little evidence toward the null hypothesis SOA 200 = SOA 300 with the evidence being 1.988631% 95%CI [-3.707585, 7.680694].

Manner #3 seems to provide evidence toward the alternative hypothesis SOA 200 > SOA 300 OR SOA 200 - SOA 300 < 0 with BF = 5.2112.

Do I find differences because #1 is two-sided while #3 is one-sided?

However, I did not manage to run #1 one-sided (with direction = "left" or "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

Or #3 two-sided ( hypothesis(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.

I am stuck, any help would be appreciated.Thanks.

How to compute contrats between levels of a parameter in bayesian mixed-effects models and produce bayes factors in R?

Answers (1)

Related Questions