library(effects) ## load the effects package
data(Cowles, package = "carData")
Model with interaction
mod.cowles <- glm(volunteer ~ sex + neuroticism*extraversion,
data=Cowles, family=binomial)
summary(mod.cowles)
##
## Call:
## glm(formula = volunteer ~ sex + neuroticism * extraversion, family = binomial,
## data = Cowles)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4749 -1.0602 -0.8934 1.2609 1.9978
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.358207 0.501320 -4.704 2.55e-06 ***
## sexmale -0.247152 0.111631 -2.214 0.02683 *
## neuroticism 0.110777 0.037648 2.942 0.00326 **
## extraversion 0.166816 0.037719 4.423 9.75e-06 ***
## neuroticism:extraversion -0.008552 0.002934 -2.915 0.00355 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1933.5 on 1420 degrees of freedom
## Residual deviance: 1897.4 on 1416 degrees of freedom
## AIC: 1907.4
##
## Number of Fisher Scoring iterations: 4
Effect plots
eff.cowles <- allEffects(mod.cowles,
xlevels=list(neuroticism=seq(0, 24, 6),
extraversion=seq(0, 24, 8)))
plot(eff.cowles, 'neuroticism:extraversion',
ylab="Prob(Volunteer)",
ticks=list(at=c(.1, .25, .5, .75, .9)),
layout=c(4,1), aspect=1)
plot(eff.cowles, 'neuroticism:extraversion',
multiline=TRUE,
ylab="Prob(Volunteer)",
key.args=list(x = .8, y = .9))
IycgLS0tDQojJyB0aXRsZTogIkNvd2xlcyBkYXRhOiBFZmZlY3QgcGxvdHMiDQojJyBhdXRob3I6ICJNaWNoYWVsIEZyaWVuZGx5Ig0KIycgZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiDQojJyBvdXRwdXQ6DQojJyAgIGh0bWxfZG9jdW1lbnQ6DQojJyAgICAgdGhlbWU6IHJlYWRhYmxlDQojJyAgICAgY29kZV9kb3dubG9hZDogdHJ1ZQ0KIycgLS0tDQoNCmxpYnJhcnkoZWZmZWN0cykgICAjIyBsb2FkIHRoZSBlZmZlY3RzIHBhY2thZ2UNCmRhdGEoQ293bGVzLCBwYWNrYWdlID0gImNhckRhdGEiKQ0KDQojJyAjIyBNb2RlbCB3aXRoIGludGVyYWN0aW9uDQptb2QuY293bGVzIDwtIGdsbSh2b2x1bnRlZXIgfiBzZXggKyBuZXVyb3RpY2lzbSpleHRyYXZlcnNpb24sIA0KICAgIGRhdGE9Q293bGVzLCBmYW1pbHk9Ymlub21pYWwpDQpzdW1tYXJ5KG1vZC5jb3dsZXMpDQoNCiMnICMjIEVmZmVjdCBwbG90cw0KZWZmLmNvd2xlcyA8LSBhbGxFZmZlY3RzKG1vZC5jb3dsZXMsIA0KCXhsZXZlbHM9bGlzdChuZXVyb3RpY2lzbT1zZXEoMCwgMjQsIDYpLCANCiAgICAgICAgICAgICAgIGV4dHJhdmVyc2lvbj1zZXEoMCwgMjQsIDgpKSkNCg0KcGxvdChlZmYuY293bGVzLCAnbmV1cm90aWNpc206ZXh0cmF2ZXJzaW9uJywgDQogICAgIHlsYWI9IlByb2IoVm9sdW50ZWVyKSIsDQogICAgdGlja3M9bGlzdChhdD1jKC4xLCAuMjUsIC41LCAuNzUsIC45KSksIA0KICAgIGxheW91dD1jKDQsMSksIGFzcGVjdD0xKQ0KDQpwbG90KGVmZi5jb3dsZXMsICduZXVyb3RpY2lzbTpleHRyYXZlcnNpb24nLCANCiAgICAgbXVsdGlsaW5lPVRSVUUsIA0KICAgIHlsYWI9IlByb2IoVm9sdW50ZWVyKSIsIA0KICAgIGtleS5hcmdzPWxpc3QoeCA9IC44LCB5ID0gLjkpKQ0KDQo=