Plotting observed data

If you want to see the individual data points plotted to better understand how the fitted lines related to the observed data, you can use the plot.points = TRUE argument.

interact_plot(fiti, pred = Illiteracy, modx = Murder, plot.points = TRUE)

For continuous moderators, as you can see, the observed data points are shaded depending on the level of the moderator variable.

It can be very enlightening, too, for categorical moderators.

interact_plot(fitiris, pred = Petal.Width, modx = Species,
              plot.points = TRUE)

Where many points are slightly overlapping as they do here, it can be useful to apply a random “jitter” to move them slightly to stop the overlap. Use the jitter argument to do this. If you provide a single number it will apply a jitter of proportional magnitude both sideways and up and down. If you provide a vector length 2, then the first is assumed to refer to the left/right jitter and the second to the up/down jitter.

If you would like to better differentiate the points with factor moderators, you can use point.shape = TRUE to give a different shape to each point. This can be especially helpful for black and white publications.

interact_plot(fitiris, pred = Petal.Width, modx = Species,
              plot.points = TRUE, jitter = 0.1, point.shape = TRUE)

If your original data are weighted, then the points will be sized based on the weight. For the purposes of our example, we’ll weight the same model we’ve been using with the population of each state.

fiti <- lm(Income ~ Illiteracy * Murder, data = states,
           weights = Population)
interact_plot(fiti, pred = Illiteracy, modx = Murder, plot.points = TRUE)

For those working with weighted data, it can be hard to use a scatterplot to explore the data unless there is some way to account for the weights. Using size is a nice middle ground.

Plotting partial residuals

In more complex regressions, plotting the observed data can sometimes be relatively uninformative because the points seem to be all over the place.

For an example, let’s take a look at this model. I am using the mpg dataset from ggplot2 and predicting the city miles per gallon (cty) based on several variables, including model year, type of car, fuel type, drive type, and an interaction between engine displacement (displ) and number of cylinders in the engine (cyl). Let’s take a look at the regression output:

library(ggplot2)
data(cars)
fitc <- lm(cty ~ year + cyl * displ + class + fl + drv, data = mpg)
summ(fitc)

## MODEL INFO:
## Observations: 234
## Dependent Variable: cty
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(16,217) = 99.73, p = 0.00
## R² = 0.88
## Adj. R² = 0.87 
## 
## Standard errors:OLS
## -------------------------------------------------------
##                            Est.    S.E.   t val.      p
## --------------------- --------- ------- -------- ------
## (Intercept)             -200.98   47.01    -4.28   0.00
## year                       0.12    0.02     5.03   0.00
## cyl                       -1.86    0.28    -6.69   0.00
## displ                     -3.56    0.66    -5.41   0.00
## classcompact              -2.60    0.93    -2.80   0.01
## classmidsize              -2.63    0.93    -2.82   0.01
## classminivan              -4.41    1.04    -4.24   0.00
## classpickup               -4.37    0.93    -4.68   0.00
## classsubcompact           -2.38    0.93    -2.56   0.01
## classsuv                  -4.27    0.87    -4.92   0.00
## fld                        6.34    1.69     3.74   0.00
## fle                       -4.57    1.66    -2.75   0.01
## flp                       -1.92    1.59    -1.21   0.23
## flr                       -0.79    1.57    -0.50   0.61
## drvf                       1.40    0.40     3.52   0.00
## drvr                       0.49    0.46     1.06   0.29
## cyl:displ                  0.36    0.08     4.56   0.00
## -------------------------------------------------------

Okay, we are explaining a lot of variance here but there is quite a bit going on. Let’s plot the interaction with the observed data.

interact_plot(fitc, pred = displ, modx = cyl, plot.points = TRUE,
              modx.values = c(4, 5, 6, 8))

Hmm, doesn’t that look…bad? You can kind of see the pattern of the interaction, but the predicted lines don’t seem to match the data very well. But I included a lot of variables besides cyl and displ in this model and they may be accounting for some of this variation. This is what partial residual plots are designed to help with. You can learn more about the technique and theory in Fox and Weisberg (2018) and another place to generate partial residual plots is in Fox’s effects package.

