[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.


1 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.

  • $\begingroup$ Thank you for your help. Yes, I meant “site” when I mentioned the variables. I have edited the question to fix this. I included variables like gender and age in the model because I hypothesised that the influence of them are existing. But I only want to know the group differences, and this is why I called "gender", "age", "IQ" and "site" as covariates. Is there a way to compare group differences when considering the existence of covariates (control variables) in emmeans? Many thanks. $\endgroup$
    – ziqian_wei
    Commented Jun 26, 2023 at 4:17
  • $\begingroup$ Yes. That's the purpose of the emmeans package. It sounds like you've done it, but don't understand what you have. Read those vignettes. $\endgroup$
    – Russ Lenth
    Commented Jun 26, 2023 at 15:26
  • $\begingroup$ Thanks. I have read the quick start and checked my results. As estimated coefficients and standard errors were the same in both model and post-hoc comparison results, I was sure these results were correct (Only p-values were different because a correction method was used in emmeans). Thank you again for your help! $\endgroup$
    – ziqian_wei
    Commented Jun 27, 2023 at 4:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.