Use normal_id_glm_lpdf in brms for fixed + random effects model?

Interfaces brms
3

The random effects can go in the intercept term I think. I think in some models we talk about these as being intercept terms. Like group level intercepts, etc.

So originally we had:

target += normal_lpdf(Y | mu, sigma);

But mu was fixed effects + random effects like here:

// initialize linear predictor term
vector[N] mu = Intercept + Xc * b;
for (n in 1:N) {
  // add more terms to the linear predictor
  mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_1_2[J_1[n]] * Z_1_2[n];
}

Now just do:

// initialize linear predictor term
vector[N] random_effects = rep_vector(Intercept, N); // Removed Xc * b
for (n in 1:N) {
  // add more terms to the linear predictor
  random_effects[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_1_2[J_1[n]] * Z_1_2[n];
}

And call:

target += normal_id_glm_lpdf(Y | Xc, random_effects, b, sigma);

And I played around with this (Edit: I didn’t play around with the glm thing, just to be clear, I’m adding a separator above to indicate this is a different topic now) a while back but I think there’s a way we might be able to make the random effects calculation more efficient.

It’s kinda weird to talk about it being more efficient cause it’s already an incredibly fast thing (just looping through and doing a little simple addition), but check here: Group orderings in regression

The motivation there is that, maybe we have 100000 people. That means we have to compute 100000 random effects. But what if the difference in most of all these people is just one thing? So we could take the nth person’s value and then just adjust it to compute the (n + 1)th person’s value and that’s cheaper? Maybe there’s a way we can reorder the calculations to do a lot less math?