Maybe this can help:

https://cran.r-project.org/web/packages/sjPlot/vignettes/plot_interactions.html

I would suggest, have the data ready, Use ggplot. You can start by plotting for one group, and if you use the ggplot framework it is possible to increase complexity.

ggplot is very powerful, and flexible. It is well worth getting used to using it.

You often end up with a number of plots in a panel which can show a 3 way interaction. The time is best on the x axis ofcourse.

Please see this link. It describes an R package that may be helpful. Also refer to D. C. Montgomery, Design of Experiments, Best, David Booth

https://www.amazon.com/Analysis-Experiments-Chapman-Statistical-Science-ebook/dp/B00UVBRX2G/ref=sr_1_1?keywords=lawson+design+experiments+r&qid=1573787208&s=books&sr=1-1

MLwiN has excellent facilities for doing this sort of plot for mixed models; it does it by simulation from the fitted model and you can place confidence intervals around the fitted lines. It can do this for population average and cluster specific predictions which can be different if nonlinear models are fitted. See the examples under customised predictions in this manual:

Article Manual supplement for MLwiN Version 2.14

It is also possible to access MLwiN's specialist facilities from within R

R2MLwiN: Running MLwiN from within R

http://www.bristol.ac.uk/cmm/software/r2mlwin/

You need both software packages on your machine.

Since sjPlot uses the ggeffects-package for marginal effects plots, you can also use ggeffects directly (documentation: https://strengejacke.github.io/ggeffects/). There is an example for mixed models (https://strengejacke.github.io/ggeffects/articles/practical_logisticmixedmodel.html) as well as for two- to four-way interactions (https://strengejacke.github.io/ggeffects/articles/ggeffects.html#two-way-three-way--and-four-way-interactions).

Thank you all for your answers!

I solved it by using get_model_data() from the sjPlot-Package to get the marginal means in table-form and then used ggplot() to draw the figure.

Maybe there are better ways to get these estimated values, but it worked!

Best,

Daniel

There are some examples here (I suggest using "ggpredict()" instead of "get_model_data"): https://strengejacke.github.io/ggeffects/articles/ggeffects.html

You can also use "plot()", which returns a ggplot-object for further modification.