Multiple comparisons with brmsfit

[vmikk]

Dear Paul and brms-users! Could you please tell me how to perform multiple comparisons with brmsfit object?

Here is a small example with dummy data:

library(brms)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
set.seed(111)

# Generate data
datt <- data.frame(
    A = rep(c("a1", "a2"), times=50),    # fixed effect 1
    B = rep(c("b1", "b2"), each=50),     # fixed effect 2
    N = rpois(100, 4),                   # dependent variable (number of successes)
    ID = rep(1:20, each = 5))            # random effect
datt$Tot <- datt$N + rpois(100, 5)       # total number of trials

# Fit the binomial model
fit <- brm(N | trials(Tot) ~ A * B + (1|ID), data = datt, family = binomial("logit"), cluster = 1, chains = 1)

# Here we can check the predicted values of the response variable:
pred.data <- data.frame(A = rep(c("a1", "a2"), times=2), B = rep(c("b1", "b2"), each=2), Tot = 100)
data.frame(pred.data,  predict(fit, newdata = pred.data, re_formula = NA))

Is it possible to test for the difference of factor combinations (e.g. a1-b1 vs a2-b2)? I'm looking something similar to the glht function from multcomp package. May I state that two combinations differ if their credible intervals (from predict) doesn't overlap? Or Estimate for one factor combination (e.g., a1-b1) doesn't belong to the interval for the other factor combination (e.g., a2-b2)? Or it is somehow possible to test it with hypothesis.brmsfit?

I really appreciate any advice you can provide.

With best regards, Vladimir

[paul-buerkner]

You can actually do that with the hypothesis method although you have to adjust your alpha-level yourself (using the alpha argument) if you want to make multiple comparisons.

First, we have to make sure we understand the coding of our factors (it's the same as for other packages). In the example, dummy coding is used. If we want to compare a1:b1 with a2:b2, we have to find out, how both of these combinations are represented. a1 and b2 are the reference levels so that a1:b2 is just the intercept. a2:b2 is represented by the intercept + the two main effects + the interaction. Accordingly, we can write:

(hyp1 <- hypothesis(fit, "Intercept + Aa2 + Bb2 + Aa2:Bb2 = Intercept"))
plot(hyp1)

Of course, the intercept cancels out, so we can remove it from the equation. To do multiple comparions, just pass the hypotheses as a charecter vector, e.g:

(hyp2 <- hypothesis(fit, c("Aa2 + Bb2 + Aa2:Bb2 = 0", "Intercept = 0"),
                    alpha = 0.025))
plot(hyp2)

It this what you had in mind?

[vmikk]

Thanks a lot for the clear answer! If I understand it correctly, I need to check the contrast matrix to find out how factor combinations are represented in the model. However, is there a simple way to extract contrast matrix from the brmsfit-object?

I tried with:

brms:::get_model_matrix(formula = fit$formula, data = fit$data)

but for binomial model it returned Error in model.matrix.default(terms, data) : model frame and formula mismatch in model.matrix(), and it seems that this function works only for models without random effects.

If the contrast coding is same in other packages and the default contrasts ("contr.treatment") are used, may I refit the model with glm or glmer and just to take the model matrix from the output?

library(lme4)
gfit <- glmer(cbind(N, Tot - N) ~ A*B + (1|ID), family = binomial, data = datt)
comb <- paste0(datt$A, datt$B)
comb <- aggregate(model.matrix(gfit) ~ comb, FUN=mean)
comb

[paul-buerkner]

The coding is the same as in other packages. If you really need to get the model.matrix to understand what contrasts you applied, you can can write:

fixed_form <- brms:::extract_effects(fit$formula)$fixed
mm <- brms:::get_model_matrix(fixed_form, data = fit$data) 

Be aware though, that internal functions may change without notice in future updates (for those two it's unlikely so it should be ok to use them).

[vmikk]

That's perfect! It is exactly what I need. Thanks!

[paul-buerkner]

I am glad to hear that! :-)

[IoSonoMarco]

Hi Paul,

I'm trying to implement the functions in the present post but they seem unreliable now. Is there any new way to get a contrast matrix and to do contrast analysis in brms right now?

[paul-buerkner]

The ways brms supports contrasts of factors is the same as in most other R packages dealing with such models. As such, brms sticks to this approach and does not implement a new one.

[IoSonoMarco]

Ok this is clear, but I have problems to call the following functions and I can't find detailed references out there :(

fixed_form <- brms:::extract_effects(fit$formula)$fixed
mm <- brms:::get_model_matrix(fixed_form, data = fit$data) 

[paul-buerkner]

This is because they are internal functions which are not designed to be applied by users.

[SaoirseConnorDesai]

Dear Paul and brms-users, Apologies I realise this is an old post, but I wondered if it is possible to obtain bayes factors for the multiple comparisons once they are coded?

[vincentott]

I wondered if it is possible to obtain bayes factors for the multiple comparisons once they are coded?

I have also wondered, and maybe this requires building standalone models for post-hoc comparisons? I.e., we only select the two cells that we want to compare and add a binary predictor that encodes cell-membership for each observation. Then, we build two models with brms. One has the additional binary predictor. Lastly, we compare the two models against each other and this gives the bayes factor?

Basically, as described here: https://www.cairn.info/revue-l-annee-psychologique-2020-1-page-73.htm (see Machiavellism post-hoc tests)

Another answer might be in here: https://vasishth.github.io/bayescogsci/book/ But I have yet to read and find out.