# 1 Interactions

• Interactions have the same meaning they did in ANOVA
• there is a synergistic effect between two variables
• However now we can examine interactions between continuous variables
• Additive Effects: $$Y = B_1X + B_2Z + B_0$$
• Interactive Effects: $$Y = B_1X + B_2Z + B_3XZ + B_0$$

## 1.1 Example of Interaction

• We need to make multivariate set of variables and then force an interaction
• DV: Reading Speed, IVs: Age + IQ
library(car)
set.seed(42)
# 250 people
n <- 250
# Uniform distrobution of Ages (5 to 15 year olds)
X <- runif(n, 5, 15)

# normal distrobution of IQ (100 +-15)
Z <- rnorm(n, 100, 15)

# Our equation to  create Y
Y <- .3*X + .5*Z + .3*X*Z +300 + rnorm(n, sd=10)
#Built our data frame

# Plot: Notes: reg.line=FALSE means remove best fit line.
# smoother=loessLine means add smoothed line called the loess
scatterplot(ReadingSpeed~Age, data= Reading.Data, reg.line=FALSE, smoother=loessLine)

scatterplot(ReadingSpeed~IQ, data= Reading.Data, reg.line=FALSE, smoother=loessLine)

scatterplot(Age~IQ, data= Reading.Data, reg.line=FALSE, smoother=loessLine)

## 1.2 Lets test an interaction

2. Model 1: Examine a main-effects model (X + Z) [not necessary but useful]
3. Model 2: Examine a main-effects + Interaction model (X + Z + X:Z)
• Note: You can simply code it as X*Z in r, as it will automatically do (X + Z + X:Z)
#Center
Reading.Data$Age.C<-scale(Reading.Data$Age, center = TRUE, scale = FALSE)[,]
Reading.Data$IQ.C<-scale(Reading.Data$IQ, center = TRUE, scale = FALSE)[,]

library(stargazer)
column.labels = c("Main Effects", "Interaction"),
intercept.bottom = FALSE,
single.row=FALSE,
notes.append = FALSE,
header=FALSE)
 Dependent variable: ReadingSpeed Main Effects Interaction (1) (2) Constant 653.207*** 653.632*** (0.980) (0.653) Age.C 29.778*** 29.851*** (0.337) (0.224) IQ.C 3.579*** 3.556*** (0.071) (0.047) Age.C:IQ.C 0.305*** (0.017) Observations 250 250 R2 0.976 0.989 Adjusted R2 0.976 0.989 Residual Std. Error 15.502 (df = 247) 10.314 (df = 246) F Statistic 5,028.365*** (df = 2; 247) 7,676.634*** (df = 3; 246) Note: p<0.1; p<0.05; p<0.01
• Also change in $$R^2$$ is significant, as we might expect
ChangeInR<-anova(Centered.Read.1,Centered.Read.2)
knitr::kable(ChangeInR, digits=4)
Res.Df RSS Df Sum of Sq F Pr(>F)
247 59357.43 NA NA NA NA
246 26169.85 1 33187.58 311.9676 0

## 1.3 How visualize this?

• Well, there are some options
• We can do a fancy-pants surface plot, but that is hard to put into a paper
• More common this is to examine slope of one factor at different levels of the other (Simple Slopes)
• What we need to decide is at which levels

### 1.3.1 Hand picking

• Manually select the what levels of each you are going to examine
• For for tracking values of interest
• here I picked -3 to 3 for age (mix and max) and -15, 0 15 for IQ (theoretical 1 SD and mean)
library(effects)
xlevels=list(Age.C=seq(-3,3, 1),
IQ.C=c(-15,0,15)))
knitr::kable(summary(Inter.1a)$effect, digits=4) -15 0 15 -3 524.445 564.0773 603.7095 -2 549.727 593.9288 638.1305 -1 575.009 623.7803 672.5515 0 600.291 653.6317 706.9725 1 625.573 683.4832 741.3935 2 650.855 713.3347 775.8144 3 676.137 743.1862 810.2354 plot(Inter.1a, multiline = TRUE) ### 1.3.2 Quantile • Examine the levels based on quantiles (bins based on probabiltiy) • We will do this into 5 equal bins based on probability cut-offs • Does not assume normality for IV IQ.Quantile<-quantile(Reading.Data$IQ.C,probs=c(0,.25,.50,.75,1))
IQ.Quantile<-round(IQ.Quantile,2)
IQ.Quantile
##     0%    25%    50%    75%   100%
## -39.89  -9.91  -0.19   9.45  37.50
Inter.1b<-effect(c("Age.C*IQ.C"), Centered.Read.2,
xlevels=list(Age.C=seq(-3,3, 1),
IQ.C=IQ.Quantile))
plot(Inter.1b, multiline = TRUE)

### 1.3.3 Based on the SD

