How to visualize interaction effects

I recently gave a workshop on data visualization. One of the topics was visualizing interaction effects (or “moderation” analyses). I think it is a great topic because it exemplifies quite well that there is not one s…

I recently gave a workshop on data visualization. One of the topics was visualizing interaction effects (or “moderation” analyses). I think it is a great topic because it exemplifies quite well that there is not one solution to all problems. In general, I would argue that trying to visualize interaction effects is great idea. Due to the conditional nature of the effects obtained from standard regression analyses that include an interaction term, it is often hard to understand the direction and size of an interaction effect only based on regression results tables. In this post, I simply want to share some code that can be used to visualize a variety of different interaction effects.

Visualizing a 2×2 experimental design

# Loading packages
library(tidyverse)   # includes 'ggplot2' which is probably the best package for data viz
library(ggExtra)     # extends ggplot2 

# Simulating the data
set.seed(100)
m1 <- sample(rep(c(0, 1), 250), 500) # Factor A
m2 <- sample(rep(c(0, 1), 250), 500) # Factor B
y <- 8 + -2*m1 + 2*m2 + 3*(m1*m2) + rnorm(500, 0, 2)
(data <- data.frame(m1 = factor(m1, 
                                label = c("A1", "A2")), 
                    m2 = factor(m2,
                                label = c("B1", "B2")),
                    y) %>%
  as_tibble)

#> # A tibble: 500 × 3
#>    m1    m2        y
#>    <fct> <fct> <dbl>
#>  1 A2    B2    11.5 
#>  2 A2    B2     9.33
#>  3 A2    B1     7.78
#>  4 A1    B2    13.2 
#>  5 A2    B2    11.2 
#>  6 A2    B2     7.93
#>  7 A1    B1     8.36
#>  8 A2    B1     6.52
#>  9 A2    B2    12.7 
#> 10 A2    B2     9.79
#> # ℹ 490 more rows

As we can see, we have a dataset that resembles data obtained from a 2×2 experimental design: m1 and m2 represent the experimental manipulations A and B resulting in four conditions (A1, A2, B1, B2). In such a case, where we are interested in the effects of both conditions on the dependent variable y and in the effect of a potential interaction, a simple bar chart with a specific type of facetting is probably the best solution.

ggplot(data, 
       aes(x = m1,
           fill = m1,  
           y = y)) +
  stat_summary(fun = mean,
               geom = "bar") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="errorbar", 
               width = 0.25) +
  ylim(0, 17) +
  facet_wrap(~m2) +
  labs(x = "Experimental Manipulations",
       y = "Dependent Variable (M)") +
  theme_minimal() +
  scale_fill_brewer(palette = "Blues") +
  theme(legend.position="none")

An alternative is a modified line plot.

# Summarizing data
data2 <- data %>% 
  group_by(m1, m2) %>% 
  summarise(y_mean = mean(y),
            y_se = psych::describe(y)$se)

# Creating the plot
data2 %>% 
  ggplot(aes(x = m1, 
             y = y_mean, 
             color = m2)) +
  geom_line(aes(group = m2)) +
  geom_point() +
  geom_errorbar(aes(ymin = y_mean-1.96*y_se, 
                      ymax = y_mean+1.96*y_se),
                width = .1) +
  ylim(0, 17) +
  labs(x = "Factor A",
       color  = "Factor B",
       y = "Dependent Variable") +
  theme_minimal() +
  scale_color_brewer(palette = "Blues")

Visualizing the moderation effect of a dichotomous variable

# Simulating the data
x <- rnorm(500, 2, 3)
m <- rep(c(0,1), 250)
y <- 5 + -4*x + 2*m + 4*(x*m) + rnorm(500, 0, 5)
(data <- data.frame(x, m = factor(m), y) %>%
  as_tibble)

#> # A tibble: 500 × 3
#>         x m          y
#>     <dbl> <fct>  <dbl>
#>  1  2.40  0     -3.95 
#>  2  6.26  1      4.15 
#>  3  2.56  0      4.38 
#>  4  2.32  1      0.932
#>  5  2.27  0     -6.47 
#>  6  2.05  1      5.34 
#>  7  0.267 0     -3.32 
#>  8 -0.157 1     15.0  
#>  9  0.843 0     -2.44 
#> 10  2.94  1     -0.957
#> # ℹ 490 more rows

The moderator m here is a factor variable. The independent variable x and the dependent variable y, however, are continuous.

When the relationship of interested is between to continuous variables and only the moderator is dichotomous (e.g., gender, an experimental manipulation, intervention vs. no intervention,…), we can try to show as much information as possible and visualize the relationships fully using a colored scatterplot with separate regression lines for each group.

# Setting up the building blocks
basic_plot <- ggplot(data,
       aes(x = x,
           y = y,
           color = m)) +
  theme_bw() +
  labs(x = "Independent variable",
       y = "Dependent variable",
       color = "Moderator")

# Colored scatterplot
basic_plot +
  geom_point()

