How to do a follow-up comparisons/contrasts for a 3-way interaction with 2 numeric predictors?

Question

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!

Answer 1

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.

Answer 2

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.