sjp.glmer: Plot estimates, predictions or effects of generalized linear mixed effects models

Description

By default, this function plots estimates (odds, risk or incidents ratios, i.e. exponentiated coefficients, depending on family and link function) with confidence intervals of either fixed effects or random effects of generalized linear mixed effects models (that have been fitted with the glmer-function of the lme4-package). Furthermore, this function also plots predicted probabilities / incidents or diagnostic plots.

Usage

sjp.glmer(fit, type = "re", vars = NULL, ri.nr = NULL,
  group.estimates = NULL, remove.estimates = NULL, emph.grp = NULL,
  sample.n = NULL, sort.est = NULL, title = NULL, legend.title = NULL,
  axis.labels = NULL, axis.title = NULL, geom.colors = "Set1",
  geom.size = NULL, show.values = TRUE, show.p = TRUE, show.ci = FALSE,
  show.legend = FALSE, show.intercept = FALSE,
  string.interc = "(Intercept)", show.scatter = TRUE, point.alpha = 0.2,
  point.color = NULL, jitter.ci = FALSE, fade.ns = FALSE,
  axis.lim = NULL, digits = 2, vline.type = 2, vline.color = "grey70",
  facet.grid = TRUE, free.scale = FALSE, y.offset = 0.1,
  prnt.plot = TRUE, ...)

Arguments

fit

A fitted model as returned by the glmer-function.

type

Type of plot. Use one of following:

"re": (default) for conditional modes (odds or incidents ratios) of random effects
"fe": for odds or incidents ratios of fixed effects
"fe.cor": for correlation matrix of fixed effects
"re.qq": for a QQ-plot of random effects (random effects quantiles against standard normal quantiles)
"ri.slope": to plot probability or incidents curves (predicted probabilities or incidents) of random intercept variances for all fixed effects coefficients. Use facet.grid to decide whether to plot each coefficient as separate plot or as integrated faceted plot. See 'Details'.
"rs.ri": for fitted probability curves (predicted probabilities) indicating the random slope-intercept pairs. Use this to visualize the random parts of random slope-intercept (or repeated measure) models. When having too many groups, use sample.n argument.
"eff": to plot marginal effects of predicted probabilities or incidents for each fixed term, where remaining co-variates are set to the mean. Use facet.grid to decide whether to plot each coefficient as separate plot or as integrated faceted plot. See 'Details'.
"pred": to plot predicted probabilities or incidents for the response, related to specific model predictors and conditioned on random effects. See 'Details'.
"pred.fe": to plot predicted probabilities or incidents for the response, related to specific model predictors, only for fixed effects. See 'Details'.
"ma": to check model assumptions. Note that only argument fit applies to this plot type. All other arguments are ignored.

vars

Numeric vector with column indices of selected variables or a character vector with variable names of selected variables from the fitted model, which should be used to plot - depending on type - estimates, fixed effects slopes or predicted values (mean, probabilities, incidents rates, ...). See 'Examples'.

ri.nr

Numeric vector. If type = "re" or type = "ri.slope", and fitted model has more than one random intercept, ri.nr indicates which random effects of which random intercept (or: which list elements of ranef) will be plotted. Default is NULL, so all random effects will be plotted.

group.estimates

Numeric or character vector, indicating a group identifier for each estimate. Dots and confidence intervals of estimates are coloured according to their group association. See 'Examples'.

remove.estimates

Character vector with coefficient names that indicate which estimates should be removed from the plot. remove.estimates = "est_name" would remove the estimate est_name. Default is NULL, i.e. all estimates are printed.

emph.grp

Numeric vector with index numbers of grouping levels (from random effect). If type = "ri.slope" and facet.grid = FALSE, an integrated plot of predicted probabilities of fixed effects resp. fixed effects slopes for each grouping level is plotted. To better find certain groups, use this argument to emphasize these groups in the plot. See 'Examples'.

sample.n

