Chapter 3 BY GEOGRAPHIC REGIONS

3.1 LIBRARIES

library(tidyverse); library(epiDisplay); library(lubridate)
library(spdep); library(sf); library(sp)
library(maptools); library(rgdal); library(INLA); 
library(MatrixModels); library(graph); library(Rgraphviz)
library(readxl)

3.2 ADMINISTRATIVE AND SPATIAL DATA

data <- read.csv("./_data/usable_data_binary_outcomes_2013-17_previous_year_corrected.csv", stringsAsFactors = F) %>%
  separate(County, c("County", "STATE"), sep = ", ")

usa.sf <- st_read("./_data/GIS/tl_2017_us_county/tl_2017_us_county.shp",  stringsAsFactors = F) %>%
  #filter(STATEFP=="01") %>%
  filter(STATEFP!="02", STATEFP!="15", STATEFP!="60", STATEFP!="66", STATEFP!="69", STATEFP!="72", STATEFP!="78") %>%
  mutate(County.Code = as.numeric(GEOID),
         STATEFP = as.numeric(STATEFP)) %>%
  left_join(read_excel("./_data/fp_codes.xlsx") %>%
              rename(STATEFP = `Numeric code`), by = "STATEFP")

## Reading layer `tl_2017_us_county' from data source `/research-home/gcarrasco/EMMERG_MAP2/_data/GIS/tl_2017_us_county/tl_2017_us_county.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3233 features and 17 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.43979
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs

3.3 GEGRAPHIC STRATA

# -----------------
# NorthEast
# -----------------
usa.N <- usa.sf %>%
  filter(loc == "N")

index.N <- usa.N %>%
  dplyr::select(County.Code) %>%
  st_set_geometry(NULL) %>%
  mutate(index = 1:nrow(usa.N),
         index2 = index)

dat.N <- data %>% 
  inner_join(index.N, by="County.Code") %>%
  mutate(fenta = as.numeric(any_fentanyl)) 

n.N<-nrow(dat.N)

nb.map.N <- poly2nb(usa.N)
nb2INLA("map.graph.N",nb.map.N)

# -----------------
# South
# -----------------
usa.S <- usa.sf %>%
  filter(loc == "S")

index.S <- usa.S %>%
  dplyr::select(County.Code) %>%
  st_set_geometry(NULL) %>%
  mutate(index = 1:nrow(usa.S),
         index2 = index)

dat.S <- data %>% 
  inner_join(index.S, by="County.Code") %>%
  mutate(fenta = as.numeric(any_fentanyl)) 

n.S<-nrow(dat.S)

nb.map.S <- poly2nb(usa.S)
nb2INLA("map.graph.S",nb.map.S)

# -----------------
# West
# -----------------
usa.W <- usa.sf %>%
  filter(loc == "W")

index.W <- usa.W %>%
  dplyr::select(County.Code) %>%
  st_set_geometry(NULL) %>%
  mutate(index = 1:nrow(usa.W),
         index2 = index)

dat.W <- data %>% 
  inner_join(index.W, by="County.Code") %>%
  mutate(fenta = as.numeric(any_fentanyl)) 

n.W<-nrow(dat.W)

nb.map.W <- poly2nb(usa.W)
nb2INLA("map.graph.W",nb.map.W)

# -----------------
# MidWest
# -----------------
usa.M <- usa.sf %>%
  filter(loc == "MW")

index.M <- usa.M %>%
  dplyr::select(County.Code) %>%
  st_set_geometry(NULL) %>%
  mutate(index = 1:nrow(usa.M),
         index2 = index)

dat.M <- data %>% 
  inner_join(index.M, by="County.Code") %>%
  mutate(fenta = as.numeric(any_fentanyl)) 

n.M<-nrow(dat.M)

nb.map.M <- poly2nb(usa.M)
nb2INLA("map.graph.M",nb.map.M)

3.4 MODEL SPECIFICATIONS

# FORMULAS

# -----------------
# NorthEast
# -----------------
N<-c(formula = fenta ~ 1,
     formula = fenta ~ 1 + f(index, model = "bym", graph = "map.graph.N") + f(Year, model = "rw1"))

# -----------------
# South
# -----------------
S<-c(formula = fenta ~ 1,
     formula = fenta ~ 1 + f(index, model = "bym", graph = "map.graph.S") + f(Year, model = "rw1"))

# -----------------
# West
# -----------------
W<-c(formula = fenta ~ 1,
     formula = fenta ~ 1 + f(index, model = "bym", graph = "map.graph.W") + f(Year, model = "rw1"))

# -----------------
# MidWest
# -----------------
M<-c(formula = fenta ~ 1,
     formula = fenta ~ 1 + f(index, model = "bym", graph = "map.graph.M") + f(Year, model = "rw1"))

3.5 Model Estimation

# -----------------
# NorthEast
# -----------------
names(N)<-N

