0
$\begingroup$

I would like to compute contrasts 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:

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.

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

$\endgroup$

1 Answer 1

1
$\begingroup$

I am going to focus on the hypothesis function, because I have used it in my own research, and therefore know generally how it works. That said, I am not a trained statistician, so my advice is going to be practical, rather than theoretical.

In order to make #3 two-sided, you can use the following code:

hyp_name <- hypothesis(brm_acc_1, "SOAsF200 = SOAsF300")

Also, it may be of use to you to know that you can call the plot() function around a hypothesis object, like this:

plot(hyp_name)

If you are interested in comparing the overlap of the posteriors for 200 and 300, you can do that with the code below. This should give you an overlap coefficient, which will tell you how similar the posteriors are

library(bayestestR)

Post_samples_200 <- posterior_samples(brm_acc_1, pars = "b_SOAsF200")
Post_samples_300 <- posterior_samples(brm_acc_1, pars = "b_SOAsF300")

overlap(Post_samples_200, Post_samples_300)

Also, some unsolicited advice, you are throwing out a lot of great information when you plot your model as you did above. Try calling the code below to show the posterior distributions of each parameter level.

library(bayesplot)
library(ggplot2)
bayesplot_theme_set(theme_default(base_size = 11, base_family = "sans"))
model_array <- as.array(brm_acc_1)
model_plot <- mcmc_areas(model_array), prob = 0.89)
model_plot
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.