> remission<-data.frame(rem=c(18,27,14,25),sex=factor(c(1,1,2,2)),treat=factor(c (1,2,1,2)),n=c(30,30,30,30)) > remission rem sex treat n 1 18 1 1 30 2 27 1 2 30 3 14 2 1 30 4 25 2 2 30
> fit1<-glm(cbind(rem,n-rem)~sex*treat,family = binomial(link = logit),
+ data=remission,subset=n>0)
> summary(fit1)
Call: glm(formula = cbind(rem, n - rem) ~ sex * treat, family = binomial(link =
logit), data = remission, subset = n > 0)
Coefficients:
Value Std. Error t value
(Intercept) 1.01964905 0.2349443 4.33996088
sex -0.28169579 0.2349443 -1.19898970
treat 0.88368219 0.2349443 3.76124133
sex:treat -0.01219754 0.2349443 -0.05191674
(Dispersion Parameter for Binomial family taken to be 1 )
Null Deviance: 18.23268 on 3 degrees of freedom
Residual Deviance: 0 on 0 degrees of freedom
Number of Fisher Scoring Iterations: 3
Correlation of Coefficients:
(Intercept) sex treat
sex -0.1532238
treat 0.3821937 -0.1419910
sex:treat -0.1419910 0.3821937 -0.1532238
> anova(fit1,test="Chisq")
Analysis of Deviance Table
Binomial model
Response: cbind(rem, n - rem)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(Chi)
NULL 3 18.23268
sex 1 1.43362 2 16.79906 0.2311748
treat 1 16.79637 1 0.00270 0.0000416
sex:treat 1 0.00270 0 0.00000 0.9585671
> fit2<-update(fit1,.~.-sex:treat)
> summary(fit2)
Call: glm(formula = cbind(rem, n - rem) ~ sex + treat, family = binomial(link =
logit), data = remission, subset = n > 0)
Deviance Residuals:
-0.02064035 0.03350603 0.02026137 -0.02719909
Coefficients:
Value Std. Error t value
(Intercept) 1.0179709 0.2323865 4.380508
sex -0.2770531 0.2169754 -1.276887
treat 0.8818647 0.2320043 3.801070
(Dispersion Parameter for Binomial family taken to be 1 )
Null Deviance: 18.23268 on 3 degrees of freedom
Residual Deviance: 0.002699 on 1 degrees of freedom
Number of Fisher Scoring Iterations: 3
Correlation of Coefficients:
(Intercept) sex
sex -0.1028273
treat 0.3668905 -0.0855042
> anova(fit2,test="Chisq")
Analysis of Deviance Table
Binomial model
Response: cbind(rem, n - rem)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(Chi)
NULL 3 18.23268
sex 1 1.43362 2 16.79906 0.2311748
treat 1 16.79637 1 0.00270 0.0000416
This analysis corresponds exactly to the derivation in McCullagh & Nelder, where the mean is the binomial proportion.
> fit3<-glm(rem/n~sex*treat,family=binomial(link=logit), data = remission, subse
t = n > 0, weight=n)
> summary(fit3)
Call: glm(formula = rem/n ~ sex * treat, family = binomial(link = logit), data =
remission, weights = n, subset = n > 0)
Coefficients:
Value Std. Error t value
(Intercept) 1.01964905 0.2349448 4.33995087
sex -0.28169579 0.2349448 -1.19898694
treat 0.88368219 0.2349448 3.76123265
sex:treat -0.01219754 0.2349448 -0.05191662
(Dispersion Parameter for Binomial family taken to be 1 )
Null Deviance: 18.23268 on 3 degrees of freedom
Residual Deviance: 0 on 0 degrees of freedom
Number of Fisher Scoring Iterations: 5
Correlation of Coefficients:
(Intercept) sex treat
sex -0.1532273
treat 0.3821965 -0.1419945
sex:treat -0.1419945 0.3821965 -0.1532273
> anova(fit3,test="Chisq")
Analysis of Deviance Table
Binomial model
Response: rem/n
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(Chi)
NULL 3 18.23268
sex 1 1.43362 2 16.79906 0.2311748
treat 1 16.79637 1 0.00270 0.0000416
sex:treat 1 0.00270 0 0.00000 0.9585671
> fit4<-update(fit3,.~.-sex:treat)
> summary(fit4)
Call: glm(formula = rem/n ~ sex + treat, family = binomial(link = logit), data =
remission, weights = n, subset = n > 0)
Deviance Residuals:
1 2 3 4
-0.02064035 0.03350603 0.02026137 -0.02719909
Coefficients:
Value Std. Error t value
(Intercept) 1.0179709 0.2323868 4.380503
sex -0.2770531 0.2169756 -1.276886
treat 0.8818647 0.2320046 3.801066
(Dispersion Parameter for Binomial family taken to be 1 )
Null Deviance: 18.23268 on 3 degrees of freedom
Residual Deviance: 0.002699 on 1 degrees of freedom
Number of Fisher Scoring Iterations: 5
Correlation of Coefficients:
(Intercept) sex
sex -0.1028275
treat 0.3668920 -0.0855044
> anova(fit4,test="Chisq")
Analysis of Deviance Table
Binomial model
Response: rem/n
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(Chi)
NULL 3 18.23268
sex 1 1.43362 2 16.79906 0.2311748
treat 1 16.79637 1 0.00270 0.0000416
> remission.expand <- data.frame(rem = c(rep(1, sum(rem)), rep(0, sum(n - rem))),
sex = c(rep(sex, rem), rep(sex, n - rem)), treat = c(rep(treat, rem),
rep(treat, n - rem)))
> remission.expand
rem sex treat
1 1 1 1
2 1 1 1
3 1 1 1
4 1 1 1
5 1 1 1
6 1 1 1
7 1 1 1
8 1 1 1
9 1 1 1
10 1 1 1
11 1 1 1
12 1 1 1
13 1 1 1
14 1 1 1
15 1 1 1
16 1 1 1
17 1 1 1
18 1 1 1
19 1 1 2
20 1 1 2
21 1 1 2
22 1 1 2
23 1 1 2
24 1 1 2
25 1 1 2
26 1 1 2
27 1 1 2
28 1 1 2
29 1 1 2
30 1 1 2
31 1 1 2
32 1 1 2
33 1 1 2
34 1 1 2
35 1 1 2
36 1 1 2
37 1 1 2
38 1 1 2
39 1 1 2
40 1 1 2
41 1 1 2
42 1 1 2
43 1 1 2
44 1 1 2
45 1 1 2
46 1 2 1
47 1 2 1
rem sex treat
48 1 2 1
49 1 2 1
50 1 2 1
51 1 2 1
52 1 2 1
53 1 2 1
54 1 2 1
55 1 2 1
56 1 2 1
57 1 2 1
58 1 2 1
59 1 2 1
60 1 2 2
61 1 2 2
62 1 2 2
63 1 2 2
64 1 2 2
65 1 2 2
66 1 2 2
67 1 2 2
68 1 2 2
69 1 2 2
70 1 2 2
71 1 2 2
72 1 2 2
73 1 2 2
74 1 2 2
75 1 2 2
76 1 2 2
77 1 2 2
78 1 2 2
79 1 2 2
80 1 2 2
81 1 2 2
82 1 2 2
83 1 2 2
84 1 2 2
85 0 1 1
86 0 1 1
87 0 1 1
88 0 1 1
89 0 1 1
90 0 1 1
91 0 1 1
92 0 1 1
93 0 1 1
94 0 1 1
rem sex treat
95 0 1 1
96 0 1 1
97 0 1 2
98 0 1 2
99 0 1 2
100 0 2 1
101 0 2 1
102 0 2 1
103 0 2 1
104 0 2 1
105 0 2 1
106 0 2 1
107 0 2 1
108 0 2 1
109 0 2 1
110 0 2 1
111 0 2 1
112 0 2 1
113 0 2 1
114 0 2 1
115 0 2 1
116 0 2 2
117 0 2 2
118 0 2 2
119 0 2 2
120 0 2 2
> remission.expand$sex<-factor(remission.expand$sex)
> remission.expand$treat<-factor(remission.expand$treat)
Note that when we fit the saturated model we get the same estimates and tests as we got from the previous fits, but there are now 116 degress of freedom for the residual deviance. Because the observed scores are either 0 or 1 but, in this example, the fitted binomial probabilities are never 0 or 1, the residual deviance can't be zero.
> fit5<-glm(rem~sex*treat,data=remission.expand,family=binomial(link=logit))
> summary(fit5)
Call: glm(formula = rem ~ sex * treat, family = binomial(link = logit), data =
remission.expand)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.145962 -1.121257 0.4590455 1.010768 1.234617
Coefficients:
Value Std. Error t value
(Intercept) 1.01964685 0.2347573 4.34340879
sex -0.28169360 0.2347573 -1.19993549
treat 0.88368000 0.2347573 3.76422822
sex:treat -0.01219535 0.2347573 -0.05194875
(Dispersion Parameter for Binomial family taken to be 1 )
Null Deviance: 146.6074 on 119 degrees of freedom
Residual Deviance: 128.3747 on 116 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
(Intercept) sex treat
sex -0.1519319
treat 0.3812090 -0.1406812
sex:treat -0.1406812 0.3812090 -0.1519319
> anova(fit5,test="Chisq")
Analysis of Deviance Table
Binomial model
Response: rem
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(Chi)
NULL 119 146.6074
sex 1 1.43362 118 145.1738 0.2311748
treat 1 16.79637 117 128.3774 0.0000416
sex:treat 1 0.00270 116 128.3747 0.9585671
> fit6<-update(fit5,.~.-sex:treat)
> summary(fit6)
Call: glm(formula = rem ~ sex + treat, family = binomial(link = logit), data =
remission.expand)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.137429 -1.118173 0.463493 1.007725 1.237822
Coefficients:
Value Std. Error t value
(Intercept) 1.0179706 0.2323177 4.381805
sex -0.2770530 0.2169430 -1.277078
treat 0.8818644 0.2319361 3.802188
(Dispersion Parameter for Binomial family taken to be 1 )
Null Deviance: 146.6074 on 119 degrees of freedom
Residual Deviance: 128.3774 on 117 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
(Intercept) sex
sex -0.1027515
treat 0.3665306 -0.0854362
> anova(fit6,test="Chisq")
Analysis of Deviance Table
Binomial model
Response: rem
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(Chi)
NULL 119 146.6074
sex 1 1.43362 118 145.1738 0.2311748
treat 1 16.79637 117 128.3774 0.0000416
This analysis gives exactly the same residual deviances as the previous analysis, but the residual from the saturated model has only 4 degrees of freedom.
> remission.group
rem count sex treat
1 1 18 1 1
2 0 12 1 1
3 1 27 1 2
4 0 3 1 2
5 1 14 2 1
6 0 16 2 1
7 1 25 2 2
8 0 5 2 2
> fit7<-glm(rem~sex*treat,data=remission.group,family=binomial(link=logit),
+ weight=count)
> summary(fit7)
Call: glm(formula = rem ~ sex * treat, family = binomial(link = logit), data =
remission.group, weights = count)
Deviance Residuals:
1 2 3 4 5 6 7 8
4.288324 -4.689454 2.38527 -3.716916 4.619515 -4.485028 3.019284 -4.232918
Coefficients:
Value Std. Error t value
(Intercept) 1.01964685 0.2347573 4.34340879
sex -0.28169360 0.2347573 -1.19993549
treat 0.88368000 0.2347573 3.76422822
sex:treat -0.01219535 0.2347573 -0.05194875
(Dispersion Parameter for Binomial family taken to be 1 )
Null Deviance: 146.6074 on 7 degrees of freedom
Residual Deviance: 128.3747 on 4 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
(Intercept) sex treat
sex -0.1519319
treat 0.3812090 -0.1406812
sex:treat -0.1406812 0.3812090 -0.1519319
> anova(fit7,test="Chisq")
Analysis of Deviance Table
Binomial model
Response: rem
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(Chi)
NULL 7 146.6074
sex 1 1.43362 6 145.1738 0.2311748
treat 1 16.79637 5 128.3774 0.0000416
sex:treat 1 0.00270 4 128.3747 0.9585671
> fit8<-update(fit7,.~.-sex:treat)
> summary(fit8)
Call: glm(formula = rem ~ sex + treat, family = binomial(link = logit), data =
remission.group, weights = count)
Deviance Residuals:
1 2 3 4 5 6 7 8
4.275416 -4.70127 2.40838 -3.702135 4.631506 -4.47269 3.000915 -4.246047
Coefficients:
Value Std. Error t value
(Intercept) 1.0179706 0.2323177 4.381805
sex -0.2770530 0.2169430 -1.277078
treat 0.8818644 0.2319361 3.802188
(Dispersion Parameter for Binomial family taken to be 1 )
Null Deviance: 146.6074 on 7 degrees of freedom
Residual Deviance: 128.3774 on 5 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
(Intercept) sex
sex -0.1027515
treat 0.3665306 -0.0854362
> anova(fit8,test="Chisq")
Analysis of Deviance Table
Binomial model
Response: rem
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(Chi)
NULL 7 146.6074
sex 1 1.43362 6 145.1738 0.2311748
treat 1 16.79637 5 128.3774 0.0000416
>