Using the argument partial.residuals = TRUE, what is plotted instead is the observed data with the effects of all the control variables accounted for. In other words, the value cty for the observed data is based only on the values of displ, cyl, and the model error. Let’s take a look.

interact_plot(fitc, pred = displ, modx = cyl, partial.residuals = TRUE,
              modx.values = c(4, 5, 6, 8))

Now we can understand how the observed data and the model relate to each other much better. One insight is how the model really underestimates the values at the low end of displacement and cylinders. You can also see how much the cylinders and displacement seem to be correlated each other, which makes it difficult to say how much we can learn from this kind of model.

Confidence bands

Another way to get a sense of the precision of the estimates is by plotting confidence bands. To get started, just set interval = TRUE. To decide how wide the confidence interval should be, express the percentile as a number, e.g., int.width = 0.8 corresponds to an 80% interval.

interact_plot(fiti, pred = Illiteracy, modx = Murder, interval = TRUE,
              int.width = 0.8)

You can also use the robust argument to plot confidence intervals based on robust standard error calculations.

Check linearity assumption

A basic assumption of linear regression is that the relationship between the predictors and response variable is linear. When you have an interaction effect, you add the assumption that relationship between your predictor and response is linear regardless of the level of the moderator.

To show a striking example of how this can work, we’ll generate two simple datasets to replicate Hainmueller et al. (2017).

The first has a focal predictor X that is normally distributed with mean 3 and standard deviation 1. It then has a dichotomous moderator W that is Bernoulli distributed with mean probability 0.5. We also generate a normally distributed error term with standard deviation 4.

set.seed(99)
x <- rnorm(n = 200, mean = 3, sd = 1)
err <- rnorm(n = 200, mean = 0, sd = 4)
w <- rbinom(n = 200, size = 1, prob = 0.5)

y_1 <- 5 - 4*x - 9*w + 3*w*x + err

We fit a linear regression model with an interaction between x and w.

model_1 <- lm(y_1 ~ x * w)
summ(model_1)

## MODEL INFO:
## Observations: 200
## Dependent Variable: y_1
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(3,196) = 25.93, p = 0.00
## R² = 0.28
## Adj. R² = 0.27 
## 
## Standard errors:OLS
## ------------------------------------------------
##                      Est.   S.E.   t val.      p
## ----------------- ------- ------ -------- ------
## (Intercept)          3.19   1.29     2.47   0.01
## x                   -3.59   0.42    -8.61   0.00
## w                   -7.77   1.90    -4.10   0.00
## x:w                  2.86   0.62     4.62   0.00
## ------------------------------------------------

In the following plot, we use linearity.check = TRUE argument to split the data by the level of the moderator W and plot predicted lines (black) and a loess line (red) within each group. The predicted lines come from the full data set. The loess line looks only at the subset of data in each pane and will be curved if the relationship is nonlinear. What we’re looking for is whether the red loess line is straight like the predicted line.

interact_plot(model_1, pred = x, modx = w, linearity.check = TRUE, 
              plot.points = TRUE)

In this case, it is. That means the assumption is satisfied.

Now we generate similar data, but this time the linearity assumption will be violated. X₂ will now be uniformly distributed across the range of -3 to 3. The true relationship between y_2 and X₂ (x_2) is non-linear.

x_2 <- runif(n = 200, min = -3, max = 3)
y_2 <- 2.5 - x_2^2 - 5*w + 2*w*(x_2^2) + err
data_2 <- as.data.frame(cbind(x_2, y_2, w))

model_2 <- lm(y_2 ~ x_2 * w, data = data_2)
summ(model_2)

## MODEL INFO:
## Observations: 200
## Dependent Variable: y_2
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(3,196) = 3.28, p = 0.02
## R² = 0.05
## Adj. R² = 0.03 
## 
## Standard errors:OLS
## ------------------------------------------------
##                      Est.   S.E.   t val.      p
## ----------------- ------- ------ -------- ------
## (Intercept)         -0.75   0.49    -1.51   0.13
## x_2                  0.53   0.30     1.77   0.08
## w                    1.63   0.72     2.28   0.02
## x_2:w               -0.25   0.42    -0.61   0.54
## ------------------------------------------------

