In LukasWallrich/GoldCoreQuants: Does not matter.

Interaction terms in linear models

Watch this video for an introduction to interaction terms in linear models:

r video_code("coyXicRccSg")

The section below essentially contains the code needed for the examples in the video, with some further annotations. Before that, just one note on terminology:

• Generally, if two variables interact, then the effect that one has on the predictor (the regression slope) depends on the value of the other. Mathematically, that relationship is symetrical - if A interacts with B, then B interacts with A.
• For interpretation, we often designate one of the variables as the moderator. That just means that we are primarily interested in the effect of the other variable on the outcome and in how the moderator changes that effect. Often, demographic variables such as age or gender serve as moderators.

Example 1: Memory and chess

pacman::p_load("sandwich")
#TK - check why required by sim_slopes (at least on ubuntu) even though robust = FALSE

The first example draws on research by Gobet & Simon (1996) but uses simulated data.

r hide("Show the code to simulate the data")

library(tidyverse)

#Generate data (errors committed) - *roughly* based on Gobet & Simon 1996
set.seed(300688) #for reproducible results
ER <- rnorm(50, 4.9,3.5) + rnorm(50, 0, 2)
NR <- rnorm(50, 15.7, 4.0) + rnorm(50, 0, 2)
EF <- rnorm(50, 21.4, 5) + rnorm(50, 0, 2)
NF <- rnorm(50, 21.8, 5) + rnorm(50, 0, 2)

obs <- data.frame(player = "expert", type = "real", errors = ER) %>% 
  rbind(data.frame(player = "novice", type = "real", errors = NR)) %>%
  rbind(data.frame(player = "expert", type = "fake", errors = EF)) %>%
  rbind(data.frame(player = "novice", type = "fake", errors = NF)) %>% 
  mutate(type=factor(type), 
         player = factor(player, levels = c("novice", "expert")))

#Adding the centrally mirrored condition
EM <- rnorm(50, 7.8, 3.5) + rnorm(50, 0, 2)
NM <- rnorm(50, 18, 3.5) + rnorm(50, 0, 2)

obs2 <- data.frame(player = "expert", type = "real", errors = ER) %>% 
  rbind(data.frame(player = "novice", type = "real", errors = NR)) %>%
  rbind(data.frame(player = "expert", type = "fake", errors = EF)) %>%
  rbind(data.frame(player = "novice", type = "fake", errors = NF)) %>% 
  rbind(data.frame(player = "expert", type = "mirrored", errors = EM)) %>% 
  rbind(data.frame(player = "novice", type = "mirrored", errors = NM)) %>% 
  mutate(type=factor(type), 
         player = factor(player, levels = c("novice", "expert")))

r unhide()

The first model tests whether player level and position type interact. That is the case, based on the very small p-value of the interaction term. Then we use a plot to understand the nature of that interaction further - because the predictor is categorical, we use the cat_plot() function.

pacman::p_load(tidyverse)
mod <- lm(errors ~ player + type + player:type, obs)
summary(mod)

pacman::p_load(interactions)
cat_plot(mod, pred="player", modx = "type", geom="line")

Next we consider a third condition - chess positions that are neither quite real nor entirely fake, but positions that are mirrored. With that, we get multiple dummy interaction terms. To test whether they are collectively significant, we need to use the Anova() function from the car package.

mod <- lm(errors ~ player + type + player:type, obs2)
summary(mod)
pacman::p_load(car)
car::Anova(mod, type=3)

An interaction plot can then help again to understand what is going on. Here it shows that two of the three conditions are very similar.

cat_plot(mod, pred="player", modx = "type", geom="line")

Example 2: link between obesity and negative emotions in the European Social Survey

In this example, I considered the link between obesity and negative emotions in Germany in the European Social Survey 2014. (Note that this relationship does not appear in the UK, which indicates that it should probably be treated as an interesting observation rather than a likely general relationship for now.)

