This document describes how to plot marginal effects of interaction terms from 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. To plot marginal effects of interaction terms, call `plot_model()`

with:

`type = "pred"`

to plot predicted values (marginal effects) for specific model terms, including interaction terms.`type = "eff"`

, which is similar to`type = "pred"`

, however, discrete predictors are held constant at their proportions (not reference level). It internally calls via .`type = "emm"`

, which is similar to`type = "eff"`

. It internally calls via .`type = "int"`

to plot marginal effects of interaction terms in a more convenient way.

`plot_model()`

supports labelled data and automatically uses variable and value labels to annotate the plot. This works with most regression modelling functions.

**Note:** For marginal effects plots, **sjPlot** calls functions from the **ggeffects-package**. If you need more flexibility when creating marginal effects plots, consider directly using the **ggeffects**-package.

*Note: To better understand the principle of plotting interaction terms, it might be helpful to read the vignette on marginal effects first.*

To plot marginal effects of interaction terms, at least two model terms need to be specified (the terms that define the interaction) in the `terms`

-argument, for which the effects are computed. To plot marginal effects for three-way-interactions, all three terms need to be specified in `terms`

.

A convenient way to automatically plot interactions is `type = "int"`

, which scans the model formula for interaction terms and then uses these as `terms`

-argument.

```
library(sjPlot)
library(sjmisc)
library(ggplot2)
data(efc)
theme_set(theme_sjplot())
# make categorical
$c161sex <- to_factor(efc$c161sex)
efc
# fit model with interaction
<- lm(neg_c_7 ~ c12hour + barthtot * c161sex, data = efc)
fit
plot_model(fit, type = "pred", terms = c("barthtot", "c161sex"))
```

For `type = "int"`

, no terms need to be specified. Note that this plot type automatically uses the first interaction term in the formula for the x-axis, while the second term is used as grouping factor. Furthermore, if continuous variables are used as second term, you can specify preset-values for this term with the `mdrt.values`

-argument, which are then used as grouping levels.

In this example, the second term is a factor with two levels (male/female), so there is no need for choosing specific values for the moderator.

`plot_model(fit, type = "int")`

To switch the terms, in this example *barthtot* and *c161sex*, simply switch the order of these terms on the `terms`

-argument and use `type = "pred"`

.

`plot_model(fit, type = "pred", terms = c("c161sex", "barthtot [0, 100]"))`

To switch the terms for plot-type `type = "int"`

, you need to re-fit the model and change the formula accordingly, i.e. using *c161sex* as first term in the interaction.

```
# fit model with interaction, switching terms in formula
<- lm(neg_c_7 ~ c12hour + c161sex * barthtot, data = efc)
fit plot_model(fit, type = "int")
```

By default, for continuous variables, the minimum and maximum values are chosen as grouping levels, which are 0 and 100 - that’s why the previous two plots are identical. You have other options as well, e.g. the mean-value and +/- 1 standard deviation (as suggested by Cohen and Cohen for continuous variables and popularized by Aiken and West 1991), which can be specified using `mdrt.values`

.

`plot_model(fit, type = "int", mdrt.values = "meansd")`

Since the `terms`

-argument accepts up to three model terms, you can also compute marginal effects for a 3-way-interaction.

```
# fit model with 3-way-interaction
<- lm(neg_c_7 ~ c12hour * barthtot * c161sex, data = efc)
fit
# select only levels 30, 50 and 70 from continuous variable Barthel-Index
plot_model(fit, type = "pred", terms = c("c12hour", "barthtot [30,50,70]", "c161sex"))
```

Again, `type = "int"`

will automatically plot the interaction terms, however, using `mdrt.values = "minmax"`

as default - in this case, the “levels” 0 and 100 from continuous variable *barthtot* are chosen by default.

```
plot_model(fit, type = "int")
#> [[1]]
```

```
#>
#> [[2]]
```

```
#>
#> [[3]]
```

```
#>
#> [[4]]
```

Aiken and West (1991). *Multiple Regression: Testing and Interpreting Interactions.*