> #------------------------------------------------------------------- > # 6. LOG LINEAR MODELS FOR POISSON DATA > #------------------------------------------------------------------- > > cars<-c(294,276,238,192) > glm<-glm(cars~1,family=poisson) > summary(glm) Call: glm(formula = cars ~ 1, family = poisson) Deviance Residuals: 1 2 3 4 2.7066 1.6171 -0.7651 -3.8259 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.52146 0.03162 174.6 <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: 25.163 on 3 degrees of freedom Residual deviance: 25.163 on 3 degrees of freedom AIC: 56.551 Number of Fisher Scoring iterations: 3 > sum((cars-fitted(glm))^2/fitted(glm)) [1] 24.48 > cbind(cars,fitted(glm),resid(glm)) cars 1 294 250 2.706637 2 276 250 1.617050 3 238 250 -0.765143 4 192 250 -3.825863 > lane<-1:4 > par(mfrow=c(2,2)) > plot(lane,cars,main="Figure 6.1 Cars vs. lane") > glm<-glm(cars~lane,family=poisson) > summary(glm) Call: glm(formula = cars ~ lane, family = poisson) Deviance Residuals: 1 2 3 4 -0.5771 0.6879 0.4897 -0.6215 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.85539 0.07353 79.634 < 2e-16 *** lane -0.13834 0.02851 -4.852 1.22e-06 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 25.1634 on 3 degrees of freedom Residual deviance: 1.4323 on 2 degrees of freedom AIC: 34.82 Number of Fisher Scoring iterations: 3 > sum((cars-fitted(glm))^2/fitted(glm)) [1] 1.432276 > cbind(cars,fitted(glm),resid(glm)) cars 1 294 304.0063 -0.5770848 2 276 264.7282 0.6879473 3 238 230.5249 0.4897084 4 192 200.7407 -0.6214794 > lines(lane,fitted(glm)) > lane<-factor(lane) > glm<-glm(cars~lane,family=poisson) > summary(glm) Call: glm(formula = cars ~ lane, family = poisson) Deviance Residuals: [1] 0 0 0 0 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.68358 0.05832 97.453 < 2e-16 *** lane2 -0.06318 0.08381 -0.754 0.4510 lane3 -0.21131 0.08720 -2.423 0.0154 * lane4 -0.42608 0.09279 -4.592 4.39e-06 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 2.5163e+01 on 3 degrees of freedom Residual deviance: -1.1990e-14 on 0 degrees of freedom AIC: 37.388 Number of Fisher Scoring iterations: 2 > sum((cars-fitted(glm))^2/fitted(glm)) [1] 2.629385e-26 > cbind(cars,fitted(glm),resid(glm)) cars 1 294 294 0 2 276 276 0 3 238 238 0 4 192 192 0 > > x<-c(1997,906,904,32) > p0<-c(0.5,0.25,0.25,0) > p1<-c(0.25,-0.25,-0.25,0.25) > ex0<-p0*sum(x) > theta<-p1*sum(x) > cbind(x,ex0,theta) x ex0 theta [1,] 1997 1919.50 959.75 [2,] 906 959.75 -959.75 [3,] 904 959.75 -959.75 [4,] 32 0.00 959.75 > glm<-glm(x~-1+theta,family=poisson(link=identity),offset=ex0) > summary(glm) Call: glm(formula = x ~ -1 + theta, family = poisson(link = identity), offset = ex0) Deviance Residuals: 1 2 3 4 0.9743 -0.6425 -0.7087 -0.3929 Coefficients: Estimate Std. Error z value Pr(>|z|) theta 0.035712 0.005838 6.117 9.52e-10 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: Inf on 4 degrees of freedom Residual deviance: 2.0187 on 3 degrees of freedom AIC: 36.057 Number of Fisher Scoring iterations: 3 > sum((x-fitted(glm))^2/fitted(glm)) [1] 2.015438 > cbind(x,fitted(glm),resid(glm)) x 1 1997 1953.77443 0.9743473 2 906 925.47557 -0.6424535 3 904 925.47557 -0.7086878 4 32 34.27443 -0.3929171 > > defect<-c(15,21,45,13, + 26,31,34, 5, + 33,17,49,20) > type<-gl(4,1,12) #type<-factor(rep(1:4,3)) > shift<-gl(3,4,12) #shift<-factor(rep(1:3,rep(4,3))) > cbind(type,shift,defect) type shift defect [1,] 1 1 15 [2,] 2 1 21 [3,] 3 1 45 [4,] 4 1 13 [5,] 1 2 26 [6,] 2 2 31 [7,] 3 2 34 [8,] 4 2 5 [9,] 1 3 33 [10,] 2 3 17 [11,] 3 3 49 [12,] 4 3 20 > glm<-glm(defect~type+shift,family=poisson) > summary(glm) Call: glm(formula = defect ~ type + shift, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -2.2406 -1.1251 0.2086 0.8537 1.9349 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.11402 0.14459 21.536 < 2e-16 *** type2 -0.06996 0.16727 -0.418 0.675785 type3 0.54797 0.14601 3.753 0.000175 *** type4 -0.66647 0.19927 -3.345 0.000824 *** shift2 0.02105 0.14502 0.145 0.884563 shift3 0.23583 0.13795 1.710 0.087346 . --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 77.409 on 11 degrees of freedom Residual deviance: 20.336 on 6 degrees of freedom AIC: 91.672 Number of Fisher Scoring iterations: 3 > sum((defect-fitted(glm))^2/fitted(glm)) [1] 19.17790 > cbind(defect,fitted(glm),resid(glm)) defect 1 15 22.51132 -1.686297122 2 21 20.99031 0.002115905 3 45 38.93850 0.947687387 4 13 11.55998 0.415169444 5 26 22.99032 0.614696133 6 31 21.43694 1.934921720 7 34 39.76703 -0.938072654 8 5 11.80595 -2.240590830 9 33 28.49838 0.822400115 10 17 26.57284 -1.989643815 11 49 49.29449 -0.041985522 12 20 14.63444 1.327803162 > > wife<-factor(c(1,1,1,2,2,2,2,3,3,3,4,4,5,5,5,5,5)) > husb<-factor(c(2,3,5,1,3,4,5,2,4,5,1,5,1,2,3,4,5)) > marr<-c(5,17,6,5,0,16,2,2,10,11,10,9,6,20,8,0,1) > cbind(wife,husb,marr) wife husb marr [1,] 1 2 5 [2,] 1 3 17 [3,] 1 5 6 [4,] 2 1 5 [5,] 2 3 0 [6,] 2 4 16 [7,] 2 5 2 [8,] 3 2 2 [9,] 3 4 10 [10,] 3 5 11 [11,] 4 1 10 [12,] 4 5 9 [13,] 5 1 6 [14,] 5 2 20 [15,] 5 3 8 [16,] 5 4 0 [17,] 5 5 1 > glm<-glm(marr~wife+husb,family=poisson) > summary(glm) Call: glm(formula = marr ~ wife + husb, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -4.22267 -1.97023 0.06253 0.14888 3.62695 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.0672 0.3471 5.955 2.59e-09 *** wife2 -0.4859 0.3075 -1.580 0.114 wife3 -0.2521 0.3065 -0.823 0.411 wife4 0.2742 0.3592 0.764 0.445 wife5 -0.3194 0.2678 -1.193 0.233 husb2 0.3108 0.3495 0.889 0.374 husb3 0.3008 0.3403 0.884 0.377 husb4 0.4400 0.3364 1.308 0.191 husb5 -0.1893 0.3008 -0.629 0.529 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 86.479 on 16 degrees of freedom Residual deviance: 76.251 on 8 degrees of freedom AIC: 150.45 Number of Fisher Scoring iterations: 4 > sum((marr-fitted(glm))^2/fitted(glm)) [1] 66.53147 > cbind(marr,fitted(glm),resid(glm)) marr 1 5 10.784044 -1.97022614 2 17 10.676186 1.78019485 3 6 6.539793 -0.21408744 4 5 4.861487 0.06252633 5 0 6.567455 -3.62421177 6 16 7.548477 2.67146522 7 2 4.022953 -1.11822641 8 2 8.380978 -2.65154407 9 10 9.536552 0.14888255 10 11 5.082496 2.26958900 11 10 10.396643 -0.12380857 12 9 8.603378 0.13420087 13 6 5.741889 0.10692360 14 20 7.835168 3.62694556 15 8 7.756803 0.08687010 16 0 8.915485 -4.22267338 17 1 4.751499 -2.09429647 > > > y<-c(2,5,4,28,35,39,38,29,22,13,9,5)*10 > opin<-gl(4,3,12) > sond<-gl(3,1,12) > glm<-glm(y~opin+sond,family=poisson) > summary(glm) Call: glm(formula = y ~ opin + sond, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -4.4288 -3.2857 -0.0338 2.2962 4.2611 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.66121 0.09930 36.869 < 2e-16 *** opin2 2.22707 0.10022 22.222 < 2e-16 *** opin3 2.09074 0.10093 20.714 < 2e-16 *** opin4 0.89794 0.11297 7.948 1.89e-15 *** sond2 -0.03774 0.05016 -0.752 0.45177 sond3 -0.14595 0.05159 -2.829 0.00467 ** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 1307.34 on 11 degrees of freedom Residual deviance: 103.99 on 6 degrees of freedom AIC: 196.41 Number of Fisher Scoring iterations: 3 > sum((y-fitted(glm))^2/fitted(glm)) [1] 101.6752 > cbind(y,fitted(glm),resid(glm)) y 1 20 38.90847 -3.3463050 2 50 37.46741 1.9467847 3 40 33.62461 1.0672012 4 280 360.78609 -4.4288402 5 350 347.42356 0.1380561 6 390 311.79043 4.2610949 7 380 314.80354 3.5575516 8 290 303.14408 -0.7604845 9 220 272.05243 -3.2654657 10 130 95.50227 3.3443225 11 90 91.96513 -0.2056537 12 50 82.53282 -3.8663082 > glm(y~opin+sond+opin:sond,family=poisson) Call: glm(formula = y ~ opin + sond + opin:sond, family = poisson) Coefficients: (Intercept) opin2 opin3 opin4 sond2 sond3 2.9957 2.6391 2.9444 1.8718 0.9163 0.6931 opin2.sond2 opin3.sond2 opin4.sond2 opin2.sond3 opin3.sond3 opin4.sond3 -0.6931 -1.1866 -1.2840 -0.3618 -1.2397 -1.6487 Degrees of Freedom: 11 Total (i.e. Null); 0 Residual Null Deviance: 1307 Residual Deviance: -7.106e-15 AIC: 104.4 > x<-codes(opin)*codes(sond) > glm<-glm(y~opin+sond+x,family=poisson) > summary(glm) Call: glm(formula = y ~ opin + sond + x, family = poisson) Deviance Residuals: 1 2 3 4 5 6 7 8 -0.3370 2.2154 -1.8053 -0.8959 -0.1888 0.9680 1.3196 -0.8725 9 10 11 12 -0.6665 -0.7165 0.4917 0.5453 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.40132 0.10867 31.300 < 2e-16 *** opin2 2.94888 0.13214 22.317 < 2e-16 *** opin3 3.46390 0.18687 18.536 < 2e-16 *** opin4 2.85310 0.24194 11.793 < 2e-16 *** sond2 0.84250 0.11021 7.644 2.10e-14 *** sond3 1.55356 0.19062 8.150 3.64e-16 *** x -0.33117 0.03586 -9.235 < 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: 1307.342 on 11 degrees of freedom Residual deviance: 14.055 on 5 degrees of freedom AIC: 108.47 Number of Fisher Scoring iterations: 3 > sum((y-fitted(glm))^2/fitted(glm)) [1] 14.41183 > cbind(y,fitted(glm),resid(glm)) y 1 20 21.54505 -0.3369678 2 50 35.92631 2.2153844 3 40 52.52866 -1.8052601 4 280 295.26061 -0.8959350 5 350 353.54471 -0.1888370 6 390 371.19467 0.9679955 7 380 354.85291 1.3196267 8 290 305.11313 -0.8725090 9 220 230.03397 -0.6664705 10 130 138.34144 -0.7165056 11 90 85.41586 0.4916681 12 50 46.24271 0.5452869 > > > > inf<-c(3,176,4,293,17,197,2,23) > clin<-gl(2,4,8) > care<-gl(2,2,8) > surv<-gl(2,1,8) > cbind(inf,clin,care,surv) inf clin care surv [1,] 3 1 1 1 [2,] 176 1 1 2 [3,] 4 1 2 1 [4,] 293 1 2 2 [5,] 17 2 1 1 [6,] 197 2 1 2 [7,] 2 2 2 1 [8,] 23 2 2 2 > glm(inf~care+surv,family=poisson) Call: glm(formula = inf ~ care + surv, family = poisson) Coefficients: (Intercept) care2 surv2 1.9665 -0.1993 3.2771 Degrees of Freedom: 7 Total (i.e. Null); 5 Residual Null Deviance: 1066 Residual Deviance: 291.5 AIC: 337.7 > glm(inf~care+surv+care:surv,family=poisson) Call: glm(formula = inf ~ care + surv + care:surv, family = poisson) Coefficients: (Intercept) care2 surv2 care2.surv2 2.303 -1.204 2.926 1.038 Degrees of Freedom: 7 Total (i.e. Null); 4 Residual Null Deviance: 1066 Residual Deviance: 285.9 AIC: 334.1 > glm(inf~care+surv+clin,family=poisson) Call: glm(formula = inf ~ care + surv + clin, family = poisson) Coefficients: (Intercept) care2 surv2 clin2 2.2528 -0.1993 3.2771 -0.6889 Degrees of Freedom: 7 Total (i.e. Null); 4 Residual Null Deviance: 1066 Residual Deviance: 211.5 AIC: 259.7 > glm(inf~care+surv+clin+care:surv,family=poisson) Call: glm(formula = inf ~ care + surv + clin + care:surv, family = poisson) Coefficients: (Intercept) care2 surv2 clin2 care2.surv2 2.589 -1.204 2.926 -0.689 1.038 Degrees of Freedom: 7 Total (i.e. Null); 3 Residual Null Deviance: 1066 Residual Deviance: 205.9 AIC: 256.1 > glm(inf~care+surv+clin+care:clin,family=poisson) Call: glm(formula = inf ~ care + surv + clin + care:clin, family = poisson) Coefficients: (Intercept) care2 surv2 clin2 care2.clin2 1.8732 0.5063 3.2771 0.1786 -2.6534 Degrees of Freedom: 7 Total (i.e. Null); 3 Residual Null Deviance: 1066 Residual Deviance: 17.83 AIC: 68.01 > glm(inf~care+surv+clin+care:clin+care:surv,family=poisson) Call: glm(formula = inf ~ care + surv + clin + care:clin + care:surv, family = poisson) Coefficients: (Intercept) care2 surv2 clin2 care2.surv2 care2.clin2 2.2093 -0.4984 2.9258 0.1786 1.0381 -2.6534 Degrees of Freedom: 7 Total (i.e. Null); 2 Residual Null Deviance: 1066 Residual Deviance: 12.22 AIC: 64.4 > glm(inf~care+surv+clin+care:clin+surv:clin,family=poisson) Call: glm(formula = inf ~ care + surv + clin + care:clin + surv:clin, family = poisson) Coefficients: (Intercept) care2 surv2 clin2 care2.clin2 surv2.clin2 0.9679 0.5063 4.2047 1.8661 -2.6534 -1.7555 Degrees of Freedom: 7 Total (i.e. Null); 2 Residual Null Deviance: 1066 Residual Deviance: 0.08229 AIC: 52.26 > glm(inf~care+surv+clin+care:clin+surv:clin+care:surv,family=poisson) Call: glm(formula = inf ~ care + surv + clin + care:clin + surv:clin + care:surv, family = poisson) Coefficients: (Intercept) care2 surv2 clin2 care2.surv2 care2.clin2 1.0343 0.3976 4.1372 1.8098 0.1104 -2.6467 surv2.clin2 -1.6991 Degrees of Freedom: 7 Total (i.e. Null); 1 Residual Null Deviance: 1066 Residual Deviance: 0.04326 AIC: 54.23 > glm<-glm(inf~care+surv+clin+care:clin+surv:clin,family=poisson) > summary(glm) Call: glm(formula = inf ~ care + surv + clin + care:clin + surv:clin, family = poisson) Deviance Residuals: 1 2 3 4 5 6 0.2216100 -0.0276932 -0.1784756 0.0214872 -0.0030436 0.0008943 7 8 0.0088945 -0.0026169 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.96788 0.38254 2.530 0.0114 * care2 0.50635 0.09462 5.351 8.74e-08 *** surv2 4.20469 0.38077 11.043 < 2e-16 *** clin2 1.86607 0.44661 4.178 2.94e-05 *** care2.clin2 -2.65345 0.23157 -11.458 < 2e-16 *** surv2.clin2 -1.75550 0.44963 -3.904 9.45e-05 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 1066.42776 on 7 degrees of freedom Residual deviance: 0.08229 on 2 degrees of freedom AIC: 52.265 Number of Fisher Scoring iterations: 3 > cbind(inf,fitted(glm),resid(glm)) inf 1 3 2.632353 0.2216100408 2 176 176.367647 -0.0276931670 3 4 4.367647 -0.1784755943 4 293 292.632353 0.0214871605 5 17 17.012552 -0.0030436313 6 197 196.987448 0.0008943334 7 2 1.987448 0.0088944543 8 23 23.012552 -0.0026168598 > >