I'm trying to visualize a three-way interaction from a rather complex linear mixed model in R (lmer function from the lme4 package; the model has a complex random-effects structure). The interaction consists of two continuous variables and one categorical variable (two experimental conditions).
So far, I have graphed the interaction via two 3D-surface plots using visreg2d from the visreg package (see image). But my reviewers found these plots confusing and asked for a different illustration, such as conditional coefficient plots (i.e., plots of the strength of coefficient 1 as coefficient 2 increases).
I've tried to find a package that allows me to create these kind of plots, but failed. The existing packages only allow coefficient plots for two-way interactions (for instance the interplot package; https://cran.r-project.org/web/packages/interplot/vignettes/interplot-vignette.html). That means I only get a conditional coefficient plot of the two-way interaction, collapsed across both levels of the categorical variable.
Is there a package for my case? If not, I probably have to manually extract fitted values from my model (e.g., using broom) and somehow plot them in ggplot2. But I don't really know how to do this, whether or not to take into account random effects (and how), etc. Any ideas would be much appreciated...