> # Handout for GLM Y. L. Tseng > # On loglinear model > # Example: Alcohol(A), Cigarette(C) and Marijuana(M) use for H.S.seniors, p. 209 > > n.table<-array(c(911, 538, 44, 456, 3, 43, 2, 279), dim = c(2, 2, 2), + dimnames = list(M.use=c("Yes", "No"), + C.use=c("Yes", "No"), A.use=c("Yes","No"))) > n.table , , A.use = Yes C.use M.use Yes No Yes 911 44 No 538 456 , , A.use = No C.use M.use Yes No Yes 3 2 No 43 279 > > > #Convert data > ACMdata<-as.data.frame(as.table(n.table)) > ACMdata # Table 7.3, p. 209 M.use C.use A.use Freq 1 Yes Yes Yes 911 2 No Yes Yes 538 3 Yes No Yes 44 4 No No Yes 456 5 Yes Yes No 3 6 No Yes No 43 7 Yes No No 2 8 No No No 279 > > #Convert data.... shorter names for easier use... > > data7.3<-data.frame(ACMdata, C=relevel(ACMdata$C.use, ref="No"), A=relevel(ACMdata$A.use, ref="No"), M=relevel(ACMdata$M.use, ref="No")) > m2<-glm(formula = Freq ~ (A + C + M)^2, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 0.3815073 Iterations - 1 Deviance = 0.3739884 Iterations - 2 Deviance = 0.3739859 Iterations - 3 > summary(m2) # Table 7.6, p. 211 Call: glm(formula = Freq ~ (A + C + M)^2, family = poisson(link = log), data = data7.3, na.action = na.exclude, control = list(epsilon = 1e-04, maxit = 50, trace = T)) Deviance Residuals: 1 2 3 4 5 6 7 0.02044 -0.02658 -0.09256 0.02890 -0.33428 0.09452 0.49134 8 -0.03690 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.63342 0.05970 94.361 < 2e-16 *** AYes 0.48772 0.07577 6.437 1.22e-10 *** CYes -1.88667 0.16270 -11.596 < 2e-16 *** MYes -5.30904 0.47505 -11.176 < 2e-16 *** AYes:CYes 2.05453 0.17406 11.803 < 2e-16 *** AYes:MYes 2.98601 0.46453 6.428 1.29e-10 *** CYes:MYes 2.84789 0.16384 17.382 < 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: 2851.46098 on 7 degrees of freedom Residual deviance: 0.37399 on 1 degrees of freedom AIC: 63.417 Number of Fisher Scoring iterations: 3 > > > mu.hatm2<-predict(object = m2, type="response") > pearson<-residuals(object = m2, type="pearson") # may needed for model checking > h<-lm.influence(model = m2)$h > standard.pearson<-pearson/sqrt(1-h) > > > m1<-glm(formula = Freq ~ A + C + M+ A*M+ C*M, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 220.1212 Iterations - 1 Deviance = 188.2610 Iterations - 2 Deviance = 187.7548 Iterations - 3 Deviance = 187.7543 Iterations - 4 > #summary(m1) > mu.hatm1<-predict(object = m1, type="response") > pearson1<-residuals(object = m1, type="pearson") > h1<-lm.influence(model = m1)$h > standard.resid1<-pearson1/sqrt(1-h1) > Table7.8<-data.frame(ACMdata, mu.hatm1 = round(mu.hatm1,1), standard.resid1=round(standard.resid1,2),mu.hatm2 = round(mu.hatm2,1),standard.pearson=round(standard.pearson,2)) > Table7.8 M.use C.use A.use Freq mu.hatm1 standard.resid1 mu.hatm2 standard.pearson 1 Yes Yes Yes 911 909.2 3.68 910.4 0.63 2 No Yes Yes 538 438.8 12.80 538.6 -0.63 3 Yes No Yes 44 45.8 -3.68 44.6 -0.63 4 No No Yes 456 555.2 -12.80 455.4 0.63 5 Yes Yes No 3 4.8 -3.70 3.6 -0.63 6 No Yes No 43 142.2 -12.81 42.4 0.63 7 Yes No No 2 0.2 3.70 1.4 0.63 8 No No No 279 179.8 12.81 279.6 -0.63 > > > m3<-glm(formula = Freq ~ A + C + M, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 2336.958 Iterations - 1 Deviance = 1422.337 Iterations - 2 Deviance = 1291.151 Iterations - 3 Deviance = 1286.034 Iterations - 4 Deviance = 1286.02 Iterations - 5 > m4<-glm(formula = Freq ~ A + C + M+C*M, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 869.765 Iterations - 1 Deviance = 561.3993 Iterations - 2 Deviance = 534.5907 Iterations - 3 Deviance = 534.2118 Iterations - 4 Deviance = 534.2117 Iterations - 5 > > m5<-glm(formula = Freq ~ A + C + M+A*M, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 1227.03 Iterations - 1 Deviance = 950.7822 Iterations - 2 Deviance = 939.5908 Iterations - 3 Deviance = 939.5626 Iterations - 4 > > m6<-glm(formula = Freq ~ A + C + M+A*C, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 1322.133 Iterations - 1 Deviance = 875.0708 Iterations - 2 Deviance = 844.0985 Iterations - 3 Deviance = 843.8267 Iterations - 4 Deviance = 843.8266 Iterations - 5 > > m7<-glm(formula = Freq ~ A + C + M+ A*C+A*M, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 676.8397 Iterations - 1 Deviance = 504.4351 Iterations - 2 Deviance = 497.3883 Iterations - 3 Deviance = 497.3693 Iterations - 4 > > m8<-glm(formula = Freq ~ A + C + M+ A*C+C*M, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 126.4209 Iterations - 1 Deviance = 93.95147 Iterations - 2 Deviance = 92.03281 Iterations - 3 Deviance = 92.01836 Iterations - 4 Deviance = 92.01836 Iterations - 5 > m9<-glm(formula = Freq ~ A + C + M+ A*M+C*M, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 220.1212 Iterations - 1 Deviance = 188.2610 Iterations - 2 Deviance = 187.7548 Iterations - 3 Deviance = 187.7543 Iterations - 4 > ms<-glm(formula = Freq ~ (A + C + M)^3, + data = data7.3, family = poisson(link = log), + na.action = na.exclude, control = list(epsilon = 0.0001, maxit = 50, trace = T)) Deviance = 3.593156e-06 Iterations - 1 Deviance = 8.987437e-13 Iterations - 2 > G2<-c(m3$deviance,m4$deviance,m5$deviance,m6$deviance,m7$deviance,m8$deviance,m9$deviance,m2$deviance,ms$deviance) > DF<-c(m3$df.residual,m4$df.residual,m5$df.residual,m6$df.residual,m7$df.residual,m8$df.residual,m9$df.residual,m2$df.residual,ms$df.residual) > Pvalue<-1-pchisq(G2,DF) # p-values for G-O-F tests > cbind(G2,DF,Pvalue) # Table 7.7, p. 212 G2 DF Pvalue [1,] 1.286020e+03 4 0.0000000 [2,] 5.342117e+02 3 0.0000000 [3,] 9.395626e+02 3 0.0000000 [4,] 8.438266e+02 3 0.0000000 [5,] 4.973693e+02 2 0.0000000 [6,] 9.201836e+01 2 0.0000000 [7,] 1.877543e+02 2 0.0000000 [8,] 3.739859e-01 1 0.5408396 [9,] 8.987437e-13 0 0.0000000 > > # Estimated odds ratios from each model, c.f. Table 7.5, p. 210 > exp(m4$coef[5]) CYes:MYes 25.13620 > exp(m5$coef[5]) AYes:MYes 61.8741 > exp(m6$coef[5]) AYes:CYes 17.703 > exp(m7$coef[5:6]) AYes:CYes AYes:MYes 17.70268 61.87151 > exp(m8$coef[5:6]) AYes:CYes CYes:MYes 17.70300 25.13620 > exp(m9$coef[5:6]) AYes:MYes CYes:MYes 61.87086 25.13620 > exp(m2$coef[5:7]) AYes:CYes AYes:MYes CYes:MYes 7.803201 19.806580 17.251329 > cbind(exp(ms$coef[5:7]),exp(ms$coef[5:7]+ms$coef[8])) [,1] [,2] AYes:CYes 7.655141 13.80304 AYes:MYes 13.460517 24.27075 CYes:MYes 9.732553 17.54883 > > # For Table 7.4, p. 210 > mu.hatm3<-round(predict(object = m3, type="response"),1) > mu.hatm4<-round(predict(object = m4, type="response"),1) > mu.hatm5<-round(predict(object = m5, type="response"),1) > mu.hatm6<-round(predict(object = m6, type="response"),1) > mu.hatm7<-round(predict(object = m7, type="response"),1) > mu.hatm8<-round(predict(object = m8, type="response"),1) > mu.hatm9<-round(predict(object = m9, type="response"),1) > mu.hatm2<-round(predict(object = m2, type="response"),1) > mu.hatms<-round(predict(object = ms, type="response"),1) > > table7.4<-data.frame(ACMdata,muhat.A.C.M=mu.hatm3,muhat.M.AC=mu.hatm6, muhat.AM.CM=mu.hatm9,muhat2way=mu.hatm2,muhats=mu.hatms, + muhat.A.CM=mu.hatm4,muhat.C.AM=mu.hatm5,muhat.AC.AM=mu.hatm7,muhat.AC.CM=mu.hatm8) > table7.4 M.use C.use A.use Freq muhat.A.C.M muhat.M.AC muhat.AM.CM muhat2way muhats 1 Yes Yes Yes 911 540.0 611.2 909.2 910.4 911 2 No Yes Yes 538 740.2 837.8 438.8 538.6 538 3 Yes No Yes 44 282.1 210.9 45.8 44.6 44 4 No No Yes 456 386.7 289.1 555.2 455.4 456 5 Yes Yes No 3 90.6 19.4 4.8 3.6 3 6 No Yes No 43 124.2 26.6 142.2 42.4 43 7 Yes No No 2 47.3 118.5 0.2 1.4 2 8 No No No 279 64.9 162.5 179.8 279.6 279 muhat.A.CM muhat.C.AM muhat.AC.AM muhat.AC.CM 1 782.7 627.3 710.0 885.9 2 497.5 652.9 739.0 563.1 3 39.4 327.7 245.0 29.4 4 629.4 341.1 255.0 470.6 5 131.3 3.3 0.7 28.1 6 83.5 211.5 45.3 17.9 7 6.6 1.7 4.3 16.6 8 105.6 110.5 276.7 264.4 >