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.03 < 2.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.4328 -13.2031 -0.0277 15.9185 21.2218
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.59346 0.06838 8.679 <2e-16 ***
## DeptB -0.05059 0.10968 -0.461 0.645
## DeptC -1.20915 0.09726 -12.432 <2e-16 ***
## DeptD -1.25833 0.10152 -12.395 <2e-16 ***
## DeptE -1.68296 0.11733 -14.343 <2e-16 ***
## DeptF -3.26911 0.16707 -19.567 <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
Main effects model
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.708 <2e-16 ***
## Gender 1 1.526 0.2167
## 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.3424 -13.0584 -0.1631 16.0167 21.3199
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.58205 0.06899 8.436 <2e-16 ***
## DeptB -0.04340 0.10984 -0.395 0.693
## DeptC -1.26260 0.10663 -11.841 <2e-16 ***
## DeptD -1.29461 0.10582 -12.234 <2e-16 ***
## DeptE -1.73931 0.12611 -13.792 <2e-16 ***
## DeptF -3.30648 0.16998 -19.452 <2e-16 ***
## GenderFemale 0.09987 0.08085 1.235 0.217
## ---
## 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.5
##
## 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'
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.5939 <2e-16 ***
## Gender 1 1.7750 0.1828
## dept1AG 1 0.0056 0.9401
## Residuals 16
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(effect('Dept:Gender', berk.mod3), multiline=TRUE)
