Spatial Demography - Homewk 6 -Spatial GLMM(s) using the INLA Approximation

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) %>% 
  mutate(Expdd=nbir*(sum(ndea)/sum(nbir))) %>% 
  subset(!state%in%c("02", "15", "60", "66", "69", "72", "78"))

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.899607279, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.2118  -0.5258  -0.0194   0.5125   7.4889  
## 
## Coefficients:
##                   Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)       -7.58730    0.04983 -152.272  < 2e-16 ***
## ruralityNon-Metro  0.26511    0.06884    3.851 0.000118 ***
## medhouvlz         -0.17240    0.03194   -5.397 6.78e-08 ***
## pctcrpopz          0.58867    0.03533   16.663  < 2e-16 ***
## ctypopdesz         0.37666    0.01060   35.544  < 2e-16 ***
## pctnhblackpr19z    0.11650    0.02286    5.095 3.49e-07 ***
## medhincz          -0.06808    0.03387   -2.010 0.044417 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(5.8996) family taken to be 1)
## 
##     Null deviance: 2283.94  on 1242  degrees of freedom
## Residual deviance:  883.55  on 1236  degrees of freedom
##   (29 observations deleted due to missingness)
## AIC: 4540.5
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  5.900 
##           Std. Err.:  0.498 
## Warning while fitting theta: alternation limit reached 
## 
##  2 x log-likelihood:  -4524.466

##       (Intercept) ruralityNon-Metro         medhouvlz         pctcrpopz 
##      0.0005068469      1.3035721968      0.8416453538      1.8015897782 
##        ctypopdesz   pctnhblackpr19z          medhincz 
##      1.4574054974      1.1235551354      0.9341852600

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.845486

# 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.899607279, link = log)
## weights: us.wt2
## 
## Moran I statistic standard deviate = 0.64628, p-value = 0.259
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.0161639772    -0.0005525404     0.0006690314

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.

Construct Spatial Relationship for INLA

mmrinla <- usctymmr  %>% 
  arrange(cofips) %>% 
  rowid_to_column("id")
#constructing a K neighborhood
nbs<-knearneigh(st_centroid(mmrinla), k = 5, longlat = T)  
nbs<-knn2nb(nbs, row.names = mmrinla$id, sym = T)
nb2INLA("cl_graph",nbs)
am_adj <-paste(getwd(),"/cl_graph",sep="")
H<-inla.read.graph(filename="cl_graph")

Model Estimation without Sptial Pattern

f1<-ndea~rurality+ medhouvlz+ pctcrpopz + ctypopdesz + pctnhblackpr19z+medhincz+ #fixed effects
  f(id, model = "iid")  #random effects

mod2<-inla(formula = f1,
           data = mmrinla,
           family = "poisson",
           E = Expdd, 
           control.compute = list(waic=T), 
           control.predictor = list(link=1),
           num.threads = 3, 
               verbose = F)

#total model summary
summary(mod2)

## 
## Call:
##    c("inla(formula = f1, family = \"poisson\", data = mmrinla, E = Expdd, 
##    ", " verbose = F, control.compute = list(waic = T), control.predictor = 
##    list(link = 1), ", " num.threads = 3)") 
## Time used:
##     Pre = 4.03, Running = 7.94, Post = 0.139, Total = 12.1 
## Fixed effects:
##                     mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)        0.134 0.048      0.039    0.135      0.228  0.135   0
## ruralityNon-Metro -0.033 0.066     -0.162   -0.033      0.096 -0.033   0
## medhouvlz         -0.150 0.034     -0.218   -0.150     -0.084 -0.150   0
## pctcrpopz          0.188 0.034      0.121    0.188      0.254  0.189   0
## ctypopdesz         0.284 0.013      0.258    0.284      0.309  0.284   0
## pctnhblackpr19z    0.169 0.023      0.123    0.169      0.214  0.169   0
## medhincz          -0.121 0.036     -0.192   -0.121     -0.051 -0.121   0
## 
## Random effects:
##   Name     Model
##     id IID model
## 
## Model hyperparameters:
##                  mean    sd 0.025quant 0.5quant 0.975quant mode
## Precision for id 3.59 0.309       3.04     3.57       4.25 3.53
## 
## Watanabe-Akaike information criterion (WAIC) ...: 6270.16
## Effective number of parameters .................: 398.25
## 
## Marginal log-Likelihood:  -3263.01 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Marginal Distributions of hyperparameters

