> admissions <- read.table("admissions.dat",header=T)
> attach(admissions, pos=1)
> admissions
Dept gender admitted count
1 a m yes 512
2 a m no 313
3 a f yes 89
4 a f no 19
5 b m yes 353
6 b m no 207
7 b f yes 17
8 b f no 8
9 c m yes 120
10 c m no 205
11 c f yes 202
12 c f no 391
13 d m yes 138
14 d m no 279
15 d f yes 131
16 d f no 244
17 e m yes 53
18 e m no 138
19 e f yes 94
20 e f no 299
21 f m yes 22
22 f m no 351
23 f f yes 24
24 f f no 317
> d <- factor(Dept); g <- factor(gender); a <- factor(admitted)
> fit <- glm(count ~ a*d + a*g + d*g, family=poisson)
> fit
Degrees of Freedom: 24 Total; 5 Residual
Residual Deviance: 20.20428
> fit.indep <- update(fit, .~. - a:g) # conditional indep. of admit and gender
> fit.indep
Degrees of Freedom: 24 Total; 6 Residual
Residual Deviance: 21.73551
> pearson <- summary.lm(fit.indep)$residuals # Pearson residuals
> hat <- lm.influence(fit.indep)$hat # hat or leverage values
> adjusted <- pearson/sqrt(1 - hat) # standardized Pearson residuals
> cbind(Dept, gender, admitted, count, fitted(fit.indep), pearson, adjusted)
Dept gender admitted count pearson adjusted
1 1 2 2 512 531.430846 -0.84290731 -4.1465776
2 1 2 1 313 293.569165 1.13391310 4.1465776
3 1 1 2 89 69.569537 2.32575149 4.1465776
4 1 1 1 19 38.431098 -3.12869525 -4.1465776
5 2 2 2 353 354.188034 -0.06312654 -0.5037048
6 2 2 1 207 205.811966 0.08281207 0.5037048
7 2 1 2 17 15.811966 0.29876748 0.5037048
8 2 1 1 8 9.188034 -0.39193581 -0.5037048
9 3 2 2 120 113.997821 0.56216078 0.8680661
10 3 2 1 205 211.002179 -0.41320484 -0.8680661
11 3 1 2 202 208.002179 -0.41617398 -0.8680661
12 3 1 1 391 384.997821 0.30590022 0.8680661
13 4 2 2 138 141.632576 -0.30523413 -0.5458732
14 4 2 1 279 275.367424 0.21890637 0.5458732
15 4 1 2 131 127.367424 0.32187369 0.5458732
16 4 1 1 244 247.632576 -0.23083985 -0.5458732
17 5 2 2 53 48.077055 0.70999484 1.0005327
18 5 2 1 138 142.922945 -0.41178810 -1.0005327
19 5 1 2 94 98.922945 -0.49496658 -1.0005327
20 5 1 1 299 294.077055 0.28707441 1.0005327
21 6 2 2 22 24.030812 -0.41427041 -0.6197504
22 6 2 1 351 348.969188 0.10871169 0.6197504
23 6 1 2 24 21.969188 0.43327252 0.6197504
24 6 1 1 317 319.030812 -0.11369817 -0.6197504
> dept.a <- c(1,1,1,1,rep(0,20))
> dept.a
[1] 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> update(fit.indep, .~. + a:g:dept.a) # add a*g assoc. only for Dept. a
Degrees of Freedom: 24 Total; 5 Residual
Residual Deviance: 2.681497
> binom.dat <- t(matrix(count, ncol=12)) # now look at equivalent logit models
# note t = transpose
> binom.dat
[,1] [,2]
[1,] 512 313
[2,] 89 19
[3,] 353 207
[4,] 17 8
[5,] 120 205
[6,] 202 391
[7,] 138 279
[8,] 131 244
[9,] 53 138
[10,] 94 299
[11,] 22 351
[12,] 24 317
> fit.b <- glm(binom.dat ~ dep + gen, family=binomial)
> fit.b
Degrees of Freedom: 12 Total; 5 Residual
Residual Deviance: 20.20428
> fit.indep.b <- update(fit.b, .~. - gen)
> fit.indep.b
Degrees of Freedom: 12 Total; 6 Residual
Residual Deviance: 21.73551
> pearson <- summary.lm(fit.indep.b)$residuals # Pearson residuals
> hat.b <- lm.influence(fit.indep.b)$hat # hat or leverage values
> pearson.b <- summary.lm(fit.indep.b)$residuals # Pearson residuals
> adjusted.b <- pearson.b/sqrt(1 - hat.b) # standardized Pearson residuals
> cbind(dep, gen, binom.dat, fitted(fit.indep.b), pearson.b, adjusted.b)
dep gen pearson.b adjusted.b
1 1 1 512 313 0.64415863 -1.4129812 -4.1530324
2 1 2 89 19 0.64415863 3.9052736 4.1530324
3 2 1 353 207 0.63247863 -0.1041288 -0.5037077
4 2 2 17 8 0.63247863 0.4928272 0.5037077
5 3 1 120 205 0.35076253 0.6976841 0.8680662
6 3 2 202 391 0.35076253 -0.5165034 -0.8680662
7 4 1 138 279 0.33964646 -0.3756167 -0.5458732
8 4 2 131 244 0.33964646 0.3960931 0.5458732
9 5 1 53 138 0.25171233 0.8207703 1.0005339
10 5 2 94 299 0.25171233 -0.5721924 -1.0005339
11 6 1 22 351 0.06442577 -0.4282971 -0.6197508
12 6 2 24 317 0.06442577 0.4479426 0.6197508