#------------------------------------------------------------------- # 6. LOG LINEAR MODELS FOR POISSON DATA #------------------------------------------------------------------- > cars <- c(294, 276, 238, 192) > glm <- glm(cars ~ 1, family = poisson) > summary(glm) Coefficients: Value Std. Error t value (Intercept) 5.521461 0.03162278 174.6039 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 25.16341 on 3 degrees of freedom Number of Fisher Scoring Iterations: 0 > 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) Coefficients: Value Std. Error t value (Intercept) 5.8553931 0.07352882 79.633986 lane -0.1383448 0.02851434 -4.851762 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 1.432349 on 2 degrees of freedom 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) > options(contrasts = c("contr.treatment", "contr.poly")) > glm <- glm(cars ~ lane, family = poisson) > summary(glm) Coefficients: Value Std. Error t value (Intercept) 5.6835798 0.05832118 97.4530924 lane2 -0.0631789 0.08381258 -0.7538117 lane3 -0.2113091 0.08719542 -2.4233967 lane4 -0.4260844 0.09278844 -4.5919989 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 0 on 0 degrees of freedom Number of Fisher Scoring Iterations: 1 > sum((cars - fitted(glm))^2/fitted(glm)) [1] 2.749881e-028 > cbind(cars, fitted(glm), resid(glm)) cars 1 294 294 2.250026e-007 2 276 276 1.632340e-007 3 238 238 -6.322027e-008 4 192 192 -1.192093e-007 > 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) Warning messages: Warning in glm.fitter(x = X, y = Y, w = w, start = start, offset = offset, family = : iterations terminated prematurely because of singularities > summary(glm) Warning messages: Warning in summary.glm(glm): This model has zero rank --- no summary is provided Call: glm(formula = x ~ -1 + theta, family = poisson(link = identity), offset = ex0) Coefficients: (1 not defined because of singularities) theta NA Degrees of Freedom: 4 Total; 4 Residual Residual Deviance: NA > sum((x - fitted(glm))^2/fitted(glm)) [1] NA > cbind(x, fitted(glm), resid(glm)) x 1 1997 NA NA 2 906 0 Inf 3 904 0 Inf 4 32 0 Inf > defect <- c(15, 21, 45, 13, 26, 31, 34, 5, 33, 17, 49, 20) > type <- factor(rep(1:4, 3)) > shift <- factor(rep(1:3, each = 4)) > 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) Coefficients: Value Std. Error t value (Intercept) 3.11401833 0.1445926 21.5364994 type2 -0.06995765 0.1672731 -0.4182242 type3 0.54796510 0.1460125 3.7528645 type4 -0.66646832 0.1992546 -3.3448072 shift2 0.02105489 0.1450153 0.1451908 shift3 0.23582879 0.1379455 1.7095792 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 20.33593 on 6 degrees of freedom Number of Fisher Scoring Iterations: 3 > sum((defect - fitted(glm))^2/fitted(glm)) [1] 19.17789 > cbind(defect, fitted(glm), resid(glm)) defect 1 15 22.51132 -1.686297053 2 21 20.99031 0.002115898 3 45 38.93850 0.947687510 4 13 11.55999 0.415166588 5 26 22.99032 0.614695720 6 31 21.43694 1.934921219 7 34 39.76704 -0.938073149 8 5 11.80596 -2.240593635 9 33 28.49838 0.822400206 10 17 26.57284 -1.989643817 11 49 49.29449 -0.041985380 12 20 14.63446 1.327799838 > 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) Coefficients: Value Std. Error t value (Intercept) 2.0672355 0.3471334 5.9551609 wife2 -0.4858911 0.3075091 -1.5800871 wife3 -0.2521026 0.3064957 -0.8225322 wife4 0.2742479 0.3592032 0.7634896 wife5 -0.3194487 0.2677991 -1.1928668 husb2 0.3108325 0.3494809 0.8894120 husb3 0.3007788 0.3403800 0.8836559 husb4 0.4399969 0.3364540 1.3077475 husb5 -0.1893293 0.3007842 -0.6294522 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 76.25068 on 8 degrees of freedom Number of Fisher Scoring Iterations: 4 > sum((marr - fitted(glm))^2/fitted(glm)) [1] 66.53164 > cbind(marr, fitted(glm), resid(glm)) marr 1 5 10.784048 -1.97022720 2 17 10.676172 1.78019960 3 6 6.539798 -0.21408934 4 5 4.861487 0.06252611 5 0 6.567434 -3.62420589 6 16 7.548444 2.67147932 7 2 4.022949 -1.11822430 8 2 8.380985 -2.65154617 9 10 9.536532 0.14888888 10 11 5.082502 2.26958571 11 10 10.396648 -0.12381009 12 9 8.603371 0.13420320 13 6 5.741881 0.10692686 14 20 7.835145 3.62695563 15 8 7.756767 0.08688325 16 0 8.915433 -4.22266091 17 1 4.751487 -2.09429179 > y <- c(2, 5, 4, 28, 35, 39, 38, 29, 22, 13, 9, 5) * 10 > opin <- rep(1:4, each = 3) > sond <- rep(1:3, 4) > Opin <- factor(opin) > Sond <- factor(sond) > glm <- glm(y ~ Opin + Sond, family = poisson) > summary(glm) Coefficients: Value Std. Error t value (Intercept) 3.66121219 0.09930144 36.8696788 Opin2 2.22707304 0.10021809 22.2222667 Opin3 2.09073658 0.10093033 20.7146504 Opin4 0.89793783 0.11297185 7.9483327 Sond2 -0.03774057 0.05015556 -0.7524704 Sond3 -0.14595395 0.05158812 -2.8292164 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 103.9902 on 6 degrees of freedom 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.90848 -3.3463055 2 50 37.46741 1.9467842 3 40 33.62461 1.0672007 4 280 360.78609 -4.4288402 5 350 347.42356 0.1380561 6 390 311.79043 4.2610948 7 380 314.80354 3.5575515 8 290 303.14408 -0.7604845 9 220 272.05243 -3.2654657 10 130 95.50227 3.3443224 11 90 91.96513 -0.2056538 12 50 82.53282 -3.8663082 > glm(y ~ Opin + Sond + Opin:Sond, family = poisson) Coefficients: (Intercept) Opin2 Opin3 Opin4 Sond2 Sond3 Opin2Sond2 2.995732 2.639057 2.944439 1.871802 0.9162907 0.6931472 -0.6931472 Opin3Sond2 Opin4Sond2 Opin2Sond3 Opin3Sond3 Opin4Sond3 -1.186581 -1.284016 -0.36179 -1.239691 -1.648659 Degrees of Freedom: 12 Total; 0 Residual Residual Deviance: 6.394885e-014 > x <- opin * sond > glm <- glm(y ~ Opin + Sond + x, family = poisson) > summary(glm) Coefficients: Value Std. Error t value (Intercept) 3.4013184 0.10866684 31.300425 Opin2 2.9488846 0.13213816 22.316677 Opin3 3.4639019 0.18686985 18.536441 Opin4 2.8530957 0.24193793 11.792676 Sond2 0.8424962 0.11021331 7.644233 Sond3 1.5535575 0.19061746 8.150132 x -0.3311723 0.03586207 -9.234613 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 14.05514 on 5 degrees of freedom 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 <- factor(rep(1:2, each = 4)) > care <- factor(rep(1:2, each = 2, length.out = 8)) > surv <- factor(rep(1:2, length.out = 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) Coefficients: (Intercept) care surv 1.96649 -0.1992522 3.277131 Degrees of Freedom: 8 Total; 5 Residual Residual Deviance: 291.5459 > glm(inf ~ care + surv + care:surv, family = poisson) Coefficients: (Intercept) care surv care:surv 2.302585 -1.203973 2.925846 1.038142 Degrees of Freedom: 8 Total; 4 Residual Residual Deviance: 285.9339 > glm(inf ~ care + surv + clin, family = poisson) Coefficients: (Intercept) care surv clin 2.252783 -0.1992562 3.277119 -0.6889495 Degrees of Freedom: 8 Total; 4 Residual Residual Deviance: 211.482 > glm(inf ~ care + surv + clin + care:surv, family = poisson) Coefficients: (Intercept) care surv clin care:surv 2.588874 -1.203979 2.92584 -0.6889504 1.038144 Degrees of Freedom: 8 Total; 3 Residual Residual Deviance: 205.8701 > glm(inf ~ care + surv + clin + care:clin, family = poisson) Coefficients: (Intercept) care surv clin care:clin 1.873201 0.5063463 3.277144 0.1785902 -2.653447 Degrees of Freedom: 8 Total; 3 Residual Residual Deviance: 17.8284 > glm(inf ~ care + surv + clin + care:clin + care:surv, family = poisson) Coefficients: (Intercept) care surv clin care:clin care:surv 2.209308 -0.4983683 2.925846 0.1785902 -2.653447 1.038137 Degrees of Freedom: 8 Total; 2 Residual Residual Deviance: 12.21642 > glm(inf ~ care + surv + clin + care:clin + surv:clin, family = poisson) Coefficients: (Intercept) care surv clin care:clin surv:clin 0.9678781 0.5063463 4.204693 1.866073 -2.653447 -1.755504 Degrees of Freedom: 8 Total; 2 Residual Residual Deviance: 0.08228918 > glm(inf ~ care + surv + clin + care:clin + surv:clin + care:surv, family = poisson) Coefficients: (Intercept) care surv clin care:clin surv:clin care:surv 1.034323 0.3976138 4.137222 1.809819 -2.646662 -1.699109 0.1103762 Degrees of Freedom: 8 Total; 1 Residual Residual Deviance: 0.04325585 > glm <- glm(inf ~ care + surv + clin + care:clin + surv:clin, family = poisson) > summary(glm) Coefficients: Value Std. Error t value (Intercept) 0.9678781 0.38254478 2.530104 care 0.5063463 0.09462343 5.351173 surv 4.2046926 0.38077153 11.042560 clin 1.8660733 0.44661044 4.178302 care:clin -2.6534465 0.23157394 -11.458312 surv:clin -1.7555041 0.44962650 -3.904361 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 0.0822892 on 2 degrees of freedom Number of Fisher Scoring Iterations: 3 > cbind(inf, fitted(glm), resid(glm)) inf 1 3 2.632353 0.2216100410 2 176 176.367647 -0.0276931670 3 4 4.367647 -0.1784755941 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