The regression output would lead us to believe there is no interaction.

Let’s do the linearity check:

interact_plot(model_2, pred = x_2, modx = w, linearity.check = TRUE, 
              plot.points = TRUE)

This is a striking example of the assumption being violated. At both levels of W, the relationship between X₂ and the response is non-linear. There really is an interaction, but it’s a non-linear one.

You could try to fit this true model with the polynomial term like this:

model_3 <- lm(y_2 ~ poly(x_2, 2) * w, data = data_2)
summ(model_3)

## MODEL INFO:
## Observations: 200
## Dependent Variable: y_2
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(5,194) = 18.32, p = 0.00
## R² = 0.32
## Adj. R² = 0.30 
## 
## Standard errors:OLS
## -----------------------------------------------------
##                           Est.   S.E.   t val.      p
## --------------------- -------- ------ -------- ------
## (Intercept)              -1.01   0.42    -2.41   0.02
## poly(x_2, 2)1             9.90   6.14     1.61   0.11
## poly(x_2, 2)2           -34.14   6.15    -5.55   0.00
## w                         1.61   0.61     2.65   0.01
## poly(x_2, 2)1:w          -6.33   8.60    -0.74   0.46
## poly(x_2, 2)2:w          75.58   8.62     8.77   0.00
## -----------------------------------------------------

interact_plot can plot these kinds of effects, too. Just provide the untransformed predictor’s name (in this case, x_2) and also include the data in the data argument. If you don’t include the data, the function will try to find the data you used but it will warn you about it and it could cause problems under some circumstances.

interact_plot(model_3, pred = x_2, modx = w, data = data_2)

Look familiar? Let’s do the linearity.check, which in this case is more like a non-linearity check:

interact_plot(model_3, pred = x_2, modx = w, data = data_2,
              linearity.check = TRUE, plot.points = TRUE)

The red loess line almost perfectly overlaps the black predicted line, indicating the assumption is satisfied. As a note of warning, in practice the lines will rarely overlap so neatly and you will have to make more difficult judgment calls.

Other options

There are a number of other options not mentioned, many relating to the appearance.

For instance, you can manually specify the axis labels, add a main title, choose a color scheme, and so on.

interact_plot(fiti, pred = Illiteracy, modx = Murder,
              x.label = "Custom X Label", y.label = "Custom Y Label",
              main.title = "Sample Plot",  legend.main = "Custom Legend Title",
              colors = "seagreen")

Because the function uses ggplot2, it can be modified and extended like any other ggplot2 object. For example, using the theme_apa function from jtools:

interact_plot(fitiris, pred = Petal.Width, modx = Species) + theme_apa()

Visualize the coefficients

Similar to what this package’s plot_coefs/plot_summs functions offer, you can save your sim_slopes output to an object and call plot on that object.

ss <- sim_slopes(fiti, pred = Illiteracy, modx = Murder, 
                 modx.values = c(0, 5, 10))
plot(ss)

Tabular output

You can also use the huxtable package to get a publication-style table from your sim_slopes output.

ss <- sim_slopes(fiti, pred = Illiteracy, modx = Murder,
                 modx.values = c(0, 5, 10))
library(huxtable)
as_huxtable(ss)

Johnson-Neyman intervals

Did you notice how I was adding the argument johnson_neyman = FALSE above? That’s because by default, sim_slopes will also calculate what is called the Johnson-Neyman interval. This tells you all the values of the moderator for which the slope of the predictor will be statistically significant. Depending on the specific analysis, it may be that all values of the moderator outside of the interval will have a significant slope for the predictor. Other times, it will only be values inside the interval—you will have to look at the output to see.

It can take a moment to interpret this correctly if you aren’t familiar with the Johnson-Neyman technique. But if you read the output carefully and take it literally, you’ll get the hang of it.

sim_slopes(fiti, pred = Illiteracy, modx = Murder, johnson_neyman = TRUE)

