Rural data exploration

There are county typology codes dating back to 1974 available here. They track the dominant economic characteristics of a county as well as other indicators such as whether there’s persistent poverty. The list of variables is below.

source("C:/Users/bean/Documents/GitHub/R-practice/my_functions.R")
librarian()

## [1] "Loaded Basic Libraries"

types <- read_csv("2015CountyTypologyCodes.csv") %>%
  clean_names()

glimpse(types)

## Rows: 3,143
## Columns: 18
## $ fip_stxt                                        <dbl> 1001, 1003, 1005, 1...
## $ state                                           <chr> "AL", "AL", "AL", "...
## $ county_name                                     <chr> "Autauga County", "...
## $ metro_nonmetro_status_2013_0_nonmetro_1_metro   <dbl> 1, 1, 0, 1, 1, 0, 0...
## $ economic_types_type_2015_update_non_overlapping <dbl> 0, 5, 3, 0, 0, 3, 0...
## $ economic_type_label                             <chr> "Nonspecialized", "...
## $ farming_2015_update                             <dbl> 0, 0, 0, 0, 0, 0, 0...
## $ mining_2015_update                              <dbl> 0, 0, 0, 0, 0, 0, 0...
## $ manufacturing_2015_update                       <dbl> 0, 0, 1, 0, 0, 1, 0...
## $ government_2015_update                          <dbl> 0, 0, 0, 0, 0, 1, 0...
## $ recreation_2015_update                          <dbl> 0, 1, 0, 0, 0, 0, 0...
## $ nonspecialized_2015_update                      <dbl> 1, 0, 0, 1, 1, 0, 1...
## $ low_education_2015_update                       <dbl> 0, 0, 1, 1, 1, 1, 0...
## $ low_employment_cnty_2008_2012_25_64             <dbl> 0, 0, 1, 1, 1, 1, 1...
## $ pop_loss_2010                                   <dbl> 0, 0, 0, 0, 0, 0, 1...
## $ retirement_dest_2015_update                     <dbl> 1, 1, 0, 0, 0, 0, 0...
## $ persistent_poverty_2013                         <dbl> 0, 0, 1, 0, 0, 1, 1...
## $ persistent_related_child_poverty_2013           <dbl> 0, 0, 1, 1, 0, 1, 1...

First we can look at the some unique policy attributes of each county.

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(albersusa)
library(leaflet)
library(leaflet.extras)

cty_sf <- counties_sf("aeqd") %>% mutate(fip_stxt= as.numeric(as.character(fips)))

epsg2163 <- leafletCRS(
  crsClass = "L.Proj.CRS",
  code = "EPSG:2163",
  proj4def = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs",
  resolutions = 2^(16:7))

map_data <- cty_sf %>%
  inner_join(types, by="fip_stxt") %>% 
  rename_with(~str_remove_all(.x, "_2015_update")) %>%
  rename_with(~str_remove_all(.x, "_2013")) %>%
  rename_with(~str_remove_all(.x, "_2010")) %>%
  rename_with(~str_remove_all(.x, "_cnty_2008_2012_25_64")) 

#low education, low employment, population loss, retirement destination, persistent poverty, persistent child-related poverty
#i'm interested to see the correlation between these scores
map_data %>%
  select(low_education:persistent_related_child_poverty) %>% 
  st_drop_geometry() %>%
  cor_matrix()

##                                  low_education low_employment pop_loss
## low_education                             1.00           0.37    -0.04
## low_employment                            0.37           1.00    -0.06
## pop_loss                                 -0.04          -0.06     1.00
## retirement_dest                          -0.07           0.05    -0.17
## persistent_poverty                        0.40           0.41     0.04
## persistent_related_child_poverty          0.40           0.49     0.02
##                                  retirement_dest persistent_poverty
## low_education                              -0.07               0.40
## low_employment                              0.05               0.41
## pop_loss                                   -0.17               0.04
## retirement_dest                             1.00              -0.07
## persistent_poverty                         -0.07               1.00
## persistent_related_child_poverty           -0.05               0.64
##                                  persistent_related_child_poverty
## low_education                                                0.40
## low_employment                                               0.49
## pop_loss                                                     0.02
## retirement_dest                                             -0.05
## persistent_poverty                                           0.64
## persistent_related_child_poverty                             1.00

