EDUC/PSY 7610

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.

- 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)
```

- Import the
`hosp`

data set into R.

`data(hosp)`

- 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))
```

- 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")
```

- 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
```

- 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? - Does this interpretation match the figure we created before? Why or why not?
- 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
```

- First, is the interaction still significant? Did it decrease or increase with the covariates added?
- Next, what does the coefficient on
`exhaust`

mean here? What about`sex`

? - Finally, interpret both
`tenure`

and`age`

.

- 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
```

- What does the coefficient on
`sex`

mean now? Compare it to the previously run model. What is different?

- 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|
```

- 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)
```