---
title: "Exploring interactions with continuous predictors in regression models"
author: "Jacob Long"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Exploring interactions with continuous predictors in regression models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r echo=FALSE}
required <- c("survey")
if (!all(sapply(required, requireNamespace, quietly = TRUE)))
  knitr::opts_chunk$set(eval = FALSE)
knitr::opts_chunk$set(message = F, warning = F, fig.width = 6, fig.height = 5)
library(jtools)
```

Understanding an interaction effect in a linear regression model is usually
difficult when using just the basic output tables and looking at the
coefficients. The `jtools` package provides several functions that can help
analysts probe more deeply.

**Categorical by categorical interactions**: All the tools described here 
require at least one variable to be continuous. A separate vignette describes
`cat_plot`, which handles the plotting of interactions in which all the focal
predictors are categorical variables. 

## Plotting interactions

First, we use example data from `state.x77` that is built into R. Let's look 
at the interaction model output with `summ` as a starting point.

```{r}
fiti <- lm(Income ~ Illiteracy * Murder, data = as.data.frame(state.x77))
summ(fiti)
```

So we see a significant main effect of `Murder` on `Income` in the presence of 
a significant interaction between `Murder` and `Illiteracy`. The positive 
estimate for `Illiteracy` does not differ significantly from zero. With that 
said, you shouldn't focus too much on the main effects of terms included in 
the interaction.

Note that if you would like to see the output if the input variables were 
mean-centered and/or standardized, `summ` can do that for you.

```{r}
summ(fiti, scale = TRUE)
```

A versatile and sometimes the most interpretable method for 
understanding interaction effects is via plotting. `jtools` provides 
`interact_plot` as a relatively pain-free method to get good-looking plots of
interactions using `ggplot2` on the backend.

```{r}
interact_plot(fiti, pred = "Illiteracy", modx = "Murder")
```

Keep in mind that the default behavior of `interact_plot` is to mean-center 
all variables not involved in the interaction so that the predicted values are 
more easily interpreted. You can disable that by adding `centered = "none"`.

By default, with a continuous moderator you get three lines---1 standard 
deviation above and below the mean and the mean itself. If you specify 
`modxvals = "plus-minus"`, the mean of the moderator is not plotted, just the
two +/- SD lines. 

```{r}
interact_plot(fiti, pred = "Illiteracy", modx = "Murder", modxvals = "plus-minus")
```

However, if your moderator is **a factor**, each level will be plotted and you 
should leave `modxvals = NULL`, the default.

```{r}
fitiris <- lm(Petal.Length ~ Petal.Width * Species, data = iris)
interact_plot(fitiris, pred = "Petal.Width", modx = "Species")
```

### 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.

```{r}
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 mdoerator variable. 

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

