sjp.glmer {sjPlot}

This document shows examples for using the sjp.glmer function of the sjPlot package.

Ressources:

(back to table of content)

Data initialization

Please refer to this document.

# load packages
library(sjPlot)
library(sjmisc)
# load sample data set.
data(efc)
sjp.setTheme(theme = "forestgrey", 
             geom.label.size = 3, 
             axis.textsize = .9, 
             axis.title.size = .9)

Customizing plot appearance

Please refer to this document

Plotting random effects of generalized linear mixed effects models

The sjp.glmer function plots effects of merMod objects, which were fitted using the glmer function of the lme4 package.

First, we need fit a sample model.

# fit model
library(lme4)
# create binary response
sleepstudy$Reaction.dicho <- 
  dicho(sleepstudy$Reaction, dich.by = "median")
# fit model
fit <- glmer(Reaction.dicho ~ Days + (Days | Subject),
             sleepstudy, family = binomial("logit"))

By default, random effects are plotted. In this example, the random effects of random intercept and random coefficient(s) are plotted as an integrated (faceted plot.)

sjp.glmer(fit, y.offset = .4)

With the sort.est argument you can specify the random intercept or any random coefficient in order to sort the effects accordingly.

# sort by predictor Days
sjp.glmer(fit, sort.est = "Days", y.offset = .4)

Since the grouping levels in a faceted plot define the x-axis for each facet, sorting can only be applied to one coefficient or the intercept. If all random intercepts and random coefficients should be sorted, turn faceting off and use sort.est = "sort.all".

# sort all predictors
sjp.glmer(fit,
          facet.grid = FALSE,
          sort.est = "sort.all",
          y.offset = .4)

Plotting fixed effects of generalized linear mixed effects models

For the next examples, we fit another model.

library(lme4)
library(sjmisc)
data(efc)
# create binary response
efc$hi_qol <- dicho(efc$quol_5)
# prepare group variable
efc$grp = as.factor(efc$e15relat)
levels(x = efc$grp) <- get_labels(efc$e15relat)
# data frame for fitted model
mydf <- data.frame(hi_qol = efc$hi_qol,
                   sex = to_factor(efc$c161sex),
                   c12hour = efc$c12hour,
                   neg_c_7 = efc$neg_c_7,
                   grp = efc$grp)
# fit glmer
fit2 <- glmer(hi_qol ~ sex + c12hour + neg_c_7 + (1 | grp),
              data = mydf, family = binomial("logit"))

With the type argument, you can decide which effects to plot. type = "fe" plots the fixed effects.

# plot fixed effects
sjp.glmer(fit2, type = "fe")

Predicted values and marginal effects

In some cases, it is easier to interprete the predicted probabilities, incidents rates or marginal effects instead of odds ratios. The sjp.glmer function has different plot typed to plot predicted values or marginal effects.

  1. type = "eff" to plot marginal effects, adjusted for all predictors.
  2. type = "pred" or type = "pred.fe" to plot predicted values against reponse, for particular model terms.
  3. type = "ri.slope" or type = "fe.slope" to plot unadjusted predicted values, i.e. the relation between model terms and response.

Marginal effect plots

Marginal effects can be plotted using type = "eff", which internally calls the effects::allEffects function. However, unlike the plot function for effects objects, the scale used in sjp.glmer is not transformed and used the complete range of the affected predictors. Hence, the effect plots may more look like a curve, instead of linear lines.

# marginal effects
sjp.glmer(fit2, type = "eff", show.ci = TRUE)

Predicted values

type = "pred" or type = "pred.fe" give predicted values against response, only fixed effects or conditional on random intercept. This plot type requires the vars argument to be specified. vars should be a character vector of length one or two, indicating the term that should be used for predictions, and a second, optional model term, which will be used as grouping level.

Predictions are based on predict(fit, type = "response", re.form = NA) resp. predict(fit, type = "response", re.form = NULL).

# predicted values for negative impact
sjp.glmer(fit2, type = "pred", vars = "neg_c_7")

The predictions can be grouped if vars is a two-element-vector.

