0

I have a linear mixed effect model with three-way interaction fitted by the following code:

m <- lmer(cog ~ PRS*poly(Age, 2, raw=T)*Score
             + gender + Edu + fam + factor(Time)
             + (1|family/DBID), 
             data = test_all, REML = F)

In this model, there is a three-way interaction between PRS, Score, and polynomial terms of age with two degrees (linear + quadratic). For this three way interaction, how can I obtain the marginal effect(slope) of one variable, conditional on the other variables? For example, what is the slope of PRS, when age = 50 and score = 1?

Second, I tried to use the following code to plot this three-way interaction:

plot <- ggpredict(m, ci.lvl=0.95, c("PRS [all]", "Age [60, 65, 70, 75, 80]", "Score[0, 0.321, 0.695, 1.492, 1.914, 3.252]"))
plot(m)

The interaction plot finally shows but R didn't give the confidence interval. The error message is Error: Confidence intervals could not be computed.

How can I plot this three-way interaction with confidence interval?

zjppdozen
  • 63
  • 5
  • You could improve your chances of finding help here by adding a [minimal reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/5963610#5963610). Adding a MRE and an example of the desired output (in code form, not tables and pictures) makes it much easier for others to find and test an answer to your question. That way you can help others to help you! P.S. Here is [a good overview on how to ask a good question](https://stackoverflow.com/help/how-to-ask) – dario Oct 22 '21 at 05:10
  • Please provide enough code so others can better understand or reproduce the problem. – Community Oct 22 '21 at 07:22
  • Have you tried `sjPlot::plot_model(m, type="int")` ? It uses ggplot2 and geom_ribbon to show a confidence band and I find this function very useful to visualize interactions. Afaik, it also works for 3-way interactions and for `lmer`-models but I may have never tried combining both. – benimwolfspelz Oct 22 '21 at 09:57

1 Answers1

0

You can use the marginaleffects package to do that (disclaimer: I am the maintainer):

library(marginaleffects)
library(lme4)

mod <- lmer(mpg ~ hp * am * vs + (1 | cyl), data = mtcars)

mfx <- marginaleffects(mod, newdata = typical(vs = 0, am = 1, cyl = 4))

summary(mfx)

## Average marginal effects 
##       type Term   Effect Std. Error z value  Pr(>|z|)    2.5 %   97.5 %
## 1 response   am  4.10167    2.13391  1.9221 0.0545884 -0.08072  8.28406
## 2 response   hp -0.03724    0.01378 -2.7016 0.0069001 -0.06426 -0.01022
## 3 response   vs -0.61237    3.52755 -0.1736 0.8621833 -7.52625  6.30151
## 
## Model type:  lmerMod 
## Prediction type:  response
Vincent
  • 15,809
  • 7
  • 37
  • 39