# How to visualize a multilevel/mixed model with two fixed effects and one random effect?

The model I want to visualize is this one:

BZ ~ 1 + SV_DMW + SV_DIFF (1 + SV_DMW| VPN)

Do you have any ideas how to do this best?

This is how I visualized the model, where just SV_DMW was included.

The single lines are 5 random chosen subjects (VPN).

Do you have any idea how to visualize the model shown above similarly, or do you maybe have entirely different recommendations how to visualize the relationship?

Further information: SV_DMW and BZ are, in general, associated positively SV_DIFF and BZ are, in general, associated negatively

• Welcome to the site. What are these variables? It helps to have context, but, in particular, are SV_Diff and SV_DMV continuous or categorical or integers or what? Are they always positive? Commented Jun 23, 2023 at 12:46
• It's easier to help you if you include a simple reproducible example with sample input and desired output that can be used to test and verify possible solutions. Commented Jun 25, 2023 at 13:07
• @PeterFlom thanks for your questions. All variables are continuous. SV_Diff has a range from 1 to 6 and is associated negatively with BZ. SV_DMW has a range from 1 to 7 and is associted positively with BZ. BZ has a range from 1 to 9. BZ = relationship satisfaction SV_DMW = Jointly averaged level of sexual desire of a couple SV_Diff = Difference of sexual desire of persons of a couple
– Bila
Commented Jun 25, 2023 at 13:16
• Instead of plotting individuals, plot the predicted values. Use two graphs. For graph 1, the x axis is SV_DIFF and the y axis is BZ. Make lines for different values of DMW (maybe the quartiles). Graph 2, just the reverse. Commented Jun 25, 2023 at 14:11
• Gelman's Data Analysis Using Regression and Multilevel/Hierarchical Models has many figures and visualizations wrt this class of models. Might be worth a look? google.com/books/edition/… Commented Jun 27, 2023 at 11:52

The sjPlot package can get you there. You can use the plot_model function.

library(lme4)
library(sjPlot)
library(ggplot2)
library(patchwork)
theme_set(theme_sjplot())

# Data creation inspired by Ben Bolker's post at
# https://stackoverflow.com/a/38296264/7050882

set.seed(101)
SV_DMW = runif(6000, 1, 24)
SV_DIFF = runif(6000, 5, 15) * -SV_DMW
VPN = as.factor(rep(seq(1:60), each = 100))
testdata <- data.frame(SV_DMW, SV_DIFF, VPN)
testdata$$BZ <- testdata$$SV_DIFF * -testdata\$SV_DMW/10 * rnorm(300, mean=30, sd=20)

# Model
fit = lmer(BZ ~ 1 + SV_DMW + SV_DIFF + (1|VPN),
data = testdata, REML = TRUE)
#> boundary (singular) fit: see help('isSingular')

plot_model(fit, type = "pred", terms = c("SV_DMW")) +
ggtitle(element_blank()) +
plot_model(fit, type = "pred", terms = c("SV_DIFF")) +
ggtitle(element_blank())
#> Warning: Ignoring unknown parameters: linewidth
#> Ignoring unknown parameters: linewidth


Created on 2023-06-27 with reprex v2.0.2

I simulated your data above. You can instead work directly on the plot_model part by replacing fit with the name of your model.