multilevel_sp.Rmd

---
output:
  pdf_document: default
---
# Spatial Multilevel Modeling of Life Expectancy

```{r message=FALSE, warning=FALSE, echo=FALSE}
require('tidycensus')
require('readstata13')
require('stringr')
require('sf')
require('dplyr')
require('spdep')
require('tibble')
require('brms')
```


## Set modeling parameters

```{r}
runs <- 20000
burnin <- 5000
thinning <- 10
cores <- 4
chains <- 4
```

## Define Moran Cluster Function

```{r}
moran_cluster <- function(data, model_residuals, listw, p_min = 0.05) {
  c <- as_tibble(localmoran(model_residuals, listw)) %>%
    bind_cols(data, .) %>%
    add_column(residuals = model_residuals, .) %>%
    rename(
      prob = 'Pr(z > 0)'
    ) %>%
    dplyr::select(-c('Ii','Z.Ii', 'Var.Ii', "E.Ii")) %>%
    mutate(
      lag_residuals = lag.listw(listw, residuals),
      type = case_when(
        residuals >= mean(residuals) & lag_residuals >= mean(lag_residuals) & prob <= p_min ~ "hh",
        residuals < mean(residuals) & lag_residuals < mean(lag_residuals) & prob <= p_min ~ "ll",
        residuals >= mean(residuals) & lag_residuals < mean(lag_residuals) & prob <= p_min ~ "hl",
        residuals < mean(residuals) & lag_residuals >= mean(lag_residuals) & prob <= p_min ~ "lh"
      ), 
      binary = case_when(
        residuals > 0 ~ "pos",
        residuals < 0 ~ "neg",
        residuals == 0 ~ "zero",
      )
    )
  return(c)
}
```

## Read Data

```{r}
options(tigris_use_cache = TRUE)

input <- read.csv('https://ftp.cdc.gov/pub/Health_Statistics/NCHS/Datasets/NVSS/USALEEP/CSV/US_A.CSV') %>%
  mutate(
    geoid=str_pad(Tract.ID, width=11, pad="0"),
    stateid=as.factor(str_pad(STATE2KX, width=2, pad="0")),
    countyid=as.factor(paste0(stateid,str_pad(CNTY2KX, width=3, pad="0"))),
    tractid=as.factor(paste0(countyid,str_pad(TRACT2KX, width=6, pad="0")))
  ) %>%
  rename(
    ex = e.0.
  ) %>%
  dplyr::select(c(geoid, stateid, countyid, tractid, ex))

# States except Wisconsin and Maine
states <- state.abb[state.abb != "WI" & state.abb != "ME"]
geom <- get_decennial(geography = "tract",
                           variables = c("P005003"),
                           # state = state.name,
                           state = states,
                           geometry = TRUE,
                           year = 2010
                           ) %>%
  rename_all(tolower) %>%
  # Conus Albers
  st_transform(5070) %>%
  dplyr::select(c(geoid, geometry)) %>%
  filter(
    !st_is_empty(geometry)
  )

merged <- merge(geom, input, 'geoid')
```

## Construct Neighbors List and Spatial Weights

```{r}
# Build neighbors list using queen contiguity condition.
nb <- poly2nb(merged, queen = TRUE)

# Remove tracts with zero cardinality (i.e., no neighbors).
merged <- subset(merged, card(nb) > 0)
nb <- subset(nb, card(nb) > 0)

# Generate spatial weights.
w_list <- nb2listw(nb, style="B")
w_mat <- nb2mat(nb, style="B")
```

## Model 1: Single-level GLM modeling.

```{r}
glm_res <- brm(
  ex ~ 1,
  data = merged,
  family = gaussian(),
  iter = runs,
  warmup = burnin,
  thin = thinning,
  chains = chains,
  cores = cores
)

# Calculate residuals.
glm_resid <- resid(glm_res)[,1]
# Calculate Widely Applicable Information Criterion.
waic(glm_res)

moran.test(glm_resid, listw = w_list)
moran.plot(glm_resid, listw = w_list)

glm_lisa <- moran_cluster(merged, glm_resid, w_list)
plot(glm_lisa[c('geometry', 'binary')], lwd=0.05)
plot(glm_lisa[c('geometry', 'type')], lwd=0.05)
```

## Model 2a: Model accounting for state hierarchy.

```{r}
glmm_st_res <- brm(
  ex ~ 1 + (1|stateid),
  data = merged,
  family = gaussian(),
  iter = runs,
  warmup = burnin,
  thin = thinning,
  chains = chains,
  cores = cores
)

glmm_st_resid <- resid(glmm_st_res)[,1]
waic(glmm_st_res)

summary(glmm_st_res)

moran.test(glmm_st_resid, listw = w_list)
moran.plot(glmm_st_resid, listw = w_list)

glmm_st_lisa <- moran_cluster(merged, glmm_st_resid, w_list)
plot(glmm_st_lisa[c('geometry', 'binary')], lwd=0.05)
plot(glmm_st_lisa[c('geometry', 'type')], lwd=0.05)
```

