All Interactions Are Wrong (?)

part 1: getting the idea

toffa et al. x psicostat

How often do we deal with interactions?

No idea, we will do a review on top journals
Probably quite often in psychology, as we frequently build on already known main effects, looking for moderators
Especially when working on individual differences here’s a typical case: Group (reference vs special) x Condition (generally one more facilitating/challenging than the baseline); both factors have a main effect on the response

Example: Poisson

This is how error counts frequently look like:

Show code

set.seed(2)
library(ggplot2)
N = 500; age = rpois(N,1.4)
ts = 26; tc = "#555555"
ggplot()+
  geom_histogram(aes(y=after_stat(density),x=age),fill="darkblue",alpha=.7,binwidth=0.5)+
  scale_x_continuous(breaks=seq(0,100,1))+
  theme(text=element_text(size=ts,color="black"))+xlab("Errors")

Example: Poisson

And this is how a linear predict may affect error count:

Show code

set.seed(1)
N = 250; age = round(runif(N,6,10),1); y = rpois(N,exp(6.5-age*.8))
d = data.frame(y,age)
ggplot(d,aes(x=age,y=y))+
  coord_cartesian(ylim=c(-0.4,max(y))) +
  geom_point(size=4,alpha=.5,color="darkblue")+
  theme(text=element_text(size=ts,color="black"))+
  scale_y_continuous(breaks=seq(0,20,2))+scale_x_continuous(breaks=seq(0,20,.5))+
  ylab("Errors")+xlab("Age (years)")

Example: Poisson

A classical linear model would be problematic

Show code

fitL = glm(y~age, data=d)

library(effects)
effL = data.frame(allEffects(fitL,xlevels=list(age=seq(min(age),max(age),.05)))$"age")
ggplot(d,aes(x=age,y=y))+
  coord_cartesian(ylim=c(-0.4,max(y))) +
  geom_point(size=4,alpha=.5,color="darkblue")+
  geom_line(data=effL,aes(y=fit),size=2,color="darkred")+
  geom_ribbon(data=effL,aes(y=fit,ymin=lower,ymax=upper),alpha=.3,fill="darkred")+
  theme(text=element_text(size=ts,color="black"))+
  scale_y_continuous(breaks=seq(0,20,2))+scale_x_continuous(breaks=seq(0,20,.5))+
  ylab("Errors")+xlab("Age (years)")

Example: Poisson

Proper model is much better

Show code

fitP = glm(y~age, family=poisson(link="log"), data=d)

effP = data.frame(allEffects(fitP,xlevels=list(age=seq(min(age),max(age),.05)))$"age")
ggplot(d,aes(x=age,y=y))+
  coord_cartesian(ylim=c(-0.4,max(y))) +
  geom_point(size=4,alpha=.5,color="darkblue")+
  geom_line(data=effP,aes(y=fit),size=2,color="blue")+
  geom_ribbon(data=effP,aes(y=fit,ymin=lower,ymax=upper),alpha=.3,fill="blue")+
  #geom_line(data=effL,aes(y=fit),size=2,color="darkred")+
  #geom_ribbon(data=effL,aes(y=fit,ymin=lower,ymax=upper),alpha=.3,fill="darkred")+
  theme(text=element_text(size=ts,color="black"))+
  scale_y_continuous(breaks=seq(0,20,2))+scale_x_continuous(breaks=seq(0,20,.5))+
  ylab("Errors")+xlab("Age (years)")

Example: Poisson

But now… let’s add another main effect

Show code

set.seed(3)
N = 250; age = round(runif(N,6,10),1); group = rbinom(N,1,.5); y = rpois(N,exp(5.5-age*.5+group*.8))
d = data.frame(id=1:length(y),y,age,group=as.factor(group))
fitP = glm(y~age*group, family=poisson(link="log"), data=d)
effP = data.frame(allEffects(fitP,xlevels=list(age=seq(min(age),max(age),.05)))$"age:group")
effP$group = as.factor(effP$group)
ggplot(d,aes(x=age,y=y,shape=group,color=group))+
  geom_point(size=4,alpha=.6)+
  scale_color_manual(values=c("darkorange2","darkgreen"))+
  scale_fill_manual(values=c("darkorange1","darkgreen"))+
  geom_line(data=effP,aes(x=age,y=fit,group=group,linetype=group),size=2)+
  geom_ribbon(data=effP,aes(x=age,y=fit,ymin=lower,ymax=upper,group=group,fill=group),alpha=.3,color=NA)+
  theme(text=element_text(size=ts,color="black"))+scale_x_continuous(breaks=seq(0,20,.5))+
  ylab("Errors")+xlab("Age (years)")

