Introduction to Interaction Effects

Testing for an Interaction Effect

The inferential question is whether the interaction effect that we are seeing in the empirical data is “real”, or whether it is an artifact of sampling error. To examine this we need a way to test whether the slopes of the two regression lines are equal. To do this, we create another predictor that is the product of the two predictors we believe interact and include that product term in the regression model along with the original predictors we used to create it (i.e., also include the constituent main-effects). In our example, we multiply the dummy-coded female predictor by the boundary-spanning work predictor to create the interaction term. Then we fit a model that includes the dummy-coded female predictor, the boundary-spanning work predictor, and the newly created interaction term.

# Create interaction term
work = work %>%
  mutate(
    bound_work_female = bound_span_work * female
    )

# View data
head(work)

# A tibble: 6 × 6
    guilt bound_span_work female authority married bound_work_female
    <dbl>           <dbl>  <dbl>     <dbl>   <dbl>             <dbl>
1 -1.64           -1.13        1     2.65        1           -1.13  
2 -0.0272          0.718       1     1.38        1            0.718 
3  0.276          -0.0164      1     1.14        1           -0.0164
4  1.65            1.29        1     1.13        1            1.29  
5 -0.149          -0.153       1     0.646       1           -0.153 
6 -1.47           -2.04        1     1.23        1           -2.04  

# Fit interaction model
lm.d = lm(guilt ~ 1 + bound_span_work + female + bound_work_female, data = work)
tidy(lm.d)

# A tibble: 4 × 5
  term              estimate std.error statistic  p.value
  <chr>                <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)         -0.325    0.0693     -4.68 4.29e- 6
2 bound_span_work      0.556    0.0718      7.74 1.60e-13
3 female               0.563    0.0890      6.33 9.09e-10
4 bound_work_female    0.161    0.0903      1.78 7.63e- 2

The fitted equation is:

\[ \hat{\mathrm{Guilt}}_i = -0.32 + 0.56(\mathrm{Boundary\mbox{-}spanning~work}_i) + 0.56(\mathrm{Female}_i) + 0.16(\mathrm{Boundary\mbox{-}spanning~work}_i)(\mathrm{Female}_i) \]

If we substitute in the appropriate dummy-coding we can obtain the fitted equations for females and non-females.

Female Fitted Equation

\[ \begin{split} \hat{\mathrm{Guilt}}_i &= -0.32 + 0.56(\mathrm{Boundary\mbox{-}spanning~work}_i) + 0.56(1) + 0.16(\mathrm{Boundary\mbox{-}spanning~work}_i)(1) \\[2ex] &= 0.24 + 0.72(\mathrm{Boundary\mbox{-}spanning~work}_i) \end{split} \]

Non-female Fitted Equation

\[ \begin{split} \hat{\mathrm{Guilt}}_i &= -0.32 + 0.56(\mathrm{Boundary\mbox{-}spanning~work}_i) + 0.56(0) + 0.16(\mathrm{Boundary\mbox{-}spanning~work}_i)(0) \\[2ex] &= -0.32 + 0.56(\mathrm{Boundary\mbox{-}spanning~work}_i) \end{split} \]

Notice that the inclusion of the interaction term allows the slopes in the two fitted equations to be different. Moreover, the female slope of \(B=0.72\) is 0.16 units higher than the non-female slope of \(B=0.56\). This difference in slopes is expressed by the estimated coefficient for the interaction term.

The hypothesis test associated with the interaction term in the tidy() output evaluates the null hypothesis that:

\[ H_0: \beta_{\mathrm{Interaction}} = 0 \]

This is a test of whether the slopes for females and non-females are the same; if the interaction term is 0, then the slopes are the same in the population. In our example, the inferential results suggest that there is some evidence that the empirical data may not be compatible with the hypothesis that the difference in slopes is 0; \(t(296)=1.78\), \(p=.076\). This suggests that the effect of boundary-spanning work on guilt for females may be different than that for non-females.

Mathematical Expression of the Interaction Model

In general, the interaction model (with two predictors) can be written as,

