3-Way Interaction – How to Do Follow-Up Comparisons/Contrasts

interactionlme4-nlmelsmeansr

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
library(lme4)
library(lmerTest)
library(emmeans)

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

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")
pairs(emt)

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!

Best Answer

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.

Related Question