population-shift-2010-2018.Rmd

---
title: "US county population shifts 2010-2018"
author: "Daniel Moul"
date: "`r Sys.Date()`"
output: 
#  html_document:
#    default
  github_document: default
---

<br>

Which US counties have seen the most change in population between the 2010 decennial census and the American Community Survey (ACS) 2014-2018 estimates?

To answer this question, I downloaded county-level population from the 2010 decennial census and corresponding estimates from the ACS 2014-2018 via the Census Bureau's Planning Database at https://www.census.gov/topics/research/guidance/planning-databases.2020.html and scraped data for the largest US cities from https://en.wikipedia.org/wiki/List_of_United_States_cities_by_population

```{r include=FALSE}
library(tidyverse)
library(readr)
library(janitor)
library(jsonlite)
library(lubridate)
library(knitr)
library(kableExtra)
library(glue)
library(scales)
library(sf)
library(tigris)
library(broom)

# remotes::install_github("hrbrmstr/hrbrthemes")
# hrbrthemes::import_roboto_condensed() 
# hrbrthemes::import_plex_sans()
library(hrbrthemes)

knitr::opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE, error=TRUE)
options(tigris_use_cache = TRUE)

# for ggplot
theme_set(theme_ipsum_ps())
my_plot_caption <- "By Daniel Moul"

```

```{r include=FALSE}
us_not_states <- tribble(
  ~statefp,                 ~statename, 
   69,                       "Commonwealth of the Northern Mariana Islands",
   66,                        "Guam",
   78,                        "United States Virgin Islands",
   60,                        "American Samoa",
   11,                        "District of Columbia",
   72,                        "Puerto Rico",
   78,                        "U.S. Virgin Islands",
   64,                        "FSM More islands",
   68,                        "MH - more islands?",
   69,                        "Northern Mariana Islands",
   70,                        "More islands",
   74,                        "U.S. Minor Outlying Islands",
   76,                        "Federated States of Micronesia"
 ) %>%
  mutate(statefp = as.character(statefp)) %>%
  arrange(statefp)

state_names <- tigris::states(cb = FALSE, class = "sf") %>% # inclues region and div when cb=FALSE
  as.data.frame() %>%
  select(region = REGION, 
         division = DIVISION, 
         statefp = STATEFP, 
         statename = NAME, 
         aland = ALAND, 
         awater = AWATER, 
         intptlat = INTPTLAT, 
         intptlon = INTPTLON) %>%
  anti_join(us_not_states, by = "statename") %>%
  arrange(statefp)

us_states48_20m <- tigris::states(cb = TRUE, resolution = "20m", class = "sf") %>%
  clean_names() %>%
  rename(statename = name) %>%
  anti_join(us_not_states, by = "statename") %>%
  arrange(statefp) %>% 
  filter( !statename %in% c("Alaska", "Hawaii") ) # continental US only

# for map
us_counties48_borders <- tigris::counties(state = state_names$statename, 
                                  cb = TRUE, 
                                  resolution = "20m", 
                                  class = "sf") %>%
  clean_names() %>%
  left_join(., 
            state_names %>% select(statefp, statename),
            by = c("statefp")) %>%
  anti_join(., us_not_states, by = "statename") %>%
  filter( !statename %in% c("Alaska", "Hawaii") )  

```