## JOHNSON-NEYMAN INTERVAL
## 
## When Murder is OUTSIDE the interval [-6.37, 8.41], the slope of Illiteracy
## is p < .05.
## 
## Note: The range of observed values of Murder is [1.40, 15.10]
## 
## SIMPLE SLOPES ANALYSIS
## 
## Slope of Illiteracy when Murder =  5.420973 (- 1 SD): 
## 
##     Est.     S.E.   t val.      p
## -------- -------- -------- ------
##   -17.43   250.08    -0.07   0.94
## 
## Slope of Illiteracy when Murder =  8.685043 (Mean): 
## 
##      Est.     S.E.   t val.      p
## --------- -------- -------- ------
##   -399.64   178.86    -2.23   0.03
## 
## Slope of Illiteracy when Murder = 11.949113 (+ 1 SD): 
## 
##      Est.     S.E.   t val.      p
## --------- -------- -------- ------
##   -781.85   189.11    -4.13   0.00

In the example above, we can see that the Johnson-Neyman interval and the simple slopes analysis agree—they always will. The benefit of the J-N interval is it will tell you exactly where the predictor’s slope becomes significant/insignificant at a specified alpha level.

You can also call the johnson_neyman function directly if you want to do something like tweak the alpha level. The johnson_neyman function will also create a plot by default — you can get them by setting jnplot = TRUE with sim_slopes.

johnson_neyman(fiti, pred = Illiteracy, modx = Murder, alpha = .05)

## JOHNSON-NEYMAN INTERVAL
## 
## When Murder is OUTSIDE the interval [-10.61, 8.16], the slope of Illiteracy
## is p < .05.
## 
## Note: The range of observed values of Murder is [1.40, 15.10]

A note on Johnson-Neyman plots: Once again, it is easy to misinterpret the meaning. Notice that the y-axis is the conditional slope of the predictor. The plot shows you where the conditional slope differs significantly from zero. In the plot above, we see that from the point Murder (the moderator) = 9.12 and greater, the slope of Illiteracy (the predictor) is significantly different from zero and in this case negative. The lower bound for this interval (about -32) is so far outside the observed data that it is not plotted. If you could have -32 as a value for Murder rate, though, that would be the other threshold before which the slope of Illiteracy would become positive.

The purpose of reminding you both within the plot and the printed output of the range of observed data is to help you put the results in context; in this case, the only justifiable interpretation is that Illiteracy has no effect on the outcome variable except when Murder is higher than 8.16. You wouldn’t interpret the lower boundary because your dataset doesn’t contain any values near it.

False discovery rate adjustment

A recent publication (Esarey & Sumner, 2017) explored ways to calculate the Johnson-Neyman interval that properly manages the Type I and II error rates. Others have noted that the alpha level implied by the Johnson-Neyman interval won’t be quite right (e.g., Bauer & Curran, 2005), but there hasn’t been any general solution that has gotten wide acceptance in the research literature just yet.

The basic problem is that the Johnson-Neyman interval is essentially doing a bunch of comparisons across all the values of the moderator, each one inflating the Type I error rate. The issue isn’t so much that you can’t possibly address it, but many solutions are far too conservative and others aren’t broadly applicable. Esarey and Sumner (2017), among other contributions, suggested an adjustment that seems to do a good job of balancing the desire to be a conservative test without missing a lot of true effects in the process. I won’t go into the details here. The implementation in johnson_neyman is based on code adapted from Esarey and Sumner’s interactionTest package, but any errors should be assumed to be from interactions, not them.

To use this feature, simply set control.fdr = TRUE in the call to johnson_neyman or sim_slopes.

sim_slopes(fiti, pred = Illiteracy, modx = Murder, johnson_neyman = TRUE,
           control.fdr = TRUE)

## JOHNSON-NEYMAN INTERVAL
## 
## When Murder is OUTSIDE the interval [-15.78, 8.84], the slope of Illiteracy
## is p < .05.
## 
## Note: The range of observed values of Murder is [1.40, 15.10]
## 
## Interval calculated using false discovery rate adjusted t = 2.35 
## 
## SIMPLE SLOPES ANALYSIS
## 
## Slope of Illiteracy when Murder =  5.420973 (- 1 SD): 
## 
##     Est.     S.E.   t val.      p
## -------- -------- -------- ------
##   -17.43   250.08    -0.07   0.94
## 
## Slope of Illiteracy when Murder =  8.685043 (Mean): 
## 
##      Est.     S.E.   t val.      p
## --------- -------- -------- ------
##   -399.64   178.86    -2.23   0.03
## 
## Slope of Illiteracy when Murder = 11.949113 (+ 1 SD): 
## 
##      Est.     S.E.   t val.      p
## --------- -------- -------- ------
##   -781.85   189.11    -4.13   0.00

