Tue Nov 23 09:48:05 2004 ======================== Some notes... Install exactLoglinTest from CRAN... library(exactLoglinTest) data(czech.dat) czech.dat[1:4,] y A B C D E F 1 44 no no no small small neg 2 40 yes no no small small neg 3 112 no yes no small small neg 4 67 yes yes no small small neg sum(czeck.dat$y) [1] 1841 # y = counts # # A = smoke? 'no' 'yes' # # B = mental work? 'no' 'yes' # # C = physical work? 'no' 'yes' # # D = blood pressure? 'small' 'large' # # E = lipoprotein level? 'small' 'large' # # F = family recall of coronary heart disease? 'neg' 'pos' # # Edwards and Havranek (1985, Bmka) ##################################################################### mod1 <- glm(y ~., family=poisson, data=czech.dat) mod0 <- glm(y ~ 1,family=poisson, data=czech.dat) mod.full <- glm(y ~ A*B*C*D*E*F, family=poisson, data=czech.dat) mod.2way <- glm(y ~ .^2, family=poisson,data=czech.dat) ##################################################################### library(MASS) result.full <- stepAIC(mod.full,direction="backward") result.2way <- stepAIC(mod.2way,direction="backward") formula(result.full) y ~ A + B + C + D + E + F + A:B + A:C + B:C + A:D + B:D + C:D + A:E + B:E + C:E + D:E + A:F + B:F + C:F + D:F + E:F + A:B:C + A:B:D + A:C:D + B:C:D + A:B:E + A:C:E + B:C:E + A:D:E + C:D:E + A:B:F + A:C:F + B:C:F + A:D:F + B:D:F + C:D:F + A:E:F + B:E:F + C:E:F + D:E:F + A:B:C:D + A:B:C:E + A:C:D:E + A:B:C:F + A:B:D:F + A:C:D:F + B:C:D:F + A:B:E:F + A:C:E:F + B:C:E:F + A:D:E:F + C:D:E:F + A:B:C:D:F + A:B:C:E:F + A:C:D:E:F formula(result.2way) y ~ A + B + C + D + E + F + A:C + A:D + A:E + B:C + B:E + B:F + C:E + D:E + E:F anova(result.2way,result.full,test="Chisq") Analysis of Deviance Table Model 1: y ~ A + B + C + D + E + F + A:C + A:D + A:E + B:C + B:E + B:F + C:E + D:E + E:F Model 2: y ~ A + B + C + D + E + F + A:B + A:C + B:C + A:D + B:D + C:D + A:E + B:E + C:E + D:E + A:F + B:F + C:F + D:F + E:F + A:B:C + A:B:D + A:C:D + B:C:D + A:B:E + A:C:E + B:C:E + A:D:E + C:D:E + A:B:F + A:C:F + B:C:F + A:D:F + B:D:F + C:D:F + A:E:F + B:E:F + C:E:F + D:E:F + A:B:C:D + A:B:C:E + A:C:D:E + A:B:C:F + A:B:D:F + A:C:D:F + B:C:D:F + A:B:E:F + A:C:E:F + B:C:E:F + A:D:E:F + C:D:E:F + A:B:C:D:F + A:B:C:E:F + A:C:D:E:F Resid. Df Resid. Dev Df Deviance 1 48 52.705 2 8 5.067 40 47.638 0.190 summary(result.2way) Call: glm(formula = y ~ A + B + C + D + E + F + A:C + A:D + A:E + B:C + B:E + B:F + C:E + D:E + E:F, family = poisson, data = czech.dat) Deviance Residuals: Min 1Q Median 3Q Max -2.56627 -0.59556 -0.06848 0.52592 1.97395 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.90128 0.08869 43.989 < 2e-16 *** Ayes -0.41447 0.08977 -4.617 3.89e-06 *** Byes 0.74165 0.09122 8.130 4.28e-16 *** Cyes 0.94610 0.09481 9.979 < 2e-16 *** Dlarge -0.28777 0.07466 -3.854 0.000116 *** Elarge -0.70690 0.12537 -5.638 1.72e-08 *** Fpos -2.01895 0.10839 -18.627 < 2e-16 *** Ayes:Cyes 0.54022 0.09534 5.666 1.46e-08 *** Ayes:Dlarge -0.35339 0.09574 -3.691 0.000223 *** Ayes:Elarge 0.48691 0.09723 5.008 5.50e-07 *** Byes:Cyes -2.78386 0.12228 -22.766 < 2e-16 *** Byes:Elarge 0.26058 0.11795 2.209 0.027152 * Byes:Fpos 0.27229 0.13492 2.018 0.043576 * Cyes:Elarge -0.29183 0.11825 -2.468 0.013593 * Dlarge:Elarge 0.37849 0.09635 3.928 8.55e-05 *** Elarge:Fpos 0.20655 0.13508 1.529 0.126230 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 2026.740 on 63 degrees of freedom Residual deviance: 52.705 on 48 degrees of freedom AIC: 374.36 Number of Fisher Scoring iterations: 4 result.2way.reduced <- update(result.2way, . ~ . - E:F) summary(result.2way.reduced) Call: glm(formula = y ~ A + B + C + D + E + F + A:C + A:D + A:E + B:C + B:E + B:F + C:E + D:E, family = poisson, data = czech.dat) Deviance Residuals: Min 1Q Median 3Q Max -2.43211 -0.62121 0.05023 0.50687 1.71741 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.89123 0.08854 43.948 < 2e-16 *** Ayes -0.41447 0.08977 -4.617 3.89e-06 *** Byes 0.73574 0.09146 8.045 8.66e-16 *** Cyes 0.94610 0.09481 9.979 < 2e-16 *** Dlarge -0.28777 0.07466 -3.854 0.000116 *** Elarge -0.68035 0.12410 -5.482 4.20e-08 *** Fpos -1.93627 0.09241 -20.954 < 2e-16 *** Ayes:Cyes 0.54022 0.09534 5.666 1.46e-08 *** Ayes:Dlarge -0.35339 0.09574 -3.691 0.000223 *** Ayes:Elarge 0.48691 0.09723 5.008 5.50e-07 *** Byes:Cyes -2.78386 0.12228 -22.766 < 2e-16 *** Byes:Elarge 0.26754 0.11782 2.271 0.023162 * Byes:Fpos 0.29251 0.13420 2.180 0.029283 * Cyes:Elarge -0.29183 0.11825 -2.468 0.013593 * Dlarge:Elarge 0.37849 0.09635 3.928 8.55e-05 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 2026.740 on 63 degrees of freedom Residual deviance: 55.034 on 49 degrees of freedom AIC: 374.69 Number of Fisher Scoring iterations: 4 par(mfrow=c(2,2)) plot(result.2way.reduced) ################################################################### decomposable.2way.reduced <- glm(y ~ A*D*E + A*B*C*E + B*F, family=poisson, data=czech.dat) plot(decomposable.2way.reduced) # What will happen when fit glm(F ~ ...,weight=y, # family=binomial,data=czech.dat)? ####################################################################