```{r}
# data from 
# * https://api.census.gov/data/2010/dec/sf1/examples.html
# * https://api.census.gov/data/2010/dec/sf1?get=P001001,NAME&for=county:*&in=state:*

# pop_county_raw <- fromJSON('./data/P001001.2010.sf1.json') %>%
#   as.data.frame()
# names(pop_county_raw) <- pop_county_raw[1, ]
# 
# pop_county <- pop_county_raw %>%
#   filter(!state == "state") %>%
#   rename(population = P001001,
#          county_name = NAME) %>%
#   mutate(population = as.numeric(population),
#          state = str_pad(state, width = 2, side = "left", pad = "0")
#          )

# Better yet, get more current county population estimates from https://www.census.gov/topics/research/guidance/planning-databases.2020.html

pop_county_raw <- read_csv('./data/data-mean-population/2020-special-pdb/pdb2020stcov2_us.csv') %>%
  clean_names()

pop_county <- pop_county_raw %>%
  filter(geog_level == "County") %>%
  select(pop2010 = tot_population_cen_2010,
         pop2018 = tot_population_acs_14_18, #ACS 2014-2018 estimates are likely to be more accurate/current
         county_name, state, county) %>%
  anti_join(., 
            us_not_states, 
            by = c("state" = "statefp")) %>%
  # setting the small number of missing values to zero is good enough for our purposes
  mutate(pop2010 = if_else(is.na(pop2010),
                           pop2018,
                           pop2010),
         pop2018 = if_else(is.na(pop2018),
                           pop2010,
                           pop2018)
         ) %>%
  replace_na(list(population = 0)) %>%
  mutate(diff_total = pop2018 - pop2010,
         pct_diff_total = if_else(pop2018 > 0 & pop2010 > 0,
                                  (pop2018 - pop2010) / pop2010,
                                  0) # or use NA_real_ however 0 is good enough for our purposes
  ) %>%
  distinct(state, county, .keep_all = TRUE) %>%
  group_by(state) %>%
  mutate(pct_state_pop = pop2010 / sum(pop2010)) %>%
  ungroup()

```

Nearly all county-level changes are within +/- 10%. Outliers are not show in the histogram below.

<br>

```{r out.width="100%"}
pop_county %>%
  filter(abs(pct_diff_total) <= 0.25) %>%
  ggplot(aes(pct_diff_total)) +
  geom_histogram(bins = 50) + 
  #scale_x_log10(labels = label_number_si(), breaks = 10^(1:7)) + #comma_format()
  scale_x_continuous(labels = percent_format()) +
  labs(title = "US counties: percent population change",
       subtitle = "Comparing 2010 decennial census and ACS 2014-2018 estimates",
       x = "Population pct change",
       y = "Number of counties",
       caption = my_plot_caption)
```

```{r include=FALSE}

us_pop_change = pop_county %>%
  summarize(diff_us = sum(pop2018 - pop2010),
            pct_diff_us = diff_us / sum(pop2010))

my_map_caption <- glue("US pop change: {comma(us_pop_change$diff_us)}",
                       " ({percent(us_pop_change$pct_diff_us, accuracy = 0.1)})",
                       "\nSource: US Census Planning Database July 2020. By Daniel Moul")

my_crs <- 2163 # for continental USA (originally used 4326)

us_counties48 <- left_join(us_counties48_borders, 
                           pop_county %>% rename(statefp = state,
                                                 countyfp = county,
                           ),
                           by = c("statefp", "countyfp")
)

```

```{r}
cities <- read_csv("./data/data-mean-population/cities.csv") %>%
  filter(!str_detect(state, "Alaska|Hawaii")) %>%
  group_by(state) %>%
  slice_max(n = 5, order_by = estimate2019) %>%
  ungroup() %>%
  st_as_sf(coords = c("long", "lat")) %>%
  st_set_crs(4269)

```

<br>

Some observations based on the plot below:

* Metropolitan area have fared well during this time period; most cities are in or surrounded by counties that have grown
* Rural areas seem to have suffered in most states
* The population changes due to the fracking boom in North Dakota are quite evident
* There is more shrinkage in the midwest and northeast and more growth in the sun belt

<br>

```{r out.width="100%"}

us_counties48 %>%
  # first deal with outliers
  mutate(pct_diff_total = case_when(
    pct_diff_total > 0.25             ~ 0.25,
    pct_diff_total < -0.25            ~ -0.25,
    TRUE                              ~ pct_diff_total),
         ) %>%
  ggplot() +
  geom_sf(aes(fill = pct_diff_total),
              color = "darkgrey", size = 0.1, alpha = 1.0) + # draw county boarders and fill colors
  geom_sf(data = st_geometry(us_states48_20m),
          fill = NA, color = "slategrey", alpha = 1.0, size = 0.3) + # draw state borders too
  geom_sf(data = cities, size = 0.75, color = "slateblue", alpha = 0.3) +
  coord_sf(crs = my_crs) +
  scale_fill_gradient2(low = "#590C25", high = "#99d8c9", mid = "white", midpoint = 0) +
  labs(title = "County population changes 2010-2018",
       subtitle = glue("Color gradient capped at +/- 25%"),
       x = "", y = "",
       fill = "% change",
       caption = my_map_caption)

```

