How to do the post hoc test in the linear mixed model if I have three predictors (two factor variables and one numeric variable)

Question

I'm using a linear mixed effects model to analyze the reaction time of learners of English as a second language. I have two factor variables - grammaticality (grammatical v.s. ungrammatical) and Finiteness (finite v.s. nonfinite); as well as one numeric variable - the proficiency score.

I'd like to analyze how their proficiency may affect their sensitivity to grammaticality and finiteness, so I built the model like this:

fitnonnative7 <- lmer(data = nonnative,
                      formula = RT ~ Grammaticality * Finiteness * Proficiency + (Grammaticality|subject) + (Grammaticality|item))

I got a significant interaction in grammaticality:finiteness:proficiency; otherwise, no other main effects or interactions were found. However, what am I supposed to do next? Which functions/code may help to understand what is happening inside these three variables?

Answer 1

You can use the marginaleffects package to estimate contrasts. (Disclaimer: I am the author.) See the contrasts vignette for a lot of details:

https://vincentarelbundock.github.io/marginaleffects/articles/contrasts.html

In this example, am and vs are binary, and wt is numeric/continuous. First, we use the comparisons() function to see what happens to the predicted outcome when vs or am move from 0 to 1, for different values of the wt variable (using the fivenum function to set the variable at Tukey’s five numbers). Then, we plot the value of the contrasts in am and vs for different values of wt using the plot_cco() function:

library(lme4)
library(marginaleffects)

mod <- lmer(mpg ~ am * vs * wt + (am | gear) + (vs | gear), data = mtcars)
#> boundary (singular) fit: see help('isSingular')

comparisons(
    mod,
    variables = c("am", "vs"),
    interaction = FALSE,
    newdata = datagrid(wt = fivenum))
#>    rowid     type term contrast comparison std.error  statistic    p.value    conf.low conf.high predicted_lo predicted_hi      am     vs gear     wt
#> 1      1 response   am    1 - 0   7.353206  3.291220  2.2341883 0.02547069   0.9025326 13.803879    23.743324  31.09652941 0.40625 0.4375    3 1.5130
#> 2      2 response   am    1 - 0   2.266254  1.678246  1.3503705 0.17689718  -1.0230479  5.555555    20.649991  22.91624453 0.40625 0.4375    3 2.5425
#> 3      3 response   am    1 - 0  -1.600225  1.765295 -0.9064916 0.36467572  -5.0601390  1.859689    18.298817  16.69859236 0.40625 0.4375    3 3.3250
#> 4      4 response   am    1 - 0  -3.206111  2.205327 -1.4538031 0.14600091  -7.5284712  1.116250    17.322291  14.11618091 0.40625 0.4375    3 3.6500
#> 5      5 response   am    1 - 0 -11.971776  5.584590 -2.1437160 0.03205565 -22.9173722 -1.026180    11.991963   0.02018736 0.40625 0.4375    3 5.4240
#> 6      1 response   vs    1 - 0   3.845119  3.632085  1.0586531 0.28975778  -3.2736381 10.963875    25.048324  28.89344284 0.40625 0.4375    3 1.5130
#> 7      2 response   vs    1 - 0   2.999110  1.827767  1.6408606 0.10082636  -0.5832466  6.581467    20.258546  23.25765591 0.40625 0.4375    3 2.5425
#> 8      3 response   vs    1 - 0   2.356078  1.771024  1.3303481 0.18340360  -1.1150648  5.827221    16.617942  18.97401990 0.40625 0.4375    3 3.3250
#> 9      4 response   vs    1 - 0   2.089004  2.213747  0.9436510 0.34534803  -2.2498596  6.427868    15.105870  17.19487395 0.40625 0.4375    3 3.6500
#> 10     5 response   vs    1 - 0   0.631191  5.801561  0.1087968 0.91336370 -10.7396589 12.002041     6.852283   7.48347423 0.40625 0.4375    3 5.4240

plot_cco(mod, effect = "am", condition = "wt")

plot_cco(mod, effect = "vs", condition = "wt")