#create bins for different poverty measures
all_labels= c("No", "Yes")

edu_bins<-c(0, 1, 2)
pal_edu <- colorBin(palette = "Blues", domain = map_data$low_education_2015_update, bins=round(edu_bins, digits = 0), na.color = "#D0D0D0")

emp_bins<-c(0, 1, 2)
pal_emp <- colorBin(palette = "Blues", domain = map_data$low_employment, bins=round(emp_bins, digits = 0), na.color = "#D0D0D0")

pop_bins<-c(0, 1, 2)
pal_pop <- colorBin(palette = "Blues", domain = map_data$pop_loss, bins=round(pop_bins, digits = 0), na.color = "#D0D0D0")

ret_bins<-c(0, 1, 2)
pal_ret <- colorBin(palette = "Blues", domain = map_data$retirement_dest, bins=round(ret_bins, digits = 0), na.color = "#D0D0D0")

pov_bins<-c(0, 1, 2)
pal_pov <- colorBin(palette = "Blues", domain = map_data$persistent_poverty, bins=round(pov_bins, digits = 0), na.color = "#D0D0D0")

cpov_bins<-c(0, 1, 2)
pal_cpov <- colorBin(palette = "Blues", domain = map_data$persistent_related_child_poverty, bins=round(cpov_bins, digits = 0), na.color = "#D0D0D0")

map_data %>%
  st_transform(crs = "+init=epsg:4326") %>%
  leaflet(width = "100%", options = leafletOptions(crs = epsg2163, zoomControl=FALSE)) %>%
  setView(-93, 37.8283, zoom=3) %>%
  addPolygons(popup = paste("<strong>County:</strong>", map_data$name, "<br>",
                            "<strong>State:</strong>", map_data$state.x, "<br>",
                            "<strong>Low Edu?:</strong>", map_data$low_education),
              group= "EDU",
              stroke = TRUE,
              smoothFactor = 1.5,
              weight = 1,
              fillOpacity = 0.9,
              opacity = 0.7,
              color = "grey",
              fillColor = ~ pal_edu(low_education)) %>%
  addPolygons(popup = paste("<strong>County:</strong>", map_data$name, "<br>",
                            "<strong>State:</strong>", map_data$state.x, "<br>",
                            "<strong>Low Emp?:</strong>", map_data$low_education),
              group= "EMP",
              stroke = TRUE,
              smoothFactor = 1.5,
              weight = 1,
              fillOpacity = 0.9,
              opacity = 0.7,
              color = "grey",
              fillColor = ~ pal_emp(low_employment)) %>%
  addPolygons(popup = paste("<strong>County:</strong>", map_data$name, "<br>",
                            "<strong>State:</strong>", map_data$state.x, "<br>",
                            "<strong>Pop loss?:</strong>", map_data$pop_loss),
              group= "POP",
              stroke = TRUE,
              smoothFactor = 1.5,
              weight = 1,
              fillOpacity = 0.9,
              opacity = 0.7,
              color = "grey",
              fillColor = ~ pal_pop(pop_loss)) %>%
  addPolygons(popup = paste("<strong>County:</strong>", map_data$name, "<br>",
                            "<strong>State:</strong>", map_data$state.x, "<br>",
                            "<strong>Retire Dest?:</strong>", map_data$retirement_dest),
              group= "RET",
              stroke = TRUE,
              smoothFactor = 1.5,
              weight = 1,
              fillOpacity = 0.9,
              opacity = 0.7,
              color = "grey",
              fillColor = ~ pal_ret(retirement_dest)) %>%
  addPolygons(popup = paste("<strong>County:</strong>", map_data$name, "<br>",
                            "<strong>State:</strong>", map_data$state.x, "<br>",
                            "<strong>Peristent Poverty?:</strong>", map_data$persistent_poverty),
              group= "POV",
              stroke = TRUE,
              smoothFactor = 1.5,
              weight = 1,
              fillOpacity = 0.9,
              opacity = 0.7,
              color = "grey",
              fillColor = ~ pal_pov(persistent_poverty)) %>%
  addPolygons(popup = paste("<strong>County:</strong>", map_data$name, "<br>",
                            "<strong>State:</strong>", map_data$state.x, "<br>",
                            "<strong>Child pov?:</strong>", map_data$persistent_related_child_poverty),
              group= "CPOV",
              stroke = TRUE,
              smoothFactor = 1.5,
              weight = 1,
              fillOpacity = 0.9,
              opacity = 0.7,
              color = "grey",
              fillColor = ~ pal_cpov(persistent_related_child_poverty)) %>%
  addLegend("bottomright", 
           pal = pal_edu, 
           values = ~ low_education,
           labFormat = function(type, cuts, p) {paste0(all_labels)},
           title = "Yes or No?",
           opacity = 1) %>%
  setMapWidgetStyle(list(background= "white")) %>%
  addLayersControl(
    overlayGroups = c("EDU", "EMP", "POP", "RET", "POV", "CPOV"),
    options = layersControlOptions(collapsed = FALSE)) %>%
  hideGroup("EMP") %>%
  hideGroup("POP") %>%
  hideGroup("RET") %>%
  hideGroup("POV") %>%
  hideGroup("CPOV")

## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.

We can also look at the economic base of each county.

unique(map_data$economic_type_label)

## [1] "Nonspecialized"           "Maufacturing"            
## [3] "Federal/State Government" "Farming"                 
## [5] "Recreation"               "Mining"

map_data <- map_data %>%
  mutate(econ_type= case_when(
    economic_type_label=="Nonspecialized" ~ 0,
    economic_type_label=="Maufacturing" ~ 1,
    economic_type_label=="Federal/State Government" ~ 2,
    economic_type_label=="Farming" ~ 3,
    economic_type_label=="Recreation" ~ 4,
    economic_type_label=="Mining" ~ 5))

econ_bins<- c(0, 1, 2, 3, 4, 5, 6)
pal_econ <- colorBin(palette = "BrBG", domain = map_data$econ_type, bins=round(econ_bins, digits = 0), na.color = "#D0D0D0")

map_data %>%
  st_transform(crs = "+init=epsg:4326") %>%
  leaflet(width = "100%", options = leafletOptions(crs = epsg2163, zoomControl=FALSE)) %>%
  setView(-93, 37.8283, zoom=3) %>%
  addPolygons(popup = paste("<strong>County:</strong>", map_data$name, "<br>",
                            "<strong>State:</strong>", map_data$state.x, "<br>",
                            "<strong>Econ Type:</strong>", map_data$economic_type_label),
              stroke = TRUE,
              smoothFactor = 1.5,
              weight = 1,
              fillOpacity = 0.9,
              opacity = 0.7,
              color = "grey",
              fillColor = ~ pal_econ(econ_type)) %>%
  addLegend("bottomright", 
           pal = pal_econ, 
           values = ~ econ_type,
           labFormat = function(type, cuts, p) {paste0(unique(map_data$economic_type_label))},
           title = "County Economic Dependence",
           opacity = 1) %>%
  setMapWidgetStyle(list(background= "white"))

There’s a lot of “nonspecialized” across the county. Not sure how accurate this data is for providing an economic overview of the county as a result (plus it’s at the county) but it’s at least helpful for identifying specific areas where ag or mining continues to be prevalent.

One last thing I’m interested in with this data is how county typologies have changed over time. Perhaps typology changes will be correlated with other outcomes such as current un/employment or vote choice. Previous codes are in a slightly different format, but it should be possible to integrate them. The codes below are from 1989.

old_types <- readxl::read_xls("typology89.xls", sheet="Data") %>% clean_names()
#this seems to not consider metro counties.. so will remove for this quick analysis

clean_types <- old_types %>%
  filter(fm != 8) %>%
  select(fips, state, county_name, fm:ns) %>%
  mutate(econ_type2= case_when(
    fm==1 ~ 3,
    mi==1 ~ 5,
    mf==1 ~ 1,
    gv==1 ~ 3,
    ts==1 ~ 4, #i dont think this is the same as recreation?
    ns==1 ~ 0))


merged_map <- map_data %>%
  left_join(clean_types, by="fips") #metro counties will be NA

comparison <- merged_map %>%
  filter(is.na(fm)==FALSE,
         econ_type != econ_type2) # a lot of mismatches.. maybe come back later