<br>

The differences are clearer when I use categories to reflect intervals of population change rates.

<br>

```{r out.width="100%"}

my_break_colors <- c("#2ca25f",
                     "#99d8c9", 
                     "#edf8fb",
                     #"#FFFEC8", 
                     #"#EED88C", 
                     "#E1B759",
                     "#EC8E43", 
                     #"#EA491D", 
                     "#BD1A2F"
                     #"#AB1130", 
                     #"#590C25"
                     )

us_counties48 %>%
  ggplot() +
  geom_sf(aes(fill = cut(pct_diff_total, c(-0.5, -0.1, -0.05, 0, 0.05, 0.1, 1.0))),
              color = "darkgrey", size = 0.1) + # draw county boarders and fill colors
  geom_sf(data = st_geometry(us_states48_20m),
          fill = NA, color = "slategrey", alpha = 1.0, size = 0.3) + # draw state borders too
  geom_sf(data = cities, size = 1, color = "slateblue", alpha = 0.3) +
  coord_sf(crs = my_crs) +
  scale_fill_manual(values = rev(my_break_colors)) +
  labs(title = "County population changes 2010-2018",
       subtitle = glue("+/- 5% and 10% or more"),
       x = "", y = "",
       fill = "% change",
       caption = my_map_caption)

```

<br>

The US population has grown `r percent(us_pop_change$pct_diff_us, accuracy = 0.1)` during this time period. Which counties have grown slower or faster than this?

<br>

```{r out.width="100%"}

my_break_colors <- c("#006d2c", #dark green
                     "#2ca25f",
                     "#99d8c9", 
                     #"#edf8fb", too grey/neurtral
                     "#FFFEC8", 
                     #"#EED88C", 
                     "#E1B759"
                     #"#EC8E43"
                     #"#EA491D", 
                    # "#BD1A2F",
                     #"#AB1130", 
                     #"#590C25"
                     )

us_counties48 %>%
  ggplot() +
  geom_sf(aes(fill = cut(pct_diff_total, c(-0.5, 0, round(us_pop_change$pct_diff_us, 3), 0.1, 0.25, 1.0))),
              color = "darkgrey", size = 0.1) + # draw county boarders and fill colors
  geom_sf(data = st_geometry(us_states48_20m),
          fill = NA, color = "slategrey", alpha = 1.0, size = 0.3) + # draw state borders too
  geom_sf(data = cities, size = 1, color = "slateblue", alpha = 0.3) +
  coord_sf(crs = my_crs) +
  scale_fill_manual(values = rev(my_break_colors)) +
  labs(title = "County population changes 2010-2018",
       subtitle = glue("Counties that lost in absolute terms, grew but slower than US,",
                       "\ngrew faster than US up to 10%, 25%, and 100%"),
       x = "", y = "",
       fill = "pct change",
       caption = my_map_caption)

```

<br>

The counties with a smaller portion of the state's population are losing people.

<br>

```{r out.width="100%"}
pop_county %>%
  ggplot(aes(pct_state_pop, pct_diff_total)) +
  geom_point(size = 1, alpha = 0.3) + 
  geom_smooth(method = "lm") +
  scale_x_log10(labels = percent_format()) +
  scale_y_continuous(labels = percent_format(), breaks = 0.1 * -3:3,
                     limits = c(-0.3, 0.3)) +
  labs(title = "US counties: percent population change in country\nby percent of state population",
       y = "Population change pct",
       x = "County's pct of state population log10",
       caption = my_plot_caption)
```

<br>

The above plot might suffer from the fact that in rural states, counties can have a high percentage of a small state population. But apparently not; the plot below shows the same trend. In general, the smaller the county population the greater the rate of loss.

<br>

```{r out.width="100%"}

pop_county %>%
  ggplot(aes(pop2010, pct_diff_total)) +
  geom_point(size = 1, alpha = 0.3) + 
  geom_smooth(method = "lm") +
  scale_x_log10(labels = label_number_si()) +
  scale_y_continuous(labels = percent_format(), breaks = 0.1 * -3:3,
                     limits = c(-0.3, 0.3)) +
  labs(title = "US counties: percent population change in country\nby county population",
       y = "Population change pct",
       x = "County population log10",
       caption = my_plot_caption)

```

<br>
<br>

(end of document)