I am having a hard time doing follow-up (post-hoc) comparisons for a 3-way interaction of 2 numeric predictors and one factor in a linear mixed model. It’s a longitudinal study where we have patients which either receive no treatment, treatment A, or treatment B (factor “group”). We measure patients at baseline, after 1 and 3 weeks (numeric predictor “time”). We have an MRI brain measure (numeric predictor “biomarker”), for which we assume it influences the efficacy of the treatment. Finally, the outcome measure is symptom severity (numeric variable “severity”).

Let’s simulate data for 60 patients (20 for each group) and run a linear mixed model:

# Simulate data 
participant <- rep(1:60, each = 3)
group <- rep(rep(0:2, each = 3), 20)
time <- rep(c(0, 1, 3), 60)
biomarker <- rnorm(180)
severity <- rnorm(180)

# Form a dataframe
df <- data.frame(participant, group, time, 
    biomarker, severity)

# Coerce factors
factors <- c("participant", "group")
df[, factors] <- lapply(df[, factors] , factor)

# Linear mixed model

model <- lmer(severity ~ group * time * 
    biomarker + (1|participant), df)

If we run anova(model) on the actual data, we get a significant main effect of time (patients generally get better), a 2-way interaction of time and group (treatments work), and a 3-way interaction of time, group, and biomarker (treatment efficacy depends on biomarker).

Now I would like to compare the interaction strength of biomarker and time between the 3 groups.

I tried:

# Group comparisons
emt <- emtrends(model, ~ biomarker*group, 
                var = "time")

I think the problem with this approach is that it only takes the average value for biomarker into account, effectively leaving its effect away and “only” comparing the effect of treatment over time. Unfortunately, emtrends() does not accept a second variable for the var-argument, which specifies the numeric predictors (emtrends(model, ~ biomarker*group, var = c("time", "biomarker")) throws an error.)

I also tried the following, but then I don’t know how to form the desired group comparisons:

emt <- emtrends(model, ~ biomarker*group, 
        at = list(biomarker = c(-1, 0, 1)), 
        var = "time")  
    # biomarker is standardized

If I plot the results of the actual data (3-way interaction for groups), I can clearly see that the biomarker only influences the outcomes for one treatment, which is perfectly in line with our hypotheses! Any help on how to get group (factor) comparisons for the interaction strength between the 2 numeric predictors would be highly appreciated!

  • $\begingroup$ This doesn't appear to be a specific programming question that's appropriate for Stack Overflow. If you have general questions about the appropriate use of various statistical methods, then you should ask such questions over at Cross Validated instead. You are more likely to get better answers there. $\endgroup$
    – MrFlick
    Mar 25, 2022 at 21:20
  • 1
    $\begingroup$ In your last emt creation, specify ~ group | biomarker,and add var = "time". This will give you separate slopes of time trends for each group and biomarker value. $\endgroup$
    – Russ Lenth
    Mar 26, 2022 at 3:24
  • 1
    $\begingroup$ Then pairs(emt) and pairs(emt, by = "group") to compare slopes by biomarker and by group $\endgroup$
    – Russ Lenth
    Mar 26, 2022 at 3:33
  • $\begingroup$ Thanks for your help @RussLenth. I tried your suggestions and but they did not really give me the comparison I was looking for. I am interested in comparing the strength of the interaction between time and biomarker between the different treatment groups. I have provided a solution I came up with below, and would be more than happy to have some feeback. $\endgroup$ Mar 29, 2022 at 18:45

2 Answers 2


The approach in @BulkySplash's answer can be done without estimating trends, because it is equivalent to looking at certain interaction contrasts among the three factors. Consider:

> EMM = emmeans(model, ~time*biomarker*group, 
+     at = list(time = c(-1,1), biomarker = c(-1,1))
+ )

This creates a grid of 2 x 2 x 3 = 12 means. Now, create interaction contrasts of these means; I'll opt for polynomial contrasts:

> contrast(EMM, interaction = "poly")
 time_poly biomarker_poly group_poly estimate    SE  df t.ratio p.value
 linear    linear         linear        0.959 0.654 168   1.466  0.1445
 linear    linear         quadratic    -2.703 1.042 167  -2.593  0.0104

Degrees-of-freedom method: kenward-roger

The first result is exactly the same as contrast(emt, custom_contrast). Because the linear contrast for group has coefficients (-1, 0, +1), it is a comparison of the linear x linear time:biomarker interaction contrasts between the 1st and 3rd group. Interestingly, the other contrast, which has group coefficients (-1, +2, -1), is more significant.

To understand all this better, it is helpful to look at the separate linear x linear contrasts for each group:

