# For Ex. on AZT use and AIDS w/ data Table 5.4 Y. L. Tseng > n.table<-array(c(14, 32, 93, 81, 11, 12, 52, 43), dim = c(2, 2, 2), + dimnames = list(AZT=c("Yes", "No"), + Symptoms=c("Yes", "No"), Race=c("White", "Black"))) > n.table , , Race = White Symptoms AZT Yes No Yes 14 93 No 32 81 , , Race = Black Symptoms AZT Yes No Yes 11 52 No 12 43 > > > #Convert data - examine the data sets at each step > azt.aids.symp<-as.data.frame(as.table(n.table)) > azt.aids.table<-xtabs(Freq ~ AZT + Race, data = azt.aids.symp) > azt.aids<-as.data.frame(azt.aids.table) > symp.yes<-azt.aids.symp[azt.aids.symp$Symptoms == "Yes",] > azt.aids2<-data.frame(azt.aids, Symptoms = symp.yes$Freq) > azt.aids2 AZT Race Freq Symptoms 1 Yes White 107 14 2 No White 113 32 3 Yes Black 63 11 4 No Black 55 12 > > > mod.fit<-glm(formula = Symptoms/Freq ~ AZT + Race, data = azt.aids2, weight = Freq, family = binomial(link = logit), + na.action = na.exclude, control = list(epsilon = 0.000001, maxit = 50, trace = T)) Deviance = 1.389124 Iterations - 1 Deviance = 1.383530 Iterations - 2 Deviance = 1.38353 Iterations - 3 > summary(mod.fit) Call: glm(formula = Symptoms/Freq ~ AZT + Race, family = binomial(link = logit), data = azt.aids2, weights = Freq, na.action = na.exclude, control = list(epsilon = 1e-06, maxit = 50, trace = T)) Deviance Residuals: 1 2 3 4 -0.5547 0.4253 0.7035 -0.6326 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.73755 0.24038 -7.228 4.89e-13 *** AZTNo 0.71946 0.27897 2.579 0.00991 ** RaceBlack -0.05548 0.28861 -0.192 0.84755 --- Signif. codes: 0 ．***・ 0.001 ．**・ 0.01 ．*・ 0.05 ．.・ 0.1 ． ・ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 8.3499 on 3 degrees of freedom Residual deviance: 1.3835 on 1 degrees of freedom AIC: 24.86 Number of Fisher Scoring iterations: 3 > data.frame(azt.aids2, obs.prop = round(azt.aids2$Symptoms/azt.aids2$Freq,4), + pi.hat = round(predict(mod.fit, azt.aids2, type = "response"),4)) AZT Race Freq Symptoms obs.prop pi.hat 1 Yes White 107 14 0.1308 0.1496 2 No White 113 32 0.2832 0.2654 3 Yes Black 63 11 0.1746 0.1427 4 No Black 55 12 0.2182 0.2547 > anova(mod.fit) Deviance = 1.42384 Iterations - 1 Deviance = 1.420614 Iterations - 2 Deviance = 1.420614 Iterations - 3 Analysis of Deviance Table Model: binomial, link: logit Response: Symptoms/Freq Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 3 8.3499 AZT 1 6.9293 2 1.4206 Race 1 0.0371 1 1.3835 > #Check what the X matrix looks like > model.matrix(mod.fit) (Intercept) AZTNo RaceBlack 1 1 0 0 2 1 1 0 3 1 0 1 4 1 1 1 attr(,"assign") [1] 0 1 2 attr(,"contrasts") attr(,"contrasts")$AZT [1] "contr.treatment" attr(,"contrasts")$Race [1] "contr.treatment" > > #Matches the estimates on p. 146 after change reference levels > azt.aids3<-data.frame(azt.aids2, AZT.new = relevel(azt.aids2$AZT, ref = "No"), Race.new = relevel(azt.aids2$Race, ref = "Black")) > mod.fit.tab4.5<-glm(formula = Symptoms/Freq ~ AZT.new + Race.new, data = azt.aids3, weight = Freq, family = binomial(link = logit), + na.action = na.exclude, control = list(epsilon = 0.000001, maxit = 50, trace = T)) Deviance = 1.389124 Iterations - 1 Deviance = 1.383530 Iterations - 2 Deviance = 1.38353 Iterations - 3 > summary(mod.fit.tab4.5) Call: glm(formula = Symptoms/Freq ~ AZT.new + Race.new, family = binomial(link = logit), data = azt.aids3, weights = Freq, na.action = na.exclude, control = list(epsilon = 1e-06, maxit = 50, trace = T)) Deviance Residuals: 1 2 3 4 -0.5547 0.4253 0.7035 -0.6326 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.07357 0.26294 -4.083 4.45e-05 *** AZT.newYes -0.71946 0.27897 -2.579 0.00991 ** Race.newWhite 0.05548 0.28861 0.192 0.84755 --- Signif. codes: 0 ．***・ 0.001 ．**・ 0.01 ．*・ 0.05 ．.・ 0.1 ． ・ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 8.3499 on 3 degrees of freedom Residual deviance: 1.3835 on 1 degrees of freedom AIC: 24.86 Number of Fisher Scoring iterations: 3 > anova(mod.fit.tab4.5) Deviance = 1.42384 Iterations - 1 Deviance = 1.420614 Iterations - 2 Deviance = 1.420614 Iterations - 3 Analysis of Deviance Table Model: binomial, link: logit Response: Symptoms/Freq Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 3 8.3499 AZT.new 1 6.9293 2 1.4206 Race.new 1 0.0371 1 1.3835 > model.matrix(mod.fit.tab4.5) (Intercept) AZT.newYes Race.newWhite 1 1 1 1 2 1 0 1 3 1 1 0 4 1 0 0 attr(,"assign") [1] 0 1 2 attr(,"contrasts") attr(,"contrasts")$AZT.new [1] "contr.treatment" attr(,"contrasts")$Race.new [1] "contr.treatment" > ################################################################# > # Change to sum to 0 option for the levels of each factor > > #Sum to zero option > options(contrasts=c("contr.sum", "contr.poly")) > > #match Agresti (1996) Table 5.6 on p. 121 (second column) > mod.fit<-glm(formula = Symptoms/Freq ~ AZT + Race, data = azt.aids2, weight = Freq, family = binomial(link = logit), + na.action = na.exclude, control = list(epsilon = 0.000001, maxit = 50, trace = T)) Deviance = 1.389124 Iterations - 1 Deviance = 1.383530 Iterations - 2 Deviance = 1.38353 Iterations - 3 > summary(mod.fit) Call: glm(formula = Symptoms/Freq ~ AZT + Race, family = binomial(link = logit), data = azt.aids2, weights = Freq, na.action = na.exclude, control = list(epsilon = 1e-06, maxit = 50, trace = T)) Deviance Residuals: 1 2 3 4 -0.5547 0.4253 0.7035 -0.6326 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.40556 0.14668 -9.582 < 2e-16 *** AZT1 -0.35973 0.13949 -2.579 0.00991 ** Race1 0.02774 0.14430 0.192 0.84755 --- Signif. codes: 0 ．***・ 0.001 ．**・ 0.01 ．*・ 0.05 ．.・ 0.1 ． ・ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 8.3499 on 3 degrees of freedom Residual deviance: 1.3835 on 1 degrees of freedom AIC: 24.86 Number of Fisher Scoring iterations: 3 > data.frame(azt.aids2, obs.prop = round(azt.aids2$Symptoms/azt.aids2$Freq,4), + pi.hat = round(predict(mod.fit, azt.aids2, type = "response"),4)) AZT Race Freq Symptoms obs.prop pi.hat 1 Yes White 107 14 0.1308 0.1496 2 No White 113 32 0.2832 0.2654 3 Yes Black 63 11 0.1746 0.1427 4 No Black 55 12 0.2182 0.2547 > > #Check what the X matrix looks like > model.matrix(mod.fit) (Intercept) AZT1 Race1 1 1 1 1 2 1 -1 1 3 1 1 -1 4 1 -1 -1 attr(,"assign") [1] 0 1 2 attr(,"contrasts") attr(,"contrasts")$AZT [1] "contr.sum" attr(,"contrasts")$Race [1] "contr.sum" > > > #First level set to 0 - change back! > options(contrasts=c("contr.treatment", "contr.poly")) > contrasts(N) 1 2 4 0 0 0 0 1 1 0 0 2 0 1 0 4 0 0 1 >