# predicted values for negative impact, grouped by gender
sjp.glmer(fit2, type = "pred", vars = c("neg_c_7", "sex"))

If you have a grouping predictor in non-faceted plots, a legend is added to see the group labels.

# predicted values for negative impact, grouped by gender
sjp.glmer(fit2, type = "pred", 
          vars = c("neg_c_7", "sex"),
          facet.grid = F)

Probabilities or incidents rates of fixed effects by grouping level

With type = "ri.slope", probability or incidents rate plots for each covariate (predicted values of coefficients) can be plotted, depending on the grouping level from the random intercept. Thus, for each covariate as many plots as grouping levels are plotted. Furthermore, with the show.ci the standard error of probabilities can be shown.

The predicted values are based on the fixed effects intercept, plus each random intercept and each specific fixed term’s estimate. All other fixed effects are set to zero (i.e. ignored), which corresponds to family(fit)$linkinv(eta = b0 + b0[r1-rn] + bi * xi) (where xi is the estimate of fixed effects, b0 is the intercept of the fixed effects and b0[r1-rn] are all random intercepts).

# plot probability curves for each covariate
# grouped by random intercepts
sjp.glmer(fit2, type = "ri.slope", show.ci = TRUE)

To select a specific predictor only, use the vars argument.

# plot probability curves for each covariate
# grouped by random intercepts
sjp.glmer(fit2,
          type = "ri.slope",
          vars = "neg_c_7",
          show.ci = TRUE)

Probabilities or incidents rates of fixed effects

The same can be done for fixed effects. With type = "fe.slope", probability or incidents rate plots for each covariate can be plotted. These predicted values are based on the fixed effects intercept, not on the random intercept. Thus, only one plot per covariate is plotted.

The predicted values are based on the fixed effects intercept’s estimate and each specific fixed term’s estimate. All other fixed effects are set to zero (i.e. ignored), which corresponds to family(fit)$linkinv(eta = b0 + bi * xi) (where xi is the estimate of fixed effects and b0 is the intercept of the fixed effects).

# plot probability curve of fixed effects
sjp.glmer(fit2, type = "fe.slope")

Predicted values in integrated (non faceted) plots

The predicted values by grouping levels can also be plotted in a single, non faceted plot with facet.grid = FALSE.

# plot probability curves for each covariate
# grouped by random intercepts
sjp.glmer(fit2,
          type = "ri.slope",
          facet.grid = FALSE)

Note that due to different x-axis ranges, the integrated, non-faceted plot cannot be applied to fixed effects only (when type = "fe.slope").

Emphasizing specific grouping levels

In non-faceted plots, grouping levels might be difficult to distinguish. However, you can emphasize specific groups with emph.grp. Note that emph.grp only works in non-faceted plots (i.e. facet.grid = FALSE)! Remaining (non-emphasized) groups have a light grey color.

# plot fixed effects depending on group levels,
# emphasize group levels 1, 2 and 5
sjp.glmer(fit2,
          type = "ri.slope",
          emph.grp = c(1, 2, 5),
          facet.grid = FALSE)

Grouping levels can also be specified via their names.

# plot fixed effects depending on group levels,
# emphasize groups "child" and "cousin",
# only for predictor "neg_c_7"
sjp.glmer(fit2,
          type = "ri.slope",
          emph.grp = c("child", "cousin"),
          facet.grid = FALSE,
          vars = "neg_c_7")

Diagnostic plots

Currently, there are two type options to plot diagnostic plots.

Correlation matrix of fixed effects

To plot a correlation matrix of the fixed effects, use type = "fe.cor"

# plot fixed effects correlation matrix
sjp.glmer(fit2, type = "fe.cor")

qq-plot of random effects

Another diagnostic plot is the qq-plot for random effects. Use type = "re.qq" to plot random against standard quantiles. The dots should be plotted along the line.

# plot qq-plot of random effects
sjp.glmer(fit2, type = "re.qq")

If you have other random effects, like random coefficients, qq-plots for these effects are plotted as well. We refer to the first model to demonstrate this.

# plot qq-plot of random effects
sjp.glmer(fit, type = "re.qq")