\[ \begin{split} Y_i = \beta_0 + \beta_1(X_{1i}) &+ \beta_2(X_{2i}) + \beta_3\left(X_{1i}\right)\left(X_{2i}\right) + \epsilon_i, \quad \mathrm{where}\\[2ex] & \epsilon_i \overset{i.i.d.}{\sim} \mathcal{N}\big(0,\sigma^2_{\epsilon}\big) \end{split} \]

First notice that if \(\beta_3\), the coefficient on the interaction term, is zero, this equation reduces to the equation for the main-effects model, namely

\[ \begin{split} Y_i = \beta_0 &+ \beta_1(X_{1i}) + \beta_2(X_{2i}) + \epsilon_i, \quad \mathrm{where}\\[2ex] & \epsilon_i \overset{i.i.d.}{\sim} \mathcal{N}\big(0,\sigma^2_{\epsilon}\big) \end{split} \]

Because of this, in practice, if the observed data are compatible with the null hypothesis that the coefficient for the interaction term is zero, it is common to drop the interaction term from the model, and re-fit the model using the main-effects model.

To further understand the terms in the interaction model, we will write out the interaction model using the context of our example. (For ease of writing, we will drop the part of the model that specifies the model’s assumptions.)

\[ \mathrm{Guilt}_i = \beta_0 + \beta_1(\mathrm{Boundary\mbox{-}spanning~work}_i) + \beta_2(\mathrm{Female}_i) + \beta_3(\mathrm{Boundary\mbox{-}spanning~work}_i)(\mathrm{Female}_i) + \epsilon_i \]

Since the predictor \(\mathrm{Female}_i\) is a dummy coded predictor (0 for non-females; 1 for females), we can use substitute those values into the model and write the regression model associated with non-females and that for females.

Regression Model for Non-females

\[ \begin{split} \mathrm{Guilt}_i &= \beta_0 + \beta_1(\mathrm{Boundary\mbox{-}spanning~work}_i) + \beta_2(0) + \beta_3(\mathrm{Boundary\mbox{-}spanning~work}_i)(0) + \epsilon_i \\[2ex] &= \beta_0 + \beta_1(\mathrm{Boundary\mbox{-}spanning~work}_i) + \epsilon_i \\ \end{split} \]

Regression Model for Females

\[ \begin{split} \mathrm{Guilt}_i &= \beta_0 + \beta_1(\mathrm{Boundary\mbox{-}spanning~work}_i) + \beta_2(1) + \beta_3(\mathrm{Boundary\mbox{-}spanning~work}_i)(1) + \epsilon_i \\[2ex] &= \beta_0 + \beta_1(\mathrm{Boundary\mbox{-}spanning~work}_i) + \beta_2 + \beta_3(\mathrm{Boundary\mbox{-}spanning~work}_i) + \epsilon_i \\[2ex] &= \bigg[\beta_0 + \beta_2\bigg] + \bigg[\beta_1 + \beta_3\bigg](\mathrm{Boundary\mbox{-}spanning~work}_i) + \epsilon_i \\ \end{split} \]

We can use these equations to understand what each term in the interaction model represents:

• The intercept term from the interaction model (\(\beta_0\)) turns out to be the intercept term for the reference group (non-females in our example).
• The slope term associated with boundary-spanning work from the interaction model (\(\beta_1\)) turns out to be the effect of boundary-spanning work for the reference group (non-females in our example).
• The slope term associated with the dummy-coded female predictor from the interaction model (\(\beta_2\)) turns out to be the difference in intercept between the group coded 1 and the reference group (difference in the female and non-female intercepts in our example).
• The slope term associated with the interaction predictor from the interaction model (\(\beta_3\)) turns out to be the difference in slopes between the group coded 1 and the reference group (difference in the female and non-female slopes in our example).

Interpreting the Fitted Model’s Coefficients

Here we will use the interaction model we fitted earlier to understand how to interpret the different coefficients in the model. Recall that the fitted equation was,

\[ \hat{\mathrm{Guilt}}_i = -0.32 + 0.56(\mathrm{Boundary\mbox{-}spanning~work}_i) + 0.56(\mathrm{Female}_i) + 0.16(\mathrm{Boundary\mbox{-}spanning~work}_i)(\mathrm{Female}_i) \]

