Output from berkeley-logit1.R

library(car)  # for Anova()
data(UCBAdmissions)
UCB.df <- as.data.frame(UCBAdmissions)
berk.mod1 <- glm(Admit == "Admitted" ~ Dept, data = UCB.df, weights = UCB.df$Freq, family = "binomial")
Anova(berk.mod1, test = "Wald")

## Analysis of Deviance Table (Type II tests)
## 
## Response: Admit == "Admitted"
##           Df Chisq Pr(>Chisq)    
## Dept       5   623     <2e-16 ***
## Residuals 18                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(berk.mod1)

## 
## Call:
## glm(formula = Admit == "Admitted" ~ Dept, family = "binomial", 
##     data = UCB.df, weights = UCB.df$Freq)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -25.433  -13.203   -0.028   15.919   21.222  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.5935     0.0684    8.68   <2e-16 ***
## DeptB        -0.0506     0.1097   -0.46     0.64    
## DeptC        -1.2091     0.0973  -12.43   <2e-16 ***
## DeptD        -1.2583     0.1015  -12.40   <2e-16 ***
## DeptE        -1.6830     0.1173  -14.34   <2e-16 ***
## DeptF        -3.2691     0.1671  -19.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6044.3  on 23  degrees of freedom
## Residual deviance: 5189.0  on 18  degrees of freedom
## AIC: 5201
## 
## Number of Fisher Scoring iterations: 6

berk.mod2 <- glm(Admit == "Admitted" ~ Dept + Gender, data = UCB.df, weights = UCB.df$Freq, 
    family = "binomial")
Anova(berk.mod2, test = "Wald")

## Analysis of Deviance Table (Type II tests)
## 
## Response: Admit == "Admitted"
##           Df  Chisq Pr(>Chisq)    
## Dept       5 534.71     <2e-16 ***
## Gender     1   1.53       0.22    
## Residuals 17                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(berk.mod2)

## 
## Call:
## glm(formula = Admit == "Admitted" ~ Dept + Gender, family = "binomial", 
##     data = UCB.df, weights = UCB.df$Freq)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -25.342  -13.058   -0.163   16.017   21.320  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.5821     0.0690    8.44   <2e-16 ***
## DeptB         -0.0434     0.1098   -0.40     0.69    
## DeptC         -1.2626     0.1066  -11.84   <2e-16 ***
## DeptD         -1.2946     0.1058  -12.23   <2e-16 ***
## DeptE         -1.7393     0.1261  -13.79   <2e-16 ***
## DeptF         -3.3065     0.1700  -19.45   <2e-16 ***
## GenderFemale   0.0999     0.0808    1.24     0.22    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6044.3  on 23  degrees of freedom
## Residual deviance: 5187.5  on 17  degrees of freedom
## AIC: 5201
## 
## Number of Fisher Scoring iterations: 6

library(effects)  ## load the effects package
berk.eff2 <- allEffects(berk.mod2)
# plot main effects
plot(berk.eff2)

# plot 'interaction' -- produces a harmless warning
plot(effect("Dept:Gender", berk.mod2), multiline = TRUE)

# include a 1 df term for dept A
UCB.df <- within(UCB.df, dept1AG <- (Dept == "A") * (Gender == "Female") * (Admit == "Admitted"))
berk.mod3 <- glm(Admit == "Admitted" ~ Dept + Gender + dept1AG, data = UCB.df, weights = UCB.df$Freq, 
    family = "binomial")
Anova(berk.mod3, test = "Wald")

## Analysis of Deviance Table (Type II tests)
## 
## Response: Admit == "Admitted"
##           Df  Chisq Pr(>Chisq)    
## Dept       5 428.59     <2e-16 ***
## Gender     1   1.78       0.18    
## dept1AG    1   0.01       0.94    
## Residuals 16                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(effect("Dept:Gender", berk.mod3), multiline = TRUE)