m2<- inla.tmarginal(
        function(x) (1/x), #invert the precision to be on variance scale
        mod2$marginals.hyperpar$`Precision for id`)
#95% credible interval for the variance
inla.hpdmarginal(.95, marginal=m2)

##                  low      high
## level:0.95 0.2333311 0.3266912

plot(m2,
     type="l",
     main=c("Posterior distibution for between county variance", "- IID model -"), 
     xlim=c(0, .5))

map1

f3<-ndea~rurality+ medhouvlz+ pctcrpopz + ctypopdesz + pctnhblackpr19z+medhincz+
  f(id, model = "bym",scale.model = T, constr = T, graph = H)

mod3<-inla(formula = f3,data = mmrinla,
           family = "poisson",
           E = Expdd, 
           control.compute = list(waic=T,return.marginals.predictor=TRUE), 
           control.predictor = list(link=1), 
           num.threads = 3, 
               verbose = F)

#total model summary
summary(mod3)

## 
## Call:
##    c("inla(formula = f3, family = \"poisson\", data = mmrinla, E = Expdd, 
##    ", " verbose = F, control.compute = list(waic = T, 
##    return.marginals.predictor = TRUE), ", " control.predictor = list(link 
##    = 1), num.threads = 3)") 
## Time used:
##     Pre = 4.01, Running = 19.5, Post = 0.37, Total = 23.8 
## Fixed effects:
##                     mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)        0.105 0.051      0.004    0.105      0.203  0.106   0
## ruralityNon-Metro -0.009 0.067     -0.139   -0.009      0.122 -0.009   0
## medhouvlz         -0.089 0.047     -0.181   -0.089      0.003 -0.089   0
## pctcrpopz          0.178 0.037      0.106    0.179      0.250  0.179   0
## ctypopdesz         0.266 0.015      0.236    0.266      0.294  0.266   0
## pctnhblackpr19z    0.113 0.033      0.048    0.113      0.177  0.113   0
## medhincz          -0.129 0.044     -0.216   -0.129     -0.043 -0.129   0
## 
## Random effects:
##   Name     Model
##     id BYM model
## 
## Model hyperparameters:
##                                      mean    sd 0.025quant 0.5quant 0.975quant
## Precision for id (iid component)     5.30 0.651       4.14     5.26       6.70
## Precision for id (spatial component) 8.25 2.386       4.70     7.86      13.99
##                                      mode
## Precision for id (iid component)     5.17
## Precision for id (spatial component) 7.14
## 
## Watanabe-Akaike information criterion (WAIC) ...: 6184.10
## Effective number of parameters .................: 371.87
## 
## Marginal log-Likelihood:  -2106.61 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

mod2$waic$waic

## [1] 6270.163

mod3$waic$waic

## [1] 6184.104

Regression effect posterior distributions

b1 <- inla.emarginal(exp, mod2$marginals.fixed[[2]])
b1 #Posterior mean

## [1] 0.9699767

b1ci <-inla.qmarginal(c(.025, .975), 
                      inla.tmarginal(exp, mod2$marginals.fixed[[2]]))
b1ci

## [1] 0.8512323 1.1001186

Hyper parameter distributions

m3a<- inla.tmarginal(
        function(x) (1/x),
        mod3$marginals.hyperpar$`Precision for id (iid component)`)
m3b<- inla.tmarginal(
        function(x) (1/x),
        mod3$marginals.hyperpar$`Precision for id (spatial component)`)

plot(m3a, type="l",
     main=c("Posterior distibution for between county variance", "- BYM model -"),
     xlim=c(0, .2)) # you may need to change this
#lines(m3b, col="red")
lines(m3b, col="green")

inla.hpdmarginal(.95,m3a)

##                  low      high
## level:0.95 0.1472058 0.2381321

inla.hpdmarginal(.95,m3b)

##                   low      high
## level:0.95 0.06573046 0.2019664

Spatial mapping of the fitted values

map2

Map of spatial random effects

map3

map4