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.
# 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!