library(vcdExtra)
data("UCBAdmissions")
structable(Dept ~ Admit + Gender, UCBAdmissions)
## Dept A B C D E F
## Admit Gender
## Admitted Male 512 353 120 138 53 22
## Female 89 17 202 131 94 24
## Rejected Male 313 207 205 279 138 351
## Female 19 8 391 244 299 317
berk.mod1 <- loglm(~Dept * (Gender + Admit), data = UCBAdmissions)
berk.mod1
## Call:
## loglm(formula = ~Dept * (Gender + Admit), data = UCBAdmissions)
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 21.74 6 0.001352
## Pearson 19.94 6 0.002840
mosaic(berk.mod1, gp = shading_Friendly)
berk.mod2 <- loglm(~(Admit + Dept + Gender)^2, data = UCBAdmissions)
berk.mod2
## Call:
## loglm(formula = ~(Admit + Dept + Gender)^2, data = UCBAdmissions)
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 20.20 5 0.001144
## Pearson 18.82 5 0.002074
mosaic(berk.mod2, gp = shading_Friendly)
anova(berk.mod1, berk.mod2, test = "Chisq")
## LR tests for hierarchical log-linear models
##
## Model 1:
## ~Dept * (Gender + Admit)
## Model 2:
## ~(Admit + Dept + Gender)^2
##
## Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1 21.74 6
## Model 2 20.20 5 1.531 1 0.21593
## Saturated 0.00 0 20.204 5 0.00114
################## same, using glm() -- need to transform the data to a data.frame
berkeley <- as.data.frame(UCBAdmissions)
berk.glm1 <- glm(Freq ~ Dept * (Gender + Admit), data = berkeley, family = "poisson")
summary(berk.glm1)
##
## Call:
## glm(formula = Freq ~ Dept * (Gender + Admit), family = "poisson",
## data = berkeley)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.478 -0.414 0.010 0.309 2.232
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.2756 0.0425 147.74 < 2e-16 ***
## DeptB -0.4057 0.0677 -5.99 2.1e-09 ***
## DeptC -1.5394 0.0831 -18.54 < 2e-16 ***
## DeptD -1.3223 0.0816 -16.21 < 2e-16 ***
## DeptE -2.4028 0.1101 -21.82 < 2e-16 ***
## DeptF -3.0962 0.1576 -19.65 < 2e-16 ***
## GenderFemale -2.0333 0.1023 -19.87 < 2e-16 ***
## AdmitRejected -0.5935 0.0684 -8.68 < 2e-16 ***
## DeptB:GenderFemale -1.0758 0.2286 -4.71 2.5e-06 ***
## DeptC:GenderFemale 2.6346 0.1234 21.35 < 2e-16 ***
## DeptD:GenderFemale 1.9271 0.1246 15.46 < 2e-16 ***
## DeptE:GenderFemale 2.7548 0.1351 20.39 < 2e-16 ***
## DeptF:GenderFemale 1.9436 0.1268 15.32 < 2e-16 ***
## DeptB:AdmitRejected 0.0506 0.1097 0.46 0.64
## DeptC:AdmitRejected 1.2091 0.0973 12.43 < 2e-16 ***
## DeptD:AdmitRejected 1.2583 0.1015 12.40 < 2e-16 ***
## DeptE:AdmitRejected 1.6830 0.1173 14.34 < 2e-16 ***
## DeptF:AdmitRejected 3.2691 0.1671 19.57 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2650.095 on 23 degrees of freedom
## Residual deviance: 21.736 on 6 degrees of freedom
## AIC: 216.8
##
## Number of Fisher Scoring iterations: 4
mosaic(berk.glm1, gp = shading_Friendly, labeling = labeling_residuals, formula = ~Admit +
Dept + Gender)
berk.glm2 <- glm(Freq ~ (Dept + Gender + Admit)^2, data = berkeley, family = "poisson")
summary(berk.glm2)
##
## Call:
## glm(formula = Freq ~ (Dept + Gender + Admit)^2, family = "poisson",
## data = berkeley)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8 9 10
## -0.7548 0.9947 1.9645 -3.1577 -0.0340 0.0445 0.1571 -0.2203 1.0127 -0.7384
## 11 12 13 14 15 16 17 18 19 20
## -0.7437 0.5490 0.0676 -0.0474 -0.0691 0.0508 1.0558 -0.6124 -0.7362 0.4268
## 21 22 23 24
## -0.2012 0.0511 0.1980 -0.0537
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.2715 0.0427 146.85 < 2e-16 ***
## DeptB -0.4032 0.0678 -5.94 2.8e-09 ***
## DeptC -1.5779 0.0895 -17.63 < 2e-16 ***
## DeptD -1.3500 0.0853 -15.83 < 2e-16 ***
## DeptE -2.4498 0.1176 -20.84 < 2e-16 ***
## DeptF -3.1379 0.1617 -19.40 < 2e-16 ***
## GenderFemale -1.9986 0.1059 -18.87 < 2e-16 ***
## AdmitRejected -0.5821 0.0690 -8.44 < 2e-16 ***
## DeptB:GenderFemale -1.0748 0.2286 -4.70 2.6e-06 ***
## DeptC:GenderFemale 2.6651 0.1261 21.14 < 2e-16 ***
## DeptD:GenderFemale 1.9583 0.1273 15.38 < 2e-16 ***
## DeptE:GenderFemale 2.7952 0.1393 20.07 < 2e-16 ***
## DeptF:GenderFemale 2.0023 0.1357 14.75 < 2e-16 ***
## DeptB:AdmitRejected 0.0434 0.1098 0.40 0.69
## DeptC:AdmitRejected 1.2626 0.1066 11.84 < 2e-16 ***
## DeptD:AdmitRejected 1.2946 0.1058 12.23 < 2e-16 ***
## DeptE:AdmitRejected 1.7393 0.1261 13.79 < 2e-16 ***
## DeptF:AdmitRejected 3.3065 0.1700 19.45 < 2e-16 ***
## GenderFemale:AdmitRejected -0.0999 0.0808 -1.24 0.22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2650.095 on 23 degrees of freedom
## Residual deviance: 20.204 on 5 degrees of freedom
## AIC: 217.3
##
## Number of Fisher Scoring iterations: 4
mosaic.glm(berk.glm2, residuals_type = "rstandard", labeling = labeling_residuals, shade = TRUE,
formula = ~Admit + Dept + Gender, main = "Model: [DeptGender][DeptAdmit][AdmitGender]")
anova(berk.glm1, berk.glm2, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: Freq ~ Dept * (Gender + Admit)
## Model 2: Freq ~ (Dept + Gender + Admit)^2
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 6 21.7
## 2 5 20.2 1 1.53 0.22
berkeley <- within(berkeley, dept1AG <- (Dept == "A") * (Gender == "Female") * (Admit == "Admitted"))
berkeley[1:6, ]
## Admit Gender Dept Freq dept1AG
## 1 Admitted Male A 512 0
## 2 Rejected Male A 313 0
## 3 Admitted Female A 89 1
## 4 Rejected Female A 19 0
## 5 Admitted Male B 353 0
## 6 Rejected Male B 207 0
berk.glm3 <- glm(Freq ~ Dept * (Gender + Admit) + dept1AG, data = berkeley, family = "poisson")
summary(berk.glm3)
##
## Call:
## glm(formula = Freq ~ Dept * (Gender + Admit) + dept1AG, family = "poisson",
## data = berkeley)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8 9 10
## 0.0000 0.0000 0.0000 0.0000 -0.0632 0.0827 0.2951 -0.4009 0.5573 -0.4152
## 11 12 13 14 15 16 17 18 19 20
## -0.4182 0.3051 -0.3066 0.2184 0.3204 -0.2314 0.6984 -0.4142 -0.4992 0.2863
## 21 22 23 24
## -0.4203 0.1086 0.4268 -0.1138
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.2383 0.0442 141.16 < 2e-16 ***
## DeptB -0.3685 0.0688 -5.36 8.5e-08 ***
## DeptC -1.5021 0.0839 -17.89 < 2e-16 ***
## DeptD -1.2851 0.0825 -15.58 < 2e-16 ***
## DeptE -2.3655 0.1108 -21.35 < 2e-16 ***
## DeptF -3.0590 0.1580 -19.36 < 2e-16 ***
## GenderFemale -2.8018 0.2363 -11.86 < 2e-16 ***
## AdmitRejected -0.4921 0.0717 -6.86 6.9e-12 ***
## dept1AG 1.0521 0.2627 4.00 6.2e-05 ***
## DeptB:GenderFemale -0.3073 0.3124 -0.98 0.33
## DeptC:GenderFemale 3.4031 0.2461 13.83 < 2e-16 ***
## DeptD:GenderFemale 2.6956 0.2468 10.92 < 2e-16 ***
## DeptE:GenderFemale 3.5233 0.2522 13.97 < 2e-16 ***
## DeptF:GenderFemale 2.7121 0.2479 10.94 < 2e-16 ***
## DeptB:AdmitRejected -0.0507 0.1118 -0.45 0.65
## DeptC:AdmitRejected 1.1078 0.0997 11.12 < 2e-16 ***
## DeptD:AdmitRejected 1.1570 0.1038 11.14 < 2e-16 ***
## DeptE:AdmitRejected 1.5816 0.1193 13.25 < 2e-16 ***
## DeptF:AdmitRejected 3.1678 0.1685 18.80 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2650.0952 on 23 degrees of freedom
## Residual deviance: 2.6815 on 5 degrees of freedom
## AIC: 199.7
##
## Number of Fisher Scoring iterations: 3
mosaic.glm(berk.glm3, residuals_type = "rstandard", labeling = labeling_residuals, shade = TRUE,
formula = ~Admit + Dept + Gender, main = "Model: [DeptGender][DeptAdmit] + DeptA*[GA]")
