Unpacking interactions in LME4 with repeated measures design

Question

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
```

Answer 1

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.