I have a category variable called
group, which has five levels indicating the types of populations (
g5). I wanted to know the effect of
group on the outcome variable
score, while other factors such as
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
df$group <- relevel(df$group, ref='g1') mod <- lmer(score ~ group + gender + age + IQ + (1|site), data=df)
Although the model results showed me how the effects of
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
g1as the reference level, so what I thought is setting other reference levels (i.e.,
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.