• The intercept term (\(\hat\beta_0=-0.32\)) indicates that non-females who experience an average amount of boundary-spanning work (boundary-spanning work = 0) have a guilt level of -0.32 standard deviations, on average.
• The slope term associated with boundary-spanning work (\(\hat\beta_1=0.56\)) indicates that, for non-females, each one-standard deviation difference in boundary-spanning work is associated with a 0.56-standard deviation difference in home-life/work guilt, on average.
• The slope term associated with the dummy-coded female predictor (\(\hat\beta_2=0.56\)) indicates that females who experience an average amount of boundary-spanning work (boundary-spanning work = 0) have a guilt level that is 0.56-standard deviations higher, on average, than non-females who experience an average amount of boundary-spanning work.
• The slope term associated with the interaction predictor (\(\hat\beta_3=0.16\)) indicates that, for females, the effect of boundary-spanning work on home-life/work guilt is 0.16-units higher, on average, than the effect of boundary-spanning work on home-life/work guilt for non-females. In other words, for females, each one-standard deviation increase in boundary-spanning work is associated with a 0.72-standard deviation difference in home-life/work guilt, on average.

In general, the easiest way to determine how to interpret the coefficients is to actually compute the partial regression equations for both groups based on the fitted interaction equation like we did earlier. Then you can just interpret the intercept and slope terms for the groups independently.

Plotting the Fitted Equation for the Interaction Model

Almost always, plotting the fitted model is a good idea once you have fitted an interaction model. To do this we use geom_abline() to plot the non-female and female fitted lines.

ggplot(data = work, aes(x = bound_span_work, y = guilt)) +
  geom_point(alpha = 0) +
  theme_bw() +
  xlab("Amount of boundary-spanning work") +
  ylab("Predicted home-life/work guilt") +
  geom_abline(intercept = 0.24, slope = 0.72, linetype = "solid", color = "#424651") +
  geom_abline(intercept = -0.32, slope = 0.56, linetype = "dashed", color = "#f5853f")

Figure 10: Fitted regression lines showing model-predicted home-life/work guilt as a function of boundary-spanning work for females (charcoal, solid line) and non-females (mango, dashed line).

In this plot, we can see some of the same effects we interpreted from the coefficients:

• Females have a higher intercept than non-females,
• The effect of boundary-spanning work for females is higher in magnitude than that for non-females.

The plot also allows us to see other things as well. For example,

• The lines do not cross at a boundary-spanning work value represented in the data. This implies that females are predicted to have more home-life/work guilt than non-females, on average, regardless of how much boundary-spanning work they experience. But, this differential diminishes for lower levels of boundary-spanning work.

In the social sciences sometimes this observed pattern of one group ALWAYS having a higher predicted level of the outcome within the range of X-values is referred to as an ordinal interaction. If the lines had crossed, that would imply that for some values of X, one group had higher predicted levels of the outcome, while for other values of X the other group had higher predicted levels of the outcome. This would be referred to as a disordinal interaction. It is generally better to show the plot and describe the effects for both groups rather than just stating that the interaction was ordinal or disordinal in nature.

Plotting the Results of the Model

To further aid our understanding of the results from fitting the interaction model, we will create several plots to visualize the different effects of interest. In these plots, we will control for job authority by setting it to the mean value of 0.83. Since marital status is also a covariate, in practice you might also set it to a fixed value and partial it out of the plot. However, here we will show this effect so you can begin to visualize differences between main-effects and interaction effects. To do this we will need to obtain four different partial fitted equations and plot them.

Non-female, non-married

\[ \begin{split} \hat{\mathrm{Guilt}}_i &= 0.05 + 0.61(\mathrm{Boundary\mbox{-}spanning~work}_i) + 0.40(0) - 0.53(0.83) + 0.23(0) + 0.12(\mathrm{Boundary\mbox{-}spanning~work}_i)(0) \\[2ex] &= -0.39 + 0.61(\mathrm{Boundary\mbox{-}spanning~work}_i) \end{split} \] Non-female, married

