This document shows examples for using the
sjp.lmer function of the sjPlot package.
(back to table of content)
Please refer to this document.
library(sjPlot) library(sjmisc) data(efc) sjp.setTheme(theme = "forestgrey", geom.label.size = 3.5, axis.textsize = .9, axis.title.size = .9)
Please refer to this document
sjp.lmer function plots effects of
merMod objects, which were fitted using the
lmer function of the
First, we need fit a sample model.
# fit model library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
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.) Note that the
y.offset argument is used to adjust the value label position. Depending on text size and screen resolution, the default position of text labels may vary.
sjp.lmer(fit, y.offset = .4)
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.lmer(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.lmer(fit, facet.grid = FALSE, sort.est = "sort.all", y.offset = .4)
For the next examples, we fit another model.
library(lme4) # for sample data and tool functions library(sjmisc) data(efc) # 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(neg_c_7 = efc$neg_c_7, sex = to_factor(efc$c161sex), c12hour = efc$c12hour, barthel = efc$barthtot, grp = efc$grp) # fit glmer fit2 <- lmer(neg_c_7 ~ sex + c12hour + barthel + (1 | grp), data = mydf)
type argument, you can decide which effects to plot.
type = "fe" plots the fixed effects. In this case, the value label for the negative estimate is not completely visible in the output, due to tight axis limits. Use
axis.lim to manually change the axis limits.
# plot fixed effects sjp.lmer(fit2, type = "fe", axis.lim = c(-5, 15))
type = "fe.std", we can plot the standardized estimates of the fixed effects.
# plot standardized fixed effects sjp.lmer(fit2, type = "fe.std")
type = "fe.slope", slopes for each fixed coefficient (i.e. each predictor regressed on the response) are plotted. Specific predictors can be selected with
vars or removed from the plot with the
# plot fixed effects slopes sjp.lmer(fit2, type = "fe.slope", vars = c("c12hour", "barthel"))
type = "eff" produces effect plots of the mixed model. The
Effect function from the effects package is used to compute the data.
# plot effects sjp.lmer(fit2, type = "eff")
type = "pred" or
type = "pred.fe" plot predicted values for response, conditional on fixed effects only or on random intercept. It’s calling
predict(fit, type = "response", re.form = NA) resp.
predict(fit, type = "response", re.form = NULL) to compute the values. This plot type requires the
vars argument to select specific terms that should be used for the x-axis and - optional - as grouping factor. Hence, vars must be a character vector with the names of one or two model predictors.
# plot effects sjp.lmer(fit2, type = "pred", vars = "c12hour")
Or grouped by gender:
# plot effects sjp.lmer(fit2, type = "pred", vars = c("c12hour", "sex"))
You can also plot predictions by groups in a non-facted plot:
# plot effects sjp.lmer(fit2, type = "pred", facet.grid = FALSE, vars = c("c12hour", "sex"))
To get a better picture of the linear relationship between fixed effects and response depending on the grouping levels (random intercepts), you can plot straight slope lines (ablines) for each coefficient with varying random intercepts.
Basically, the formula is
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).
type = "ri.slope" for this kind of plots.
# random intercepts ranef(fit2)
## $grp ## (Intercept) ## spouse/partner 0.62318882 ## child 0.42390010 ## sibling -0.05435509 ## daughter or son -in-law 0.05759517 ## ancle/aunt -0.10111187 ## nephew/niece -0.55060420 ## cousin -0.11383605 ## other, specify -0.28477688
# fixed effects fixef(fit2)
## (Intercept) sex2 c12hour barthel ## 14.135947301 0.478500190 0.003365667 -0.047946653
# plot fixed effects depending on group levels sjp.lmer(fit2, type = "ri.slope")
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.lmer(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 barthel sjp.lmer(fit2, type = "ri.slope", emph.grp = c("child", "cousin"), facet.grid = FALSE, vars = "barthel")
Use this plot type to visualize the random parts of random slope-intercept (or repeated measure) models. When having too many groups, use the
sample.n argument to randomly select a specific amount of subjects.
# plot random-slope-intercept sjp.lmer(fit, type = "rs.ri", sample.n = 15)
sample.n is a vector of length greater than one, the specifc “subjects” indicated by
sample.n are plotted.
# plot random-slope-intercept, plot subjects 1, 5 and 7. sjp.lmer(fit, type = "rs.ri", sample.n = c(1, 5, 7), show.legend = TRUE)
If you want to plot only selected coefficients for each random intercept, use the
vars argument and specify the coefficient’s name.
# plot fixed effect "c12hour" only, # depending on group levels sjp.lmer(fit2, type = "ri.slope", vars = "c12hour")
Currently, there are two
type options to plot diagnostic plots.
To plot a correlation matrix of the fixed effects, use
type = "fe.cor"
# plot fixed effects correlation matrix sjp.lmer(fit2, type = "fe.cor")
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.lmer(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.lmer(fit, type = "re.qq")