Marginal likelihood for linear model with random effects to do Bayesian model comparison

Suppose I have behavioral data from multiple participants to four different conditions (four observations per participant per condition). The conditions can be characterized in terms of two fixed effects, x1 and x2 (i.e., the experiment is a 2-by-2 design). That is to say, condition 1 has properties such that x1 is 0 and x2 is 0, condition 2 has properties such that x1 is 1 and x2 is 0, condition 3 has properties such that x1 is 0 and x2 is 1, and condition 4 has properties such that x1 is 1 and x2 is 1.

Here is some sample data for two participants:

data <- data.frame(
Participant = c(rep("1", 16), rep("2", 16)),
ItemNumber = as.character(rep(1:16, 2)),
x1 = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
x2 = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,
0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
y = c(0.604638544849808, -0.550491809490124, 0.604638544849808,
1.18220372201977, 0.0270733676798422, 0.604638544849808,
1.75976889918974, 0.604638544849808, 0.0270733676798422,
0.604638544849808, 0.604638544849808, -0.550491809490124,
0.0270733676798422, 0.604638544849808, 1.18220372201977,
0.604638544849808, -0.341514434548145, 1.29775485128295,
0.204908660728887, -0.341514434548145, -0.887937529825177,
0.751331756005919, 1.29775485128295, -0.341514434548145,
0.204908660728887, 0.751331756005919, 0.204908660728887,
0.751331756005919, -0.887937529825177, -1.43436062510221,
-1.43436062510221, 1.29775485128295))


I would like to determine, for each participant, whether their data was more likely to have been generated by model 1 vs. model 2, where model 1 only has fixed effects for x1 and x2, but model 2 additionally has a fixed effect for the interaction between x1 and x2.

I was thinking I could do this using Bayesian model comparison. That is, assuming model 1 and model 2 are equally likely, I could calculate the marginal likelihoods, $$P(X_i\mid m_1)$$ and $$P(X_i\mid m_2)$$ for each participant $$i$$, and then just divide one by the other to get a Bayes Factor for each participant.

I have two questions:

1. Is this a reasonable approach to determining whether a given participant's data was more likely to have been generated by model 1 vs. model 2?

2. If so, how can I go about doing this?

Assuming this is a reasonable approach, I'm not sure what to do to calculate the marginal likelihoods. Since the observations are not IID (i.e., there are multiple observations per participant per condition), each model should have a random effect for participant and for item number. Let's say that there are only random intercepts (one for participant and one for item). Then model 1 would look like this:

$$y \sim \mathcal{N}(\alpha_{\text{participant}} + \alpha_{\text{item}} + \beta_0 + \beta_1x_1 + \beta_2x_2, \sigma^2)$$

and model 2 would look like this:

$$y \sim \mathcal{N}(\alpha_{\text{participant}} + \alpha_{\text{item}} + \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2, \sigma^2)$$

If there weren't random effects, I think I could just analytically calculate the following integral using the normal-gamma distribution as the conjugate prior to the normal distribution with unknown mean and variance:

$$P(X_i\mid m ) = \int_{\theta} P(X_i\mid \theta)P(\theta\mid m) d\theta$$

However, since the observations aren't IID and random effects are needed in the two models, I'm not sure how to calculate the marginal likelihoods. Is there an analytical solution available, or is this something that must instead be estimated? If it must be estimated, how can I go about doing this?

Yes, this is a reasonable thing to do. As for actually computing the marginal likelihoods, the easiest thing would probably be to fit your model in brms/stan, and then use the bridgesampling package to compute a bayes factor.

For your example, it might look something like this:

library(bridgesampling)
library(brms)

m1 <- brm(y ~ (1 | participant) + (1 | item) + x1 + x2, data=data,
family=gaussian(), save_pars=save_pars(all=True))
m2 <- brm(y ~ (1 | participant) + (1 | item) + x1 + x2 + x1 * x2, data=data,
family=gaussian(), save_pars=save_pars(all=True))

bayes_factor(m1, m2)


The bridgesampling algorithm uses an iterative procedure to approximate the marginal likelihood, via samples from the posterior. If you want to learn more about the algorithm, a good reference can be found here.