\[ \begin{split} \hat{\mathrm{Guilt}}_i &= 0.05 + 0.61(\mathrm{Boundary\mbox{-}spanning~work}_i) + 0.40(0) - 0.53(0.83) + 0.23(1) + 0.12(\mathrm{Boundary\mbox{-}spanning~work}_i)(0) \\[2ex] &= -0.17 + 0.61(\mathrm{Boundary\mbox{-}spanning~work}_i) \end{split} \] Female, non-married

\[ \begin{split} \hat{\mathrm{Guilt}}_i &= 0.05 + 0.61(\mathrm{Boundary\mbox{-}spanning~work}_i) + 0.40(1) - 0.53(0.83) + 0.23(0) + 0.12(\mathrm{Boundary\mbox{-}spanning~work}_i)(1) \\[2ex] &= 0.01 + 0.73(\mathrm{Boundary\mbox{-}spanning~work}_i) \end{split} \]

Female, married

\[ \begin{split} \hat{\mathrm{Guilt}}_i &= 0.05 + 0.61(\mathrm{Boundary\mbox{-}spanning~work}_i) + 0.40(1) - 0.53(0.83) + 0.23(1) + 0.12(\mathrm{Boundary\mbox{-}spanning~work}_i)(1) \\[2ex] &= 0.24 + 0.73(\mathrm{Boundary\mbox{-}spanning~work}_i) \end{split} \]

# Female
ggplot(data = work, aes(x = bound_span_work, y = guilt)) +
  geom_point(alpha = 0) +
  theme_bw() +
  xlab("Amount of boundary-spanning work") +
  ylab("Predicted home-life/work guilt") +
  geom_abline(intercept = 0.24, slope = 0.73, linetype = "dashed", color = "#424651") +
  geom_abline(intercept = 0.01, slope = 0.73, linetype = "solid", color = "#424651") +
  ggtitle("Female")

# Non-female
ggplot(data = work, aes(x = bound_span_work, y = guilt)) +
  geom_point(alpha = 0) +
  theme_bw() +
  xlab("Amount of boundary-spanning work") +
  ylab("Predicted home-life/work guilt") +
  geom_abline(intercept = -0.39, slope = 0.61, linetype = "solid", color = "#f5853f") +
  geom_abline(intercept = -0.17, slope = 0.61, linetype = "dashed", color = "#f5853f") +
  ggtitle("Non-female")

Figure 7: Fitted regression lines showing model-predicted home-life/work guilt as a function of boundary-spanning work for married (dashed lines) and non-married (solid lines) employees conditioned on gender. Job authority is controlled for by setting this to its mean value.

Notice that the slope of the lines represents the effect of boundary-spanning work and that this differs by gender; the lines across gender are not parallel. Across marital status, however, the slopes are the same (at least within gender); parallel lines. This is because marital status was only included in the model as a main effect, whereas gender was included as part of an interaction effect.

We could also have facetted on marital status rather than gender. In this plot, the partial fitted lines within each panel of theplot are not parallel, but the set of lines are parallel across panels.

# Non-married
ggplot(data = work, aes(x = bound_span_work, y = guilt)) +
  geom_point(alpha = 0) +
  theme_bw() +
  xlab("Amount of boundary-spanning work") +
  ylab("Predicted home-life/work guilt") +
  geom_abline(intercept = 0.01, slope = 0.73, linetype = "solid", color = "#424651") +
  geom_abline(intercept = -0.39, slope = 0.61, linetype = "dashed", color = "#f5853f") +
  ggtitle("Non-married")

# Married
ggplot(data = work, aes(x = bound_span_work, y = guilt)) +
  geom_point(alpha = 0) +
  theme_bw() +
  xlab("Amount of boundary-spanning work") +
  ylab("Predicted home-life/work guilt") +
  geom_abline(intercept = 0.24, slope = 0.73, linetype = "solid", color = "#424651") +
  geom_abline(intercept = -0.17, slope = 0.61, linetype = "dashed", color = "#f5853f") +
  ggtitle("Married")

Figure 14: Fitted regression lines showing model-predicted home-life/work guilt as a function of boundary-spanning work for females (charcoal, solid line) and non-females (mango, dashed line) conditioned on marital status. Job authority is controlled for by setting this to its mean value.