Unpacking interactions in LME4 with repeated measures design

I am conducting a study where I look at the interaction of 3 categorical variables and 1 continuous variable. However, I want to be able to see all the possible comparisons of these 4 variables. In the past, I have used emmeans but I noticed that emmeans only takes the lowest and highest value of the continuous variable which does not make sense in repeated measures since it basically takes the lowest participant compared to the highest participant.


Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: RT ~ Domain * ShiftType * TrialType * VA_k + (1 | Probe) + (1 |      Story_order) + (1 | Subject)
Data: bsmu[bsmu\$ACC == 1, ]

REML criterion at convergence: 79528.9

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.2523 -0.4384 -0.1779  0.1482 12.8338

Random effects:
Groups      Name        Variance Std.Dev.
Probe       (Intercept)  169631   411.9
Subject     (Intercept)  545749   738.7
Story_order (Intercept)  101769   319.0
Residual                3042405  1744.2
Number of obs: 4472, groups:  Probe, 380; Subject, 60; Story_order, 8

Fixed effects:
Estimate Std. Error      df t value Pr(>|t|)
(Intercept)                                      2254.31     228.51   80.13   9.865 1.74e-15 ***
DomainL2                                          268.94     188.33 4375.18   1.428   0.1534
ShiftTypeNo Shift                                 -30.19     200.03 1992.88  -0.151   0.8800
ShiftTypeUnchanged                                244.29     201.03 1959.52   1.215   0.2244
TrialTypeSpace                                    482.29     206.74 2100.11   2.333   0.0198 *
VA_k                                              150.41     105.24  170.10   1.429   0.1548
DomainL2:ShiftTypeNo Shift                       -132.44     264.37 4376.51  -0.501   0.6164
DomainL2:ShiftTypeUnchanged                      -193.74     265.03 4372.25  -0.731   0.4648
DomainL2:TrialTypeSpace                          -300.61     276.90 4372.10  -1.086   0.2777
ShiftTypeNo Shift:TrialTypeSpace                  -78.89     289.34 2148.38  -0.273   0.7851
ShiftTypeUnchanged:TrialTypeSpace                -444.28     289.46 2117.47  -1.535   0.1250
DomainL2:VA_k                                     -97.90     101.58 4222.47  -0.964   0.3352
ShiftTypeNo Shift:VA_k                             56.42     101.62 4242.86   0.555   0.5788
ShiftTypeUnchanged:VA_k                           -82.56     100.25 4236.84  -0.824   0.4103
TrialTypeSpace:VA_k                                39.49     104.60 4268.38   0.377   0.7058
DomainL2:ShiftTypeNo Shift:TrialTypeSpace         198.60     385.45 4374.32   0.515   0.6064
DomainL2:ShiftTypeUnchanged:TrialTypeSpace        446.45     385.71 4379.42   1.157   0.2471
DomainL2:ShiftTypeNo Shift:VA_k                    53.32     142.70 4224.51   0.374   0.7087
DomainL2:ShiftTypeUnchanged:VA_k                  292.02     142.67 4209.58   2.047   0.0407 *
DomainL2:TrialTypeSpace:VA_k                      164.48     148.53 4217.14   1.107   0.2682
ShiftTypeNo Shift:TrialTypeSpace:VA_k            -101.93     148.10 4262.76  -0.688   0.4913
ShiftTypeUnchanged:TrialTypeSpace:VA_k             55.99     145.95 4256.12   0.384   0.7013
DomainL2:ShiftTypeNo Shift:TrialTypeSpace:VA_k   -152.50     208.23 4232.78  -0.732   0.4640
DomainL2:ShiftTypeUnchanged:TrialTypeSpace:VA_k  -422.48     206.97 4219.86  -2.041   0.0413 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 24 > 12.
Use print(x, correlation=TRUE)  or
vcov(x)        if you need it
$$$$


I wasn't entirely sure about what you mean by all possible comparisons, but one thing you could do would be to compare the continuous predictor's slopes between different levels of your categorical predictors. For this, you can use emmeans emtrends:

(Va_k seems to be the continuous predictor. I'll use the Shifttype as the example categorical predictor):

library(emmeans)
em<-emtrends(model, "ShiftType", var="Va_k")
pairs(em)


Edited to add: for 3-way interactions you can similarly use

em2<-emtrends(model, ~ShiftType|TrialType, var="Va_k")
pairs(em2)
`

but I suspect there may be some better way to do the latter (better than just running all possible contrasts). For instance this and this response here might be useful.

• This is exactly what I needed! Thank you so much! Commented Apr 16, 2023 at 1:06
• Nice, I'm glad it helped! Commented Apr 16, 2023 at 9:40