• We select 3 values: $$M-1SD$$, $$M$$, $$M+1SD$$
• Not it does not have to be 1 SD, it can be 1.5,2 or 3
• Assumes normality for IV
IQ.SD<-c(mean(Reading.Data$IQ.C)-sd(Reading.Data$IQ.C),
mean(Reading.Data$IQ.C), mean(Reading.Data$IQ.C)+sd(Reading.Data$IQ.C)) IQ.SD<-round(IQ.SD,2) IQ.SD ## [1] -13.89 0.00 13.89 Inter.1c<-effect(c("Age.C*IQ.C"), Centered.Read.2, xlevels=list(Age.C=seq(-3,3, 1), IQ.C=IQ.SD)) plot(Inter.1c, multiline = TRUE) ## 1.4 Why center your IVs? • The center of each IV represents the mean of that IV • Thus when you interact them $$X*Z$$, that means $$0*0 = 0$$ • Also this can reduce multicollineatity issues • This is because if $$X*Z$$ creates a line, it means you have added a new predictor (XZ) that strongly correlates with X and Z X<-c(0,2,4,6,8,10) Z<-c(0,2,4,6,8,10) XZ<-X*Z plot(XZ) c(cor(X,Z), cor(X,XZ), cor(Z,XZ)) ## [1] 1.0000000 0.9598833 0.9598833 -But if you center them, now you will make a U-shaped interaction which will be orthogonal to X and Z! X.C<-scale(X, center = TRUE, scale = FALSE)[,] Z.C<-scale(Z, center = TRUE, scale = FALSE)[,] XZ.C<-X.C*Z.C plot(XZ.C) c(cor(X.C,Z.C), cor(X.C,XZ.C), cor(Z.C,XZ.C)) ## [1] 1 0 0 • Lets apply this to our class example • For example, see below where I manually multiply the values and you can see a very strong positive slope Reading.Data$Age.X.IQ<-Reading.Data$Age*Reading.Data$IQ
scatterplot(ReadingSpeed~Age.X.IQ, data= Reading.Data, reg.line=FALSE, smoother=loessLine)

• But when you center them, things change
Reading.Data$Age.X.IQ.C<-Reading.Data$Age.C*Reading.Data$IQ.C scatterplot(ReadingSpeed~Age.X.IQ.C, data= Reading.Data, reg.line=FALSE, smoother=loessLine) ### 1.4.1 Lets run an uncentered regression • as you can see the the terms have changed from centered models Uncentered.Read.1<-lm(ReadingSpeed~IQ+Age,Reading.Data) Uncentered.Read.2<-lm(ReadingSpeed~IQ*Age,Reading.Data) stargazer(Uncentered.Read.1,Uncentered.Read.2,type="html", column.labels = c("Main Effects", "Interaction"), intercept.bottom = FALSE, single.row=FALSE, notes.append = FALSE, header=FALSE)  Dependent variable: ReadingSpeed Main Effects Interaction (1) (2) Constant -3.036 304.494*** (7.979) (18.203) IQ 3.579*** 0.482*** (0.071) (0.182) Age 29.778*** -0.427 (0.337) (1.725) IQ:Age 0.305*** (0.017) Observations 250 250 R2 0.976 0.989 Adjusted R2 0.976 0.989 Residual Std. Error 15.502 (df = 247) 10.314 (df = 246) F Statistic 5,028.365*** (df = 2; 247) 7,676.634*** (df = 3; 246) Note: p<0.1; p<0.05; p<0.01 • lets check the multicollinearity of the uncentered model vif(Uncentered.Read.2) ## IQ Age IQ:Age ## 14.89057 59.14549 71.06919 • We lost our main effect of Age as the variances got all inflated • We did not have this problem in the centered model vif(Centered.Read.2) ## Age.C IQ.C Age.C:IQ.C ## 1.001540 1.001961 1.001144 • Note: Its more standard to plot the uncentered version, but run the analysis on the centered data # 2 Testing Simple Slopes • The goal here it figure out when the slope at a given level of another variable is different from zero • We chop up the interaction at specific places as we did with the interactions plots (-1 SD, M, +1 SD) on the moderating variable (a third variable that affects the strength of the relationship between a dependent and independent variable) • Next we test the slopes on the predictor variable (IV of main interest) • Example: Performance on a task, predicted by Level of training, moderated by how hard the challange is set.seed(42) # 250 people n <- 250 #Training (low to high) X <- runif(n, -10, 10) # normal distrobution of Challenge Difficulty Z <- rnorm(n, 0, 15) # Our equation to create Y Y <- 2.2*X + 2.46*Z + .75*X*Z + 16 + rnorm(n, sd=50) #Built our data frame Skill.Data<-data.frame(Performance=Y,Training=X,Challenge=Z) #run our regression Skill.Model.1<-lm(Performance~Training+Challenge,Skill.Data) Skill.Model.2<-lm(Performance~Training*Challenge,Skill.Data) stargazer(Skill.Model.1,Skill.Model.2,type="html", column.labels = c("Main Effects", "Interaction"), intercept.bottom = FALSE, single.row=FALSE, notes.append = FALSE, header=FALSE)  Dependent variable: Performance Main Effects Interaction (1) (2) Constant 14.238*** 16.173*** (4.909) (3.268) Training 0.895 1.542*** (0.843) (0.562) Challenge 2.855*** 2.600*** (0.354) (0.236) Training:Challenge 0.762*** (0.043) Observations 250 250 R2 0.210 0.652 Adjusted R2 0.204 0.648 Residual Std. Error 77.510 (df = 247) 51.571 (df = 246) F Statistic 32.867*** (df = 2; 247) 153.486*** (df = 3; 246) Note: p<0.1; p<0.05; p<0.01 ## 2.1 Ploting • We will use the rockchalk library to visualize our interaction and test our simple slopes library(rockchalk) ## ## Attaching package: 'rockchalk' ## The following object is masked from 'package:plyr': ## ## summarize ## The following object is masked from 'package:MASS': ## ## mvrnorm m1ps <- plotSlopes(Skill.Model.2, modx = "Challenge", plotx = "Training", n=3, modxVals="std.dev") m1psts <- testSlopes(m1ps) ## Values of Challenge OUTSIDE this interval: ## lo hi ## -3.5050807 -0.5730182 ## cause the slope of (b1 + b2*Challenge)Training to be statistically significant round(m1psts$hypotests,4)
##        "Challenge"   slope Std. Error  t value Pr(>|t|)
## (m-sd)      -14.50 -9.5013     0.8131 -11.6848   0.0000
## (m)          -0.61  1.0771     0.5611   1.9196   0.0561
## (m+sd)       13.28 11.6554     0.8282  14.0737   0.0000
• We see the slope at the mean is not significant,
• The +1SD and -1SD are significant, this is what is driving the interaction
• Well trained people at low levels of challenge get worse with more training, while high levels of challenge with more training get better
• The bonus using rockchalk is that you can change the number of lines (change n=3) you want to see or you can test quantiles as well (change modxVals=“quantile”)
m1pQ <- plotSlopes(Skill.Model.2, modx = "Challenge", plotx = "Training", n=5, modxVals="quantile")

