> "GunReg" <- + structure(.Data = list( + "pop" = c(4089, 2372, 30380, 3291, 598, 13277, 1135, 2795, 11543, 5996, + 4860, 9368, 4432, 5158, 6737, 635, 7760, 18058, 10939, 11961, 1004, + 3560, 4953, 17349, 1770, 5018, 570, 3750, 3377, 680, 6623, + 1039, 5610, 2495, 3713, 4252, 1235, 2592, 808, 1593, 1105, + 1548, 1284, 3175, 2922, 703, 6286, 567, 4955, 1801, 460.), + "area" = c(52.4, 53.2, 163.7, 5.5, 0.1, 65.8, 10.9, 56.3, 57.9, 10.6, + 12.4, 96.8, 86.9, 69.7, 53.8, 70.7, 8.7, 54.5, 44.8, 46.1, 1.5, 32, + 42.1, 268.6, 84.9, 71.3, 656.4, 114, 104.1, 2.5, 59.4, 83.6, 36.4, + 82.3, 40.4, 51.8, 35.4, 48.4, 147, 77.4, 9.4, 121.6, 110.6, 69.9, + 98.4, 77.1, 42.8, 9.6, 65.5, 24.2, 97.8), + "urban" = c(60, 54, 93, 79, 100, 85, 89, 61, 85, 84, 81, 70, 71, 53, + 50, 53, 89, 84, 74, 69, 86, 55, 61, 80, 87, 76, 68, 88, 82, + 73, 63, 57, 65, 69, 52, 68, 45, 47, 53, 66, 51, 73, 88, + 68, 71, 50, 69, 32, 66, 36, 65.), + "poverty" = c(19, 18.4, 14.2, 5.8, 19.2, 14.1, 10, 10.1, 13.3, 10.2, + 9.3, 13.9, 12, 13.6, 13.2, 13.5, 9, 14.1, 11.8, 10.8, 8.2, 16.5, + 16.9, 16.8, 9.8, 26.2, 11.2, 14.2, 12.1, 8.1, 16, 13.7, 14.1, 11.1, + 17.4, 22, 12.5, 23.8, 15.8, 10.9, 7.1, 20.9, 10.7, 15.8, 11.3, 13.5, + 10.6, 7.1, 9.2, 17.2, 10.6), + "gunreg" = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.), + "homicides" = c(410, 240, 3710, 170, 489, 1300, 44, 62, 1270, 200, 540, + 1020, 100, 550, 730, 11, 350, 2550, 760, 740, 38, 350, 470, 2660, + 43, 220, 56, 290, 155, 32, 720, 21, 380, 150, 260, 760, + 23, 370, 29, 43, 32, 160, 135, 220, 120, 9, 550, 24, 240, + 135, 20.)), + names = c("pop", "area", "urban", "poverty", "gunreg", "homicides"), + row.names = c("AL", "AR", "CA", "CT", "DC", "FL", "HI", "IA", "IL", "MA", + "MD", "MI", "MN", "MO", "NC", "ND", "NJ", "NY", "OH", "PA", "RI", "SC", + "TN", "TX", "UT", "WA", "AK", "AZ", "CO", "DE", "GA", "ID", "IN", "KS", + "KY", "LA", "ME", "MS", "MT", "NE", "NH", "NM", "NV", "OK", "OR", "SD", + "VA", "VT", "WI", "WV", "WY"), class = "data.frame") > summary(GunReg) pop area urban Min. : 460 Min. : 0.10 Min. : 32.00 1st Qu.: 1260 1st Qu.: 35.90 1st Qu.: 56.00 Median : 3377 Median : 56.30 Median : 69.00 Mean : 4945 Mean : 74.26 Mean : 68.51 3rd Qu.: 5803 3rd Qu.: 84.25 3rd Qu.: 81.50 Max. :30380 Max. :656.40 Max. :100.00 poverty gunreg homicides Min. : 5.80 Min. :0.0000 Min. : 9.0 1st Qu.:10.60 1st Qu.:0.0000 1st Qu.: 50.0 Median :13.30 Median :1.0000 Median : 220.0 Mean :13.47 Mean :0.5098 Mean : 469.8 3rd Qu.:15.90 3rd Qu.:1.0000 3rd Qu.: 545.0 Max. :26.20 Max. :1.0000 Max. :3710.0 > head(GunReg) pop area urban poverty gunreg homicides AL 4089 52.4 60 19.0 1 410 AR 2372 53.2 54 18.4 1 240 CA 30380 163.7 93 14.2 1 3710 CT 3291 5.5 79 5.8 1 170 DC 598 0.1 100 19.2 1 489 FL 13277 65.8 85 14.1 1 1300 > GunReg[order(GunReg$homicides),] pop area urban poverty gunreg homicides SD 703 77.1 50 13.5 0 9 ND 635 70.7 53 13.5 1 11 WY 460 97.8 65 10.6 0 20 ID 1039 83.6 57 13.7 0 21 ME 1235 35.4 45 12.5 0 23 VT 567 9.6 32 7.1 0 24 MT 808 147.0 53 15.8 0 29 DE 680 2.5 73 8.1 0 32 NH 1105 9.4 51 7.1 0 32 RI 1004 1.5 86 8.2 1 38 UT 1770 84.9 87 9.8 1 43 NE 1593 77.4 66 10.9 0 43 HI 1135 10.9 89 10.0 1 44 AK 570 656.4 68 11.2 0 56 IA 2795 56.3 61 10.1 1 62 MN 4432 86.9 71 12.0 1 100 OR 2922 98.4 71 11.3 0 120 NV 1284 110.6 88 10.7 0 135 WV 1801 24.2 36 17.2 0 135 KS 2495 82.3 69 11.1 0 150 CO 3377 104.1 82 12.1 0 155 NM 1548 121.6 73 20.9 0 160 CT 3291 5.5 79 5.8 1 170 MA 5996 10.6 84 10.2 1 200 WA 5018 71.3 76 26.2 1 220 OK 3175 69.9 68 15.8 0 220 AR 2372 53.2 54 18.4 1 240 WI 4955 65.5 66 9.2 0 240 KY 3713 40.4 52 17.4 0 260 AZ 3750 114.0 88 14.2 0 290 NJ 7760 8.7 89 9.0 1 350 SC 3560 32.0 55 16.5 1 350 MS 2592 48.4 47 23.8 0 370 IN 5610 36.4 65 14.1 0 380 AL 4089 52.4 60 19.0 1 410 TN 4953 42.1 61 16.9 1 470 DC 598 0.1 100 19.2 1 489 MD 4860 12.4 81 9.3 1 540 MO 5158 69.7 53 13.6 1 550 VA 6286 42.8 69 10.6 0 550 GA 6623 59.4 63 16.0 0 720 NC 6737 53.8 50 13.2 1 730 PA 11961 46.1 69 10.8 1 740 OH 10939 44.8 74 11.8 1 760 LA 4252 51.8 68 22.0 0 760 MI 9368 96.8 70 13.9 1 1020 IL 11543 57.9 85 13.3 1 1270 FL 13277 65.8 85 14.1 1 1300 NY 18058 54.5 84 14.1 1 2550 TX 17349 268.6 80 16.8 1 2660 CA 30380 163.7 93 14.2 1 3710 > glm(homicides ~ gunreg, family="poisson", data=GunReg) Call: glm(formula = homicides ~ gunreg, family = "poisson", data = GunReg) Coefficients: (Intercept) gunreg 5.285 1.310 Degrees of Freedom: 50 Total (i.e. Null); 49 Residual Null Deviance: 35580 Residual Deviance: 27270 AIC: 27640 > summary(glm(homicides ~ gunreg, family="poisson", data=GunReg)) Call: glm(formula = homicides ~ gunreg, family = "poisson", data = GunReg) Deviance Residuals: Min 1Q Median 3Q Max -36.732 -15.944 -7.440 3.592 78.027 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.28503 0.01424 371.24 <2e-16 *** gunreg 1.31049 0.01598 82.03 <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: 35577 on 50 degrees of freedom Residual deviance: 27272 on 49 degrees of freedom AIC: 27637 Number of Fisher Scoring iterations: 5 > exp(glm(homicides ~ gunreg, family="poisson", data=GunReg)$coef) (Intercept) gunreg 197.360000 3.707984 > GunReg[order(GunReg$pop),] pop area urban poverty gunreg homicides WY 460 97.8 65 10.6 0 20 VT 567 9.6 32 7.1 0 24 AK 570 656.4 68 11.2 0 56 DC 598 0.1 100 19.2 1 489 ND 635 70.7 53 13.5 1 11 DE 680 2.5 73 8.1 0 32 SD 703 77.1 50 13.5 0 9 MT 808 147.0 53 15.8 0 29 RI 1004 1.5 86 8.2 1 38 ID 1039 83.6 57 13.7 0 21 NH 1105 9.4 51 7.1 0 32 HI 1135 10.9 89 10.0 1 44 ME 1235 35.4 45 12.5 0 23 NV 1284 110.6 88 10.7 0 135 NM 1548 121.6 73 20.9 0 160 NE 1593 77.4 66 10.9 0 43 UT 1770 84.9 87 9.8 1 43 WV 1801 24.2 36 17.2 0 135 AR 2372 53.2 54 18.4 1 240 KS 2495 82.3 69 11.1 0 150 MS 2592 48.4 47 23.8 0 370 IA 2795 56.3 61 10.1 1 62 OR 2922 98.4 71 11.3 0 120 OK 3175 69.9 68 15.8 0 220 CT 3291 5.5 79 5.8 1 170 CO 3377 104.1 82 12.1 0 155 SC 3560 32.0 55 16.5 1 350 KY 3713 40.4 52 17.4 0 260 AZ 3750 114.0 88 14.2 0 290 AL 4089 52.4 60 19.0 1 410 LA 4252 51.8 68 22.0 0 760 MN 4432 86.9 71 12.0 1 100 MD 4860 12.4 81 9.3 1 540 TN 4953 42.1 61 16.9 1 470 WI 4955 65.5 66 9.2 0 240 WA 5018 71.3 76 26.2 1 220 MO 5158 69.7 53 13.6 1 550 IN 5610 36.4 65 14.1 0 380 MA 5996 10.6 84 10.2 1 200 VA 6286 42.8 69 10.6 0 550 GA 6623 59.4 63 16.0 0 720 NC 6737 53.8 50 13.2 1 730 NJ 7760 8.7 89 9.0 1 350 MI 9368 96.8 70 13.9 1 1020 OH 10939 44.8 74 11.8 1 760 IL 11543 57.9 85 13.3 1 1270 PA 11961 46.1 69 10.8 1 740 FL 13277 65.8 85 14.1 1 1300 TX 17349 268.6 80 16.8 1 2660 NY 18058 54.5 84 14.1 1 2550 CA 30380 163.7 93 14.2 1 3710 > GunReg$rate <- GunReg$homicides/GunReg$pop > summary(GunReg$rate) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.01280 0.03992 0.06929 0.08619 0.10420 0.81770 > GunReg[order(GunReg$rate),] pop area urban poverty gunreg homicides rate SD 703 77.1 50 13.5 0 9 0.01280228 ND 635 70.7 53 13.5 1 11 0.01732283 ME 1235 35.4 45 12.5 0 23 0.01862348 ID 1039 83.6 57 13.7 0 21 0.02021174 IA 2795 56.3 61 10.1 1 62 0.02218247 MN 4432 86.9 71 12.0 1 100 0.02256318 UT 1770 84.9 87 9.8 1 43 0.02429379 NE 1593 77.4 66 10.9 0 43 0.02699309 NH 1105 9.4 51 7.1 0 32 0.02895928 MA 5996 10.6 84 10.2 1 200 0.03335557 MT 808 147.0 53 15.8 0 29 0.03589109 RI 1004 1.5 86 8.2 1 38 0.03784861 HI 1135 10.9 89 10.0 1 44 0.03876652 OR 2922 98.4 71 11.3 0 120 0.04106776 VT 567 9.6 32 7.1 0 24 0.04232804 WY 460 97.8 65 10.6 0 20 0.04347826 WA 5018 71.3 76 26.2 1 220 0.04384217 NJ 7760 8.7 89 9.0 1 350 0.04510309 CO 3377 104.1 82 12.1 0 155 0.04589873 DE 680 2.5 73 8.1 0 32 0.04705882 WI 4955 65.5 66 9.2 0 240 0.04843592 CT 3291 5.5 79 5.8 1 170 0.05165603 KS 2495 82.3 69 11.1 0 150 0.06012024 PA 11961 46.1 69 10.8 1 740 0.06186774 IN 5610 36.4 65 14.1 0 380 0.06773619 OK 3175 69.9 68 15.8 0 220 0.06929134 OH 10939 44.8 74 11.8 1 760 0.06947619 KY 3713 40.4 52 17.4 0 260 0.07002424 WV 1801 24.2 36 17.2 0 135 0.07495836 AZ 3750 114.0 88 14.2 0 290 0.07733333 VA 6286 42.8 69 10.6 0 550 0.08749602 TN 4953 42.1 61 16.9 1 470 0.09489198 FL 13277 65.8 85 14.1 1 1300 0.09791369 AK 570 656.4 68 11.2 0 56 0.09824561 SC 3560 32.0 55 16.5 1 350 0.09831461 AL 4089 52.4 60 19.0 1 410 0.10026901 AR 2372 53.2 54 18.4 1 240 0.10118044 NM 1548 121.6 73 20.9 0 160 0.10335917 NV 1284 110.6 88 10.7 0 135 0.10514019 MO 5158 69.7 53 13.6 1 550 0.10663048 NC 6737 53.8 50 13.2 1 730 0.10835684 GA 6623 59.4 63 16.0 0 720 0.10871206 MI 9368 96.8 70 13.9 1 1020 0.10888130 IL 11543 57.9 85 13.3 1 1270 0.11002339 MD 4860 12.4 81 9.3 1 540 0.11111111 CA 30380 163.7 93 14.2 1 3710 0.12211982 NY 18058 54.5 84 14.1 1 2550 0.14121165 MS 2592 48.4 47 23.8 0 370 0.14274691 TX 17349 268.6 80 16.8 1 2660 0.15332296 LA 4252 51.8 68 22.0 0 760 0.17873942 DC 598 0.1 100 19.2 1 489 0.81772575 > hist(GunReg$rate, breaks=50, col="blue", xlab="Homicide Rate per 1000 person-years") > Y <- aggregate(homicides ~ gunreg, data = GunReg, sum) > Y gunreg homicides 1 0 4934 2 1 19027 > T <- aggregate(pop ~ gunreg, data = GunReg, sum) > T gunreg pop 1 0 63143 2 1 189038 > cbind(Y,T) gunreg homicides gunreg pop 1 0 4934 0 63143 2 1 19027 1 189038 > by.strata <- cbind(Y,T)[,-3] > by.strata$homrate <- by.strata$homicides/by.strata$pop > by.strata gunreg homicides pop homrate 1 0 4934 63143 0.07814009 2 1 19027 189038 0.10065172 > lm(rate ~ gunreg, data=GunReg) Call: lm(formula = rate ~ gunreg, data = GunReg) Coefficients: (Intercept) gunreg 0.06623 0.03917 > plot(lm(rate ~ gunreg, data=GunReg)) Hit to see next plot: Hit to see next plot: Hit to see next plot: Hit to see next plot: > sort(fitted(lm(rate ~ gunreg, data=GunReg))) AK GA NE NH NM 0.06622606 0.06622606 0.06622606 0.06622606 0.06622606 NV AZ CO DE ID 0.06622606 0.06622606 0.06622606 0.06622606 0.06622606 IN KS KY LA ME 0.06622606 0.06622606 0.06622606 0.06622606 0.06622606 MS MT OK OR SD 0.06622606 0.06622606 0.06622606 0.06622606 0.06622606 VA VT WI WV WY 0.06622606 0.06622606 0.06622606 0.06622606 0.06622606 AR IA MN ND UT 0.10539351 0.10539351 0.10539351 0.10539351 0.10539351 CA CT FL HI IL 0.10539351 0.10539351 0.10539351 0.10539351 0.10539351 MA MD MI MO NC 0.10539351 0.10539351 0.10539351 0.10539351 0.10539351 NJ NY OH PA RI 0.10539351 0.10539351 0.10539351 0.10539351 0.10539351 SC TN TX WA DC 0.10539351 0.10539351 0.10539351 0.10539351 0.10539351 AL 0.10539351 > boxplot(GunReg$rate ~ GunReg$gunreg, xlab="Gun registration laws", ylab="Homicide Rate per 1000 person-years") > summary(glm(homicides ~ gunreg, family="poisson", data=GunReg) + ) Call: glm(formula = homicides ~ gunreg, family = "poisson", data = GunReg) Deviance Residuals: Min 1Q Median 3Q Max -36.732 -15.944 -7.440 3.592 78.027 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.28503 0.01424 371.24 <2e-16 *** gunreg 1.31049 0.01598 82.03 <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: 35577 on 50 degrees of freedom Residual deviance: 27272 on 49 degrees of freedom AIC: 27637 Number of Fisher Scoring iterations: 5 > mod1 <- glm(homicides ~ gunreg, family=poisson, offset=log(pop), data=GunReg) > summary(mod1) Call: glm(formula = homicides ~ gunreg, family = poisson, data = GunReg, offset = log(pop)) Deviance Residuals: Min 1Q Median 3Q Max -19.827 -8.012 -2.853 2.114 34.513 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.54925 0.01424 -179.07 <2e-16 *** gunreg 0.25316 0.01598 15.85 <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: 5742.1 on 50 degrees of freedom Residual deviance: 5478.5 on 49 degrees of freedom AIC: 5843.6 Number of Fisher Scoring iterations: 4 > exp(mod1$coef) (Intercept) gunreg 0.07814009 1.28809315 > glm(rate ~ gunreg, family="poisson", data=GunReg) Call: glm(formula = rate ~ gunreg, family = "poisson", data = GunReg) Coefficients: (Intercept) gunreg -2.7147 0.4646 Degrees of Freedom: 50 Total (i.e. Null); 49 Residual Null Deviance: 3.487 Residual Deviance: 3.257 AIC: Inf There were 50 or more warnings (use warnings() to see the first 50) > pchisq(summary(mod1)$deviance, df=49, lower.tail=FALSE) [1] 0 > summary(mod1)$deviance [1] 5478.541 > GunReg$x <- rnorm(51, 40, 5) > summary(glm(homicides ~ gunreg+x, family=poisson, offset=log(pop), data=GunReg)) Call: glm(formula = homicides ~ gunreg + x, family = poisson, data = GunReg, offset = log(pop)) Deviance Residuals: Min 1Q Median 3Q Max -19.953 -8.107 -2.628 2.561 35.042 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.933165 0.043883 -66.840 <2e-16 *** gunreg 0.220231 0.016383 13.443 <2e-16 *** x 0.009892 0.001066 9.278 <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: 5742.1 on 50 degrees of freedom Residual deviance: 5392.9 on 48 degrees of freedom AIC: 5759.9 Number of Fisher Scoring iterations: 4 > mod2 <- glm(homicides ~ gunreg+poverty+urban, family=poisson, offset=log(pop), data=GunReg) > summary(mod2) Call: glm(formula = homicides ~ gunreg + poverty + urban, family = poisson, data = GunReg, offset = log(pop)) Deviance Residuals: Min 1Q Median 3Q Max -31.521 -5.950 -2.002 2.657 26.876 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.2546654 0.0508790 -83.623 < 2e-16 *** gunreg 0.1427145 0.0175846 8.116 4.82e-16 *** poverty 0.0640389 0.0016936 37.812 < 2e-16 *** urban 0.0116133 0.0005619 20.668 < 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: 5742.1 on 50 degrees of freedom Residual deviance: 3919.4 on 47 degrees of freedom AIC: 4288.4 Number of Fisher Scoring iterations: 5 > sort(residuals(mod2, type="deviance")) WA MN MA NJ UT -31.5214929 -16.1075129 -15.2999021 -12.3952427 -10.2610340 IA CO ID NE SD -10.0270924 -7.5947120 -6.7605900 -6.2420521 -6.2303196 HI ND ME OR OH -6.2273917 -6.1676780 -6.1075265 -5.7930802 -4.8192722 RI PA MT AZ FL -4.7735853 -4.5792935 -4.2653284 -4.1887716 -3.7337517 OK NM AL WI KY -3.3280417 -2.6314556 -2.0902829 -2.0340375 -2.0150354 NH IN CT WY DE -2.0016726 -1.8873902 -1.8700140 -1.4842768 -0.9787958 KS TN AR CA VT -0.8535794 -0.7354904 0.2420688 1.2024746 1.2480878 WV SC IL AK NV 1.6371827 1.7857251 2.3278148 2.9861384 3.2672072 MS MI GA VA MD 4.4039747 5.9233166 7.1624536 7.5121385 8.2736975 MO LA TX NC NY 8.5291348 8.7293578 11.5116719 11.6294032 13.3240860 DC 26.8763999 > plot(rate ~ poverty, data=GunReg) > points(predict(mod2, type="response")/GunReg$pop ~ GunReg$poverty, pch=3, col="red") > legend("topleft", c("Observed Rate", "Fitted Rate"), pch=c(1,3), col=c(1,"red")) >