How to calculate marginal effects in mixed models

library(lme4)
m1 <- lmer(resp ~ pred1 + pred2 + pred3 + (1|site), data = Xs)
mod.predictions1 <- predict(m1, re.form = NA)
plot(pred1, mod.predictions1) # plot 1

Xs$$pred2 <- mean(Xs$$pred2)
Xs$$pred3 <- mean(Xs$$pred3)
mod.predictions2 <- predict(m1, newdata = Xs, re.form = NA)
plot(pred1, mod.predictions2) # plot 2


Between plot1 and plot2, what is the difference? If I want to show how does my resp change with pred1, should I show plot1 or plot2?

If you want to show how your outcome resp is related to pred1, you should go for the second approach, because in the first one the predictions you obtained are affected by changing values of pred2 and pred3.

In general, you can something like this

plot_data <- with(Xs, expand.grid(
pred1 = seq(min(pred1), max(pred1), length = 100),
pred2 = mean(pred2),
pred3 = mean(pred3)
))

preds <- predict(m1, newdata = plot_data, re.form = NA)

• Since re.form = NA, that means that the random effect is set to zero in the fiited model when using the predict() function, so we are dealing with a "typical" (or "average") site. Then the predict() function can be used as explained by Dimitris for the second approach to determine, for example, how the mean value of resp changes as a function of pred1 for a "typical" site for which pred2 and pred3 are known to have average values. – Isabella Ghement Oct 8 '18 at 15:22

You could use the ggeffects-package to compute (and plot) marginal effects:

library(lme4)
library(ggeffects)
m1 <- lmer(resp ~ pred1 + pred2 + pred3 + (1|site), data = Xs)
pr <- ggpredict(resp, "pred1")
pr
plot(pr)


There's a vignette showing the basics (here) and one that shows the specialties of mixed model (here). You can find an overview of all vignettes here.