Introduction

Chapter 13 and 14 cover interactions (also known as moderation), including their use, interpretation, and limitations. The use in R is straightforward as we’ll highlight below.

Interactions (Moderation)

  1. Let’s start by loading the tidyverse package, the furniture package, and the educ7610 package.
library(tidyverse)
library(furniture)
library(educ7610)
library(jtools)
library(interactions)
  1. Import the hosp data set into R.
data(hosp)
  1. Because when a data set was imported from SPSS (as this one originally was), there are labelled variable types, we need to adjust these. In this data set, all variables can be re-classified as numeric so we can use the code below to do this.
hosp <- hosp %>%
  map_df(~as.numeric(.x))
  1. Interactions are often either theory-driven or data-driven. Of the two, theory-driven is often the most respected reason to test and report on an interaction. Herein, we’ll explore some relationships with visualizations that we’ll then test (since we don’t have any good theories about this data set).
hosp %>%
  ggplot(aes(exhaust, safety, group = sex, color = sex)) +
    geom_point(alpha = .6) +
    geom_smooth(method = "lm")

  1. There looks like there could be an interaction between exhaust and sex on safety. To test this, let’s run a linear regression model with the cross product of safety and sex. R is smart enough that if we put exhaust * sex in the formula, it will include all the right effects for us.
hosp %>%
  lm(safety ~ exhaust * sex, data = .) %>%
  summary()
## 
## Call:
## lm(formula = safety ~ exhaust * sex, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.23439 -0.61719  0.00054  0.74033  2.61092 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.32847    0.27543  12.085  < 2e-16 ***
## exhaust      0.13674    0.07956   1.719  0.08673 .  
## sex         -0.45969    0.40175  -1.144  0.25345    
## exhaust:sex  0.31860    0.11456   2.781  0.00576 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9816 on 296 degrees of freedom
## Multiple R-squared:  0.1812, Adjusted R-squared:  0.1729 
## F-statistic: 21.83 on 3 and 296 DF,  p-value: 8.46e-13

Interpretation

  1. If exhaust is how exhausted the individual feels, safety is the amount the individual works around safety protocols (avoids performing the full safety protocol) and sex is the biological sex of the individual (0 = female, 1 = male), what is the interpretation of each coefficient?
  2. Does this interpretation match the figure we created before? Why or why not?
  3. Let’s include some covariates in the model; namely, tenure (the job experience the individual has) and age (the individual’s age in years).
hosp %>%
  lm(safety ~ exhaust * sex + tenure + age, data = .) %>%
  summary()
## 
## Call:
## lm(formula = safety ~ exhaust * sex + tenure + age, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.26295 -0.60779 -0.00117  0.66936  2.17351 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.377140   0.504865  10.651  < 2e-16 ***
## exhaust      0.067240   0.077642   0.866 0.387184    
## sex         -1.293346   0.422770  -3.059 0.002424 ** 
## tenure      -0.096794   0.027183  -3.561 0.000431 ***
## age         -0.022167   0.008041  -2.757 0.006201 ** 
## exhaust:sex  0.401353   0.111166   3.610 0.000359 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9423 on 294 degrees of freedom
## Multiple R-squared:  0.2505, Adjusted R-squared:  0.2377 
## F-statistic: 19.65 on 5 and 294 DF,  p-value: < 2.2e-16
  1. First, is the interaction still significant? Did it decrease or increase with the covariates added?
  2. Next, what does the coefficient on exhaust mean here? What about sex?
  3. Finally, interpret both tenure and age.

Conditional Effects

  1. Given that the coefficient in the previous model on sex was interpretted when exhaust = 0 (something not actually seen in the data), let’s look at the difference between males and females when exhaust is at its mean value.
hosp %>%
  mutate(exhaust = exhaust - mean(exhaust, na.rm=TRUE)) %>%
  lm(safety ~ exhaust * sex + tenure + age, data = .) %>%
  summary()
## 
## Call:
## lm(formula = safety ~ exhaust * sex + tenure + age, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.26295 -0.60779 -0.00117  0.66936  2.17351 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.603221   0.391249  14.321  < 2e-16 ***
## exhaust      0.067240   0.077642   0.866 0.387184    
## sex          0.056135   0.171171   0.328 0.743183    
## tenure      -0.096794   0.027183  -3.561 0.000431 ***
## age         -0.022167   0.008041  -2.757 0.006201 ** 
## exhaust:sex  0.401353   0.111166   3.610 0.000359 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9423 on 294 degrees of freedom
## Multiple R-squared:  0.2505, Adjusted R-squared:  0.2377 
## F-statistic: 19.65 on 5 and 294 DF,  p-value: < 2.2e-16
  1. What does the coefficient on sex mean now? Compare it to the previously run model. What is different?

Probing Interactions

  1. Most aspects of probing an interaction are not necessary the need for assessing the interaction in more detail. We will use the Johnson-Neyman Technique for this as it applies in several situations. Here, we are using the book’s example of exhaust * tenure. A new function jtools::summ() gives some nice summary information of the results.
interaction_model <- lm(safety ~ exhaust * tenure + age + sex, data = hosp)
jtools::summ(interaction_model)
## MODEL INFO:
## Observations: 300
## Dependent Variable: safety
## Type: OLS linear regression 
## 
## MODEL FIT:
## F(5,294) = 20.01, p = 0.00
## R² = 0.25
## Adj. R² = 0.24 
## 
## Standard errors: OLS
## 
## |               |  Est.| S.E.| t val.|    p|
## |:--------------|-----:|----:|------:|----:|
## |(Intercept)    |  3.51| 0.52|   6.77| 0.00|
## |exhaust        |  0.62| 0.11|   5.70| 0.00|
## |tenure         |  0.10| 0.06|   1.73| 0.08|
## |age            | -0.02| 0.01|  -2.53| 0.01|
## |sex            |  0.02| 0.17|   0.09| 0.93|
## |exhaust:tenure | -0.06| 0.02|  -3.80| 0.00|
  1. We can then look at the difference in the slope of exhaust based on different levels of tenure using interactions::sim_slopes() and interactions::interact_plot(). Both will give us similar, useful information. (See https://interactions.jacob-long.com/index.html for lots of examples on how to use the interactions package.)
interactions::sim_slopes(interaction_model, pred = exhaust, modx = tenure)
## JOHNSON-NEYMAN INTERVAL 
## 
## When tenure is OUTSIDE the interval [7.97, 15.27], the slope of
## exhaust is p < .05.
## 
## Note: The range of observed values of tenure is [1.00, 18.00]
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of exhaust when tenure = 8.79 (+ 1 SD): 
## 
##  Est.   S.E.   t val.      p
## -----  -----  -------  -----
##  0.08   0.07     1.07   0.29
## 
## Slope of exhaust when tenure = 5.51 (Mean): 
## 
##  Est.   S.E.   t val.      p
## -----  -----  -------  -----
##  0.28   0.06     5.08   0.00
## 
## Slope of exhaust when tenure = 2.24 (- 1 SD): 
## 
##  Est.   S.E.   t val.      p
## -----  -----  -------  -----
##  0.48   0.08     6.05   0.00
interactions::interact_plot(interaction_model, 
                            pred = exhaust, 
                            modx = tenure, 
                            plot.points = TRUE, 
                            interval = TRUE)