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

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:

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

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