INLA:::inla.dynload.workaround()
m.N <- N %>% purrr::map(~inla.batch.safe(formula = ., dat1 = dat.N))

N.s <- m.N %>%
  purrr::map(~Rsq.batch.safe(model = ., dic.null = m.N[[1]]$dic, n = n.N)) %>%
  bind_rows(.id = "formula") %>% mutate(id = row_number())

## [1] 0
## [1] 0.6336362

N.s %>% plot_score()

# -----------------
# South
# -----------------
names(S)<-S

INLA:::inla.dynload.workaround()
m.S <- S %>% purrr::map(~inla.batch.safe(formula = ., dat1 = dat.S))

S.s <- m.S %>%
  purrr::map(~Rsq.batch.safe(model = ., dic.null = m.S[[1]]$dic, n = n.S)) %>%
  bind_rows(.id = "formula") %>% mutate(id = row_number())

## [1] 0
## [1] 0.2958469

S.s %>% plot_score()

# -----------------
# West
# -----------------
names(W)<-W

INLA:::inla.dynload.workaround()
m.W <- W %>% purrr::map(~inla.batch.safe(formula = ., dat1 = dat.W))

W.s <- m.W %>%
  purrr::map(~Rsq.batch.safe(model = ., dic.null = m.W[[1]]$dic, n = n.W)) %>%
  bind_rows(.id = "formula") %>% mutate(id = row_number())

## [1] 0
## [1] 0.1556095

W.s %>% plot_score()

# -----------------
# MidWest
# -----------------
names(M)<-M

INLA:::inla.dynload.workaround()
m.M <- M %>% purrr::map(~inla.batch.safe(formula = ., dat1 = dat.M))

M.s <- m.M %>%
  purrr::map(~Rsq.batch.safe(model = ., dic.null = m.M[[1]]$dic, n = n.M)) %>%
  bind_rows(.id = "formula") %>% mutate(id = row_number())

## [1] 0
## [1] 0.2660458

M.s %>% plot_score()

3.6 Output

3.6.1 Year Random Effects

# -----------------
# NorthEast
# -----------------
m.N[[2]] %>% plot_random2(y_lab = "Odd Ratio")

summary(m.N[[2]])

## 
## Call:
##    c("inla(formula = formula, family = \"binomial\", data = dat1, verbose 
##    = F, ", " control.compute = list(config = T, dic = T, cpo = T, waic = 
##    T), ", " control.predictor = list(link = 1, compute = TRUE), 
##    control.fixed = list(correlation.matrix = T))" ) 
## Time used:
##     Pre = 1.57, Running = 6.87, Post = 4.75, Total = 13.2 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -5.595 0.713     -7.132   -5.522     -4.414 -5.485   0
## 
## Linear combinations (derived):
##             ID   mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept)  1 -5.595 0.713      -7.12   -5.539       -4.4 -5.57   0
## 
## Random effects:
##   Name     Model
##     index BYM model
##    Year RW1 model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for index (iid component)     1863.210 1836.900    126.232 1321.073
## Precision for index (spatial component)    0.017    0.005      0.010    0.017
## Precision for Year                         0.114    0.078      0.023    0.096
##                                         0.975quant    mode
## Precision for index (iid component)       6721.569 344.814
## Precision for index (spatial component)      0.029   0.016
## Precision for Year                           0.314   0.060
## 
## Expected number of effective parameters(stdev): 161.05(5.42)
## Number of equivalent replicates : 8.08 
## 
## Deviance Information Criterion (DIC) ...............: Inf
## Deviance Information Criterion (DIC, saturated) ....: Inf
## Effective number of parameters .....................: Inf
## 
## Watanabe-Akaike information criterion (WAIC) ...: 443.50
## Effective number of parameters .................: 85.78
## 
## Marginal log-Likelihood:  -323.28 
## CPO and PIT are computed
## 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

# -----------------
# South
# -----------------
m.S[[2]] %>% plot_random2(y_lab = "Odd Ratio")

summary(m.S[[2]])

## 
## Call:
##    c("inla(formula = formula, family = \"binomial\", data = dat1, verbose 
##    = F, ", " control.compute = list(config = T, dic = T, cpo = T, waic = 
##    T), ", " control.predictor = list(link = 1, compute = TRUE), 
##    control.fixed = list(correlation.matrix = T))" ) 
## Time used:
##     Pre = 2.09, Running = 47.9, Post = 29, Total = 79.1 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -14.379 1.269    -17.047  -14.159    -12.294 -13.938   0
## 
## Linear combinations (derived):
##             ID    mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  1 -14.382 1.269    -17.045  -14.173     -12.29 -13.982   0
## 
## Random effects:
##   Name     Model
##     index BYM model
##    Year RW1 model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for index (iid component)     1800.285 1796.984    121.839 1268.088
## Precision for index (spatial component)    0.013    0.002      0.009    0.013
## Precision for Year                         0.187    0.127      0.036    0.158
##                                         0.975quant    mode
## Precision for index (iid component)       6571.712 332.931
## Precision for index (spatial component)      0.018   0.013
## Precision for Year                           0.513   0.097
## 
## Expected number of effective parameters(stdev): 470.88(9.64)
## Number of equivalent replicates : 18.12 
## 
## Deviance Information Criterion (DIC) ...............: Inf
## Deviance Information Criterion (DIC, saturated) ....: Inf
## Effective number of parameters .....................: Inf
## 
## Watanabe-Akaike information criterion (WAIC) ...: 950.71
## Effective number of parameters .................: 173.96
## 
## Marginal log-Likelihood:  -624.83 
## CPO and PIT are computed
## 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

