```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## 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 `rlm` package.
```{r}
library(tidyverse)
library(furniture)
library(rlm)
```
2. Import the `hosp` data set into R.
```{r}
data(hosp)
```
3. Because when a data set was imported from SPSS (as this one originally was), there are `labelled` variable types, we need to adjust these using a little loop based on whether that variable is labelled.
```{r}
hosp[] <- hosp %>%
map_if(labelled::is.labelled, ~as_factor(.x))
```
4. 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).
```{r}
hosp %>%
ggplot(aes(exhaust, safety, group = sex, color = sex)) +
geom_point(alpha = .6) +
geom_smooth(method = "lm")
```
5. 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.
```{r}
hosp %>%
lm(safety ~ exhaust * sex, data = .) %>%
summary()
```
### Interpretation
6. 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?
7. Does this interpretation match the figure we created before? Why or why not?
8. Let's include some covariates in the model; namely, `tenure` (the job experience the individual has) and `age` (the individual's age in years).
```{r}
hosp %>%
lm(safety ~ exhaust * sex + tenure + age, data = .) %>%
summary()
```
9. First, is the interaction still significant? Did it decrease or increase with the covariates added?
10. Next, what does the coefficient on `exhaust` mean here? What about `sex`?
11. Finally, interpret both `tenure` and `age`.
### Conditional Effects
12. 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.
```{r}
hosp %>%
mutate(exhaust = exhaust - mean(exhaust, na.rm=TRUE)) %>%
lm(safety ~ exhaust * sex + tenure + age, data = .) %>%
summary()
```
13. What does the coefficient on `sex` mean now? Compare it to the previously run model. What is different?
### Probing Interactions
14. Most aspects of probing an interaction are not necessary but you may come across reviewers that are... let's say... less than helpful that may ask for some probing. We will use the Johnson-Neyman Technique for this as it should apply in several situations. The `jn()` function that `rlm` re-exports for us can do this. For this method to be useful, we need an interaction between two continuous variables. Here, we are using the book's example of `exhaust * tenure`. Note that we use this same interaction model for the average marginal effects below.
```{r}
interaction_model <- lm(safety ~ exhaust * tenure + age + sex, data = hosp)
jn(interaction_model, dv = "safety", iv = "exhaust", mod = "tenure")
```
15. For what values of `tenure` is `exhaust` significantly different than zero? Where does the effect of `exhaust` turn negative?
16. Do these results agree with the following figure? Why or Why not?
```{r}
hosp %>%
mutate(tenure_groups = cut_number(tenure, 6)) %>%
ggplot(aes(exhaust, safety, group = tenure_groups, color = tenure_groups)) +
geom_point(alpha = .6) +
geom_smooth(method = "lm",
se = FALSE,
fullrange = TRUE) +
labs(color = "Tenure") +
scale_color_hue() +
theme_bw()
```
17. What are some of the shortcomings of this plot? Is there any way to work around it?
### Average Marginal Effects
18. Obtaining the average marginal effect of the exhaust variable is straightforward using the `margins` package.
```{r}
#install.packages("margins")
library(margins)
```
19. Let's use the `margins()` function with the last model we ran to see what the average effect of `exhaust` is. Using `margins(lm_obj) %>% summary()` we obtain the effects, their SEs, p-values, and 95% confidence intervals.
```{r}
margins(interaction_model) %>%
summary()
```
20. Is the overall, average effect of exhaust significant? Interpret the AME on exhaust.
## Conclusion
Interactions are an essential element of understanding real-world phenomena. There are many ways to organize, understand, and communicate their meaning. This walk-through helped demonstrate some of the basic ways to understand these interesting relationships.