Purpose

Create a function that plots three-way interactions with estimated marginal means.

Motivation

When dealing with experimental data, you may need to create a plot that shows a three-way interaction. If one or more of the variables in your data are collected, not experimentally manipulated, it is recommended that you use estimated marginal means instead of simple averages. I was unable to find an easy way to plot these three-ways interactions with estimated marginal means, so I created a function to do so.

Example

Using the function plot_threeways_cont and about 8-lines of code, you will be able to plot a graph that looks like this:

plot_threeways_cont(df = raw, 
              y = raw$y_var, 
              x1 = raw$x_var1, 
              x2 =  raw$x_var2, 
              x3 = raw$x_var3, 
              title = "Graph Title", 
              xlabel = "X Variable Label", 
              ylabel = "Y Variable Label")

Background

To build the function, here are the packages you will need:

library(ggplot2)     # for plotting 
library(dplyr)       # for manipulating the data
library(emmeans)     # for calculating the estimated marginal means

If you don’t have them, here’s the code to install them:

install.packages(c('ggplot2', 'dplyr', 'emmeans'))

The Function

graph_threeway <- function (df, dependent_var, x1, x2, x3, title, xlabel, ylabel, output=TRUE) { 
  
  model_output <- lm(data = df,
               dependent_var ~ x1 * x2 * x3) #Estimate the regression model
  
  # Remove Quotes from Strings: 
  title <- as.symbol(title) 
  xlabel <- as.symbol(xlabel)
  ylabel <- as.symbol(ylabel)

  # Use emmip function to build the graph: 
  
    plot_test <- emmip(model_output,  x1 ~ x2 | x3, CI=TRUE) + 
    labs(title= deparse(substitute(title)), 
         x=(deparse(substitute(xlabel))),
         y=(deparse(substitute(ylabel)))) +
    theme_classic()  + 
    theme(legend.title=element_blank()) #Remove legend name
    
    if(output==TRUE) {

      ggsave(paste(deparse(substitute(title)), ".pdf", sep=""),
             plot = last_plot(),
             width = 8, height = 6,
             units = "in")
          return(plot_test)

    } 
    
  else {
    return(plot_test)
  }
}

The function takes the following arguments:

Note: Make sure that x2, and x3 are both coded as afactor. This function will not work if you try to use it on numeric variables.

Feel free to leave comments with questions or suggestions on how to make this function better. If you know of a function that already accomplishes the same purpose, let me know. I would love to check it out.