6

In the past, I had used the sjp.glmer from the package sjPlot to visualize the different slopes from a generalized mixed effects model. However, with the new package, I can't figure out how to plot the individual slopes, as in the figure for the probabilities of fixed effects by (random) group level, located here

Here is the code that, I think, should allow for the production of the figure. I just can't seem to get it in the new version of sjPlot.

library(lme4)
library(sjPlot)
data(efc)
# create binary response
efc$hi_qol = 0
efc$hi_qol[efc$quol_5 > mean(efc$quol_5,na.rm=T)] = 1
# prepare group variable
efc$grp = as.factor(efc$e15relat)
# data frame for 2nd fitted model
mydf <- na.omit(data.frame(hi_qol = as.factor(efc$hi_qol),
                           sex = as.factor(efc$c161sex),
                           c12hour = as.numeric(efc$c12hour),
                           neg_c_7 = as.numeric(efc$neg_c_7),
                           grp = efc$grp))
# fit 2nd model
fit2 <- glmer(hi_qol ~ sex + c12hour + neg_c_7 + (1|grp),
              data = mydf,
              family = binomial("logit"))

I have tried to graph the model using the following code.

plot_model(fit2,type="re")
plot_model(fit2,type="prob")
plot_model(fit2,type="eff")

I think that I may be missing a flag, but after reading through the documentation, I can't find out what that flag may be.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
user44796
  • 1,053
  • 1
  • 11
  • 20
  • just to clarify: you're talking about `type = "ri.prob"` from the linked blog post? And can you specify what you've tried & what's not working? – Ben Bolker Dec 19 '18 at 16:57
  • I am talking about the plot type = "ri.prob". I will edit my question to show what I have tried. – user44796 Dec 19 '18 at 17:02

2 Answers2

7

Looks like this might do what you want:

(pp <- plot_model(fit2,type="pred",
       terms=c("c12hour","grp"),pred.type="re"))
  • type="pred": plot predicted values
  • terms=c("c12hour", "grp"): include c12hour (as the x-axis variable) and grp in the predictions
  • pred.type="re": random effects

I haven't been able to get confidence-interval ribbons yet (tried ci.lvl=0.9, but no luck ...)

pp+facet_wrap(~group) comes closer to the plot shown in the linked blog post (each random-effects level gets its own facet ...)

enter image description here

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • 1
    I'm using your code here: http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#predictions-andor-confidence-or-prediction-intervals-on-predictions to compute conf. int. or pred. int. for mixed models, however, I'm not sure how this might work when plotting specific levels of random effects, to get the CI / PI for the random effects. – Daniel Dec 22 '18 at 14:35
3

Ben already posted the correct answer. sjPlot uses the ggeffects-package for marginal effects plot, so an alternative would be using ggeffects directly:

ggpredict(fit2, terms = c("c12hour", "grp"), type="re") %>% plot()

There's a new vignette describing how to get marginal effects for mixed models / random effects. However, confidence intervals are currently not available for this plot-type.

The type = "ri.prob" option in the linked blog-post did not adjust for covariates, that's why I first removed that option and later re-implemented it (correctly) in ggeffects / sjPlot. The confidence intervals shown in the linked blog-post are not correct, either. Once I figure out a way how to obtain CI or prediction intervals, I'll add this option as well.

Daniel
  • 7,252
  • 6
  • 26
  • 38
  • Is it possible to exp() the predicted values and the confidence intervals so instead of getting lines, get the curves ? I have been trying to use the transform inside the ggpredict but I keep having non parallel lines – Rosa Maria Mar 20 '22 at 15:46
  • No, it's only possible to plot predictions on the response scale. However, you may try [`modelbased::estimate_link()`](https://easystats.github.io/modelbased/reference/estimate_expectation.html) to plot predictions on the link-scale. – Daniel Mar 22 '22 at 10:13