# Colored scatterplot and regression lines
basic_plot +
  geom_point(alpha = .3, 
             size = .9) +
  geom_smooth(method = "lm")

ggsave("test2.png",
       height = 4,
       width = 5)

# Only regression lines
basic_plot + 
  geom_smooth(method = "lm",
              se = F)

Note: Many papers only show the last plot. However, I would always prefer the second plot which shows the actual scatterplots for each groups, the regression lines and confidence bounds.

Visualizing the moderation effect of metric variable

# Simulating the data
n <- 500
x <- runif(n)
m <- rnorm(n, 0.5, 1)
y <- 2*x + -2*m + -6*(x*m) + rnorm(500, 0, 5.88)
(data <- data.frame(x, m, y) %>%
  as_tibble)

#> # A tibble: 500 × 3
#>        x      m        y
#>    <dbl>  <dbl>    <dbl>
#>  1 0.595  0.524   3.46  
#>  2 0.679  1.23    0.224 
#>  3 0.676  0.200   2.95  
#>  4 0.904  1.90  -10.5   
#>  5 0.174 -1.05   13.9   
#>  6 0.335  0.873   0.136 
#>  7 0.674  1.43  -11.0   
#>  8 0.283 -1.47    6.16  
#>  9 0.252  2.03  -11.2   
#> 10 0.294  0.902   0.0474
#> # ℹ 490 more rows

If the moderator (in this case m) is also continuous, visualizing interaction effects becomes more complex. There are generally three ways of presenting results from such a moderation analysis:

1. Estimating the regression for different values of the moderator (e.g., -SD, M, +SD) and reporting all three results in a table
2. Visualizing three regression lines for three (arbitrary) groups, also know as the “pick-a-point” approach
3. Plotting the conditional effects for all values of the moderator (Conditional plot)

Before we build these plots in R, here is an interactive version of approaches 2 and 3. Change the size of the interaction and the sample size, and drag the moderator to see where the effect of x becomes significant:

Interactive: simple slopes and the conditional effect, estimated on simulated data in your browser.

Pick-a-point approach

For this approach, we have to manually recode the moderator variable into three meaningful subgroups. Cutting points could be -SD, M, +SD or certain quantiles. In this case, we simply use the cut function to create subgroups¹

data$m_groups <- cut(data$m, breaks = 3) %>% 
  factor(., labels = c("small", "average", "large"))
table(data$m_groups)

#> 
#>   small average   large 
#>      71     310     119

Next, we simply repeat the procedure that we already used for visualizing interaction effects with a dichotomous moderator. We use the newly computed variable to color the points in the scatterplot and create separate regression lines.

ggplot(data,
       aes(x = x,
           y = y,
           color = m_groups)) +
  geom_point(size = .9,
             alpha = .3) +
  geom_smooth(method = "lm") +
  theme_bw() +
  scale_color_brewer(type = "qual", 
                     palette = 3) +
  labs(x = "Independent variable",
       y = "Dependent variable",
       color = "Moderator")

Conditional effect plot

A problem with the pick-a-point approach is, however, that we only see the conditional effects for three (somewhat arbitrary) groups. A conditional plot which visualizes the conditional effects for all possible values of the moderator is much more valuable.

For this approach, we first need to specify the interaction model itself.

model <- lm(y ~ x + m + x:m, data)
summary(model)

#> 
#> Call:
#> lm(formula = y ~ x + m + x:m, data = data)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -18.2856  -3.9007   0.2177   4.0035  18.0900 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -0.09295    0.61490  -0.151  0.87991    
#> x            2.69892    1.11273   2.425  0.01564 *  
#> m           -1.80938    0.56348  -3.211  0.00141 ** 
#> x:m         -6.67974    0.99706  -6.699 5.69e-11 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 6.148 on 496 degrees of freedom
#> Multiple R-squared:  0.423,  Adjusted R-squared:  0.4195 
#> F-statistic: 121.2 on 3 and 496 DF,  p-value: < 2.2e-16

Now, we have to create a function that uses the output of this model and creates a data frame that includes the conditional effects for several values of the moderator.

# Function that estimates the conditinal effects for certain values of the moderator
conditional.effects <-  function(model, x, m, quantiles = 10){
  interact = paste0(x,':',m)
  beta.hat = coef(model) 
  covs = vcov(model)
  z0  = quantile(model$model[,m], seq(0 , 1, 1/quantiles))
  dy.dx = beta.hat[x] + beta.hat[interact]*z0
  se.dy.dx = sqrt(covs[x, x] + z0^2*covs[interact, interact] + 2*z0*covs[x, interact])
  upr = dy.dx+1.96*se.dy.dx
  lwr = dy.dx-1.96*se.dy.dx
  data.frame(m = z0, b = dy.dx, lwr, upr)
}

(con_effects <- conditional.effects(model, 
                                   x = "x", 
                                   m = "m", 
                                   quantiles = 500) %>%
  as_tibble)