```{r}
interact_plot(fitiris, pred = "Petal.Width", modx = "Species", plot.points = 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.

```{r}
fiti <- lm(Income ~ Illiteracy * Murder, data = as.data.frame(state.x77),
           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.

### 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.

```{r}
interact_plot(fiti, pred = "Illiteracy", modx = "Murder", interval = TRUE, int.width = 0.8)
```

### 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.

```{r}
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",
              color.class = "Oranges")
```

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`:

```{r}
interact_plot(fitiris, pred = "Petal.Width", modx = "Species") + theme_apa()
```

## Simple slopes analysis and Johnson-Neyman intervals

Simple slopes analysis gives researchers a way to express the interaction 
effect in terms that are easy to understand to those who know how to interpret
direct effects in regression models. This method is designed for **continuous
variable by continuous variable** interactions, but can work when the 
moderator is binary.

In simple slopes analysis, researchers are interested in the *conditional*
slope of the focal predictor; that is, what is the slope of the predictor when
the moderator is held at some particular value? The regression output we get
when including the interaction term tells us what the slope is when the
moderator is held at zero, which is often not a practically/theoretically
meaningful value. To better understand the nature of the interaction, simple
slopes analysis allows the researcher to specify meaningful values at which to
hold the moderator value.

While the computation behind doing so isn't exactly rocket science, it is 
inconvenient and prone to error. The `sim_slopes` function from `jtools` 
accepts a regression model (with an interaction term) as an input and 
automates the simple slopes procedure. The function will, by default, do the 
following:

* Mean-center all non-focal predictors (so that setting them to zero means
setting them to their mean)
* For continuous moderators, it will choose the mean as well as the mean 
plus/minus 1 standard deviation as values at which to find the slope of the
focal predictor.
* For categorical/binary moderators, it will find the slope of the focal
predictor at each level of the moderator.

In its most basic use case, `sim_slopes` needs three arguments: a linear model 
(with support for `svyglm` models), the name of the focal predictor in 
quotations as the argument for `pred =`, and the name of the moderator in 
quotations as the argument for `modx =`. Let's go through an example.

Don't sweat the change in the significance of the lower-order focal
variables---they  aren't directly interpretable, anyway.

Now let's do the most basic simple slopes analysis:

```{r}
sim_slopes(fiti, pred = Illiteracy, modx = Murder, johnson_neyman = FALSE)
```

So what we see in this example is that when the value of `Murder` is high, the
slope of `Illiteracy` is negative and significantly different from zero. The 
value for `Illiteracy` when `Murder` is high is in the opposite direction from 
its coefficient estimate for the first version of the model fit with `lm` but 
this result makes sense considering the interaction coefficient was negative; 
it means that as one of the variables goes up, the other goes down. Now we 
know the effect of `Illiteracy` only exists when `Murder` is high.

You may also choose the values of the moderator yourself with the
`modxvals =` argument.

```{r}
sim_slopes(fiti, pred = Illiteracy, modx = Murder, modxvals = c(0, 5, 10),
           johnson_neyman = FALSE)
```

### 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.

```{r}
sim_slopes(fiti, pred = Illiteracy, modx = Murder, johnson_neyman = TRUE)
```

So 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. 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`.

```{r}
johnson_neyman(fiti, pred = Illiteracy, modx = Murder, alpha = 0.01)
```

One 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 9.12. 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 `jtools`, not them.

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

```{r}
sim_slopes(fiti, pred = Illiteracy, modx = Murder, johnson_neyman = TRUE,
           control.fdr = TRUE)
```

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

#### Conditional intercepts

Sometimes it is informative to know the conditional intercepts in addition to 
the slopes. It might be interesting to you that individuals low on the 
moderator have a positive slope and individuals high on it don't, but that 
doesn't mean that individuals low on the moderator will have higher values of 
the dependent variable. You would only know that if you know the conditional 
intercept.

You can print the conditional intercepts with the `cond.int = TRUE` argument.

```{r}
sim_slopes(fiti, pred = Illiteracy, modx = Murder, cond.int = TRUE)
```

This example shows you that while the slope associated with `Illiteracy` is 
negative when `Murder` is high, the conditional intercept is also high when 
`Murder` is high. That tells you that increases in `Illiteracy` for 
high-`Murder` observations will tend towards being equal on `Income` to 
observations with lower values of `Murder`.

#### Robust standard errors

Certain models require heteroskedasticity-robust standard errors. To be 
consistent with the reporting of heteroskedasticity-robust standard errors 
offered by `summ`, `sim_slopes` will do the same with the use of the
`robust = TRUE` option so you can consistently report standard errors across
models.

```{r}
sim_slopes(fiti, pred = Illiteracy, modx = Murder, robust = TRUE)
```

These data are a relatively rare case in which the robust standard errors are
even more efficient than typical OLS standard errors. Note that you must have 
the `sandwich` package installed to use this feature.

#### Centering and standardizing variables

By default, all non-focal variables are mean-centered. You can additionally
specify that the focal predictors are centered with `centered = "all"`. You 
may also request that no variables be centered with `centered = "none"`. You 
may also request specific variables to center by providing a vector of quoted 
variable names---no others will be centered in that case.

Note that the moderator is centered around the specified values. Factor 
variables are ignored in the centering process.

You can standardize the variables as well using the `scale = TRUE` 
argument.

```{r}
sim_slopes(fiti, pred = Illiteracy, modx = Murder, scale = TRUE, 
           centered = "all")
```

Standardization is only meaningful if applied to the focal predictors. If you 
use `scale = TRUE` when `centered = NULL`, the output won't change since 
the scale of the other predictors has no direct effect on the others.

## Do simple slopes and plotting with one function

You won't have to use these functions long before you may find yourself using
both of them for each model you want to explore. To streamline the process,
this package offers `probe_interaction()` as a convenience function that calls
both `sim_slopes()` and `interact_plot()`, taking advantage of their 
overlapping syntax.

```{r}
library(survey)
data(api)
dstrat <- svydesign(id = ~1, strata = ~stype, weights = ~pw,
                    data = apistrat, fpc = ~fpc)
regmodel <- svyglm(api00 ~ avg.ed * growth, design = dstrat)

probe_interaction(regmodel, pred = growth, modx = avg.ed, cond.int = TRUE,
                  interval = TRUE,  jnplot = TRUE)
```

Note in the above example that you can provide arguments that only apply to one
function and they will be used appropriately. On the other hand, you cannot
apply their overlapping functions selectively. That is, you can't have one 
`scale = TRUE` and the other `scale = FALSE`. If you want that 
level of control, just call each function separately.

Also, the above example comes from the survey package as a means to show that,
yes, these tools can be used with `svyglm` objects, though should only be 
applied to linear models.

It returns an object with each of the two functions' return objects:

```{r}
out <- probe_interaction(regmodel, pred = growth, modx = avg.ed,
                         cond.int = TRUE, interval = TRUE, jnplot = TRUE)
names(out)
```

## 3-way interactions

If 2-way interactions can be hard to grasp by looking at regular regression 
output, then 3-way interactions are outright inscrutable. The aforementioned 
functions also support 3-way interactions, however. Plotting these effects is 
particularly helpful.

Note that Johnson-Neyman intervals are still provided, but only insofar as you
get the intervals for chosen levels of the second moderator. This does go 
against some of the distinctiveness of the J-N technique, which for 2-way 
interactions is a way to avoid having to choose points of the moderator to 
check whether the predictor has a significant slope.

```{r}
fita3 <- lm(rating ~ privileges * critical * learning, data = attitude)
probe_interaction(fita3, pred = critical, modx = learning, mod2 = privileges)
```

You can change the labels for each plot via the `mod2.labels` argument.

`interact_plot` has a bit more flexibility here than `sim_slopes`, allowing
for factor moderators with greater than 2 levels. And don't forget that you can
use `theme_apa` to format for publications or just to make more economical
use of space.

```{r}
mtcars$cyl <- factor(mtcars$cyl,
                     labels = c("4 cylinder", "6 cylinder", "8 cylinder"))
fitc3 <- lm(mpg ~ hp * wt * cyl, data = mtcars)
interact_plot(fitc3, pred = hp, modx = wt, mod2 = cyl) + 
  theme_apa(legend.pos = "bottomright")
```


You can get Johnson-Neyman plots for 3-way interactions as well, but keep in
mind what I mentioned earlier in this section about the J-N technique for 3-way
interactions. You will also need the `cowplot` package, which is used on the 
backend to mush together the separate J-N plots.

```{r fig.height = 8}
regmodel3 <- svyglm(api00 ~ avg.ed * growth * enroll, design = dstrat)
sim_slopes(regmodel3, pred = growth, modx = avg.ed, mod2 = enroll,
          jnplot = TRUE)
```

Notice that at one of the three values of the second moderator, there were no
Johnson-Neyman interval values so it wasn't plotted. The more levels of the 
second moderator you want to plot, the more likely that the resulting plot
will be unwieldy and hard to read. You can resize your window to help, though.

## Generalized linear models, mixed models, et al.

`interact_plot` is designed to be as general as possible and has been tested
with `glm`, `svyglm`, and `merMod` models. When dealing with generalized
linear models, it can be immensely useful to get a look at the predicted 
values on their response scale (e.g., the probability scale for logit models).

To give an example of how such a plot might look, I'll generate some example 
data.

```{r}
set.seed(5)
x <- rnorm(100)
m <- rnorm(100)
prob <- boot::inv.logit(.25 + .3*x + .3*m + -.5*(x*m) + rnorm(100))
y <- rep(0, 100)
y[prob >= .5] <- 1
logit_fit <- glm(y ~ x * m, family = binomial)
```

Here's some summary output, for reference:

```{r}
summ(logit_fit)
```

Now let's plot our logit model's interaction:

```{r}
interact_plot(logit_fit, pred = x, modx = m)
```

A beautiful transverse interaction with the familiar logistic curve.

## 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](https://doi.org/10.1177/0010414017730080)

Hainmueller, J., Mummolo, J., & Xu, Y. (2016). 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](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.

