> Berk.grp = data.frame(Sex = rep(c("Male","Female"),each=6), + Program = rep(c("A","B","C","D","E","F"),2), + Admit = c(511,352,120,137,53,22, + 89,17,202,132,95,24), + Deny = c(314,208,205,270,138,351, + 19,8,391,243,298,317)) > Berk.grp Sex Program Admit Deny 1 Male A 511 314 2 Male B 352 208 3 Male C 120 205 4 Male D 137 270 5 Male E 53 138 6 Male F 22 351 7 Female A 89 19 8 Female B 17 8 9 Female C 202 391 10 Female D 132 243 11 Female E 95 298 12 Female F 24 317 > (511*19)/(314*89) [1] 0.34742 > mod4 <- glm(cbind(Admit,Deny) ~ Sex*Program, family=binomial, data=Berk.grp) > summary(mod4) Call: glm(formula = cbind(Admit, Deny) ~ Sex * Program, family = binomial, data = Berk.grp) Deviance Residuals: [1] 0 0 0 0 0 0 0 0 0 0 0 0 Coefficients: Estimate Std. Error z value (Intercept) 1.5442 0.2527 6.110 SexMale -1.0572 0.2627 -4.025 ProgramB -0.7904 0.4977 -1.588 ProgramC -2.2046 0.2672 -8.252 ProgramD -2.1545 0.2749 -7.838 ProgramE -2.6874 0.2788 -9.638 ProgramF -4.1250 0.3297 -12.512 SexMale:ProgramB 0.8295 0.5104 1.625 SexMale:ProgramC 1.1821 0.2995 3.946 SexMale:ProgramD 0.9890 0.3028 3.266 SexMale:ProgramE 1.2435 0.3302 3.766 SexMale:ProgramF 0.8683 0.4027 2.156 Pr(>|z|) (Intercept) 9.94e-10 *** SexMale 5.71e-05 *** ProgramB 0.112241 ProgramC < 2e-16 *** ProgramD 4.58e-15 *** ProgramE < 2e-16 *** ProgramF < 2e-16 *** SexMale:ProgramB 0.104085 SexMale:ProgramC 7.93e-05 *** SexMale:ProgramD 0.001091 ** SexMale:ProgramE 0.000166 *** SexMale:ProgramF 0.031046 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 8.7026e+02 on 11 degrees of freedom Residual deviance: -3.7526e-14 on 0 degrees of freedom AIC: 92.938 Number of Fisher Scoring iterations: 3 > mod1 <- glm(cbind(Admit,Deny) ~ Sex, family=binomial, data=Berk.grp) > mod2 <- glm(cbind(Admit,Deny) ~ Program, family=binomial, data=Berk.grp) > mod3 <- glm(cbind(Admit,Deny) ~ Sex+Program, family=binomial, data=Berk.grp) > anova(mod4, test="LRT") Analysis of Deviance Table Model: binomial, link: logit Response: cbind(Admit, Deny) Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 11 870.26 Sex 1 92.51 10 777.75 < 2.2e-16 *** Program 5 757.52 5 20.23 < 2.2e-16 *** Sex:Program 5 20.23 0 0.00 0.001131 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > summary(mod3) Call: glm(formula = cbind(Admit, Deny) ~ Sex + Program, family = binomial, data = Berk.grp) Deviance Residuals: 1 2 3 4 5 6 7 -1.2574 -0.0600 1.2485 0.1427 1.1625 -0.2096 3.7442 8 9 10 11 12 0.2900 -0.9208 -0.1466 -0.8097 0.2072 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.67647 0.09909 6.827 8.69e-12 *** SexMale -0.09900 0.08087 -1.224 0.221 ProgramB -0.04613 0.10975 -0.420 0.674 ProgramC -1.25745 0.10660 -11.795 < 2e-16 *** ProgramD -1.27089 0.10610 -11.978 < 2e-16 *** ProgramE -1.72506 0.12593 -13.699 < 2e-16 *** ProgramF -3.30147 0.16996 -19.425 < 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: 870.26 on 11 degrees of freedom Residual deviance: 20.23 on 5 degrees of freedom AIC: 103.17 Number of Fisher Scoring iterations: 4 > deviance(mod3) [1] 20.23002 > sum( residuals(mod3, type="deviance")^2 ) [1] 20.23002 > sum( residuals(mod3, type="pearson")^2 ) [1] 18.82189 > plot(residuals(mod3, type="deviance") ~ fitted(mod3), ylab="Deviance Residuals", xlab="Fitted Probabilities") > abline(h=0, type=2, col="red") Warning message: In int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : graphical parameter "type" is obsolete > identify(residuals(mod3, type="deviance") ~ fitted(mod3)) ## Click on points then hit Esc [1] 3 5 7 Warning message: graphical parameter "type" is obsolete > Berk.grp[7,] Sex Program Admit Deny 7 Female A 89 19 > 89/(89+19) [1] 0.8240741 > fitted(mod3)[7] 7 0.6629512 > mod3 Call: glm(formula = cbind(Admit, Deny) ~ Sex + Program, family = binomial, data = Berk.grp) Coefficients: (Intercept) SexMale ProgramB ProgramC 0.67647 -0.09900 -0.04613 -1.25745 ProgramD ProgramE ProgramF -1.27089 -1.72506 -3.30147 Degrees of Freedom: 11 Total (i.e. Null); 5 Residual Null Deviance: 870.3 Residual Deviance: 20.23 AIC: 103.2 > beta.hat <- coef(mod3) > beta.hat (Intercept) SexMale ProgramB ProgramC 0.67647380 -0.09899799 -0.04613349 -1.25745307 ProgramD ProgramE ProgramF -1.27089482 -1.72505949 -3.30146910 > V.hat <- vcov(mod3) > V.hat (Intercept) SexMale ProgramB (Intercept) 0.009819515 -0.0058028642 -0.0042580822 SexMale -0.005802864 0.0065400235 -0.0004650577 ProgramB -0.004258082 -0.0004650577 0.0120441161 ProgramC -0.007805033 0.0035324758 0.0044195283 ProgramD -0.006844001 0.0024493596 0.0044965482 ProgramE -0.007984295 0.0037345093 0.0044051618 ProgramF -0.006913046 0.0025271758 0.0044910147 ProgramC ProgramD ProgramE (Intercept) -0.007805033 -0.006844001 -0.007984295 SexMale 0.003532476 0.002449360 0.003734509 ProgramB 0.004419528 0.004496548 0.004405162 ProgramC 0.011364587 0.005993698 0.006687849 ProgramD 0.005993698 0.011257973 0.006069364 ProgramE 0.006687849 0.006069364 0.015857180 ProgramF 0.006035729 0.005617195 0.006113798 ProgramF (Intercept) -0.006913046 SexMale 0.002527176 ProgramB 0.004491015 ProgramC 0.006035729 ProgramD 0.005617195 ProgramE 0.006113798 ProgramF 0.028886835 > A = matrix(c(0,0,1,0,-1,0,0),nrow=1) > A [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0 0 1 0 -1 0 0 > A %*% beta.hat # On log scale [,1] [1,] 1.224761 > exp(A %*% beta.hat) # On odds scale [,1] [1,] 3.403354 > sqrt( A %*% V.hat %*% t(A) ) [,1] [1,] 0.1196202 > exp( A%*%beta.hat + c(-1,1)*1.96*sqrt( A %*% V.hat %*% t(A) )) [1] 2.692057 4.302589 > source("http://www.math.montana.edu/shancock/courses/stat539/r/GillenRFunctions.R") > linContr.glm(contr.names=c("ProgramB","ProgramD"), contr.coef=c(1, -1), model=mod3) Test of H_0: exp( 1*ProgramB + -1*ProgramD ) = 1 : exp( Est ) se.est zStat pVal ci95.lo ci95.hi 1 3.403 0.12 10.239 0 2.692 4.303 >