#> # A tibble: 501 × 4
#>        m     b   lwr   upr
#>    <dbl> <dbl> <dbl> <dbl>
#>  1 -2.32  18.2 12.2   24.2
#>  2 -2.27  17.9 12.0   23.8
#>  3 -2.04  16.3 10.9   21.8
#>  4 -1.75  14.4  9.48  19.3
#>  5 -1.74  14.3  9.42  19.2
#>  6 -1.67  13.9  9.08  18.6
#>  7 -1.57  13.2  8.60  17.8
#>  8 -1.53  12.9  8.41  17.5
#>  9 -1.47  12.5  8.10  16.9
#> 10 -1.44  12.3  7.96  16.7
#> # ℹ 491 more rows

The dataframe includes the values of the moderator (m), the respective conditional effects (b) and their 95% confidence intervals. We can now use this data frame to create a new plot.

(mod_plot <- ggplot(con_effects, 
       aes(x = m, 
           y = b, 
           ymin = lwr,
           ymax = upr)) +
  geom_smooth(stat = "identity",
              color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  theme_bw() +
  labs(x = "Moderator",
       y = "Conditional effect of x on y"))

This plot now shows the conditional effects for each value of the moderator. We can see that the effect of x on y is positive for lower values of m (the moderator) and turns negative for higher values of m. We can also asses at which values of m, the conditional effect of x on y becomes non-significant (Johnson-Neyman technique: see also https://rpubs.com/bachl/jn-plot).

I would suggest to also plot a histogram of the moderator so that one can assess how many values provide information for the information of the conditional effects.

# ggExtra now needs a point layer, so we add invisible points (alpha = 0)
ggMarginal(mod_plot + geom_point(alpha = 0), # function of the package `ggExtra`
           margins = c("x"),
           type = "histogram",
           color = "white",
           fill = "darkgrey",
           size = 3)

Hope this is helpful to some people! 😉

Note5 comments from the original post

These comments were left on the earlier WordPress version of this post.

Rosie · 6 October 2020

Thank you so much!

Daniel · 8 June 2021

This is awesome! Is there also the code for visualizing poisson interaction effects?

Jordan · 8 September 2022

Whether or not you can visualize poisson effects might depend on the type of variables you use (ie, categorical or continuous) for both x and m. For two continuous variables the ‘Conditional effect plot’ model should work.

I successfully created the ‘Conditional effect plot’ visual with a glm model specified as a ‘poisson’ model, ie.:

model <-glm(y ~ x + m + x:m, family = “poisson”, data = dataframe) [1]

Whether or not the model is true to the actual effects of the moderator, I am unsure; I would think the ‘conditional.effects’ function creates conditional effects which can be used without inferential errors for a ‘poisson’ specified glm model with an interaction effect similar to a non-specified linear model (lm) with an interaction effect because I believe the variables requested by the ‘conditional.effects’ function come from the same objects within the model, ie, coef(model). If both the glm and lm models store the requested variables in the same locations, the conditional effects (represented in the ‘con_effects’ data frame) should provide accurate estimations; that is, unless of course, users must ‘correct’ the variables requested in the ‘conditional.effects’ for ‘poisson’ specified generalized linear models. I’m curious to see if the author will respond to this point.

Assuming users do not need to correct for any mathematical or to be clear, inferential, based errors, the ‘family’ argument set to ‘poisson’ within the glm function should run your model under the ‘poisson’ parametric. You can use the code Philipp provides under the ‘# Function that estimates the conditi(o)nal effects for certain values of the moderator’ with the poisson specified model to get the conditional effect model.

For more information on the different types of models you can specify within a generalized linear model (glm) I recommend:

The page provides examples of the functions with different specifications (similar to my examples ([1]), above). Just make sure to ignore the link part of the argument within the argument (As I did, and as they do in their examples! :)).

Lucas Leonhard · 7 March 2022

I really love the pick a point approach! However, is it possible to use the estimates of a previous estimated model? How would the code look like, if we would call the previous estimates model “model1”?

Jordan · 8 September 2022

A quick word of advice: when running the conditional effects function (created under the “# Function that estimates the conditi(o)nal effects for certain values of the moderator” notation) you’ll want to be sure both your x and m variable are numerics within your original dataset. If both your x and m variables are not numerics (ie, you may run one of them ‘as a numeric’ with the ‘as.numeric’ function within your model) the function will not work. To ensure your x and m variables are numeric you can simply run:

dataframe$x <- as.numeric(dataframe$x)
dataframe$m <- as.numeric(dataframe$m)

This will ensure both your x and m values are coded as numbers within your dataset (it should also ensure you don’t have to run the ‘as.numeric’ function within your model (ie, model <- lm(y ~ x + m + x:m)) which is helpful! 🙂

Footnotes

1. For this function, we have to specify the breaks. When breaks is specified as a single number, the range of the data is divided into breaks pieces of equal length, and then the outer limits are moved away by 0.1% of the range to ensure that the extreme values both fall within the break intervals.