To begin, I just mention that I am totally new to bayesian statistic (and I don’t have a very solid statistical background).

I ran the following model

```
prior1 <- c(
prior(normal(0, 50), class = b),
prior(exponential(0.1), class = sd),
prior(exponential(0.1), class = sigma))
BMvpa <- brm (
RT ~ 1 + GroupC*PrimeC*CongC*EmoC*SexC
+ (1 + CongC *PrimeC*EmoC || ID)
+ (1 + GroupC*PrimeC*CongC*EmoC*SexC || Target),
data = df1,
family = exgaussian(),
prior = prior1,
warmup = 2000, iter = 5000,
chains = 3,
cores= parallel::detectCores(),
sample_prior = TRUE
)
```

Her is part of my output:

```
Predictors Estimate SE Lower Upper Rhat BF01
1 Intercept 700.175 13.470 674.330 726.617 1.000 <NA>
2 Group 49.027 33.542 -17.261 115.122 1.000 0.51
3 Prime -12.197 2.816 -17.799 -6.655 1.000 0.017
4 Congruency -15.879 2.798 -21.507 -10.435 1.001 <0.001
5 Emotion 17.092 6.860 3.740 30.373 1.000 0.32
6 Sex 5.362 24.381 -42.347 52.871 1.001 2.031
...
21 Group x Prime x Emotion 22.339 8.509 5.464 39.136 1.001 0.24
...
```

Each variable (Group, Sex, Prime, Congruency , Emotion) is dichotomous and coded -0.5 (TD, M, LSF, ICG, Joy), +0.5 (ASD, F, HSF, CG, Anger).

I would like to go more in details in the interaction Group x Prime x Emotion (represented on the figure that i posted here) and would like to know the posterior distribution regarding the effect of Prime for each group and each emotion.

I thought about 2 strategies.

1/ **First using emmeans** :

```
BMGPE_emm <- emmeans(BMvpa, ~ GroupC:PrimeC:EmoC)
BMGPE_fac <- update(BMGPE_emm, levels =list(GroupC= c("TD","ASD"), PrimeC= c("LSF","HSF"), EmoC=c("joy","anger")))
contRT1 <- as.data.frame(contrast(BMGPE_fac, method = "pairwise", by = c("GroupC","EmoC")))
```

output:

```
1 LSF - HSF TD joy 8.706978 -0.3446188 18.33661
2 LSF - HSF ASD joy 20.487280 10.2622944 30.46115
3 LSF - HSF TD anger 15.029082 6.2702713 24.67623
4 LSF - HSF ASD anger 4.412052 -5.6749680 14.60393
```

I am not sure about this because I would have expected only negative estimates (HSF primes reducing response time). Additionally, is there a possibility to compute a Bayes factor here?

2/ **Second, using the hypothesis function** (I read the blogpost of Matti Vuorre). I think it would be the best but the result I got are even more strange and I think I probably made a mistake (I was expected only negative estimates) :

```
> hypothesis(BMvpa,c(qAJ = "PrimeC + 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 qAJ 3.77 17.35 -30.51 37.79 3.56 0.78
---
> hypothesis(BMvpa,c(qAA = "PrimeC + 0.5*GroupC + 0.5*EmoC = 0"))
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 qAA 20.86 17.46 -13.54 55.14 1.63 0.62
---
> hypothesis(BMvpa,c(qTJ = "PrimeC - 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 qTJ -45.26 17.34 -78.93 -10.78 0.15 0.13 *
---
> hypothesis(BMvpa,c(qTA = "PrimeC - 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 qTA -45.26 17.34 -78.93 -10.78 0.15 0.13 *
---
```

I am sure I am wrong on what I did but I don’t manage to find what is wrong (I check the coding of my variables but I think they are right).

Could somebody help me?

(I add a small part of my data)

df1.csv (598.1 KB)