[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.
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.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
andg5
), 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.