# Emmeans and contrasts after fitting a linear mixed effects model (repeated measures) - A complex model

I am facing a really complex model and tried several models and post-hoc tests -with a great help from StackExchange- and would really appreciate your opinion. I need to fit a linear mixed effects model as I have repeated measures. Specifically, I have several biomarkers (continuous) being measured twice (Timepoint is: T1 and T2). I need to see the differences in each biomarker both within each timepoint as well as across timepoints. These timepoints-related comparisons need to be done between and across each sex group (female, male) and another group (group1, group2). Therefore, I fit the following model:

model <- lme(bio1 ~ Sex * Group * Timepoint, data=bio12, method="ML", random = ~ 1 | Subject);

summary(model);

Then, I'm doing the post-hoc tests using emmeans and contrasts:

emms.bio1.timepoint <- emmeans(bio1, ~ Group*Sex | Timepoint);

tukey.bio1.timepoint <- contrast(emms.bio1.timepoint, method="pairwise");

summary(tukey.bio1.timepoint);

Timepoint = T1:

contrast estimate SE df t.ratio p.value

G1 1 - G2 1 1.365 2.13 407 0.640 0.9189

G1 1 - G1 2 -12.651 2.07 407 -6.118 <.0001

G1 1 - G2 2 -13.004 2.03 407 -6.391 <.0001

G2 1 - G1 2 -14.016 2.05 407 -6.837 <.0001

G2 1 - G2 2 -14.370 2.02 407 -7.125 <.0001

G1 2 - G2 2 -0.353 1.95 407 -0.181 0.9979

Timepoint = T2:

contrast estimate SE df t.ratio p.value

G1 1 - G2 1 -0.258 2.13 407 -0.121 0.9994

G1 1 - G1 2 -12.946 2.07 407 -6.261 <.0001

G1 1 - G2 2 -15.356 2.03 407 -7.547 <.0001

G2 1 - G1 2 -12.688 2.05 407 -6.189 <.0001

G2 1 - G2 2 -15.097 2.02 407 -7.486 <.0001

G1 2 - G2 2 -2.410 1.95 407 -1.236 0.6040

*1 = male, 2 = female

Then I repeat for the other pairwise comparisons I want to see:

emms.bio1.sex <- emmeans(bio1, ~ Group*Timepoint| Sex);

tukey.bio1.sex <- contrast(emms.bio1.sex, method="pairwise");

summary(tukey.bio1.sex);

emms.bio1.group<- emmeans(bio1, ~ Sex*Timepoint| Group);

tukey.bio1.group<- contrast(emms.bio1.group, method="pairwise");

summary(tukey.bio1.group);

So my questions are:

1. Do you think this is the right approach?
2. If yes, are the comparisons from the contrasts I get reasonable? Since my model is for repeated measures, is it reasonable that I get the within-timepoint comparisons?

Thank you very much in advance and any help would be greatly appreciated.

This could be the right approach, but I suggest doing some model diagnostics and some exploration before plunging into post hoc tests. Yeah, I know, looking at plots, etc. biases statistical tests; but not looking and just turning a crank can be dangerous. You need to know if the model fits the data. And the decision as to what post hoc comparisons are needed depends on what interactions are important.

So what I'd do is to at least plot residuals versus predicted, and make sure you don't have a "horn of plenty" pattern or some such. If it looks OK, I suggest doing car::Anova(model) and perhaps consider re-fitting the model with the really less significant interactions (if any) excluded.

For post hoc, I also suggest re-fitting with method = "REML". THE ML estimates are good for model comparisons, but REML estimates are less biased and better for follow-up testing.

Get an idea, graphically, of what you have overall; for example, via

library(emmeans)
emmip(model, Sex ~ Timepoint | Group)


Are there lots of parallel trends, or do things go all over the place? It matters which. Usually, you do not need separate sets of means if you can marginalize (average over one or more factors) if you're not averaging over substantial interaction effects (non parallel trends in the plot).

Finally, you have a lot more code than needed, even for those results. Only need one emmeans() call is needed; look at the results of all the emm... objects you created and notice that all have the same estimates. you could get your tukey.bio.sex results, for example, via contrast(emms.bio1.timepoint, method="pairwise", by="sex"). On the other hand, if there are no big interactions between sex and the other two factors, it is not necessary to separate by sex so you could marginalize via EMM.TG <- emmeans(emms.bio1.timepoint, ~Timepoint*Group) and contrast(EMM.TG, "pairwise"). (That's just a for-instance; obviously, I don't know which factors interact.)

I hope this gives the flavor of how I'd approach it. Don't just follow a prescription; look at what you have and what provides the most meaningful descriptions. Try to simplify where complexity is unneeded.

• Hi, I really liked your answer, regarding the emmeans, I read the vignettes but it is still confusing for me, is this function estimating the significance in the mean between two interactions? Apr 2 at 19:53