# When should grouping variables interact in a mixed-effects model?

I was reading this post which is relevant to a research project I'm working on now. I think that I understand the difference between crossed and nested random effects, e.g. as described here.

The first post I linked asks about differences between the random effects structure Ticks~(1|Year)+(1|Site) vs. Ticks~(1|Year)+(1|Site)+(1|Year:Site), i.e. allowing the number of ticks at each Site to vary across years, in addition to varying across years and across sites. Since there are multiple observations of each site during each year, this model is estimable and more flexible than the model that omits the term (1|Year:Site).

While the excellent GLMM FAQ does not mention interactions between grouping variables, this is mentioned in another writing by Ben Bolker here, which also describes this term as allowing the effect of each year to vary by site.

I'm trying to understand more about this interaction between the grouping variables and when it might be necessary, and whether this differs between terms with a 1 on the LHS or a variable on the LHS of the grouping term. Assuming such a model is estimable, under what circumstances would one want to use the model y ~ (x|f) + (x|g) instead of y ~ (x|f) + (x|g) + (x|f:g)?

From a lot of what I've read, it seems that the former model is often recommended for fully crossed designs. Does this justify leaving out the extra term? Why?

Am I just overthinking this problem, and the inclusion of this term should be based on domain knowledge and the specific data one has?

• These models are nested within one another. Have you tried using a likelihood ratio test (anova() in R) to determine whether the added complexity of the (1|f:g) intercept fits the data better than the less complex model without it? If you have substantive reasons to prefer the more complex model, then so be it. But you can also use the LRT to help you if you don't have such reasons. Commented Dec 19, 2023 at 0:02
• @ErikRuzek Yes, and it appears to benefit the model. But as we know, just because an LRT or AIC says to include the term doesn't mean it's always best from an inferential standpoint, so I was curious if something about the study design should affect this decision. Commented Dec 19, 2023 at 17:13

I am going to work from the simpler model in which you do not have a random (or varying) slope for x and instead just have a set of random (or varying) intercepts:

m1 <- lmer(dv ~ 1 + (1|f) + (1|g) + (1|f:g), d)


The random (or varying) intercepts can be interpreted as such:

1. (1|f) - the random intercept for f is shared across all groups g for a given f.
2. (1|g) - the random intercept for g is shared across all groups f for a given g.
3. (1|f:g) - the random intercept for f:g is shared among the unique groups comprised by combinations of f and g.

Keep in mind that in order to validly estimate #3, you need to have many cases in each combination. If you only have one case per combination, then this if completely confounded with the residual term.

So how do you decide whether to include (1|f:g) in the model? Substantively, you can ask yourself whether you think the combinations of f and g are distinct enough from their fellow f and/or g members that they need their own unique variance term. One may say, "of course they are!" But the problem is that one might not have enough data to reliably estimate this additional random parameter. That could mean that you get a model convergence error, a weird parameter estimate, etc.

All that said, model fit criteria is a perfectly fine justification one can use to determine whether the added complexity of the additional interaction intercept improves the fit of the model relative to the simpler additive model. This answers the question of whether the more complex model does a better job of accounting for the unexplained variance in the outcome than the simpler model.