Pairwise comparisons in emmeans and brms

Modeling
1

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

2

@Dan-Zapata hello, I haven’t tried the ‘emmeans’ methods much for brms models but I suspect that this will fulfil what you’re looking for (they are the posterior mean and highest posterior density intervals, for the difference in the population predicted value of the response).

In this case I’d say that you can also reach this manually via the fitted() function by passing in newdata, setting summary = FALSE (so that you are returned the samples), which will generate predictions corresponding to the emmeans output. This is how the output of conditional_effects() for example is generated. Then you can obtain the contrasts by calculating the difference between the predictions across the samples. The interaction terms are all implicitly taken into account.