Cowles data: Effect plots

library(effects)   ## load the effects package
data(Cowles)

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

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))