m1pQts <- testSlopes(m1pQ)
## Values of Challenge OUTSIDE this interval:
##         lo         hi
## -3.5050807 -0.5730182
## cause the slope of (b1 + b2*Challenge)Training to be statistically significant
round(m1pQts\$hypotests,4)
##      "Challenge"    slope Std. Error  t value Pr(>|t|)
## 0%      -40.4989 -29.3016     1.7993 -16.2847   0.0000
## 25%     -10.5190  -6.4694     0.6990  -9.2555   0.0000
## 50%      -0.7985   0.9335     0.5610   1.6640   0.0974
## 75%       8.8382   8.2726     0.6994  11.8278   0.0000
## 100%     36.8939  29.6394     1.7214  17.2183   0.0000

### 2.1.1 Notes

• You often don’t need to run the t-tests on the simple slopes, but they can useful to test very specific hypotheses
• You can run simple slopes analysis on three-way interactions, but lets leave that aside for now as you would have to use a different R-package, again often the plot tells the story

# 3 Non-Linear interactions

• You can have non-linear terms interacting with other linear and non-linear terms
• Example: Quit smoking, X = fear of you health, Z = moderated by Self-Efficacy for quiting
set.seed(42)
# 250 people
n <- 250
#Fear of Health
X <- runif(n, 1, 9)
#Centered Self-Efficacy
Z <- runif(n, 0, 15)
# Our equation to  create Y
Y <- .178*X - .052*X^2 + .164*X*Z + .164*(X^2)*Z+ 3.651 + rnorm(n, sd=50)
#Built our data frame
Smoke.Data<-data.frame(Smoking=Y,Health=X,SelfEfficacy=Z)
#run our regression
Smoke.Model.1<-lm(Smoking~poly(Health,2)+SelfEfficacy,Smoke.Data)
Smoke.Model.2<-lm(Smoking~poly(Health,2)*SelfEfficacy,Smoke.Data)
stargazer(Smoke.Model.1,Smoke.Model.2,type="html",
column.labels = c("Main Effects", "Interaction"),
intercept.bottom = FALSE,
single.row=FALSE,
notes.append = FALSE,
header=FALSE)
 Dependent variable: Smoking Main Effects Interaction (1) (2) Constant 1.121 0.888 (6.187) (5.857) poly(Health, 2)1 464.121*** 49.056 (51.816) (90.966) poly(Health, 2)2 118.649** 50.228 (51.721) (93.073) SelfEfficacy 5.713*** 5.931*** (0.745) (0.706) poly(Health, 2)1:SelfEfficacy 60.253*** (11.088) poly(Health, 2)2:SelfEfficacy 10.879 (11.441) Observations 250 250 R2 0.353 0.425 Adjusted R2 0.345 0.413 Residual Std. Error 51.685 (df = 246) 48.903 (df = 244) F Statistic 44.648*** (df = 3; 246) 36.081*** (df = 5; 244) Note: p<0.1; p<0.05; p<0.01