Numeric vector. only applies, if type = "rs.ri". If plot has many random intercepts (grouping levels), overplotting of regression lines may occur. In this case, consider random sampling of grouping levels. If sample.n is of length 1, a random sample of sample.n observation is selected to plot random intercepts. If sample.n is of length > 1, random effects indicated by the values in sample.n are selected to plot random effects. Use the latter option to always select a fixed, identical set of random effects for plotting (useful when ecomparing multiple models).

sort.est

Determines in which way estimates are sorted in the plot:

If NULL (default), no sorting is done and estimates are sorted in order of model coefficients.
If sort.est = "sort.all", estimates are re-sorted for each coefficient (only applies if type = "re" and facet.grid = FALSE), i.e. the estimates of the random effects for each predictor are sorted and plotted to an own plot.
If type = "fe" or type = "fe.std", TRUE will sort estimates
If type = "re", specify a predictor's / coefficient's name to sort estimates according to this coefficient.

See 'Examples'.

title

Character vector with one or more labels that are used as plot title.

legend.title

Character vector, used as title for the plot legend. Note that only some plot types have legends (e.g. type = "pred" or when grouping estimates with group.estimates).

axis.labels

Character vector with labels for the model terms, used as axis labels. For mixed models, should either be vector of fixed effects variable labels (if type = "fe" or type = "fe.std") or a vector of group (value) labels from the random intercept's categories (if type = "re").

axis.title

Character vector of length one or two (depending on the plot function and type), used as title(s) for the x and y axis. If not specified, a default labelling is chosen. To set multiple axis titles (e.g. with type = "eff" for many predictors), axis.title must be a character vector of same length of plots that are printed. In this case, each plot gets an own axis title (applying, for instance, to the y-axis for type = "eff"). Note: Some plot types do not support this argument. In such cases, use the return value and add axis titles manually with labs, e.g.: $plot.list[[1]] + labs(x = ...)

geom.colors

User defined color palette for geoms. If group.estimates is not specified, must either be vector with two color values or a specific color palette code (see 'Details' in sjp.grpfrq). Else, if group.estimates is specified, geom.colors must be a vector of same length as groups. See 'Examples'.

geom.size

size resp. width of the geoms (bar width, line thickness or point size, depending on plot type and function). Note that bar and bin widths mostly need smaller values than dot sizes.

show.values

Logical, whether values should be plotted or not.

show.p

Logical, adds significance levels to values, or value and variable labels.

show.ci

Logical, if TRUE, depending on type, a confidence interval or region is added to the plot. For frequency plots, the confidence interval for the relative frequencies are shown.

show.legend

logical, if TRUE, and depending on plot type and function, a legend is added to the plot.

show.intercept

Logical, if TRUE, the intercept of the fitted model is also plotted. Default is FALSE. For glm's, please note that due to exponential transformation of estimates, the intercept in some cases can not be calculated, thus the function call is interrupted and no plot printed.

string.interc

String, axis label of intercept estimate. Only applies, if show.intercept = TRUE and axis.labels is not NULL.

show.scatter

Logical, if TRUE (default), adds a scatter plot of data points to the plot. Only applies for slope-type or predictions plots. For most plot types, dots are jittered to avoid overplotting, hence the points don't reflect exact values in the data.

point.alpha

Alpha value of point-geoms in the scatter plots. Only applies, if show.scatter = TRUE.

point.color

Color of of point-geoms in the scatter plots. Only applies, if show.scatter = TRUE.

jitter.ci

Logical, if TRUE and show.ci = TRUE and confidence bands are displayed as error bars, adds jittering to lines and error bars to avoid overlapping.

fade.ns

Logical, if TRUE, non significant estimates will be printed in slightly faded colors.

axis.lim

Numeric vector of length 2, defining the range of the plot axis. Depending on plot type, may effect either x- or y-axis, or both. For multiple plot outputs (e.g., from type = "eff" or type = "slope" in sjp.glm), axis.lim may also be a list of vectors of length 2, defining axis limits for each plot (only if non-faceted).

digits

Numeric, amount of digits after decimal point when rounding estimates and values.