# -----------------
# West
# -----------------
m.W[[2]] %>% plot_random2(y_lab = "Odd Ratio")

summary(m.W[[2]])

## 
## Call:
##    c("inla(formula = formula, family = \"binomial\", data = dat1, verbose 
##    = F, ", " control.compute = list(config = T, dic = T, cpo = T, waic = 
##    T), ", " control.predictor = list(link = 1, compute = TRUE), 
##    control.fixed = list(correlation.matrix = T))" ) 
## Time used:
##     Pre = 1.97, Running = 12.8, Post = 12.1, Total = 26.8 
## Fixed effects:
##              mean   sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -9.42 1.35    -12.446    -9.23     -7.284 -9.096   0
## 
## Linear combinations (derived):
##             ID   mean   sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)  1 -9.422 1.35    -12.455    -9.25     -7.275 -9.187   0
## 
## Random effects:
##   Name     Model
##     index BYM model
##    Year RW1 model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for index (iid component)     1326.825 1170.834     72.581  988.791
## Precision for index (spatial component)    0.032    0.012      0.015    0.030
## Precision for Year                         1.081    0.881      0.174    0.845
##                                         0.975quant    mode
## Precision for index (iid component)        4290.94 182.735
## Precision for index (spatial component)       0.06   0.027
## Precision for Year                            3.42   0.466
## 
## Expected number of effective parameters(stdev): 99.01(6.79)
## Number of equivalent replicates : 25.09 
## 
## Deviance Information Criterion (DIC) ...............: Inf
## Deviance Information Criterion (DIC, saturated) ....: Inf
## Effective number of parameters .....................: Inf
## 
## Watanabe-Akaike information criterion (WAIC) ...: 267.75
## Effective number of parameters .................: 38.99
## 
## Marginal log-Likelihood:  -145.37 
## CPO and PIT are computed
## 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

# -----------------
# MidWest
# -----------------
m.M[[2]] %>% plot_random2(y_lab = "Odd Ratio")

summary(m.M[[2]])

## 
## Call:
##    c("inla(formula = formula, family = \"binomial\", data = dat1, verbose 
##    = F, ", " control.compute = list(config = T, dic = T, cpo = T, waic = 
##    T), ", " control.predictor = list(link = 1, compute = TRUE), 
##    control.fixed = list(correlation.matrix = T))" ) 
## Time used:
##     Pre = 2.55, Running = 34.2, Post = 20.7, Total = 57.4 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -12.445 1.156    -14.865  -12.335    -10.493 -12.169   0
## 
## Linear combinations (derived):
##             ID    mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  1 -12.445 1.156    -14.865  -12.346     -10.48 -12.245   0
## 
## Random effects:
##   Name     Model
##     index BYM model
##    Year RW1 model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for index (iid component)     1899.686 1859.968    127.271 1351.493
## Precision for index (spatial component)    0.022    0.004      0.014    0.021
## Precision for Year                         0.287    0.195      0.060    0.240
##                                         0.975quant    mode
## Precision for index (iid component)        6.8e+03 346.524
## Precision for index (spatial component)    3.1e-02   0.020
## Precision for Year                         7.9e-01   0.155
## 
## Expected number of effective parameters(stdev): 273.83(8.98)
## Number of equivalent replicates : 23.09 
## 
## Deviance Information Criterion (DIC) ...............: Inf
## Deviance Information Criterion (DIC, saturated) ....: Inf
## Effective number of parameters .....................: Inf
## 
## Watanabe-Akaike information criterion (WAIC) ...: 696.84
## Effective number of parameters .................: 115.55
## 
## Marginal log-Likelihood:  -414.44 
## CPO and PIT are computed
## 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

3.6.2 Spatial Random Effects

# -----------------
# NorthEast
# -----------------
m.N[[2]] %>% plot_random_map2(map = usa.N, id = County.Code, col1 = NA)

# -----------------
# South
# -----------------
m.S[[2]] %>% plot_random_map2(map = usa.S, id = County.Code, col1 = NA)

# -----------------
# West
# -----------------
m.W[[2]] %>% plot_random_map2(map = usa.W, id = County.Code, col1 = NA)

# -----------------
# MidWest
# -----------------
m.M[[2]] %>% plot_random_map2(map = usa.M, id = County.Code, col1 = NA)