How to Create an APA Style Interaction Graph in R (Part 2)

When running regression analyses, analyzing the interactions between predictor variables is an essential step. Examining p-values and coefficients is not enough, though, we also need to plot the lines to really understand what is happening between an independent variable and a moderator.

Following up on Part 1 of this topic, the code below allows you to create a simple APA-style interaction graph in R with two predictors (a 2-way interaction). This example works with a continuous DV, a continuous IV, and a two-level categorical moderator (factor). Please see Part 1 if you have two continuous predictor variables.

##
# Generate an APA compatible interaction plot with a:
# * continous DV
# * continous IV
# * categorical moderator 
##

library(tidyverse)
library(ggeffects)
library(jtools)
library(forcats)

# Declare data frame
# Here we reduce species to two levels
df <- iris |> 
  filter(Species %in% c("setosa", "versicolor")) |> 
  droplevels()

# Declare names of variables and levels in data frame
dv_var <- "Sepal.Length"
iv_var <- "Petal.Length"
mod_var <- "Species"     # name of mod factor variable
mod_lvl_1 <- "setosa"    # declare the first mod level
mod_lvl_2 <- "versicolor"# declare the second mod level

# Declare plot labels -- this is how they will appear in the plot
dv_label <- "Sepal Length"
iv_label <- "Petal Length"
mod_label <- "Species"
mod_lvl_label_1 <- "Setosa"
mod_lvl_label_2 <- "Versicolor"

# Declare location of legend
legend_pos <- "bottomright"

##
# The rest of the code should not be modified unless you want
# to customize
##

# Rename variables in df
df <- df |>
  rename(dv = all_of(dv_var),
         iv = all_of(iv_var),
         modv = all_of(mod_var))

# Generate your model, using the new var names
model <- lm(dv ~ iv * modv, data = df)

# Calculate one standard deviation above and below
# the predictor and moderator
cond1 <- mean(df$iv) - sd(df$iv)
cond2 <- mean(df$iv) + sd(df$iv)

mods <- c(mod_lvl_1, mod_lvl_2)
conds <- c(cond1, cond2)

# Generate the dependent variable predictions. These
# will be our y-axis values
predictions <- ggpredict(model,
                         terms = c("iv [conds]",
                                   "modv [mods]"))

# Update the data frame with the variables needed
# to build the plot
plot_vars <- predictions |>
  mutate(mod_order = ifelse(group == mod_lvl_1, 1, 2)) |>
  mutate(iv1_order = ifelse(x == cond1, 1, 2)) |>
  mutate(mod = ifelse(group == mod_lvl_1,
                      paste(mod_lvl_label_1),
                      paste(mod_lvl_label_2))) |>
  mutate(iv1 = ifelse(x == cond1,
                      paste("Low", iv_label),
                      paste("High", iv_label)))

# Generate the plot
p <- ggplot(plot_vars,
            aes(x = fct_reorder(iv1, iv1_order),
                y = predicted,
                shape = fct_reorder(mod, mod_order),
                group = fct_reorder(mod, mod_order))) +
  geom_point(size = 2) +
  geom_line(aes(linetype = fct_reorder(mod, mod_order))) +
  labs(x = '', y = dv_label) +
  ylim(0,NA) +
  theme_apa(legend.pos = all_of(legend_pos)) +
  scale_colour_grey() +
  theme(legend.background = element_rect(fill = "white",
                                         colour = "black"),
        axis.text = element_text(size=12, colour = "black"))

# Output the plot
p

For a simple 2-way interaction graph, the above code only requires declaring:

1. the data frame,
2. variable names,
3. moderator levels,
4. labels to be used in the plot, and optionally,
5. the legend position.

Of course, there are other combinations of IVs and moderators that you may need to plot. Hopefully, this post and Part 1 will give you a head start to create journal ready APA interaction plots.