##Load the data year is numeric, but we want to treat it as a factor

data(Arrests, package = "carData")
Arrests$year <- as.factor(Arrests$year)

Main effects model, with two interactions

arrests.mod <- glm(released ~ employed + citizen + checks +
                   colour*year + colour*age,
                   family=binomial, data=Arrests)

Anova(arrests.mod)
## Analysis of Deviance Table (Type II tests)
## 
## Response: released
##             LR Chisq Df Pr(>Chisq)    
## employed      72.673  1  < 2.2e-16 ***
## citizen       25.783  1  3.820e-07 ***
## checks       205.211  1  < 2.2e-16 ***
## colour        19.572  1  9.687e-06 ***
## year           6.087  5  0.2978477    
## age            0.459  1  0.4982736    
## colour:year   21.720  5  0.0005917 ***
## colour:age    13.886  1  0.0001942 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot all effects

arrests.effects <- allEffects(arrests.mod, 
                              xlevels=list(age=seq(15,45,5)))
plot(arrests.effects, 
     ylab="Probability(released)")

Plot 3-way effect (not in model)

plot(effect("colour:year:age", arrests.mod, 
            xlevels=list(age=15:45)),
     multiline=TRUE, 
     ylab="Probability(released)", 
     rug=FALSE)

colour x year interaction

plot(effect("colour:year", arrests.mod),
      multiline=TRUE, ylab="Probability(released)")

colour x age interaction

plot(effect("colour:age", arrests.mod),
      multiline=FALSE, ylab="Probability(released)")

IycgLS0tDQojJyB0aXRsZTogIkFycmVzdHMgZGF0YTogRWZmZWN0cyBwbG90cyINCiMnIGF1dGhvcjogIk1pY2hhZWwgRnJpZW5kbHkiDQojJyBkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCINCiMnIG91dHB1dDoNCiMnICAgaHRtbF9kb2N1bWVudDoNCiMnICAgICB0aGVtZTogcmVhZGFibGUNCiMnICAgICBjb2RlX2Rvd25sb2FkOiB0cnVlDQojJyAtLS0NCg0KIysgZWNobz1GQUxTRQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KA0KICB3YXJuaW5nID0gRkFMU0UsICAgIyBhdm9pZCB3YXJuaW5ncyBhbmQgbWVzc2FnZXMgaW4gdGhlIG91dHB1dA0KICBtZXNzYWdlID0gRkFMU0UNCikNCg0KDQojIEV4YW1wbGVzIGZyb20gSm9obiBGb3gncyBKU1MgcGFwZXINCg0KbGlicmFyeShlZmZlY3RzKSAgICAjIGVmZmVjdCBwbG90cw0KbGlicmFyeShjYXIpICAgICAgICAjIGZvciBBbm92YSgpDQoNCiMnICMjTG9hZCB0aGUgZGF0YQ0KIycgYHllYXJgIGlzIG51bWVyaWMsIGJ1dCB3ZSB3YW50IHRvIHRyZWF0IGl0IGFzIGEgZmFjdG9yDQpkYXRhKEFycmVzdHMsIHBhY2thZ2UgPSAiY2FyRGF0YSIpDQpBcnJlc3RzJHllYXIgPC0gYXMuZmFjdG9yKEFycmVzdHMkeWVhcikNCg0KIycgIyMgTWFpbiBlZmZlY3RzIG1vZGVsLCB3aXRoIHR3byBpbnRlcmFjdGlvbnMNCmFycmVzdHMubW9kIDwtIGdsbShyZWxlYXNlZCB+IGVtcGxveWVkICsgY2l0aXplbiArIGNoZWNrcyArDQogICAgICAgICAgICAgICAgICAgY29sb3VyKnllYXIgKyBjb2xvdXIqYWdlLA0KICAgICAgICAgICAgICAgICAgIGZhbWlseT1iaW5vbWlhbCwgZGF0YT1BcnJlc3RzKQ0KDQpBbm92YShhcnJlc3RzLm1vZCkNCg0KIycgIyMjIHBsb3QgYWxsIGVmZmVjdHMNCiMrIGZpZy53aWR0aD05LCBmaWcuaGVpZ2h0PTQNCmFycmVzdHMuZWZmZWN0cyA8LSBhbGxFZmZlY3RzKGFycmVzdHMubW9kLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHhsZXZlbHM9bGlzdChhZ2U9c2VxKDE1LDQ1LDUpKSkNCnBsb3QoYXJyZXN0cy5lZmZlY3RzLCANCiAgICAgeWxhYj0iUHJvYmFiaWxpdHkocmVsZWFzZWQpIikNCg0KIycgIyMjIFBsb3QgMy13YXkgZWZmZWN0IChub3QgaW4gbW9kZWwpDQpwbG90KGVmZmVjdCgiY29sb3VyOnllYXI6YWdlIiwgYXJyZXN0cy5tb2QsIA0KICAgICAgICAgICAgeGxldmVscz1saXN0KGFnZT0xNTo0NSkpLA0KICAgICBtdWx0aWxpbmU9VFJVRSwgDQogICAgIHlsYWI9IlByb2JhYmlsaXR5KHJlbGVhc2VkKSIsIA0KICAgICBydWc9RkFMU0UpDQoNCiMnICMjIyBjb2xvdXIgeCB5ZWFyIGludGVyYWN0aW9uDQpwbG90KGVmZmVjdCgiY29sb3VyOnllYXIiLCBhcnJlc3RzLm1vZCksDQogICAgICBtdWx0aWxpbmU9VFJVRSwgeWxhYj0iUHJvYmFiaWxpdHkocmVsZWFzZWQpIikNCg0KIycgIyMjIGNvbG91ciB4IGFnZSBpbnRlcmFjdGlvbg0KcGxvdChlZmZlY3QoImNvbG91cjphZ2UiLCBhcnJlc3RzLm1vZCksDQogICAgICBtdWx0aWxpbmU9RkFMU0UsIHlsYWI9IlByb2JhYmlsaXR5KHJlbGVhc2VkKSIpDQoNCiAgICAg