In this case, you can see that the interval is just a little bit wider. The output reports the adjusted test statistic, which is 2.33, not much different than the (approximately) 2 that would be used otherwise. In other cases it may be quite a bit larger.

Additional options

sim_slopes(fiti, pred = Illiteracy, modx = Murder, cond.int = TRUE)

## JOHNSON-NEYMAN INTERVAL
## 
## When Murder is OUTSIDE the interval [-6.37, 8.41], the slope of Illiteracy
## is p < .05.
## 
## Note: The range of observed values of Murder is [1.40, 15.10]
## 
## SIMPLE SLOPES ANALYSIS
## 
## When Murder =  5.420973 (- 1 SD): 
## 
##                                  Est.     S.E.   t val.      p
## --------------------------- --------- -------- -------- ------
## Slope of Illiteracy            -17.43   250.08    -0.07   0.94
## Conditional intercept         4618.50   229.76    20.10   0.00
## 
## When Murder =  8.685043 (Mean): 
## 
##                                  Est.     S.E.   t val.      p
## --------------------------- --------- -------- -------- ------
## Slope of Illiteracy           -399.64   178.86    -2.23   0.03
## Conditional intercept         5097.72   239.69    21.27   0.00
## 
## When Murder = 11.949113 (+ 1 SD): 
## 
##                                  Est.     S.E.   t val.      p
## --------------------------- --------- -------- -------- ------
## Slope of Illiteracy           -781.85   189.11    -4.13   0.00
## Conditional intercept         5576.95   340.70    16.37   0.00

sim_slopes(fiti, pred = Illiteracy, modx = Murder, robust = "HC3")

## JOHNSON-NEYMAN INTERVAL
## 
## When Murder is OUTSIDE the interval [-3.34, 8.14], the slope of Illiteracy
## is p < .05.
## 
## Note: The range of observed values of Murder is [1.40, 15.10]
## 
## SIMPLE SLOPES ANALYSIS
## 
## Slope of Illiteracy when Murder =  5.420973 (- 1 SD): 
## 
##     Est.     S.E.   t val.      p
## -------- -------- -------- ------
##   -17.43   227.37    -0.08   0.94
## 
## Slope of Illiteracy when Murder =  8.685043 (Mean): 
## 
##      Est.     S.E.   t val.      p
## --------- -------- -------- ------
##   -399.64   158.77    -2.52   0.02
## 
## Slope of Illiteracy when Murder = 11.949113 (+ 1 SD): 
## 
##      Est.     S.E.   t val.      p
## --------- -------- -------- ------
##   -781.85   156.96    -4.98   0.00

References

Bauer, D. J., & Curran, P. J. (2005). Probing interactions in fixed and multilevel regression: Inferential and graphical techniques. Multivariate Behavioral Research, 40, 373–400.

Esarey, J., & Sumner, J. L. (2017). Marginal effects in interaction models: Determining and controlling the false positive rate. Comparative Political Studies, 1–33. https://doi.org/10.1177/0010414017730080

Fox, J., & Weisberg, S. (2018). Visualizing fit and lack of fit in complex regression models with predictor effect plots and partial residuals. Journal of Statistical Software, 87, 1–27. https://doi.org/10.18637/jss.v087.i09

Hainmueller, J., Mummolo, J., & Xu, Y. (2017). How much should we trust estimates from multiplicative interaction models? Simple tools to improve empirical practice. SSRN Electronic Journal. https://doi.org/10.2139/ssrn.2739221

Johnson, P. O., & Fay, L. C. (1950). The Johnson-Neyman technique, its theory and application. Psychometrika, 15, 349–367.