vline.type

Linetype of the vertical "zero point" line. Default is 2 (dashed line).

vline.color

Color of the vertical "zero point" line. Default value is "grey70".

facet.grid

TRUE to arrange the lay out of of multiple plots in a grid of an integrated single plot. This argument calls facet_wrap or facet_grid to arrange plots. Use plot_grid to plot multiple plot-objects as an arranged grid with grid.arrange.

free.scale

Logical, if TRUE and facet.grid = TRUE, each facet grid gets its own fitted scale. If free.scale = FALSE, each facet in the grid has the same scale range.

y.offset

numeric, offset for text labels when their alignment is adjusted to the top/bottom of the geom (see hjust and vjust).

prnt.plot

logical, if TRUE (default), plots the results as graph. Use FALSE if you don't want to plot any graphs. In either case, the ggplot-object will be returned as value.

...

Other arguments passed down to further functions. Currently, following arguments are supported:

?effects::effect: Any arguments accepted by the effect resp. allEffects function, for type = "eff".
width: The width-argument for error bars.
alpha: The alpha-argument for confidence bands.
level: The level-argument confidence bands.

Value

(Insisibily) returns, depending on the plot type

The ggplot-object (plot). For multiple plots and if facet.grid = FALSE) a plot.list is returned.
A data frame data with the data used to build the ggplot-object(s), or a list of data frames (data.list).

Details

type = "re": plots the conditional modes of the random effects, inclduing predicion intervals. It basically does the same as dotplot(exp(ranef(fit, condVar = TRUE)[[i]]), where i denotes the random effect index.
type = "fe.slope": 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; the inverse link-function is used). This plot type may give similar results as type = "pred", however, type = "fe.slope" does not adjust for other predictors.
type = "eff": plots the marginal effects of model predictors. Unlike type = "fe.slope", the predicted values computed by type = "eff" are adjusted for all co-variates, which are set to the mean, as returned by the allEffects function. You can pass further arguments down to allEffects for flexible function call via the ...-argument.
type = "ri.slope": 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).
type = "rs.ri": the predicted values are based on the fixed effects intercept, plus each random intercept and random slope. This plot type is intended to plot the random part, i.e. the predicted probabilities or incident rates of each random slope for each random intercept. Since the random intercept specifies the deviance from the gloabl intercept, the global intercept is always included. In case of overplotting, use the sample.n argument to randomly sample a limited amount of groups.
type = "coef": forest plot of joint fixed and random effect coefficients, as retrieved by coef.merMod, it's simply ranef + fixef.
type = "pred" or type = "pred.fe": predicted values against response, only fixed effects or conditional 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. See 'Examples'.

Examples

Run this code

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

# simple plot
sjp.glmer(fit)

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

# }
# NOT RUN {
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 = to_factor(efc$hi_qol),
                   sex = to_factor(efc$c161sex),
                   education = to_factor(efc$c172code),
                   c12hour = efc$c12hour,
                   neg_c_7 = efc$neg_c_7,
                   grp = efc$grp)

# fit glmer, with categorical predictor with more than 2 levels
fit <- glmer(hi_qol ~ sex + education + c12hour + neg_c_7 + (1|grp),
             data = mydf, family = binomial("logit"))

# plot and sort fixed effects, axis labels automatically retrieved
sjp.glmer(fit, type = "fe", sort.est = TRUE)

# plot probability curves (predicted probabilities)
# for each covariate, grouped by random intercepts
# in integrated plots, emphasizing groups 1 and 4
sjp.glmer(fit, type = "ri.slope", emph.grp = c(1, 4), facet.grid = FALSE)

# plot predicted probabilities for response,
# non faceted, with ci
sjp.glmer(fit, type = "pred.fe", vars = c("neg_c_7", "education"),
          show.ci = TRUE, facet.grid = FALSE)

# predictions by gender and education
sjp.glmer(fit, type = "pred.fe", vars = c("neg_c_7", "sex", "education"))
# }
# NOT RUN {
# }