knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "png", fig.width = 7, fig.height = 5, message = FALSE, warning = FALSE)

This document describes how to plot estimates as forest plots (or dot whisker plots) of various regression models, using the `plot_model()`

function. `plot_model()`

is a generic plot-function, which accepts many model-objects, like `lm`

, `glm`

, `lme`

, `lmerMod`

etc.

`plot_model()`

allows to create various plot tyes, which can be defined via the `type`

-argument. The default is `type = "fe"`

, which means that fixed effects (model coefficients) are plotted. For mixed effects models, only fixed effects are plotted by default as well.

library(sjPlot) library(sjlabelled) library(sjmisc) library(ggplot2) data(efc) theme_set(theme_sjplot())

First, we fit a model that will be used in the following examples. The examples work in the same way for any other model as well.

# create binary response y <- ifelse(efc$neg_c_7 < median(na.omit(efc$neg_c_7)), 0, 1) # create data frame for fitting model df <- data.frame( y = to_factor(y), sex = to_factor(efc$c161sex), dep = to_factor(efc$e42dep), barthel = efc$barthtot, education = to_factor(efc$c172code) ) # set variable label for response set_label(df$y) <- "High Negative Impact" # fit model m1 <- glm(y ~., data = df, family = binomial(link = "logit"))

The simplest function call is just passing the model object as argument. By default, estimates are sorted in descending order, with the highest effect at the top.

plot_model(m1)

The "neutral" line, i.e. the vertical intercept that indicates no effect (x-axis position 1 for most glm's and position 0 for most linear models), is drawn slightly thicker than the other grid lines. You can change the line color with the `vline.color`

-argument.

plot_model(m1, vline.color = "red")

By default, the estimates are sorted in the same order as they were introduced into the model. Use `sort.est = TRUE`

to sort estimates in descending order, from highest to lowest value.

plot_model(m1, sort.est = TRUE)

Another way to sort estimates is to use the `order.terms`

-argument. This is a numeric vector, indicating the order of estimates in the plot. In the summary, we see that "sex2" is the first term, followed by the three dependency-categories (position 2-4), the Barthel-Index (5) and two levels for intermediate and high level of education (6 and 7).

```
summary(m1)
```

Now we want the educational levels (6 and 7) first, than gender (1), followed by dependency (2-4)and finally the Barthel-Index (5). Use this order as numeric vector for the `order.terms`

-argument.

plot_model(m1, order.terms = c(6, 7, 1, 2, 3, 4, 5))

By default, `plot_model()`

automatically exponentiates coefficients, if appropriate (e.g. for models with log or logit link). You can explicitley prevent transformation by setting the `transform`

-argument to `NULL`

, or apply any transformation by using a character vector with the function name.

plot_model(m1, transform = NULL) plot_model(m1, transform = "plogis")

By default, just the dots and error bars are plotted. Use `show.values = TRUE`

to show the value labels with the estimates values, and use `show.p = FALSE`

to suppress the asterisks that indicate the significance level of the p-values. Use `value.offset`

to adjust the relative positioning of value labels to the dots and lines.

plot_model(m1, show.values = TRUE, value.offset = .3)

As seen in the above examples, by default, the plotting-functions of **sjPlot** retrieve value and variable labels if the data is *labelled*, using the sjlabelled-package. If the data is not labelled, the variable names are used. In such cases, use the arguments `title`

, `axis.labels`

and `axis.title`

to annotate the plot title and axes. If you want variable names instead of labels, even for labelled data, use `""`

as argument-value, e.g. `axis.labels = ""`

, or set `auto.label`

to `FALSE`

.

Furthermore, `plot_model()`

applies case-conversion to all labels by default, using the snakecase-package. This converts labels into human-readable versions. Use `case = NULL`

to turn case-conversion off, or refer to the package-vignette of the **snakecase**-package for further options.

data(iris) m2 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Species, data = iris) # variable names as labels, but made "human readable" # separating dots are removed plot_model(m2) # to use variable names even for labelled data plot_model(m1, axis.labels = "", title = "my own title")

Use `terms`

resp. `rm.terms`

to select specific terms that should (not) be plotted.

# keep only coefficients sex2, dep2 and dep3 plot_model(m1, terms = c("sex2", "dep2", "dep3")) # remove coefficients sex2, dep2 and dep3 plot_model(m1, rm.terms = c("sex2", "dep2", "dep3"))

For linear models, you can also plot standardized beta coefficients, using `type = "std"`

or `type = "std2"`

. These two options differ in the way how coefficients are standardized. `type = "std2"`

plots standardized beta values, however, standardization follows Gelman's (2008) suggestion, rescaling the estimates by dividing them by two standard deviations instead of just one.

plot_model(m2, type = "std")

`plot_model()`

also supports stan-models fitted with the **rstanarm** or **brms** packages. However, there are a few differences compared to the previous plot examples.

First, of course, there are no *confidence intervals*, but *uncertainty intervals* - high density intervals, to be precise.

Second, there's not just one interval range, but an *inner* and *outer* probability. By default, the inner probability is fixed to `.5`

(50%), while the outer probability is specified via `ci.lvl`

(which defaults to `.89`

(89%) for Bayesian models). However, you can also use the arguments `prob.inner`

and `prob.outer`

to define the intervals boundaries.

Third, the point estimate is by default the *median*, but can also be another value, like mean. This can be specified with the `bpe`

-argument.

if (require("rstanarm", quietly = TRUE)) { # make sure we apply a nice theme library(ggplot2) theme_set(theme_sjplot()) data(mtcars) m <- stan_glm(mpg ~ wt + am + cyl + gear, data = mtcars, chains = 1) # default model plot_model(m) # same model, with mean point estimate, dot-style for point estimate # and different inner/outer probabilities of the HDI plot_model( m, bpe = "mean", bpe.style = "dot", prob.inner = .4, prob.outer = .8 ) }

There are several options to customize the plot appearance:

- The
`colors`

-argument either takes the name of a valid colorbrewer palette (see also the related vignette),`"bw"`

or`"gs"`

for black/white or greyscaled colors, or a string with a color name. `value.offset`

and`value.size`

adjust the positioning and size of value labels, if shown.`dot.size`

and`line.size`

change the size of dots and error bars.`vline.color`

changes the neutral "intercept" line.`width`

,`alpha`

and`scale`

are passed down to certain ggplot-geoms, like`geom_errorbar()`

or`geom_density_ridges()`

.

plot_model( m1, colors = "Accent", show.values = TRUE, value.offset = .4, value.size = 4, dot.size = 3, line.size = 1.5, vline.color = "blue", width = 1.5 )

Gelman A (2008) *Scaling regression inputs by dividing by two standard deviations.* Statistics in Medicine 27: 2865–2873.

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.