## 3.1 Plotting the Simple Slopes when there are Curves

• These can be done with the ‘effects’ package as I showed you above and last week
• or we can simply use the plotCurves function from the ‘rockchalk’ package
SmokeInterPlot <- plotCurves(Smoke.Model.2, modx = "SelfEfficacy", plotx = "Health",n=3,modxVals="std.dev")

• Notice that the second order term was not significant in the main effect model
• This means you could have missed it. As the model 1 plot would have looked like this:
plotCurves(Smoke.Model.1, plotx = "Health")

• In other words, you cannot see the higher order effects as they are hidden by the moderating variable
• How to find these hidden things?
• You have to test for interactions and changes in $$R^2$$
• You have try to add higher order terms, but they should be motivated by some theory
• What if you did not know the second order term was in the model above.
Smoke.Model.1.L<-lm(Smoking~Health+SelfEfficacy,Smoke.Data)
Smoke.Model.2.L<-lm(Smoking~Health*SelfEfficacy,Smoke.Data)
stargazer(Smoke.Model.1.L,Smoke.Model.2.L,type="html",
column.labels = c("Main Effects", "Interaction"),
intercept.bottom = FALSE,
single.row=FALSE,
notes.append = FALSE,
header=FALSE)
 Dependent variable: Smoking Main Effects Interaction (1) (2) Constant -62.385*** -5.746 (9.811) (14.105) Health 12.606*** 1.370 (1.420) (2.497) SelfEfficacy 5.649*** -2.384 (0.751) (1.664) Health:SelfEfficacy 1.627*** (0.304) Observations 250 250 R2 0.339 0.407 Adjusted R2 0.333 0.400 Residual Std. Error 52.129 (df = 247) 49.444 (df = 246) F Statistic 63.249*** (df = 2; 247) 56.387*** (df = 3; 246) Note: p<0.1; p<0.05; p<0.01
ChangeInR.Smoke<-anova(Smoke.Model.1.L,Smoke.Model.2.L)
knitr::kable(ChangeInR.Smoke, digits=4)
Res.Df RSS Df Sum of Sq F Pr(>F)
247 671215.5 NA NA NA NA
246 601410.1 1 69805.32 28.5531 0
• Reploting as linear by linear interaction
library(rockchalk)
SmokeInterPlot.Linear <- plotSlopes(Smoke.Model.2.L, modx = "SelfEfficacy", plotx = "Health", n=3, modxVals="std.dev")

• lets check the change in $$R^2$$ from Linear interaction model to poly interaction model now
stargazer(Smoke.Model.2.L,Smoke.Model.2,type="html",
column.labels = c("Linear", "Poly"),
intercept.bottom = FALSE,
single.row=FALSE,
notes.append = FALSE,
header=FALSE)
 Dependent variable: Smoking Linear Poly (1) (2) Constant -5.746 0.888 (14.105) (5.857) Health 1.370 (2.497) poly(Health, 2)1 49.056 (90.966) poly(Health, 2)2 50.228 (93.073) SelfEfficacy -2.384 5.931*** (1.664) (0.706) Health:SelfEfficacy 1.627*** (0.304) poly(Health, 2)1:SelfEfficacy 60.253*** (11.088) poly(Health, 2)2:SelfEfficacy 10.879 (11.441) Observations 250 250 R2 0.407 0.425 Adjusted R2 0.400 0.413 Residual Std. Error 49.444 (df = 246) 48.903 (df = 244) F Statistic 56.387*** (df = 3; 246) 36.081*** (df = 5; 244) Note: p<0.1; p<0.05; p<0.01
ChangeInR.Smoke.2<-anova(Smoke.Model.2.L,Smoke.Model.2)
knitr::kable(ChangeInR.Smoke.2, digits=4)
Res.Df RSS Df Sum of Sq F Pr(>F)
246 601410.1 NA NA NA NA
244 583529.7 2 17880.49 3.7383 0.0252

## 3.2 Notes

• So model with poly is a significantly better fit, but really its not a huge change
• Would it change the story by much?
• Most people would ignore the higher order term, but you the choice depends on your theory!
• Ignoring the higher order term can make it easier to test simple slopes
• But ignoring the higher order terms can voilate your assumptions, lets compare our linear vs poly models residuals
plot(Smoke.Model.2.L, which=1)

plot(Smoke.Model.2, which=1)