Hello,

I am currently running the following model in which I aim to summarize the response of males and females to three treatments (D, F, M), with Pair ID (a pair being a couple) and Individual as random effects.

```
# Set prior
prior <- set_prior("normal(0,4)", class = "b")
# Run model
fit <- brm(Response ~ 1 + Treatment + Sex + Treatment:Sex + (1|Pair.ID/Individual),
data = Dmod,
warmup = 2500, iter = 10000,
cores = 2, chains = 4, prior = prior, sample_prior = T)
```

I am interested in getting pairwise comparisons for each sex in each treatment (in the same way as frequentists perform a post-hoc Tukey after running an ANOVA), but I do not know exactly how to do it in `brms`

.

I have tried using the `emmeans`

package for that:

```
emmeans(fit , specs = pairwise ~ Treatment:Sex)
```

Which seems to be giving what I am looking for:

```
$contrasts
contrast estimate lower.HPD upper.HPD
D Female - F Female -0.531 -1.0061 -0.0509
D Female - M Female 1.175 0.6762 1.6507
D Female - D Male 1.393 0.9143 1.8812
D Female - F Male -0.290 -0.7796 0.2056
D Female - M Male 1.811 1.3284 2.3208
F Female - M Female 1.704 1.2303 2.1950
F Female - D Male 1.923 1.4480 2.4118
F Female - F Male 0.241 -0.2591 0.7167
F Female - M Male 2.340 1.8366 2.8298
M Female - D Male 0.219 -0.2842 0.7079
M Female - F Male -1.465 -1.9528 -0.9596
M Female - M Male 0.637 0.1292 1.1362
D Male - F Male -1.683 -2.1655 -1.2128
D Male - M Male 0.418 -0.0766 0.8919
F Male - M Male 2.101 1.6098 2.5859
Point estimate displayed: median
HPD interval probability: 0.95
```

However, I do not know if this approach is correct. Would there be an analog of this in `brms`

?

I checked Matti Vuorre [Sometimes I R: How to calculate contrasts from a fitted brms model)] post on the use of hypothesis to do contrasts in `brms`

, but I am not sure if I understood it right. If I am not wrong, for each comparison I should write:

```
q1 <- c("Intercept = TreatmentF")
q2 <- c("Intercept = TreatmentM")
q3 <- c("TreatmentF + Intercept = TreatmentM + Intercept")
... # And so on until q15
# Run contrasts
hypothesis(fit, c(q1, q2, q3, ..., q15))
```

If this is correct, how should I specify the interaction terms? Maybe…

```
q7 <- c("TreatmentF + Intercept = TreatmentF:SexMale + Intercept")
```

If not, I would really appreciate is someone could let me know (1) if the use of `emmeans`

in a `brms`

model is reliable, and (2) what is the approach to perform contrasts in `brms`

?

Thanks in advance,

D