> "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") > > GunReg$lrate <- log(GunReg$homicides/GunReg$pop) > mod.lm = lm(lrate ~ gunreg, data=GunReg) > summary(mod.lm) ## Note low r-squared --> gun registration Call: lm(formula = lrate ~ gunreg, data = GunReg) Residuals: Min 1Q Median 3Q Max -1.4570 -0.4689 0.2090 0.4318 2.4428 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.9011 0.1476 -19.659 <2e-16 gunreg 0.2571 0.2067 1.244 0.219 (Intercept) *** gunreg --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7379 on 49 degrees of freedom Multiple R-squared: 0.03061, Adjusted R-squared: 0.01083 F-statistic: 1.547 on 1 and 49 DF, p-value: 0.2195 > 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 (Intercept) -4.2546654 0.0508790 -83.623 gunreg 0.1427145 0.0175846 8.116 poverty 0.0640389 0.0016936 37.812 urban 0.0116133 0.0005619 20.668 Pr(>|z|) (Intercept) < 2e-16 *** gunreg 4.82e-16 *** poverty < 2e-16 *** urban < 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 > plot(residuals(mod2, type="pearson")^2 ~ fitted(mod2), xlab="Fitted Homicide Counts", ylab="Squared Pearson Residuals") > sfit <- supsmu(fitted(mod2), residuals(mod2, type="pearson")^2) > lines(sfit$x[ order(sfit$x) ], sfit$y[ order(sfit$x) ], col="red", lwd=2) > abline(h=1, lty=2, col="blue", lwd=2) > identify(fitted(mod2), residuals(mod2, type="pearson")^2, label=row.names(GunReg)) [1] 5 10 13 18 24 26 > presid2 <- residuals(mod2,type="pearson")^2 > phihat <- sum(presid2)/mod2$df.residual > phihat [1] 85.98062 > mod2.quasi <- glm(homicides ~ gunreg+poverty+urban, family=quasipoisson, offset=log(pop), data=GunReg) > summary(mod2.quasi) Call: glm(formula = homicides ~ gunreg + poverty + urban, family = quasipoisson, 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 t value Pr(>|t|) (Intercept) -4.25467 0.47178 -9.018 8.06e-12 gunreg 0.14271 0.16306 0.875 0.385886 poverty 0.06404 0.01570 4.078 0.000174 urban 0.01161 0.00521 2.229 0.030634 (Intercept) *** gunreg poverty *** urban * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 85.98062) Null deviance: 5742.1 on 50 degrees of freedom Residual deviance: 3919.4 on 47 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 > summary(mod2.quasi) Call: glm(formula = homicides ~ gunreg + poverty + urban, family = quasipoisson, 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 t value Pr(>|t|) (Intercept) -4.25467 0.47178 -9.018 8.06e-12 *** gunreg 0.14271 0.16306 0.875 0.385886 poverty 0.06404 0.01570 4.078 0.000174 *** urban 0.01161 0.00521 2.229 0.030634 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 85.98062) Null deviance: 5742.1 on 50 degrees of freedom Residual deviance: 3919.4 on 47 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 > summary(mod2)$ttable NULL > summary(mod2)$coef Estimate Std. Error z value Pr(>|z|) (Intercept) -4.25466540 0.050878993 -83.623224 0.000000e+00 gunreg 0.14271449 0.017584616 8.115872 4.823104e-16 poverty 0.06403887 0.001693599 37.812292 0.000000e+00 urban 0.01161333 0.000561889 20.668362 6.674260e-95 > summary(mod2.quasi)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) -4.25466540 0.471779202 -9.0183403 8.057210e-12 gunreg 0.14271449 0.163054645 0.8752556 3.858860e-01 poverty 0.06403887 0.015704024 4.0778638 1.742239e-04 urban 0.01161333 0.005210157 2.2289780 3.063425e-02 > 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 > summary(mod2.quasi) Call: glm(formula = homicides ~ gunreg + poverty + urban, family = quasipoisson, 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 t value Pr(>|t|) (Intercept) -4.25467 0.47178 -9.018 8.06e-12 *** gunreg 0.14271 0.16306 0.875 0.385886 poverty 0.06404 0.01570 4.078 0.000174 *** urban 0.01161 0.00521 2.229 0.030634 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 85.98062) Null deviance: 5742.1 on 50 degrees of freedom Residual deviance: 3919.4 on 47 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 > anova(mod2) Analysis of Deviance Table Model: poisson, link: log Response: homicides Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 50 5742.1 gunreg 1 263.56 49 5478.5 poverty 1 1117.54 48 4361.0 urban 1 441.62 47 3919.4 > anova(mod2, test="LRT") Analysis of Deviance Table Model: poisson, link: log Response: homicides Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 50 5742.1 gunreg 1 263.56 49 5478.5 < 2.2e-16 *** poverty 1 1117.54 48 4361.0 < 2.2e-16 *** urban 1 441.62 47 3919.4 < 2.2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(mod2.quasi, test="LRT") Analysis of Deviance Table Model: quasipoisson, link: log Response: homicides Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 50 5742.1 gunreg 1 263.56 49 5478.5 0.0799789 . poverty 1 1117.54 48 4361.0 0.0003119 *** urban 1 441.62 47 3919.4 0.0234316 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > lrtest(mod2.quasi) Error: could not find function "lrtest" > lrTest(mod2.quasi) Error: could not find function "lrTest" > ?lrtest No documentation for ‘lrtest’ in specified packages and libraries: you could try ‘??lrtest’ > ??lrtest starting httpd help server ... done > library(lmtest) Loading required package: zoo Attaching package: ‘zoo’ The following objects are masked from ‘package:base’: as.Date, as.Date.numeric > lrtest(mod2.quasi) Likelihood ratio test Model 1: homicides ~ gunreg + poverty + urban Model 2: homicides ~ 1 #Df LogLik Df Chisq Pr(>Chisq) 1 4 2 1 -3 > anova(mod2.quasi, test="deviance") Error in match.arg(test) : 'arg' should be one of “Rao”, “LRT”, “Chisq”, “F”, “Cp” > anova(mod2.quasi, test="Chisq") Analysis of Deviance Table Model: quasipoisson, link: log Response: homicides Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 50 5742.1 gunreg 1 263.56 49 5478.5 0.0799789 . poverty 1 1117.54 48 4361.0 0.0003119 *** urban 1 441.62 47 3919.4 0.0234316 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(mod2.quasi, test="LRT") Analysis of Deviance Table Model: quasipoisson, link: log Response: homicides Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 50 5742.1 gunreg 1 263.56 49 5478.5 0.0799789 . poverty 1 1117.54 48 4361.0 0.0003119 *** urban 1 441.62 47 3919.4 0.0234316 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(mod2.quasi, test="LRT") Analysis of Deviance Table Model: quasipoisson, link: log Response: homicides Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 50 5742.1 gunreg 1 263.56 49 5478.5 0.0799789 . poverty 1 1117.54 48 4361.0 0.0003119 *** urban 1 441.62 47 3919.4 0.0234316 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(update(mod2.quasi, .~.-gunreg), mod2.quasi) Analysis of Deviance Table Model 1: homicides ~ poverty + urban Model 2: homicides ~ gunreg + poverty + urban Resid. Df Resid. Dev Df Deviance 1 48 3986.4 2 47 3919.4 1 67.048 > anova(update(mod2.quasi, .~.-gunreg), mod2.quasi, test="Chisq") Analysis of Deviance Table Model 1: homicides ~ poverty + urban Model 2: homicides ~ gunreg + poverty + urban Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 48 3986.4 2 47 3919.4 1 67.048 0.3772 > anova(update(mod2.quasi, .~.-gunreg), mod2.quasi, test="LRT") Analysis of Deviance Table Model 1: homicides ~ poverty + urban Model 2: homicides ~ gunreg + poverty + urban Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 48 3986.4 2 47 3919.4 1 67.048 0.3772 > ts <- 67.048/85.98062 > ts [1] 0.7798036 > pchisq(ts, df=1, lower.tail=FALSE) [1] 0.3772012 > library(MASS) > ?glm.nb # Uses log link by default > mod2.nb <- glm.nb(homicides ~ gunreg+poverty+urban+offset(log(pop)), data=GunReg) > summary(mod2.nb) Call: glm.nb(formula = homicides ~ gunreg + poverty + urban + offset(log(pop)), data = GunReg, init.theta = 3.418144423, link = log) Deviance Residuals: Min 1Q Median 3Q Max -2.78734 -0.84712 -0.06245 0.28497 2.73334 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -5.262017 0.478178 -11.004 < 2e-16 *** gunreg 0.158352 0.166324 0.952 0.34106 poverty 0.100595 0.018321 5.491 4.01e-08 *** urban 0.017749 0.005608 3.165 0.00155 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for Negative Binomial(3.4181) family taken to be 1) Null deviance: 106.179 on 50 degrees of freedom Residual deviance: 53.877 on 47 degrees of freedom AIC: 634.09 Number of Fisher Scoring iterations: 1 Theta: 3.418 Std. Err.: 0.677 2 x log-likelihood: -624.088 > 1/3.418 [1] 0.2925688 > ?dnbinom >