Multiple group comparisons in a linear mixed-effect model

Question

[The background]

I have a category variable called group, which has five levels indicating the types of populations (g1, g2, g3, g4 and g5). I wanted to know the effect of group on the outcome variable score, while other factors such as gender, age, IQ and site were included as covariates (i.e., the influence of these variables are existing, but it is not what I focus on). And I hypothesised site works as a random factor.

As far as I know, Linear Mixed-effect Models (LMMs) should be a suitable approach in this case. I constructed a mixed-effect model in R and I defined g1 as the reference level for the variable group.

df$group <- relevel(df$group, ref='g1')
mod <- lmer(score ~ group + gender + age + IQ + (1|site), data=df)

[The question]

Although the model results showed me how the effects of g2 to g5 on score differ from g1 (i.e., the reference level), I wanted to further know the differences between other levels (e.g., "g2 vs. g3"). I have tried the two approaches as follows.

1. Conducting post-hoc comparisons using emmeans(mod, pairwise~group). But this raised two issues. First, it did not include covariate variables (e.g., age). Second, it utilised correction methods to handle the multiple comparison problems. Therefore, results such as "g1 vs. g2" were different between the model results and the post-hoc comparisons. So it was a bit confusing when interpreting the findings.

2. Setting different reference levels and running many models. I have already used g1 as the reference level, so what I thought is setting other reference levels (i.e., g2, g3, g4 and g5), running new models, and combining all the results into a table for reporting. However, I was not sure if this is appropriate and I did not find out relevant statistical papers/textbooks to support my analysis plan. Also, the table would be very large as I have five levels and 10 comparisons (from "g1 vs. g2" to "g4 vs. g5").

How should I analyse group differences in this case? Any comments are appreciated.

Answer 1

You did not get gender- or age-specific estimates because you didn't ask for them. The call emmeans(mod, pairwise ~ group asks for marginal means for each group, and pairwise comparisons thereof. Those group means are marginal means averaged over gender and age (unless age is numeric, in which case the average age is used). And region is excluded because it is not used as a predictor in your model. Or did you mean to say site?

The emmeans package has a whole lot of vignettes with discussions and examples. I suggest you start by looking at the quick start one and the basics one.