## Model 2b: Model accounting for county hierarchy.

```{r}
glmm_cnty_res <- brm(
  ex ~ 1 + (1|countyid),
  data = merged,
  family = gaussian(),
  iter = runs,
  warmup = burnin,
  thin = thinning,
  chains = chains,
  cores = cores
)

glmm_cnty_resid <- resid(glmm_cnty_res)[,1]
waic(glmm_cnty_res)

summary(glmm_cnty_res)

moran.test(glmm_cnty_resid, listw = w_list)
moran.plot(glmm_cnty_resid, listw = w_list)

glmm_cnty_lisa <- moran_cluster(merged, glmm_cnty_resid, w_list)
plot(glmm_cnty_lisa[c('geometry', 'binary')], lwd=0.05)
plot(glmm_cnty_lisa[c('geometry', 'type')], lwd=0.05)
```

## Model 2c: Model accounting for nested state and county hierarchy.

```{r}
glmm_nest_res <- brm(
  ex ~ 1 + (1|stateid/countyid),
  data = merged,
  family = gaussian(),
  iter = runs,
  warmup = burnin,
  thin = thinning,
  chains = chains,
  cores = cores
)

glmm_nest_resid <- resid(glmm_nest_res)[,1]
waic(glmm_nest_res)

summary(glmm_nest_res)

moran.test(glmm_nest_resid, listw = w_list)
moran.plot(glmm_nest_resid, listw = w_list)

glmm_nest_lisa <- moran_cluster(merged, glmm_nest_resid, w_list)
plot(glmm_nest_lisa[c('geometry', 'binary')], lwd=0.05)
plot(glmm_nest_lisa[c('geometry', 'type')], lwd=0.05)

```

# Model 3: Two-level model accounting for spatial clustering.

```{r}
sp_res <- brm(
  ex ~ 1,
  data = merged,
  family = gaussian(),
  iter = runs,
  warmup = burnin,
  thin = thinning,
  chains = chains,
  cores = cores,
  autocor = cor_car(w_mat)
)

sp_resid <- resid(sp_res)[,1]
waic(sp_res)

moran.test(sp_resid, listw = w_list)
moran.plot(sp_resid, listw = w_list)

sp_lisa <- moran_cluster(merged, sp_resid, w_list, p_min = 0.1)
plot(sp_lisa[c('geometry', 'binary')], lwd=0.05)
plot(sp_lisa[c('geometry', 'type')], lwd=0.05)
```

# Model 4a: Model accounting for state hierarchy and space.

```{r}
sp_st_res <- brm(
  ex ~ 1 + (1|stateid),
  data = merged,
  family = gaussian(),
  iter = runs,
  warmup = burnin,
  thin = thinning,
  chains = chains,
  cores = cores,
  autocor = cor_car(w_mat)
)

sp_st_resid <- resid(sp_st_res)[,1]
waic(sp_st_res)

moran.test(sp_st_resid, listw = w_list)
moran.plot(sp_st_resid, listw = w_list)

sp_st_lisa <- moran_cluster(merged, sp_st_resid, w_list)
plot(sp_st_lisa[c('geometry', 'binary')], lwd=0.05)
plot(sp_st_lisa[c('geometry', 'type')], lwd=0.05)

```

# Model 4b: Model accounting for county hierarchy and space.

```{r}
sp_cnty_res <- brm(
  ex ~ 1 + (1|countyid),
  data = merged,
  family = gaussian(),
  iter = runs,
  warmup = burnin,
  thin = thinning,
  chains = chains,
  cores = cores,
  autocor = cor_car(w_mat)
)

sp_cnty_resid <- resid(sp_cnty_res)[,1]
waic(sp_cnty_res)

moran.test(sp_cnty_resid, listw = w_list)
moran.plot(sp_cnty_resid, listw = w_list)

sp_cnty_lisa <- moran_cluster(merged, sp_cnty_resid, w_list)
plot(sp_cnty_lisa[c('geometry', 'binary')], lwd=0.05)
plot(sp_cnty_lisa[c('geometry', 'type')], lwd=0.05)
```

# Model 4c: Model accounting for state-county hierarchy and space.

```{r}
sp_nest_res <- brm(
  ex ~ 1 + (1|stateid/countyid),
  data = merged,
  family = gaussian(),
  iter = runs,
  warmup = burnin,
  thin = thinning,
  chains = chains,
  cores = cores,
  autocor = cor_car(w_mat)
)

sp_nest_resid <- resid(sp_nest_res)[,1]
waic(sp_nest_res)

moran.test(sp_nest_resid, listw = w_list)
moran.plot(sp_nest_resid, listw = w_list)

sp_nest_lisa <- moran_cluster(merged, sp_nest_resid, w_list)
plot(sp_nest_lisa[c('geometry', 'binary')], lwd=0.05)
plot(sp_nest_lisa[c('geometry', 'type')], lwd=0.05)
```