Spatial Demography - Homewk 4 -Generalized Linear Models for Spatial Count data

Data Manipulations

# Get US county socio-economic variables from Area resource file 2019-2020

arf2020<-import("ahrf2020.sas7bdat") 
 arf2020<-arf2020%>%
  mutate(cofips=as.factor(f00004), 
         coname=f00010,
         state = f00011,
         medhouvl=f1461314,
         pctcrpop= round(100*(f1492010/f0453010),2),
         medhinc= f1434514,
         pctctyunemp=round(100*(f1451214/f1451014),2),
         ctypopdes= (f0453010/f1387410),  
         rucc= as.numeric(f0002013),
         pctperpov= f1332118,
         obgyn10_pc= 1000*(f1168410/ f0453010) )%>%
  dplyr::select(medhouvl, pctcrpop,medhinc, pctctyunemp,state, cofips, coname, ctypopdes, rucc, obgyn10_pc,pctperpov)%>%
    mutate(rurality= car::Recode(rucc, recodes ="1:3 ='Metro';4:9='Non-Metro'")) 


# Get other  US county socioeconomic variables from the 5 year (2015-2019) ASC data source.

 US_county <- get_acs(geography = 'county',variables = c(totPop19 = "B03003_001", 
                                                   hispanic ="B03003_003", 
                                                   afrAm = "B03002_004",
                                                   white="B03002_003",
                                                   gini="B19083_001"), 
                    year = 2019, geometry = TRUE, shift_geo = TRUE) %>% 
            dplyr::select(GEOID, NAME, variable, estimate) %>% 
            spread(variable, estimate) %>% 
            mutate(pcthispPr19  = round((hispanic/totPop19)*100),1, 
                   pctnhblackpr19 = round((afrAm/totPop19)*100),1 , pctnhpwhite19=round((white/totPop19)*100),1, cofips=as.factor(GEOID)) %>%
            dplyr::select(cofips,totPop19,pcthispPr19, pctnhblackpr19, pctnhpwhite19,gini)


#Merge the socioeconomic variables from the Area resource file 2019-2020 with the other socioeconomic variables from the 5 year (2015-2019) ACS estimates.

 sevar <- geo_join( US_county, arf2020, by_sp="cofips", by_df="cofips",how="left") 

 #load and calculate maternal mortality rates for counties in US for  a 5year period (2015-2019)
alldat<-readRDS("~/OneDrive - University of Texas at San Antonio/maternalmortality/aggregate/alldat.rds") 
 countyMMR <-  alldat %>% 
  group_by(cofips)%>%
  summarise(nbir=sum(nbirths, na.rm = T), ndea = sum(ndeaths, na.rm = T))  %>% 
  mutate(mmrate =round(100000*(ndea/nbir),2))
 
  #Merge the socioeconomic variables with US counties maternal mortality rates for a 5- years period(2015-2019)
MMRdata<- geo_join( sevar, countyMMR, by_sp="cofips", by_df="cofips",how="left") %>% 
  mutate(rmmrate=ifelse(mmrate>=1,mmrate,NA)) %>% 
  filter(!is.na(ndea), mmrate!=Inf) 

Descriptive Maps

Defining Variables

The denominator/ offset variable for this analysis is Total number of births in US county. Five socioeconomic variables were used as predictors in the models–percentage of county population that is rural;the median household income in the county; the county’s median house value; the population density per square mile of the county; and ercentage of non Hispanic black population in the counties. All variables were Z scored

Estimating Generalized Linear Model

## 
## Call:
## glm.nb(formula = ndea ~ offset(log(nbir)) + rurality + medhouvlz + 
##     pctcrpopz + ctypopdesz + pctnhblackpr19z + medhincz, data = usctymmrr, 
##     init.theta = 5.874759265, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.1469  -0.5259  -0.0197   0.5194   7.5038  
## 
## Coefficients:
##                   Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)       -7.59014    0.04983 -152.335  < 2e-16 ***
## ruralityNon-Metro  0.26739    0.06872    3.891 9.97e-05 ***
## medhouvlz         -0.15762    0.03133   -5.030 4.89e-07 ***
## pctcrpopz          0.58810    0.03537   16.629  < 2e-16 ***
## ctypopdesz         0.37213    0.01049   35.473  < 2e-16 ***
## pctnhblackpr19z    0.11609    0.02286    5.078 3.82e-07 ***
## medhincz          -0.07698    0.03367   -2.286   0.0222 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(5.8748) family taken to be 1)
## 
##     Null deviance: 2286.23  on 1245  degrees of freedom
## Residual deviance:  889.32  on 1239  degrees of freedom
##   (30 observations deleted due to missingness)
## AIC: 4557.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  5.875 
##           Std. Err.:  0.496 
## Warning while fitting theta: alternation limit reached 
## 
##  2 x log-likelihood:  -4541.753

##       (Intercept) ruralityNon-Metro         medhouvlz         pctcrpopz 
##      0.0005054108      1.3065501892      0.8541755812      1.8005727263 
##        ctypopdesz   pctnhblackpr19z          medhincz 
##      1.4508233576      1.1230981105      0.9259127913

The model above shows a 30% increase in relative risk for maternal death in non metropolitan counties compared to metropolitan counties. Also, for every one unit increase in county’s median house value and median household income, maternal death decrease by 15% and 7.5% respectively. One unit increase in county’s percentage of rural population, County’s population density, and percentage of Non-Hispanic black increase maternal mortality significantly.

Testing for Overdispersion

# Testing for overdispersion  using the ratio of residual deviance and residual degree of freedom from the poisson model. 
scale<-sqrt(fit_nb$deviance/fit_nb$df.residual)
scale

## [1] 0.8472132

# Testing for overdispersion using the goodness of fit statistic.
1-pchisq(fit_nb$deviance, df = fit_nb$df.residual)

## [1] 1

The overdispersion test shows that the model fits the data properly. The scale test is approximately equals to 1 showing that the mean and the variance of the model is equal to 1. Likewise the goodness of fit statistic gave a large p value of 1. Hence, the model fit the data.

Testing model for Residual Autocorrelation

nbs<-knearneigh(coordinates(as(usctymmrr, "Spatial")), k = 2, longlat = T)
nbs<-knn2nb(nbs, sym = T)
us.wt2<-nb2listw(nbs, style = "W")
lm.morantest(fit_nb, listw = us.wt2)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: glm.nb(formula = ndea ~ offset(log(nbir)) + rurality + medhouvlz
## + pctcrpopz + ctypopdesz + pctnhblackpr19z + medhincz, data =
## usctymmrr, init.theta = 5.874759265, link = log)
## weights: us.wt2
## 
## Moran I statistic standard deviate = 0.55844, p-value = 0.2883
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.0136988488    -0.0006232626     0.0006577452

Given the result from the above model; testing for autocorrelation, it appears there is no autocorrelation in our data. The global Moran I is 0.013 with a p-value that is larger than 0.05. At this point, I am not sure on how to move forward with the rest of the assignment.