> contrast(EMM, interaction = "poly", by = "group")
group = 0:
 time_poly biomarker_poly estimate    SE  df t.ratio p.value
 linear    linear           -0.734 0.443 166  -1.658  0.0992

group = 1:
 time_poly biomarker_poly estimate    SE  df t.ratio p.value
 linear    linear            1.097 0.406 166   2.702  0.0076

group = 2:
 time_poly biomarker_poly estimate    SE  df t.ratio p.value
 linear    linear            0.225 0.481 168   0.467  0.6410

Degrees-of-freedom method: kenward-roger 

Here we see that the interaction is strongest in group 1.

  • $\begingroup$ Thanks so much for your clarifications and response, Russ! contrast(EMM, interaction = "poly", by = "group") is very interesting and insightful. I am wondering if some family-wise p-value correction for these 3 tests would be adequate or maybe is already implemented. I tried .contrast(EMM, interaction = "poly", by = "group", adjust = "tukey") and contrast(EMM, interaction = "poly", by = "group", adjust = "none"), but that seemed to make no difference? $\endgroup$ Mar 31, 2022 at 19:14
  • $\begingroup$ Yeah, it's annoying but the by variable gets in the way if you want to adjust for everything. What you have is NOT adjusted no matter what because there is only one test in each group. If you want to adjust the family of 3 tests, you have to do something like update(contrast(..., by = "group"), by = NULL, adjust = "bonf"). But a Tukey adjustment will not be allowed because you can use that only for a single complete set of pairwise comparisons, as required by its theoretical underpinnings. $\endgroup$
    – Russ Lenth
    Mar 31, 2022 at 20:33
  • $\begingroup$ I see. Good to know that by-variables are preventing the correction. I asked because I think reporting the interaction in each group separately might be a more insightful follow-up test than group comparisons... Anyways, thank you very much, Russ Lenth, for your time to answer here, but also - unless I am mistaken - for the creation of the emmeans package!!! The vignettes have been educating me. :-) $\endgroup$ Apr 1, 2022 at 12:30

I am the poster of the question and here is the solution I found so far.

I am interested in comparing the interaction strength of biomarker and time between groups. In other words, I want to know if the influence of the biomarker on symptom severity is different for different groups. More specifically, I assume that the biomarker has an influence on treatment efficacy for treatment B, but not for treatment A or placebo. I could not find a "standard" way of getting this comparison. So I used a custom contrast, based on the following emtrends object (equal to the last one in the question):

emt <- emtrends(model, ~ biomarker*group, 
        at = list(biomarker = c(-1, 0, 1)),  # biomarker is standardized
        var = "time")  

I derived the custom contrast of interest, which is an interaction of the group effect (treatment B vs placebo) and the effect of the biomarker. So, I coded the effect for group and biomarker and then formed the custom contrast as the product of these two.

Factor group Factor biomarker Effect group Effect biomarker Custom contrast
placebo -1 -1 -1 1
placebo 0 -1 0 0
placebo 1 -1 1 -1
treatment A -1 0 -1 0
treatment A 0 0 0 0
treatment A 1 0 1 0
treatment B -1 1 -1 -1
treatment B 0 1 0 0
treatment B 1 1 1 1

When emt is called, one can see the order of the possible condition combinations and it is in this order that contrast coefficients have to be provided. The order is equal to the one I've used in the table. Thus, I copied the contrast from the table and called the contrast() function on the emtrends object:

custom_contrast <- list(
   treatmentB_vs_placebo = c(1, 0, -1, 0, 0, 0, -1, 0, 1)

contrast(emt, custom_contrast)

This is just an example, but with the table, it is easy to form other contrasts (such as treatment A vs placebo, or treatment A vs treatment B).

I am a bit unsure about the time aspect, but the emtrends object provides a separate slope for time trend for each combination of factors, so I guess that deals with that. I am also unsure about the coding of the biomarker, but since it is standardized to a z-distribution (mean 0, standard deviation 1), I guess the values -1, 0, and 1 should be adequate and enough to recognize interactions with other factors. I'd be more than happy to have some input on these 2 points.

However, these issues aside, I believe this custom contrast is correct and providing me with the result I was looking for. I also think this contrast is a much more specific test for my hypothesis and I am considering to edit the manuscript to only contain contrast analyses, without reporting the originally posted 3-way interaction.

  • $\begingroup$ Note: Your custom contrast is actually an interaction contrast, and you can get the same result via contrast(emt[c(1,3,7,9)], interaction = "consec"). $\endgroup$
    – Russ Lenth
    Mar 30, 2022 at 16:07

Your Answer

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

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