Example: Poisson

library(glmmTMB)

fitP_0 = glmmTMB(y ~ age + group + (1|id), family=poisson(link="log"), data=d)
fitP_1 = glmmTMB(y ~ age * group + (1|id), family=poisson(link="log"), data=d)
anova(fitP_0, fitP_1)

Data: d
Models:
fitP_0: y ~ age + group + (1 | id), zi=~0, disp=~1
fitP_1: y ~ age * group + (1 | id), zi=~0, disp=~1
       Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
fitP_0  4 1210.6 1224.7 -601.31   1202.6                        
fitP_1  5 1212.5 1230.1 -601.25   1202.5 0.116      1     0.7334

Example: Poisson

But an interaction would emerge with a classical linear model

fitL_0 = glmmTMB(y ~ age + group + (1|id), data=d)
fitL_1 = glmmTMB(y ~ age * group + (1|id), data=d)
anova(fitL_0, fitL_1)

Data: d
Models:
fitL_0: y ~ age + group + (1 | id), zi=~0, disp=~1
fitL_1: y ~ age * group + (1 | id), zi=~0, disp=~1
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
fitL_0  5 1367.2 1384.8 -678.62   1357.2                             
fitL_1  6 1320.3 1341.5 -654.17   1308.3 48.899      1  2.695e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Example: Poisson

But an interaction would emerge with a classical linear model

Show code

effL = data.frame(allEffects(fitL_1,xlevels=list(age=seq(min(age),max(age),.05)))$"age:group")
ggplot(d,aes(x=age,y=y,shape=group,color=group))+
  geom_point(size=4,alpha=.6)+
  scale_color_manual(values=c("darkorange2","darkgreen"))+
  scale_fill_manual(values=c("brown1","darkslategray"))+
  geom_line(data=effL,aes(x=age,y=fit,group=group,linetype=group),size=2,color=rep(c("brown1","darkslategray"),each=nrow(effL)/2))+
  geom_ribbon(data=effL,aes(x=age,y=fit,ymin=lower,ymax=upper,group=group,fill=group),alpha=.3,color=NA)+
  theme(text=element_text(size=ts,color=tc))+scale_x_continuous(breaks=seq(0,20,.5))+
  ylab("Errors")+xlab("Age (years)")

Example: Poisson

family=gaussian(link="log")

actually… issue with interaction is link function rather than family

fitL_log_0 = glmmTMB(y~age+group+(1|id), family=gaussian(link="log"), data=d, start = c(1, 0, 0))
fitL_log_1 = glmmTMB(y~age*group+(1|id), family=gaussian(link="log"), data=d, start = c(1, 0, 0, 0))
anova(fitL_log_0, fitL_log_1)

Data: d
Models:
fitL_log_0: y ~ age + group + (1 | id), zi=~0, disp=~1
fitL_log_1: y ~ age * group + (1 | id), zi=~0, disp=~1
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
fitL_log_0  5 1228.8 1246.4 -609.41   1218.8                         
fitL_log_1  6 1230.3 1251.5 -609.16   1218.3 0.4946      1     0.4819

Example: Poisson

family=gaussian(link="log")

Show code

fitL_log_1 = glm(y~age*group, family=gaussian(link="log"), data=d, start = c(1, 0, 0, 0))
effL = data.frame(allEffects(fitL_log_1,xlevels=list(age=seq(min(age),max(age),.05)))$"age:group")
ggplot(d,aes(x=age,y=y,shape=group,color=group))+
  geom_point(size=4,alpha=.6)+
  scale_color_manual(values=c("darkorange2","darkgreen"))+
  scale_fill_manual(values=c("brown1","darkslategray"))+
  geom_line(data=effL,aes(x=age,y=fit,group=group,linetype=group),size=2,color=rep(c("brown1","darkslategray"),each=nrow(effL)/2))+
  geom_ribbon(data=effL,aes(x=age,y=fit,ymin=lower,ymax=upper,group=group,fill=group),alpha=.3,color=NA)+
  theme(text=element_text(size=ts,color=tc))+scale_x_continuous(breaks=seq(0,20,.5))+
  ylab("Errors")+xlab("Age (years)")