Homework 6
Mahmuda Sultana
3/14/2022
1) Define a count outcome for the dataset of your choosing, the Area Resource File used in class provides many options here
##a. State a research question about your outcome:
My Count outcome is Infant Mortality Rate. Research queation is what is the prevalence of infant mortality rate at county level in USA?
temp <- tempfile()
#Download the SAS dataset as a ZIP compressed archive
download.file("https://data.hrsa.gov/DataDownload/AHRF/AHRF_2019-2020_SAS.zip", temp)
#Read SAS data into R
ahrf<-haven::read_sas(unz(temp,
filename = "ahrf2020.sas7bdat"))
rm(temp)library(tidyverse)
ahrf2<-ahrf%>%
mutate(cofips = f00004,
coname = f00010,
state = f00011,
popn = f1198414,
births1618 = f1254616,
infantmor1618 = f1194816,
fampov14 = f1443214,
infantmortrate1618 = 1000*(f1194816/f1254616), #Rate per 1000 births
hpsa16= case_when(.$f0978716 == 0 ~ 'no shortage',
.$f0978716 == 1 ~ 'whole county shortage',
.$f0978716 == 2 ~ 'partial county shortage'),
rucc = as.factor(f0002013) )%>%
mutate(rucc = droplevels(rucc, ""))%>%
dplyr::select(births1618,
infantmor1618,
infantmortrate1618,
state,
cofips,
coname,
popn,
fampov14,
hpsa16,
rucc)%>%
filter(complete.cases(.))%>%
as.data.frame()options(tigris_class="sf")
library(tigris)## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.library(sf)## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUEusco<-counties(cb = T, year= 2016)##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%usco$cofips<-usco$GEOID
sts<-states(cb = T, year = 2016)##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%sts<-st_boundary(sts)%>%
filter(!STATEFP %in% c("02", "15", "60", "66", "69", "72", "78"))%>%
st_transform(crs = 2163)
ahrf_m<-left_join(usco, ahrf2,
by = "cofips")%>%
filter(is.na(infantmortrate1618 )==F,
!STATEFP %in% c("02", "15", "60", "66", "69", "72", "78"))%>%
st_transform(crs = 2163)
glimpse(ahrf_m)## Rows: 452
## Columns: 20
## $ STATEFP <chr> "01", "01", "01", "01", "04", "04", "05", "05", "06~
## $ COUNTYFP <chr> "003", "073", "097", "101", "001", "013", "131", "1~
## $ COUNTYNS <chr> "00161527", "00161562", "00161575", "00161577", "00~
## $ AFFGEOID <chr> "0500000US01003", "0500000US01073", "0500000US01097~
## $ GEOID <chr> "01003", "01073", "01097", "01101", "04001", "04013~
## $ NAME <chr> "Baldwin", "Jefferson", "Mobile", "Montgomery", "Ap~
## $ LSAD <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06~
## $ ALAND <dbl> 4117584019, 2878224117, 3184131335, 2031491982, 290~
## $ AWATER <dbl> 1133130502, 32446239, 1073840930, 40270068, 5417416~
## $ cofips <chr> "01003", "01073", "01097", "01101", "04001", "04013~
## $ births1618 <dbl> 2286, 8599, 5553, 3150, 959, 52871, 1737, 918, 1471~
## $ infantmor1618 <dbl> 12, 83, 48, 28, 10, 293, 11, 10, 91, 475, 26, 26, 1~
## $ infantmortrate1618 <dbl> 5.249344, 9.652285, 8.643976, 8.888889, 10.427529, ~
## $ state <chr> "01", "01", "01", "01", "04", "04", "05", "05", "06~
## $ coname <chr> "Baldwin", "Jefferson", "Mobile", "Montgomery", "Ap~
## $ popn <dbl> 200111, 660793, 415123, 226189, 71828, 4087191, 126~
## $ fampov14 <dbl> 7.3, 13.2, 15.0, 16.3, 28.2, 10.6, 14.9, 12.8, 19.7~
## $ hpsa16 <chr> "partial county shortage", "partial county shortage~
## $ rucc <fct> 03, 01, 02, 02, 06, 01, 02, 04, 02, 01, 02, 02, 02,~
## $ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((1157264 -15..., MULTIP~Here is a ggplot() histogram of the infant mortality rate for US counties.
ahrf_m%>%
ggplot()+
geom_histogram(aes(x = infantmortrate1618))+
labs(title = "Distribution of Infant Mortality Rates in US Counties",
subtitle = "2016 - 2018")+
xlab("Rate per 1,000 Live Births")+
ylab ("Frequency")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
library(tmap)
tm_shape(ahrf_m)+
tm_polygons(col = "infantmortrate1618",
border.col = NULL,
title="Infant Mortality Rt",
palette="Blues",
style="quantile",
n=5,
showNA=T, colorNA = "grey50")+
tm_format(format= "World",
main.title="US Infant Mortality Rate by County",
legend.position = c("left", "bottom"),
main.title.position =c("center"))+
tm_scale_bar(position = c(.1,0))+
tm_compass()+
tm_shape(sts)+
tm_lines( col = "black")
glm1<- glm(infantmortrate1618 ~ fampov14 + rucc + hpsa16,
data = ahrf_m,
family =gaussian)
glm1<-glm1%>%
tbl_regression()
summary(glm1)## Length Class Mode
## table_body 24 broom.helpers list
## table_styling 7 -none- list
## N 1 -none- numeric
## n 1 -none- numeric
## model_obj 30 glm list
## inputs 10 -none- list
## call_list 15 -none- listglmb<- glm(cbind(infantmor1618, births1618-infantmor1618) ~ rucc + fampov14 + hpsa16,
data = ahrf_m,
family = binomial)
glmb%>%
tbl_regression()| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| rucc | |||
| 01 | — | — | |
| 02 | 0.10 | 0.06, 0.14 | <0.001 |
| 03 | 0.23 | 0.17, 0.30 | <0.001 |
| 04 | 0.36 | 0.15, 0.55 | <0.001 |
| 05 | 0.52 | -0.18, 1.1 | 0.10 |
| 06 | 0.28 | -0.42, 0.85 | 0.4 |
| % Families Below Poverty Level 2014-18 | 0.02 | 0.02, 0.03 | <0.001 |
| hpsa16 | |||
| no shortage | — | — | |
| partial county shortage | -0.02 | -0.09, 0.05 | 0.5 |
| whole county shortage | -0.03 | -0.23, 0.17 | 0.8 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
summary(glmb)##
## Call:
## glm(formula = cbind(infantmor1618, births1618 - infantmor1618) ~
## rucc + fampov14 + hpsa16, family = binomial, data = ahrf_m)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.4343 -0.5974 0.1788 0.9634 5.6783
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.419224 0.036735 -147.522 < 2e-16 ***
## rucc02 0.100698 0.018421 5.466 4.59e-08 ***
## rucc03 0.234864 0.031839 7.377 1.62e-13 ***
## rucc04 0.356595 0.103241 3.454 0.000552 ***
## rucc05 0.516925 0.318248 1.624 0.104315
## rucc06 0.276874 0.319680 0.866 0.386437
## fampov14 0.021709 0.001827 11.882 < 2e-16 ***
## hpsa16partial county shortage -0.022653 0.036479 -0.621 0.534606
## hpsa16whole county shortage -0.027322 0.102078 -0.268 0.788962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1465.7 on 451 degrees of freedom
## Residual deviance: 1189.7 on 443 degrees of freedom
## AIC: 3483.5
##
## Number of Fisher Scoring iterations: 4glmb%>%
tbl_regression(exponentiate=TRUE)| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| rucc | |||
| 01 | — | — | |
| 02 | 1.11 | 1.07, 1.15 | <0.001 |
| 03 | 1.26 | 1.19, 1.35 | <0.001 |
| 04 | 1.43 | 1.16, 1.74 | <0.001 |
| 05 | 1.68 | 0.84, 2.95 | 0.10 |
| 06 | 1.32 | 0.66, 2.33 | 0.4 |
| % Families Below Poverty Level 2014-18 | 1.02 | 1.02, 1.03 | <0.001 |
| hpsa16 | |||
| no shortage | — | — | |
| partial county shortage | 0.98 | 0.91, 1.05 | 0.5 |
| whole county shortage | 0.97 | 0.79, 1.18 | 0.8 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
b. Is an offset term necessary? why or why not? If an offset term is need, what is the appropriate offset?
summary(ahrf_m$infantmor1618)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.00 13.00 20.00 35.88 38.00 475.00ahrf_m$infantmor1618<- ahrf_m$infantmor1618 +1Yes an offset term is neccesary. One of the conditions of using the poisson model is that each area or person has the same risk or population size. Since this is rear in demography and not the case for my population of study, poisson second type modelling( rate model) will be used.
The offset term is how unequal population size is incorporated. Possion is a strictly integer count so we can’t use the rate/probability so we log instead.
After doing the summary of the data, the outcome variable( infant mortality rat) contain zero. This would bring an error if done directly. So I added 1 to the variable(population). I used 1 because they are individuals and we can’t have 0.9 or 0.1 individual.
2) Consider a Poisson regression model for the outcome
glmp_s <- glm(infantmor1618 ~ offset(log(births1618)) + fampov14 + rucc + hpsa16,
data=ahrf_m,
family=poisson)
glmp_s%>%
tbl_regression(exp = TRUE)| Characteristic | IRR1 | 95% CI1 | p-value |
|---|---|---|---|
| % Families Below Poverty Level 2014-18 | 1.02 | 1.02, 1.02 | <0.001 |
| rucc | |||
| 01 | — | — | |
| 02 | 1.13 | 1.09, 1.17 | <0.001 |
| 03 | 1.33 | 1.25, 1.41 | <0.001 |
| 04 | 1.53 | 1.26, 1.85 | <0.001 |
| 05 | 1.82 | 0.95, 3.12 | 0.047 |
| 06 | 1.46 | 0.75, 2.51 | 0.2 |
| hpsa16 | |||
| no shortage | — | — | |
| partial county shortage | 0.95 | 0.89, 1.02 | 0.2 |
| whole county shortage | 0.94 | 0.77, 1.14 | 0.5 |
|
1
IRR = Incidence Rate Ratio, CI = Confidence Interval
|
|||
a. Evaluate the level of dispersion in the outcome
scale<-sqrt(glmp_s$deviance/glmp_s$df.residual)
scale## [1] 1.651197Here the scale factor is only 1.65, so it doesn’t show so much overdispersion and not that bad.
1-pchisq(glmp_s$deviance,
df = glmp_s$df.residual)## [1] 0b. Is the Poisson model a good choice?
No it is not a good choice since the goodness of fit results to zero. Since its zero it means the poisson distribution is a wrong choice
3) Consider a Negative binomial model
library(MASS) ##
## Attaching package: 'MASS'## The following object is masked from 'package:gtsummary':
##
## select## The following object is masked from 'package:srvyr':
##
## select## The following object is masked from 'package:dplyr':
##
## selectglmnb<- glm.nb(infantmor1618 ~ offset(log(births1618)) + fampov14 + rucc + hpsa16,
data=ahrf_m)
glmnb%>%
tbl_regression( exp= T)| Characteristic | IRR1 | 95% CI1 | p-value |
|---|---|---|---|
| % Families Below Poverty Level 2014-18 | 1.03 | 1.02, 1.04 | <0.001 |
| rucc | |||
| 01 | — | — | |
| 02 | 1.11 | 1.04, 1.17 | <0.001 |
| 03 | 1.27 | 1.16, 1.38 | <0.001 |
| 04 | 1.39 | 1.09, 1.77 | 0.008 |
| 05 | 1.65 | 0.78, 3.26 | 0.2 |
| 06 | 1.21 | 0.57, 2.42 | 0.6 |
| hpsa16 | |||
| no shortage | — | — | |
| partial county shortage | 0.92 | 0.84, 1.02 | 0.10 |
| whole county shortage | 0.92 | 0.72, 1.18 | 0.5 |
|
1
IRR = Incidence Rate Ratio, CI = Confidence Interval
|
|||
summary(glmnb)##
## Call:
## glm.nb(formula = infantmor1618 ~ offset(log(births1618)) + fampov14 +
## rucc + hpsa16, data = ahrf_m, init.theta = 25.76070655, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2105 -0.4910 0.0722 0.6014 2.0557
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.362080 0.050038 -107.161 < 2e-16 ***
## fampov14 0.027745 0.003279 8.462 < 2e-16 ***
## rucc02 0.101056 0.029872 3.383 0.000717 ***
## rucc03 0.235238 0.042458 5.540 3.02e-08 ***
## rucc04 0.329877 0.123602 2.669 0.007611 **
## rucc05 0.498916 0.361543 1.380 0.167598
## rucc06 0.190596 0.365942 0.521 0.602481
## hpsa16partial county shortage -0.078926 0.048441 -1.629 0.103243
## hpsa16whole county shortage -0.082220 0.127529 -0.645 0.519111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(25.7607) family taken to be 1)
##
## Null deviance: 536.11 on 451 degrees of freedom
## Residual deviance: 385.65 on 443 degrees of freedom
## AIC: 3049.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 25.76
## Std. Err.: 2.82
##
## 2 x log-likelihood: -3029.254library(gamlss)## Warning: package 'gamlss' was built under R version 4.1.3## Loading required package: gamlss.data##
## Attaching package: 'gamlss.data'## The following object is masked from 'package:datasets':
##
## sleep## Loading required package: gamlss.dist## Warning: package 'gamlss.dist' was built under R version 4.1.3## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## Loading required package: parallel## ********** GAMLSS Version 5.4-1 **********## For more on GAMLSS look at https://www.gamlss.com/## Type gamlssNews() to see new features/changes/bug fixes.##
## Attaching package: 'gamlss'## The following object is masked from 'package:caret':
##
## calibrationglmnb2<-gamlss(infantmor1618 ~ offset(log(births1618)) + fampov14 + rucc + hpsa16,
family = NBII,
data=ahrf_m)## GAMLSS-RS iteration 1: Global Deviance = 3947.202
## GAMLSS-RS iteration 2: Global Deviance = 3464.864
## GAMLSS-RS iteration 3: Global Deviance = 3226.917
## GAMLSS-RS iteration 4: Global Deviance = 3216.66
## GAMLSS-RS iteration 5: Global Deviance = 3216.633
## GAMLSS-RS iteration 6: Global Deviance = 3216.632summary(glmnb2)## ******************************************************************
## Family: c("NBII", "Negative Binomial type II")
##
## Call: gamlss(formula = infantmor1618 ~ offset(log(births1618)) +
## fampov14 + rucc + hpsa16, family = NBII, data = ahrf_m)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.334033 0.055923 -95.382 < 2e-16 ***
## fampov14 0.017416 0.002804 6.210 1.22e-09 ***
## rucc02 0.144745 0.029057 4.981 9.07e-07 ***
## rucc03 0.336701 0.048260 6.977 1.11e-11 ***
## rucc04 0.518595 0.151743 3.418 0.00069 ***
## rucc05 0.698215 0.455533 1.533 0.12605
## rucc06 0.506609 0.457940 1.106 0.26921
## hpsa16partial county shortage -0.062074 0.055714 -1.114 0.26581
## hpsa16whole county shortage -0.091592 0.155558 -0.589 0.55630
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4780 0.1084 4.41 1.3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 452
## Degrees of Freedom for the fit: 10
## Residual Deg. of Freedom: 442
## at cycle: 6
##
## Global Deviance: 3216.632
## AIC: 3236.632
## SBC: 3277.769
## ******************************************************************#4) Compare the model fits of the alternative models using AIC
AIC(glm1, glmb, glmp_b, glmnb, glmnb2)
AIC(glmp_s, glmb,glmnb, glmnb2)Among these four models, the poisson model fits the worst and the NB1 model, estimated by glm.nb() fits the best, with the lowest AIC among these four models.
LS0tDQp0aXRsZTogIkhvbWV3b3JrIDYiDQphdXRob3I6ICJNYWhtdWRhIFN1bHRhbmEiDQpkYXRlOiAiMy8xNC8yMDIyIg0Kb3V0cHV0OiANCiAgICAgIGh0bWxfZG9jdW1lbnQ6DQogICAgICAgIGRmX3ByaW50OiBwYWdlZA0KICAgICAgICBmaWdfaGVpZ2h0OiA3DQogICAgICAgIGZpZ193aWR0aDogNw0KICAgICAgICB0b2M6IHllcw0KICAgICAgICB0b2NfZmxvYXQ6IHllcw0KICAgICAgICBjb2RlX2Rvd25sb2FkOiB5ZXMNCg0KLS0tDQojIDEpIERlZmluZSBhIGNvdW50IG91dGNvbWUgZm9yIHRoZSBkYXRhc2V0IG9mIHlvdXIgY2hvb3NpbmcsIHRoZSBBcmVhIFJlc291cmNlIEZpbGUgdXNlZCBpbiBjbGFzcyBwcm92aWRlcyBtYW55IG9wdGlvbnMgaGVyZQ0KIyNhLiBTdGF0ZSBhIHJlc2VhcmNoIHF1ZXN0aW9uIGFib3V0IHlvdXIgb3V0Y29tZToNCiAgIA0KIE15IENvdW50IG91dGNvbWUgaXMgSW5mYW50IE1vcnRhbGl0eSBSYXRlLiANCiBSZXNlYXJjaCBxdWVhdGlvbiBpcyB3aGF0IGlzIHRoZSBwcmV2YWxlbmNlIG9mIGluZmFudCBtb3J0YWxpdHkgcmF0ZSBhdCBjb3VudHkgbGV2ZWwgaW4gVVNBPw0KDQpgYGB7ciBpbmNsdWRlPUZBTFNFfQ0KbGlicmFyeShzdGFyZ2F6ZXIsIHF1aWV0bHkgPSBUKQ0KbGlicmFyeShzdXJ2ZXksIHF1aWV0bHkgPSBUKQ0KbGlicmFyeShjYXIsIHF1aWV0bHkgPSBUKQ0KbGlicmFyeShxdWVzdGlvbnIsIHF1aWV0bHkgPSBUKQ0KbGlicmFyeShkcGx5ciwgcXVpZXRseSA9IFQpDQpsaWJyYXJ5KGZvcmNhdHMsIHF1aWV0bHkgPSBUKQ0KbGlicmFyeSh0aWR5dmVyc2UsIHF1aWV0bHkgPSBUKQ0KbGlicmFyeShzcnZ5ciwgcXVpZXRseSA9IFQpDQpsaWJyYXJ5KCBndHN1bW1hcnksIHF1aWV0bHkgPSBUKQ0KbGlicmFyeShjYXJldCwgcXVpZXRseSA9IFQpDQpsaWJyYXJ5KFZHQU0sIHF1aWV0bHkgPSBUKQ0KbGlicmFyeShnZ3Bsb3QyLCBxdWlldGx5ID0gVCkNCmxpYnJhcnkoc3Z5VkdBTSwgcXVpZXRseSA9IFQpDQpgYGANCg0KDQpgYGB7cn0NCnRlbXAgPC0gdGVtcGZpbGUoKQ0KDQojRG93bmxvYWQgdGhlIFNBUyBkYXRhc2V0IGFzIGEgWklQIGNvbXByZXNzZWQgYXJjaGl2ZQ0KZG93bmxvYWQuZmlsZSgiaHR0cHM6Ly9kYXRhLmhyc2EuZ292L0RhdGFEb3dubG9hZC9BSFJGL0FIUkZfMjAxOS0yMDIwX1NBUy56aXAiLCB0ZW1wKQ0KDQojUmVhZCBTQVMgZGF0YSBpbnRvIFINCmFocmY8LWhhdmVuOjpyZWFkX3Nhcyh1bnoodGVtcCwNCiAgICAgICAgICAgICAgICAgICAgICAgICAgZmlsZW5hbWUgPSAiYWhyZjIwMjAuc2FzN2JkYXQiKSkNCg0Kcm0odGVtcCkNCg0KDQpgYGANCg0KYGBge3J9DQoNCiANCmxpYnJhcnkodGlkeXZlcnNlKQ0KDQphaHJmMjwtYWhyZiU+JQ0KICBtdXRhdGUoY29maXBzID0gZjAwMDA0LA0KICAgICAgICAgY29uYW1lID0gZjAwMDEwLA0KICAgICAgICAgc3RhdGUgPSBmMDAwMTEsDQogICAgICAgICBwb3BuID0gIGYxMTk4NDE0LA0KICAgICAgICAgYmlydGhzMTYxOCA9IGYxMjU0NjE2LA0KICAgICAgICAgaW5mYW50bW9yMTYxOCA9IGYxMTk0ODE2LA0KICAgICAgICAgZmFtcG92MTQgPSBmMTQ0MzIxNCwNCiAgICAgICAgIGluZmFudG1vcnRyYXRlMTYxOCA9IDEwMDAqKGYxMTk0ODE2L2YxMjU0NjE2KSwgI1JhdGUgcGVyIDEwMDAgYmlydGhzDQogICAgICAgICBocHNhMTY9IGNhc2Vfd2hlbiguJGYwOTc4NzE2ID09IDAgfiAnbm8gc2hvcnRhZ2UnLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgLiRmMDk3ODcxNiA9PSAxIH4gJ3dob2xlIGNvdW50eSBzaG9ydGFnZScsDQogICAgICAgICAgICAgICAgICAgICAgICAgICAuJGYwOTc4NzE2ID09IDIgfiAncGFydGlhbCBjb3VudHkgc2hvcnRhZ2UnKSwNCiAgICAgICAgIA0KICAgICAgICAgcnVjYyA9IGFzLmZhY3RvcihmMDAwMjAxMykgKSU+JQ0KICAgICAgICAgDQogICAgICAgDQogIG11dGF0ZShydWNjID0gZHJvcGxldmVscyhydWNjLCAiIikpJT4lDQogIGRwbHlyOjpzZWxlY3QoYmlydGhzMTYxOCwNCiAgICAgICAgICAgICAgICBpbmZhbnRtb3IxNjE4LA0KICAgICAgICAgICAgICAgaW5mYW50bW9ydHJhdGUxNjE4LA0KICAgICAgICAgICAgICAgIHN0YXRlLA0KICAgICAgICAgICAgICAgIGNvZmlwcywNCiAgICAgICAgICAgICAgICBjb25hbWUsDQogICAgICAgICAgICAgICAgcG9wbiwNCiAgICAgICAgICAgICAgICBmYW1wb3YxNCwNCiAgICAgICAgICAgICAgICBocHNhMTYsDQogICAgICAgICAgICAgICAgcnVjYyklPiUNCiAgZmlsdGVyKGNvbXBsZXRlLmNhc2VzKC4pKSU+JQ0KICBhcy5kYXRhLmZyYW1lKCkNCg0KYGBgDQoNCg0KYGBge3J9DQpvcHRpb25zKHRpZ3Jpc19jbGFzcz0ic2YiKQ0KbGlicmFyeSh0aWdyaXMpDQpsaWJyYXJ5KHNmKQ0KdXNjbzwtY291bnRpZXMoY2IgPSBULCB5ZWFyPSAyMDE2KQ0KDQp1c2NvJGNvZmlwczwtdXNjbyRHRU9JRA0KDQpzdHM8LXN0YXRlcyhjYiA9IFQsIHllYXIgPSAyMDE2KQ0KDQpzdHM8LXN0X2JvdW5kYXJ5KHN0cyklPiUNCiAgZmlsdGVyKCFTVEFURUZQICVpbiUgYygiMDIiLCAiMTUiLCAiNjAiLCAiNjYiLCAiNjkiLCAiNzIiLCAiNzgiKSklPiUNCiAgc3RfdHJhbnNmb3JtKGNycyA9IDIxNjMpDQoNCmFocmZfbTwtbGVmdF9qb2luKHVzY28sIGFocmYyLA0KICAgICAgICAgICAgICAgICAgICBieSA9ICJjb2ZpcHMiKSU+JQ0KICBmaWx0ZXIoaXMubmEoaW5mYW50bW9ydHJhdGUxNjE4ICk9PUYsDQogICAgICAgICAhU1RBVEVGUCAlaW4lIGMoIjAyIiwgIjE1IiwgIjYwIiwgIjY2IiwgIjY5IiwgIjcyIiwgIjc4IikpJT4lDQogIHN0X3RyYW5zZm9ybShjcnMgPSAyMTYzKQ0KDQpnbGltcHNlKGFocmZfbSkNCg0KDQpgYGANCg0KDQpIZXJlIGlzIGEgZ2dwbG90KCkgaGlzdG9ncmFtIG9mIHRoZSBpbmZhbnQgbW9ydGFsaXR5IHJhdGUgZm9yIFVTIGNvdW50aWVzLg0KIA0KYGBge3J9DQoNCmFocmZfbSU+JQ0KICBnZ3Bsb3QoKSsNCiAgZ2VvbV9oaXN0b2dyYW0oYWVzKHggPSBpbmZhbnRtb3J0cmF0ZTE2MTgpKSsNCiAgbGFicyh0aXRsZSA9ICJEaXN0cmlidXRpb24gb2YgSW5mYW50IE1vcnRhbGl0eSBSYXRlcyBpbiBVUyBDb3VudGllcyIsDQogICAgICAgc3VidGl0bGUgPSAiMjAxNiAtIDIwMTgiKSsNCiAgICAgICB4bGFiKCJSYXRlIHBlciAxLDAwMCBMaXZlIEJpcnRocyIpKw0KICB5bGFiICgiRnJlcXVlbmN5IikNCg0KYGBgDQoNCmBgYHtyfQ0KbGlicmFyeSh0bWFwKQ0KDQp0bV9zaGFwZShhaHJmX20pKw0KICB0bV9wb2x5Z29ucyhjb2wgPSAiaW5mYW50bW9ydHJhdGUxNjE4IiwNCiAgICAgICAgICAgICAgYm9yZGVyLmNvbCA9IE5VTEwsDQogICAgICAgICAgICAgIHRpdGxlPSJJbmZhbnQgTW9ydGFsaXR5IFJ0IiwNCiAgICAgICAgICAgICAgcGFsZXR0ZT0iQmx1ZXMiLA0KICAgICAgICAgICAgICBzdHlsZT0icXVhbnRpbGUiLA0KICAgICAgICAgICAgICBuPTUsDQogICAgICAgICAgICAgIHNob3dOQT1ULCBjb2xvck5BID0gImdyZXk1MCIpKw0KICAgdG1fZm9ybWF0KGZvcm1hdD0gIldvcmxkIiwNCiAgICAgICAgICAgICBtYWluLnRpdGxlPSJVUyBJbmZhbnQgTW9ydGFsaXR5IFJhdGUgYnkgQ291bnR5IiwNCiAgICAgICAgICAgIGxlZ2VuZC5wb3NpdGlvbiA9ICBjKCJsZWZ0IiwgImJvdHRvbSIpLA0KICAgICAgICAgICAgbWFpbi50aXRsZS5wb3NpdGlvbiA9YygiY2VudGVyIikpKw0KICB0bV9zY2FsZV9iYXIocG9zaXRpb24gPSBjKC4xLDApKSsNCiAgdG1fY29tcGFzcygpKw0KdG1fc2hhcGUoc3RzKSsNCiAgdG1fbGluZXMoIGNvbCA9ICJibGFjayIpDQoNCmBgYA0KDQoNCmBgYHtyfQ0KZ2xtMTwtIGdsbShpbmZhbnRtb3J0cmF0ZTE2MTggfiBmYW1wb3YxNCArICBydWNjICsgaHBzYTE2LA0KICAgICAgICAgIGRhdGEgPSBhaHJmX20sDQogICAgICAgICAgZmFtaWx5ID1nYXVzc2lhbikNCg0KZ2xtMTwtZ2xtMSU+JQ0KICB0YmxfcmVncmVzc2lvbigpDQoNCnN1bW1hcnkoZ2xtMSkNCmBgYA0KDQpgYGB7cn0NCiANCmdsbWI8LSBnbG0oY2JpbmQoaW5mYW50bW9yMTYxOCwgYmlydGhzMTYxOC1pbmZhbnRtb3IxNjE4KSB+ICBydWNjICsgZmFtcG92MTQgKyBocHNhMTYsDQogICAgICAgICAgZGF0YSA9IGFocmZfbSwNCiAgICAgICAgICBmYW1pbHkgPSBiaW5vbWlhbCkNCg0KZ2xtYiU+JQ0KICB0YmxfcmVncmVzc2lvbigpDQoNCnN1bW1hcnkoZ2xtYikNCg0KYGBgDQpgYGB7cn0NCiANCmdsbWIlPiUNCiAgdGJsX3JlZ3Jlc3Npb24oZXhwb25lbnRpYXRlPVRSVUUpDQoNCmBgYA0KDQojIGIuIElzIGFuIG9mZnNldCB0ZXJtIG5lY2Vzc2FyeT8gd2h5IG9yIHdoeSBub3Q/IElmIGFuIG9mZnNldCB0ZXJtIGlzIG5lZWQsIHdoYXQgaXMgdGhlIGFwcHJvcHJpYXRlIG9mZnNldD8NCg0KDQpgYGB7cn0NCnN1bW1hcnkoYWhyZl9tJGluZmFudG1vcjE2MTgpDQpgYGANCg0KYGBge3J9DQphaHJmX20kaW5mYW50bW9yMTYxODwtIGFocmZfbSRpbmZhbnRtb3IxNjE4ICsxDQpgYGANCg0KIA0KIFllcyBhbiBvZmZzZXQgdGVybSBpcyBuZWNjZXNhcnkuIE9uZSBvZiB0aGUgY29uZGl0aW9ucyBvZiB1c2luZyB0aGUgcG9pc3NvbiAgbW9kZWwgaXMgdGhhdCBlYWNoIGFyZWEgb3IgcGVyc29uIGhhcyB0aGUgc2FtZSByaXNrIG9yIHBvcHVsYXRpb24gc2l6ZS4gU2luY2UgdGhpcyBpcyByZWFyIGluIGRlbW9ncmFwaHkgYW5kIG5vdCB0aGUgY2FzZSBmb3IgbXkgcG9wdWxhdGlvbiBvZiBzdHVkeSwgcG9pc3NvbiBzZWNvbmQgdHlwZSBtb2RlbGxpbmcoIHJhdGUgbW9kZWwpIHdpbGwgYmUgdXNlZC4NCiANCiBUaGUgb2Zmc2V0IHRlcm0gaXMgaG93IHVuZXF1YWwgcG9wdWxhdGlvbiBzaXplIGlzIGluY29ycG9yYXRlZC4gUG9zc2lvbiBpcyBhIHN0cmljdGx5IGludGVnZXIgY291bnQgc28gd2UgY2FuJ3QgdXNlIHRoZSByYXRlL3Byb2JhYmlsaXR5IHNvIHdlIGxvZyBpbnN0ZWFkLg0KIA0KIEFmdGVyIGRvaW5nIHRoZSBzdW1tYXJ5IG9mIHRoZSBkYXRhLCB0aGUgb3V0Y29tZSB2YXJpYWJsZSggaW5mYW50IG1vcnRhbGl0eSByYXQpIGNvbnRhaW4gemVyby4gVGhpcyB3b3VsZCBicmluZyBhbiBlcnJvciBpZiBkb25lIGRpcmVjdGx5LiBTbyBJIGFkZGVkIDEgdG8gdGhlIHZhcmlhYmxlKHBvcHVsYXRpb24pLiBJIHVzZWQgMSBiZWNhdXNlIHRoZXkgYXJlIGluZGl2aWR1YWxzIGFuZCB3ZSBjYW4ndCBoYXZlIDAuOSBvciAwLjEgaW5kaXZpZHVhbC4NCiANCiANCiMgMikgQ29uc2lkZXIgYSBQb2lzc29uIHJlZ3Jlc3Npb24gbW9kZWwgZm9yIHRoZSBvdXRjb21lDQogDQpgYGB7cn0NCmdsbXBfcyA8LSBnbG0oaW5mYW50bW9yMTYxOCB+IG9mZnNldChsb2coYmlydGhzMTYxOCkpICsgZmFtcG92MTQgKyBydWNjICsgaHBzYTE2LA0KICAgICAgICAgICAgICBkYXRhPWFocmZfbSwNCiAgICAgICAgICAgICAgZmFtaWx5PXBvaXNzb24pDQoNCmdsbXBfcyU+JQ0KICB0YmxfcmVncmVzc2lvbihleHAgPSBUUlVFKQ0KDQoNCmBgYA0KDQogIA0KICANCiMgYS4gRXZhbHVhdGUgdGhlIGxldmVsIG9mIGRpc3BlcnNpb24gaW4gdGhlIG91dGNvbWUNCg0KYGBge3J9DQpzY2FsZTwtc3FydChnbG1wX3MkZGV2aWFuY2UvZ2xtcF9zJGRmLnJlc2lkdWFsKQ0Kc2NhbGUNCg0KYGBgDQpIZXJlIHRoZSBzY2FsZSBmYWN0b3IgaXMgb25seSAxLjY1LCBzbyBpdCBkb2Vzbid0IHNob3cgc28gbXVjaCBvdmVyZGlzcGVyc2lvbiBhbmQgbm90IHRoYXQgYmFkLg0KDQogDQoNCmBgYHtyfQ0KMS1wY2hpc3EoZ2xtcF9zJGRldmlhbmNlLA0KICAgICAgICAgZGYgPSBnbG1wX3MkZGYucmVzaWR1YWwpDQoNCmBgYA0KDQojIGIuIElzIHRoZSBQb2lzc29uIG1vZGVsIGEgZ29vZCBjaG9pY2U/DQoNCk5vIGl0IGlzIG5vdCBhIGdvb2QgY2hvaWNlIHNpbmNlICB0aGUgZ29vZG5lc3Mgb2YgZml0IHJlc3VsdHMgdG8gemVyby4gU2luY2UgaXRzIHplcm8gaXQgbWVhbnMgdGhlIHBvaXNzb24gZGlzdHJpYnV0aW9uIGlzIGEgd3JvbmcgY2hvaWNlDQoNCiMgMykgQ29uc2lkZXIgYSBOZWdhdGl2ZSBiaW5vbWlhbCBtb2RlbA0KDQpgYGB7cn0NCmxpYnJhcnkoTUFTUykgDQogDQogDQpnbG1uYjwtIGdsbS5uYihpbmZhbnRtb3IxNjE4ICB+IG9mZnNldChsb2coYmlydGhzMTYxOCkpICsgZmFtcG92MTQgKyBydWNjICsgaHBzYTE2LA0KICAgICAgICAgICAgICBkYXRhPWFocmZfbSkNCg0KDQoNCg0KZ2xtbmIlPiUNCiAgdGJsX3JlZ3Jlc3Npb24oIGV4cD0gVCkNCg0KDQpzdW1tYXJ5KGdsbW5iKQ0KYGBgDQoNCg0KDQpgYGB7cn0NCmxpYnJhcnkoZ2FtbHNzKQ0KDQpnbG1uYjI8LWdhbWxzcyhpbmZhbnRtb3IxNjE4IH4gb2Zmc2V0KGxvZyhiaXJ0aHMxNjE4KSkgKyBmYW1wb3YxNCArIHJ1Y2MgKyBocHNhMTYsDQogICAgICAgICAgICAgICBmYW1pbHkgPSBOQklJLA0KICAgICAgICAgICAgICAgZGF0YT1haHJmX20pDQoNCnN1bW1hcnkoZ2xtbmIyKQ0KYGBgDQoNCiM0KSBDb21wYXJlIHRoZSBtb2RlbCBmaXRzIG9mIHRoZSBhbHRlcm5hdGl2ZSBtb2RlbHMgdXNpbmcgQUlDDQoNCkFJQyhnbG0xLCBnbG1iLCBnbG1wX2IsIGdsbW5iLCBnbG1uYjIpDQoNCmBgYHtyfQ0KDQpBSUMoZ2xtcF9zLCBnbG1iLGdsbW5iLCBnbG1uYjIpDQogDQoNCg0KYGBgDQoNCkFtb25nIHRoZXNlIGZvdXIgbW9kZWxzLCB0aGUgcG9pc3NvbiBtb2RlbCBmaXRzIHRoZSB3b3JzdCBhbmQgdGhlIE5CMSBtb2RlbCwgZXN0aW1hdGVkIGJ5IGdsbS5uYigpIGZpdHMgdGhlIGJlc3QsIHdpdGggdGhlIGxvd2VzdCBBSUMgYW1vbmcgdGhlc2UgZm91ciBtb2RlbHMuDQo=