r hide("Click here to see the code that loads and prepares the data")

ess <- read_rds(url("http://empower-training.de/Gold/round7.RDS"))

pacman::p_load(psych) #Used to create a scale

a <- ess %>% select(cldgng, fltsd, enjlf, fltlnl, wrhpp, slprl, flteeff, fltdpr) %>% 
  haven::zap_labels() %>% 
  #This is sometimes needed when SPSS files have been imported with haven 
  #and errors related to data classes appear.
  psych::alpha(check.keys = TRUE)

ess$depr <- a$scores

ess$bmi <- ess$weight/((ess$height/100)^2)

#Filter for Germans with realistic BMI who are not underweight 
#and reported their gender
essDE <- ess %>% filter(bmi < 60, bmi>=19, gndr != "No answer", cntry=="DE")

#Data prep for example 3 - strip unnecessary labels
ess$stflife <- haven::zap_labels(ess$stflife)

r unhide()

The lm() output below shows that there is a significant interaction between gender and BMI, with women having a stronger relationship between BMI and the frequency of experiencing negative emotions. This is again shown in an interaction plot - as the predictor variable is continuous, we now use the interact_plot() function.

pacman::p_load(interactions)
mod <- lm(depr ~ bmi + gndr + bmi:gndr, essDE)
summary(mod)

(ref:int plot) interact_plot() shows different slopes for men and women.

interact_plot(mod, pred="bmi", modx = "gndr")

Note that ggplot2 automatically includes an interaction when fitting regression lines as soon as a categorical variable is mapped to the colour (or linetype) aesthetic - so in simple cases, geom_smooth(method = "lm") can be used as an alternative to interact_plot()

ggplot(essDE, aes(x=bmi, y=depr, colour=gndr)) + geom_smooth(method="lm", se=FALSE)

Example 3: link between working hours, income and life satisfaction

This example considers whether working hours (wkhtot) affect the link between income (hinctnta) and life satisfaction - i.e. does working very long hours make income less valuable? This time, the effect appeared in the UK data in the 2014 European Social Survey - again, the question might be whether that is just an incident of spurious data mining, or whether it reveals a broader relationship.

In any case, let's have a look at the interaction. Note that * in the lm() formula is an abbreviation for + and :, so that the command below could also be written as lm(stflife ~ wkhtot + hinctnta + wkhtot:hinctnta)

essUK <- filter(ess, cntry == "GB")
mod <- lm(stflife ~ wkhtot*hinctnta, essUK) 
summary(mod)

We can see that while income positively predicts life satisfaction, this effects appears to be weaker when both income and working hours are high. However, we need to be careful with interpretation here - technically, the coefficients for income and working hours now reflect the impact of that variable when the other variable is 0, and are thus unlikely to be meaningful. Therefore, it is better to look at the interaction plot and at simple slopes analyses. Both can now be done in two ways, depending on which variable we see as the primary predictor/variable of interest.

interact_plot(mod, pred="hinctnta", modx = "wkhtot")

interact_plot(mod, pred="wkhtot", modx = "hinctnta")

Simple slopes analyses (sim_slopes()) offer similar information together with significance tests and thus help to decide which of the slopes should be interpreted. They contain the Johnson-Neyman interval, which is the range of values of the moderator for which the predictor has a significant effect on the outcome.

pacman::p_load(interactions)
sim_slopes(mod, pred="hinctnta", modx = "wkhtot")
sim_slopes(mod, pred="wkhtot", modx = "hinctnta")

Finally, the johnson_neyman() function creates a plot showing the slope of one variable depending on the value of the other, while highlighting the regions of significance. This can help with an understanding of the relationship, but is not (yet?) widely used.

johnson_neyman(mod, pred="wkhtot", modx = "hinctnta")

LukasWallrich/GoldCoreQuants documentation built on Nov. 27, 2021, 1:58 a.m.