nc-pop-change.qmd

# NC population and population change; under/over 18 years old

While family size in the USA has been shrinking since before the start of this dataset^[<https://www.statista.com/statistics/183657/average-size-of-a-family-in-the-us/> ], the decline was not uniform across NC counties. Net population change hides interesting dynamics. The following uses decennial US Census results 1970 to 2020 provided by IPUMS NHGIS.

Growth rates are over ten years.

:::{.callout-note collapse="true"}
## Note on 2020 census differential privacy

The `tidycensus` package offers the following helpful warning when returning results from 2020, which is relevant to this endeavor:

>2020 decennial Census data use differential privacy, a technique that introduces errors into data to preserve respondent confidentiality. Small counts should be interpreted with caution. See <https://www.census.gov/library/fact-sheets/2021/protecting-the-confidentiality-of-the-2020-census-redistricting-data.html> for additional guidance.
:::

<br>

```{r}
#| label: setup
#| echo: false
#| warning: false
#| message: false

source("./scripts/my-setup.R") # use for all explorationXX.qmd files
source("./scripts/constants2.R")

# TODO: create functions2.R
source("./scripts/functions2.R") 

# library(geofacet)
# library(ggridges)
# library(ggrepel)

#library(elevatr)
library(terra)
library(stars)
library(rcartocolor)
library(fs)

library(english) #for ordinal() to spell out words
library(geomtextpath)

# install.packages("Rcpp", repos="https://rcppcore.github.io/drat")
# to address https://github.com/rspatial/terra/issues/30 until the next Rccp release
# "Error in x$.self$finalize() : attempt to apply non-function"
# which I only see when running 'quarto render' from the command line 2022-11-25


```

```{r}
nc_city_boundaries <- tar_read(nc_city_boundaries)

county_region_df <- tar_read(nc_county_region_mapping)

nc_county_boundaries <- tar_read(nc_county_boundaries) |>
  inner_join(
    county_region_df,
    by = "county"
  )

urban_crescent_boundary <- nc_county_boundaries |>
  filter(county %in% urban_crescent) |>
  summarize(boundary = st_union(geometry))

coastal_boundary <-  nc_county_boundaries |>
  filter(county %in% coastal) |>
  summarize(boundary = st_union(geometry))

mountain_boundary <- nc_county_boundaries |>
  filter(county %in% mountain) |>
  summarize(boundary = st_union(geometry))

manufacturing_boundary <- nc_county_boundaries |>
  filter(county %in% manufacturing) |>
  summarize(boundary = st_union(geometry))

# agriculture boundary is everything else

nc_state_boundary <- st_union(nc_county_boundaries)
nc_state_boundary_sf <- st_as_sf(nc_state_boundary)

```

```{r}
#| label: define-pop_u18
#| message: false

# TODO: put this in a function and pre-compute using target

years_in_scope <- c("1970", "1980", "1990", "2000", "2010", "2020")
# can use 2020 for these; otherwise need to use 2015-2019:
# Table 1: (A00) Total Population
#         A00AA:       Persons: Total
# Table 2: (D08) Persons by Age [2]: Children and Adults
#         D08AA:       Persons: Under 18 years
#         D08AB:       Persons: 18 years and over

pop_u18 <- read_csv("./data/ipums/nhgis0008_csv/nhgis0008_ts_nominal_county.csv") |>
  filter(STATE == "North Carolina",
         YEAR %in% years_in_scope) |>
  select(GISJOIN, YEAR, COUNTY = NAME, 
        pop  = A00AA, # Persons: Total,
        pop_u18 = D08AA, # Persons: Under 18 years
        pop_18p = D08AB  # Persons: 18 years and over
         ) |>
  clean_names() |>
  mutate(county = str_to_title(str_extract(county, "[A-z]+")),
         county = case_when(
           county == "New"         ~ "New Hanover",
           county == "Mc"          ~ "McDowell", # prob in 1970 data
           county == "Mcdowell"    ~ "McDowell",
           TRUE                    ~ county
           ),
         pop_pct_u18 = pop_u18 / pop,
         pop_pct_18p = pop_18p / pop) |>
  left_join(county_region_df,
            by = "county") |>
  mutate(region = factor(region, levels = region_colors$region))

nc_total <- pop_u18 |>
  group_by(year) |>
  summarize(
    pop = sum(pop),
    pop_u18 = sum(pop_u18),
    pop_18p = sum(pop_18p)
  ) |>
  ungroup() |>
  mutate(
    pop_pct_u18 = pop_u18 / pop,
    pop_pct_18p = pop_18p / pop,
    county = "NC total"
  ) |>
  mutate(change_since_1970 = pop_pct_u18 / first(pop_pct_u18),
         region = "NC")

nc_total_long <- nc_total |>
  pivot_longer(
    cols = starts_with("pop"),
    names_to = "variable",
    values_to = "value"
  ) |>
  mutate(sign_value = sign(value)) |>
  arrange(region, variable, year) |>
  group_by(region, variable) |>
  # growth_rate for any row refers to the growth between this an the next census
  mutate(
    pct_growth_naive = abs(value / lag(value)) - 1,
    pct_growth_factor = case_when(
      sign_value + lag(sign_value) == 2     ~ 1,
      sign_value + lag(sign_value) == 0     ~ -1,
      sign_value + lag(sign_value) == -2    ~ -1,
      TRUE                                   ~ 0 # error
    ),
    pct_growth = pct_growth_naive * pct_growth_factor,
         # pct_growth_original = lead(value) / value - 1,
         # pct_growth = lead(value_xformed) / value_xformed - 1,
         positive_growth = case_when(
           pct_growth > 0          ~ "positive",
           pct_growth < 0          ~ "negative",
           abs(pct_growth) < 1e15  ~ "zero",
           is.na(pct_growth)       ~ "zz info not available",
           TRUE                    ~ "error in case_when() input"
         )
         ) |>
  ungroup()

pop_u18_summary <- bind_rows(
  pop_u18,
  nc_total) |>
  group_by(region, year) |>
  summarize(pop = sum(pop),
            pop_u18 = sum(pop_u18),
            pop_18p = sum(pop_18p),
            pop_pct_u18 = pop_u18 / pop,
            pop_pct_18p = pop_18p / pop) |>
  ungroup() |>
  mutate(region = factor(region, levels = region_colors$region))

d_pop_growth_u18_long <- pop_u18 |>
  pivot_longer(
    cols = starts_with("pop"),
    names_to = "variable",
    values_to = "value"
  ) |>
  mutate(sign_value = sign(value)) |>
  arrange(county, variable, year) |>
  group_by(county, variable) |>
  # growth_rate for any row refers to the growth between this an the next census
  mutate(
    pct_growth_naive = abs(value / lag(value)) - 1,
    pct_growth_factor = case_when(
      sign_value + lag(sign_value) == 2     ~ 1,
      sign_value + lag(sign_value) == 0     ~ -1,
      sign_value + lag(sign_value) == -2    ~ -1,
      TRUE                                   ~ 0 # error
    ),
    pct_growth = pct_growth_naive * pct_growth_factor,
         # pct_growth_original = lead(value) / value - 1,
         # pct_growth = lead(value_xformed) / value_xformed - 1,
         positive_growth = case_when(
           pct_growth > 0          ~ "positive",
           pct_growth < 0          ~ "negative",
           abs(pct_growth) < 1e15  ~ "zero",
           is.na(pct_growth)       ~ "zz info not available",
           TRUE                    ~ "error in case_when() input"
         )
         ) |>
  ungroup() #|>
  #filter(year != "1970") # TODO: make sure consuming sections do this if needed

d_pop_growth_u18_summary_long <- pop_u18_summary |>
  pivot_longer(
    cols = starts_with("pop"),
    names_to = "variable",
    values_to = "value"
  ) |>
  mutate(sign_value = sign(value)) |>
  arrange(region, variable, year) |>
  group_by(region, variable) |>
  # growth_rate for any row refers to the growth between this an the next census
  mutate(
    pct_growth_naive = abs(value / lag(value)) - 1,
    pct_growth_factor = case_when(
      sign_value + lag(sign_value) == 2     ~ 1,
      sign_value + lag(sign_value) == 0     ~ -1,
      sign_value + lag(sign_value) == -2    ~ -1,
      TRUE                                   ~ 0 # error
    ),
    pct_growth = pct_growth_naive * pct_growth_factor,
         # pct_growth_original = lead(value) / value - 1,
         # pct_growth = lead(value_xformed) / value_xformed - 1,
         positive_growth = case_when(
           pct_growth > 0          ~ "positive",
           pct_growth < 0          ~ "negative",
           abs(pct_growth) < 1e15  ~ "zero",
           is.na(pct_growth)       ~ "zz info not available",
           TRUE                    ~ "error in case_when() input"
         )
         ) |>
  ungroup() #|>
  #filter(year != "1970") # TODO: make sure consuming sections do this if needed

```

## By region

Using the regions defined earlier, let's look at the total population, 18+ and under 18 for each region. Population is on a linear scale with each row having its own scale to maximize the detail within each row. In all regions the 18+ group has increased, while the under 18 group stayed about the same in regions other than the urban crescent and the coastal counties.

<br>

```{r fig.height=8, fig.width=8}
#| label: fig-summary_population-facet-grid-1
#| fig-cap: "NC regional population (1980 - 2022). Each row provides comparison across the variables for a region"
#| fig-height: 8
#| fig-width: 8
#| column: screen-inset-right

region_levels_new <- rev(region_colors$region)

region_color_levels_new = rev(region_colors$color)

d_pop_growth_u18_summary_long |>
  filter(!str_detect(variable, "pct")) |>
  mutate(region = factor(region, levels = region_levels_new)) |> 
  ggplot(aes(year, value, color = region)) +
  geom_point(size = 2, alpha = 0.8) +
  geom_line() +
  facet_grid(region ~ variable, scales = "free_y") +
  scale_x_continuous(breaks = c(1970, 1990, 2010)) +
  scale_y_continuous(labels = label_number(scale_cut = cut_short_scale())) +
  scale_color_manual(values = region_color_levels_new) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "none"
        ) +
  labs(title = "NC population (1970 - 2022)",
       subtitle = "US decennial census counts",
       x = NULL,
       y = "Population (scale varies by row)",
       caption = my_caption_nhgis)

```

<br>

The plot below allows comparison of total population, 18+ and under 18 across the regions. The uniqueness of the urban crescent stands out.

<br>

```{r fig.height=8, fig.width=8}
#| label: fig-summary_population-facet-grid-2
#| fig-cap: "NC regional population (1980 - 2022). Each row provides comparison across the regions for each variable."
#| fig-height: 8
#| fig-width: 8
#| column: screen-inset-right

d_pop_growth_u18_summary_long |>
  filter(!str_detect(variable, "pct")) |>
  mutate(region = factor(region, levels = region_levels_new)) |>
  ggplot(aes(year, value, color = region)) +
  geom_point(size = 2, alpha = 0.8) +
  geom_line() +
  facet_grid(variable ~ region, scales = "free_y") +
  scale_x_continuous(breaks = c(1970, 1990, 2010)) +
  scale_y_continuous(labels = label_number(scale_cut = cut_short_scale())) +
  scale_color_manual(values = region_color_levels_new) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "none"
        ) +
  labs(title = "NC population (1970 - 2022)",
       subtitle = "US decennial census counts",
       x = NULL,
       y = "Population (scale varies by row)",
       caption = my_caption_nhgis)

```

<br>

## Ranked by county

We can rank counties by their population. On the left is a simple ranking; on the right population is plotted on the y axis.

<br>

```{r fig.height=14, fig.width=10}
#| label: fig-nc-counties-total-pop-rank
#| fig-cap: "County rank: total population"
#| fig-width: 10
#| fig-height: 14
#| column: screen-inset-right

plot_rank_population(my_title = "NC counties ranked by total population"
          )

```

<br>

## Growth rates

The plot below allows comparison of growth rates for total population, 18+ and under 18 *for each region*. In the ten years to 2020 the 18+ population grew in most regions. In contrast, the under 18 population was about net flat; growth in the urban crescent was offset by declines in the other regions.

<br>

```{r fig.height=8, fig.width=8}
#| label: fig-summary-population-growth-rates-facet-grid-1
#| fig-cap: "NC regional population growth rates (1980 - 2020). Each row provides comparison across the variables for a region."
#| fig-height: 8
#| fig-width: 8
#| column: screen-inset-right

region_levels_new <- rev(region_colors$region)

region_color_levels_new = rev(region_colors$color)

d_pop_growth_u18_summary_long |>
  filter(!str_detect(variable, "pct"),
         year != 1970) |>
  mutate(region = factor(region, levels = region_levels_new)) |> 
  ggplot(aes(year, pct_growth, color = region)) +
  geom_hline(yintercept = 0, lty = 2, linewidth = 0.25, alpha = 0.2) +
  geom_point(size = 2, alpha = 0.8) +
  geom_line() +
  scale_x_continuous(breaks = c(1970, 1990, 2010)) +
  scale_y_continuous(labels = percent_format()) +
  scale_color_manual(values = region_color_levels_new) +
  facet_grid(region ~ variable) + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "none"
        ) +
  labs(title = "NC population growth rate (1980 - 2022)",
       subtitle = "US decennial census counts",
       x = NULL,
       y = "Popluation growth rate",
       caption = my_caption_nhgis)
```

<br>

The plot below allows comparison of growth rates for total population, 18+ and under 18 *across the regions*.

<br>

```{r fig.height=8, fig.width=8}
#| label: fig-summary-population-growth-rates-facet-grid-2
#| fig-cap: "NC regional population growth rates (1980 - 2020). Each row provides comparison across the regions for a variable."
#| fig-height: 8
#| fig-width: 8
#| column: screen-inset-right

d_pop_growth_u18_summary_long |>
  filter(!str_detect(variable, "pct"),
         year != 1970) |>
  mutate(region = factor(region, levels = region_levels_new)) |> 
  ggplot(aes(year, pct_growth, color = region)) +
  geom_hline(yintercept = 0, lty = 2, linewidth = 0.25, alpha = 0.2) +
  geom_point(size = 2, alpha = 0.8) +
  geom_line() +
  facet_grid(variable ~ region) + 
  scale_x_continuous(breaks = c(1970, 1990, 2010)) +
  scale_y_continuous(labels = percent_format()) +
  scale_color_manual(values = region_color_levels_new) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "none"
        ) +
  labs(title = "NC population growth rate (1980 - 2020)",
       subtitle = "US decennial census counts",
       x = NULL,
       y = "Popluation growth rate",
       caption = my_caption_nhgis)

```

<br>

In recent censuses it's generally true that the counties with larger population grew while those with smaller, shrank. This wasn't the case in the 10 years leading to the 1980 census. By 2000 it was quite pronounced in the under 18 population, and since then the pattern is visible in all three categories: total population, 18+, and under 18.

The grow rate slowed a lot in the 10 years to 2020.

The regional dynamics are visible when plotting county growth rates on a map.  For rural counties in the east and west, the 2000 census offered reason for hope for the future, since the number of children increased; since then the under 18 population has declined again.

<br>

```{r fig.height=6, fig.width=10}
#| label: fig-nc-map-fact-u18-pct
#| fig-cap: "NC county growth rates 1990-2020"
#| fig-height: 8
#| fig-width: 14
#| column: screen-inset-right

data_for_plot <- d_pop_growth_u18_long |>
  filter(year != 1970) |>
  group_by(variable, county) |>
  mutate(order = row_number()) |>
  fill(c(pct_growth, positive_growth), .direction = "down") |>
  ungroup() |>
  mutate(order = factor(as.character(order), levels = c("1", "2", "3", "4")))

nc_county_boundaries |>
  left_join(data_for_plot,
          by = "county",
          multiple = "all") |> 
  filter(variable %in% c("pop", "pop_18p", "pop_u18")) |> #"pop_pct_18p", 
  ggplot() +
  geom_sf(aes(fill = pct_growth), color = "white", size = 0.1) +
  facet_grid(year ~ variable) +
  scale_fill_gradient2(low = "firebrick", high = "lightskyblue", mid = "grey80", 
                       midpoint = 0, labels = percent) +
  theme_map() +
  theme(plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "right") +
  labs(title = "NC county growth rates",
       subtitle = glue("Total, 18+, and under 18 population growth in the 10 years to 1990-2020"),
       fill = "Population growth",
       caption = my_caption_nhgis)

```

<br>

## Components of population change

Population changes are due to (1) births and deaths ("natural increase"); and (2) migration in and out ("net migration").  In the case of North Carolina as a whole, most of the increase is due to net migration, as seen in @fig-estimates-state-rates-natural-migration and @fig-estimates-five-year-four-vars-state-summary. 

In @fig-estimates-state-rates-natural-migration each thin line shows one county's rate of natural increase or rate of net migration. The thick line is the weighted mean.

<br>

```{r}
#| label: define-multi-year-pop-change-components
#| include: false

# state-level
my_vars <- c("BIRTHS", "DEATHS", "NATURALINC", 
                "RBIRTH", "RDEATH", "RNATURALINC", 
                "DOMESTICMIG", "INTERNATIONALMIG", "NETMIG",
                "RDOMESTICMIG", "RINTERNATIONALMIG", "RNETMIG")

get_nc_pop <- function(my_year) {
  get_estimates(
  geography = "state",
  state = "NC",
  variables = "POP",
  year = my_year
  ) |>
  mutate(value = if_else(str_detect(variable, "^R"),
                         value / 1000,
                         value),
         datayear = my_year)
}

multi_year_nc_pop <- map_df(2015:2019, get_nc_pop) |>
  mutate(NAME = "North Carolina") |>
  rename(state = NAME)

######

get_nc_estimates <- function(my_year) {
  get_estimates(
  geography = "state",
  state = "NC",
  variables = my_vars,
  year = my_year
  ) |>
  mutate(value = if_else(str_detect(variable, "^R"), 
                         value / 1000, 
                         value),
         datayear = my_year,
         variable = factor(variable, levels = my_vars))
}

multi_year_nc <- map_df(2015:2019, get_nc_estimates) |>
  mutate(NAME = "North Carolina") |>
  rename(state = NAME) |>
  inner_join(multi_year_nc_pop |>
               rename(pop_total = value) |>
               select(-c(variable, GEOID, state)),
            by = "datayear")
# The Population Estimates API is not available in tidycensus for years prior to 2015


# county-level

get_nc_county_pop <- function(my_year) {
  get_estimates(
  geography = "county",
  state = "NC",
  variables = "POP",
  year = my_year
  ) |>
  mutate(value = if_else(str_detect(variable, "^R"),
                         value / 1000,
                         value),
         datayear = my_year)
}

multi_year_nc_county_pop <- map_df(2015:2019, get_nc_county_pop) |>
  #mutate(NAME = "North Carolina") |>
  rename(county = NAME,
         pop_total = value) |>
  mutate(county = str_remove(county, " County.*")) |>
  select(-variable)

######

get_nc_county_estimates <- function(my_year) {
  get_estimates(
  geography = "county",
  state = "NC",
  variables = my_vars,
  year = my_year
  ) |>
  mutate(value = if_else(str_detect(variable, "^R"), 
                         value / 1000, 
                         value),
         datayear = my_year,
         variable = factor(variable, levels = my_vars))
}

multi_year_nc_counties <- map_df(2015:2019, get_nc_county_estimates) |>
  mutate(NAME = str_remove(NAME, " County.*")) |>
  rename(county = NAME)
# The Population Estimates API is not available in tidycensus for years prior to 2015

```

```{r fig.height=6, fig.width=10}
#| label: fig-estimates-state-rates-natural-migration
#| fig-cap: "Rates of population change in North Carolina (2015-2019)"
#| fig-height: 6
#| fig-width: 10
#| column: screen-inset-right

data_for_mean_plot <- multi_year_nc_counties |>
  #filter(str_detect(variable, "^R")) |>
  filter(str_detect(variable, "NETMIG|NATURALINC")) |>
  pivot_wider(id_cols = c(county, GEOID, datayear), names_from = "variable", values_from = "value") |>
  group_by(datayear) |>
  reframe(mean_rnaturalinc = weighted.mean(RNATURALINC, w = abs(NATURALINC)),
         mean_rnetmig = weighted.mean(RNETMIG, w = abs(NETMIG))) |>
  ungroup() |>
  #filter(variable %in% c("RNETMIG", "RNATURALINC")) |>
  pivot_longer(cols = starts_with("mean_"), names_to = "variable", values_to = "value") |>
  mutate(variable = if_else(variable == "mean_rnetmig", "Rate of net migration", "Rate of natural increase"))

multi_year_nc_counties |>
  filter(str_detect(variable, "^R")) |>
  filter(variable %in% c("RNETMIG", "RNATURALINC")) |>
  mutate(variable = if_else(variable == "RNETMIG", "Rate of net migration", "Rate of natural increase")) |>
  ggplot(aes(x = datayear, y = value, color = variable)) +
  geom_hline(yintercept = 0, lty = 2, linewidth = 1, color = "firebrick", alpha = 0.25) +
  geom_line(aes(group = county), alpha = 0.1, show.legend = FALSE) +
  #geom_smooth(method = "lm", formula = 'y ~ x', se = FALSE, show.legend = FALSE) +
  geom_line(data = data_for_mean_plot,
            aes(x = datayear, y = value, color = variable),
            linewidth = 1, alpha = 0.8) +
  scale_y_continuous(labels = label_percent()) +
  scale_color_viridis_d(end = 0.8) +
  coord_cartesian(ylim = c(-0.02, 0.02)) +
  facet_wrap(~ variable, ncol = 2) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "top") +
  guides(color = guide_legend(override.aes = list(linewidth = 3))) +
  labs(
    title = "NC population increase 2015-2019",
    subtitle = glue("Roughtly one third was due to natural increase; two thirds due to net migration"),
    x = NULL,
    y = "Yearly rate",
    color = NULL,
    caption = my_caption
  )

```

<!-- br>< -->

<!-- In this time period the state-level pattern repeats each year: there are more births than deaths, and domestic and international migration are net positive, resulting in a population growth of about one percent each year. -->

<!-- <br> -->

```{r fig.height=6, fig.width=8}
#| label: fig-estimates-five-year-four-vars-state-summary
#| fig-cap: "Population change in NC counties over the five years 2015-2019"
#| fig-height: 6
#| fig-width: 8
#| column: screen-inset-right
#| eval: false


#factor_county_stretch <- 0 # Not needed for state summary
factor_datayear_stretch <- 8
factor_variable_stretch <- 2

data_for_plot_tmp <- multi_year_nc |>
  filter(!str_detect(variable, "^R"),
         !variable %in% c("NETMIG", "NATURALINC")) |>
  mutate(variable = fct_drop(variable),
         value = if_else(variable == "DEATHS", 
                         -1 * value,
                         value),
         value_pct_pop = value / pop_total
         ) |>
  arrange(datayear, variable) |>
  mutate(
    xend = cumsum(value_pct_pop),
    xstart = lag(xend, default = 0),
    row_id = row_number()
    ) |>
  arrange(datayear, variable)

datayear_levels <- data_for_plot_tmp  |>
  distinct(datayear) |>
  arrange(datayear) |>
  mutate(datayear_id = (row_number() - 1) *  factor_datayear_stretch)

variable_levels <- data_for_plot_tmp  |>
  distinct(variable) |>
  arrange(variable) |>
  mutate(variable_id = (row_number() - 1) *  factor_variable_stretch)

data_for_plot <- data_for_plot_tmp |>
  inner_join(datayear_levels,
             by = "datayear") |>
  inner_join(variable_levels,
             by = "variable") |>
  mutate(
    y_level = datayear_id + variable_id,
    group = case_when(
      variable %in% c("BIRTHS", "DEATHS", "NATURALINC")      ~ "natural",
      str_detect(variable, "MIG")                            ~ "migration",
      TRUE                                                   ~ "ERROR"
    )
  ) |>
  mutate(
    label_datayear = if_else(variable_id == 2,
                           as.character(datayear),
                           "")
  )

data_for_plot |>
  ggplot() +
  geom_vline(xintercept = 0, lty = 2, linewidth = 0.5, color = "firebrick", alpha = 0.25) +
  geom_hline(aes(yintercept = if_else(variable_id == 0,
                                      datayear_id - 1,
                                      NA_real_)),
             na.rm = TRUE, linewidth = 0.2, alpha = 0.5) +
  geom_segment(aes(x = xstart, xend = xend, y = y_level, yend = y_level, color = variable),
               linewidth = 3) +
  geom_text(aes(x = -0.005, #x_value_label,
                y = y_level,
                label = label_datayear),
            hjust = 1, na.rm = TRUE, size = 4) +
  scale_x_continuous(labels = label_percent()) +
  scale_y_continuous(expand = expansion(mult = c(0, .025))) +
  scale_color_viridis_d(option = "H", end = 0.8) +
  coord_cartesian(xlim = c(-0.01, NA)) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "top") +
  guides(color = guide_legend(override.aes = list(linewidth = 3))) +
  
  labs(
    title = "Population change in NC in the five years 2015-2019",
    subtitle = "Percent of state population",
    x = "Change in state population since 2014",
    y = NULL,
    color = NULL,
    caption = my_caption
  )

```

<br>

As is hinted in the "RNETMIG" panel in @fig-estimates-state-rates-natural-migration, county-level dynamics vary considerably. First, let's look at this by region.

Observations

* The coastal counties and those in the urban crescent are the only regions that are all growing
* Domestic in- and out-migration are the biggest contributors to county population growth or decline

<br>

```{r fig.height=10, fig.width=12}
#| label: fig-estimates-five-year-four-vars-6
#| fig-cap: "Population change in NC counties over the five years 2015-2019"
#| fig-height: 10
#| fig-width: 12
#| column: screen-inset-right


factor_county_stretch <- 40
factor_datayear_stretch <- 8
factor_variable_stretch <- 2

data_for_plot_tmp <- multi_year_nc_counties |>
  filter(!str_detect(variable, "^R"),
         !variable %in% c("NETMIG", "NATURALINC")) |>
  inner_join(multi_year_nc_county_pop,
               by = c("GEOID", "county", "datayear")) |>
  mutate(variable = fct_drop(variable),
         value = if_else(variable == "DEATHS", 
                         -1 * value,
                         value),
         value_pct_pop = value / pop_total
         ) |>
  arrange(datayear, variable) |>
  inner_join(county_region_df,
             by = "county") |>
  group_by(GEOID, county, region) |>
  mutate(
    xend = cumsum(value_pct_pop),
    xstart = lag(xend, default = 0),
    row_id = row_number()
    ) |>
  ungroup() |>
  arrange(GEOID, county, datayear, variable) 

county_levels <- data_for_plot_tmp  |>
  group_by(region, county) |>
  mutate(end_xend = xend[max(row_id)]) |>
  ungroup() |>
  distinct(county, region, end_xend) |>
  group_by(region) |>
  arrange(desc(end_xend)) |>
  mutate(county_id = (row_number() - 1) *  factor_county_stretch) |>
  ungroup()

datayear_levels <- data_for_plot_tmp  |>
  distinct(datayear) |>
  arrange(datayear) |>
  mutate(datayear_id = (row_number() - 1) *  factor_datayear_stretch)

variable_levels <- data_for_plot_tmp  |>
  distinct(variable) |>
  arrange(variable) |>
  mutate(variable_id = (row_number() - 1) *  factor_variable_stretch)

data_for_plot <- data_for_plot_tmp |>
  inner_join(county_levels,
             by = c("county", "region")) |>
  inner_join(datayear_levels,
             by = "datayear") |>
  inner_join(variable_levels,
             by = "variable") |>
  mutate(
    y_level = county_id + datayear_id + variable_id,
    group = case_when(
      variable %in% c("BIRTHS", "DEATHS", "NATURALINC")      ~ "natural",
      str_detect(variable, "MIG")                            ~ "migration",
      TRUE                                                   ~ "ERROR"
    )
  ) |>
  mutate(
    label_county = if_else(datayear_id == 16 & variable_id == 0, # TODO: generalize: 16 is 3rd of 5 datayear values
                           county,
                           "")
  ) |>
  group_by(region) |>
  mutate(x_value_label = min(pmin(xstart, xend)) - 0.01) |>
  ungroup()

data_for_plot |>
  ggplot() +
  geom_vline(xintercept = 0, lty = 2, linewidth = 0.5, color = "firebrick", alpha = 0.25) +
  geom_segment(aes(x = xstart, xend = xend, y = y_level, yend = y_level, color = variable),
               linewidth = 0.5) +
  geom_text(aes(x = x_value_label,
                y = y_level,
                label = label_county),
            hjust = 1, na.rm = TRUE, size = 2.5) +
  facet_wrap( ~ region, nrow = 1) +
  scale_x_continuous(labels = label_percent()) +
  scale_y_reverse() +
  scale_color_viridis_d(option = "H", end = 0.8) +
  coord_cartesian(xlim = c(-0.20, 0.17)) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "top") +
  guides(color = guide_legend(override.aes = list(linewidth = 3))) +
  
  labs(
    title = "Population changes in NC counties in the five years 2015-2019",
    subtitle = "Percent of county population; sorted within region by net population chage in this period",
    x = "Change in county population since 2014",
    y = NULL,
    color = NULL,
    caption = my_caption
  )

```

<br>

Here too some patterns are easier to see when the counties are put in approximate geographical position (@fig-estimates-five-year-four-vars-geofacet).

Observations

* Most of agricultural counties growing the fastest border a county containing of the larger cities: (1) Chatham, Franklin, Harnett, and Johnston surround Wake; (2) Iredell (spill-over from Charlotte); (3) Henderson and Haywood (near Asheville / Buncombe county)

* Exceptions: (2) Hoke (spill-over from Fayetteville in Cumberland); (2) Watauga (contains App State University and the town of Boone); (3) Jackson (not obvious why it's growing; perhaps benefiting from Buncombe-Haywood growth?)

* Onslow, Hoke and Cumberland experienced largest net natural increase:

```{r}
data_for_table <- multi_year_nc_counties |>
  filter(str_detect(variable, "RNATURALINC")) |>
  group_by(GEOID, county) |>
  summarize(avg_yearly_rate = mean(value)) |>
  ungroup() |>
  arrange(desc(avg_yearly_rate)) |>
  select(-GEOID) |>
  left_join(
    county_region_df,
    by = "county"
  ) |>
  select(county, region, avg_yearly_rate)

bind_rows(
  data_for_table |> head(10),
  data_for_table |> tail(10)
) |>
  gt() |>
  tab_header(md(glue("***Ten counties with greatest positive and negative<br>percent natural increase 2015-2019***",
                     "<br>Simple average of natural increase rate"))) |>
  fmt_percent(columns = avg_yearly_rate,
              decimals = 1)
  
```

<br>

* Brunswick, Chatham, and Pender experienced largest percent increases due net migration:

```{r}
data_for_table <- multi_year_nc_counties |>
  filter(str_detect(variable, "^RNETMIG")) |>
  group_by(GEOID, county) |>
  summarize(avg_yearly_rate = mean(value)) |>
  ungroup() |>
  arrange(desc(avg_yearly_rate)) |>
  select(-GEOID) |>
  left_join(
    county_region_df,
    by = "county"
  ) |>
  select(county, region, avg_yearly_rate)

bind_rows(
  data_for_table |> head(10),
  data_for_table |> tail(10)
) |>
  gt() |>
  tab_header(md(glue("***Ten counties with greatest positive and negative<br>percent net migration 2015-2019***",
                     "<br>Simple average of net migration rate"))) |>
  fmt_percent(columns = avg_yearly_rate,
              decimals = 1)
  
```

<br>

```{r fig.height=10, fig.width=24}
#| label: fig-estimates-five-year-four-vars-geofacet
#| fig-cap: "Population change in NC counties over the five years 2015-2019 (approximate geographical position)"
#| fig-height: 10
#| fig-width: 24
#| column: screen-inset-right

factor_county_stretch <- 40
factor_datayear_stretch <- 8
factor_variable_stretch <- 2

data_for_plot_tmp <- multi_year_nc_counties |>
  filter(!str_detect(variable, "^R"),
         !variable %in% c("NETMIG", "NATURALINC")) |>
  inner_join(multi_year_nc_county_pop,
               by = c("GEOID", "county", "datayear")) |>
  mutate(variable = fct_drop(variable),
         value = if_else(variable == "DEATHS", 
                         -1 * value,
                         value),
         value_pct_pop = value / pop_total
         ) |>
  arrange(datayear, variable) |>
  inner_join(county_region_df,
             by = "county") |>
  group_by(GEOID, county, region) |>
  mutate(
    xend = cumsum(value_pct_pop),
    xstart = lag(xend, default = 0),
    row_id = row_number()
    ) |>
  ungroup() |>
  arrange(GEOID, county, datayear, variable) 

county_levels <- data_for_plot_tmp  |>
  distinct(county) |>
  mutate(county_id = 0)  # since each county is its own facet

datayear_levels <- data_for_plot_tmp  |>
  distinct(datayear) |>
  arrange(datayear) |>
  mutate(datayear_id = (row_number() - 1) *  factor_datayear_stretch)

variable_levels <- data_for_plot_tmp  |>
  distinct(variable) |>
  arrange(variable) |>
  mutate(variable_id = (row_number() - 1) *  factor_variable_stretch)

data_for_plot <- data_for_plot_tmp |>
  inner_join(county_levels,
             by = c("county")) |>
  inner_join(datayear_levels,
             by = "datayear") |>
  inner_join(variable_levels,
             by = "variable") |>
  mutate(
    y_level = county_id + datayear_id + variable_id,
    group = case_when(
      variable %in% c("BIRTHS", "DEATHS", "NATURALINC")      ~ "natural",
      str_detect(variable, "MIG")                            ~ "migration",
      TRUE                                                   ~ "ERROR"
    )
  ) |>
  mutate(
    label_datayear = if_else(variable_id == 2,
                           as.character(datayear),
                           "")
  ) |>
  # group_by(region) |>
  # mutate(x_value_label = min(pmin(xstart, xend)) - 0.02) |>
  # ungroup() |>
  mutate(county = if_else(county == "McDowell", "Mcdowell", county))


data_for_plot |>
  ggplot() +
  geom_vline(xintercept = 0, lty = 2, linewidth = 0.5, color = "firebrick", alpha = 0.25) +
  geom_hline(aes(yintercept = if_else(variable_id == 0,
                                      datayear_id - 1,
                                      NA_real_)),
             na.rm = TRUE, linewidth = 0.1, alpha = 0.5) +
  geom_segment(aes(x = xstart, xend = xend, y = y_level, yend = y_level, color = variable),
               linewidth = 2) +
  geom_text(aes(x = -0.10,  #x_value_label,
                y = y_level,
                label = label_datayear),
            hjust = 1, na.rm = TRUE, size = 2) +
  scale_x_continuous(labels = label_percent()) +
  scale_y_continuous(expand = expansion(mult = c(0, .02))) +
  coord_cartesian(xlim = c(-0.16, NA)) +
  scale_color_viridis_d(option = "H", end = 0.8) +
  facet_geo( ~ county,
             scales = "free_y",
             grid = "us_nc_counties_grid1") +
  theme(panel.grid = element_blank(),
        axis.title.y = element_text(hjust = 1),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.title = element_text(size = rel(3.0), face = "bold"),
        legend.position = "top") +
  labs(title = "NC county population change 2015-2019",
       subtitle = "Percent of county population change from births, deaths, domestic migration, and international migration",
       x = NULL,
       y = NULL,
       color = NULL,
       caption = my_caption)

```

<br>

For ease of reference, below are the same plots made larger and placed in alphabetical order.

<br>

```{r fig.height=16, fig.width=20}
#| label: fig-estimates-five-year-four-vars-facet
#| fig-cap: "Population change in NC counties over the five years 2015-2019"
#| fig-height: 16
#| fig-width: 20
#| column: screen-inset-right


data_for_plot |>
  ggplot() +
  geom_vline(xintercept = 0, lty = 2, linewidth = 0.5, color = "firebrick", alpha = 0.25) +
  geom_hline(aes(yintercept = if_else(variable_id == 0,
                                      datayear_id - 1,
                                      NA_real_)),
             na.rm = TRUE, linewidth = 0.1, alpha = 0.5) +
  geom_segment(aes(x = xstart, xend = xend, y = y_level, yend = y_level, color = variable),
               linewidth = 2) +
  geom_text(aes(x = -0.09,
                y = y_level,
                label = label_datayear),
            hjust = 1, na.rm = TRUE, size = 2) +
  scale_x_continuous(labels = label_percent()) +
  scale_y_continuous(expand = expansion(mult = c(0, .02))) +
  coord_cartesian(xlim = c(-0.12, NA)) +
  scale_color_viridis_d(option = "H", end = 0.8) +
  facet_wrap( ~ county, ncol = 10,) +
  theme(panel.grid = element_blank(),
        axis.title.y = element_text(hjust = 1),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.title = element_text(size = rel(3.0), face = "bold"),
        legend.position = "top") +
  labs(title = "NC county population change 2015-2019",
       subtitle = "Percent of county population change from births, deaths, domestic migration, and international migration",
       x = NULL,
       y = NULL,
       color = NULL,
       caption = my_caption)

```

<br>

## Population under 18

From 1970 to 1990 there was a slight association between a larger portion under 18 and smaller counties. Starting in 2000 there was a strong move in the other direction, and the association has only gotten stronger since: counties with a larger proportion of children tend to have larger populations.

<br>

```{r fig.width=8, fig.height=12}
#| label: fig-nc-counties-pop-by-growth
#| fig-cap: "Under 18 population has become increasingly associated with larger county population"
#| fig-width: 8
#| fig-height: 12
#| column: screen-inset-right
#| warning: false

pop_u18 |>
  rename(pop_total = pop,
         pct_u18 = pop_pct_u18,
         pct_18p = pop_pct_18p) |>
  pivot_longer(cols = starts_with("pct_"),
               names_to = "group", 
               values_to = "value") |> 
  mutate(year = factor(year),
         group = fct_reorder(group, value, sum)) |>
  mutate(name_label = as.character(county)) |>
  ggplot(aes(value, pop_total, color = region)) + #
  geom_hline(yintercept = 0, lty = 2, linewidth = 0.5, alpha = 0.5) +
  geom_point(alpha = 0.6, size = 2, na.rm = TRUE) +
  scale_y_log10(labels = label_number(scale_cut = cut_short_scale())) +
  scale_x_continuous(labels = percent_format()) +
  scale_color_manual(values = region_colors$color) +
  geom_smooth(formula = y ~ x, method = "lm", se = FALSE, color = "black", linewidth = 0.5, na.rm = TRUE) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = "top") +
  facet_grid(year ~ group, scales = "free_x") +
  labs(
    title = "Percent under 18 has become strongly associated\nwith larger county population",
    subtitle = glue("Portion of the population under 18 and 18+ (1970-2020)"),
    y = "Population (log10 scale)",
    x = "Portion of the population",
    color = "Region",
    caption = my_caption_nhgis
  )
```

<br>

While growth rates are interesting, they hide changes in the mix (the proportion of one group compared to the others). Since the 18 and over group has been growing faster than the under 18 group as seen in @fig-summary_population-facet-grid-1, the proportion of youth has been declining; as noted earlier, this is a national trend as families size declined.

In every NC county, the *proportion* of the population under 18 declined since 1970--and except for counties the urban crescent and at the coast, the *count* declined in most counties in the 10 years before 2020 (see @fig-summary-population-growth-rates-facet-grid-2). In most counties the biggest drops in proportion occurred in the 10 years to 1980 and 1990 (and maybe in earlier decades not included in this dataset). In the ~5% of cases in which the proportion of the population under 18 grew between two decennial censuses, arrows in the plot below point to the right, completing a loop.

<br>

```{r fig.height=14, fig.width=10}
#| label: fig-nc-counties-portion-u18-change
#| fig-cap: "Portion of NC county population under 18 (1970 - 2020)"
#| fig-width: 10
#| fig-height: 14
#| column: screen-inset-right

xx_path <- d_pop_growth_u18_long |>
  filter(variable == "pop_pct_u18") |>
  mutate(name_label = as.character(county)) |>
  select(year, county, variable, value, positive_growth) |>
  arrange(county, year) |>
  mutate(year_to = lead(year),
         value_to = lead(value),
         positive_growth_to = lead(positive_growth)) |>
  filter(year_to != 1970) |>
  mutate(year = factor(year))

county_levels <- pop_u18 |>
  filter(year == 2020) |>
  arrange(pop_pct_u18) |>
  pull(county)

region_region_levels <- pop_u18 |>
  filter(year == 2020) |>
  arrange(pop_pct_u18) |>
  pull(region) |>
  unique() |>
  as.character()

region_color_levels <- pop_u18 |>
  filter(year == 2020) |>
  arrange(pop_pct_u18) |>
  left_join(region_colors,
            by = "region") |>
  pull(color) |>
  unique() |>
  as.character()

pop_u18_for_plot <- pop_u18 |>
  mutate(region = factor(region, levels = region_region_levels),
         ) |>
  group_by(county) |>
  mutate(change_in_portion = last(pop_pct_u18) - first(pop_pct_u18),
         change_since_1970 = pop_pct_u18 / first(pop_pct_u18)) |>
  ungroup() |>
  mutate(county = factor(county, levels = county_levels),
         curve = even_or_odd(year / 10)
  ) |>
  arrange(county) |>
  left_join(region_colors,
            by = "region") |>
  mutate(region = factor(as.character(region), levels = region_region_levels),
         color = factor(color, levels = region_color_levels)
  ) |>
  mutate(year = factor(year))

nudge <- pop_u18_for_plot |>
  arrange(county) |>
  distinct(county, .keep_all = TRUE) |>
  mutate(row = row_number(),
         nudge = if_else(even_or_odd(row) == 1, .05, 0)) |>
  pull(nudge)

change_in_pct_u18 <- pop_u18_for_plot |>
  mutate(year = as.numeric(as.character(year))) |>
  filter(year %in% c(1970, 2020)) |>
  group_by(year) |>
  summarize(ave_pop_pct_u18 = mean(pop_pct_u18)) |>
  ungroup()

pop_u18_for_plot |>
  ggplot() +
  geom_point(data = pop_u18_for_plot |> filter(year == 2020),
             aes(pop_pct_u18, county, size = pop_u18, color = year),
             alpha = 0.3, show.legend = FALSE) +
  geom_point(data = pop_u18_for_plot |> filter(year %in% c(1970, 2020)),
             aes(pop_pct_u18, county, color = year),
             alpha = 0.8, shape = "|", size = 6, show.legend = FALSE) +
  geom_curve(data = xx_path,
             aes(x = value, xend = value_to, y = county, yend = county, color = year), curvature = 0.1, #curve * 0.2
            arrow = arrow(length = unit(1, "mm"), type = "closed"), linewidth = 0.5) +
  geom_segment(data = pop_u18_for_plot |>
                 #add_row(nc_total) |>
                 filter(as.character(year) == "2020") |>
                 group_by(county) |>
                 mutate(xend = min(pop_pct_u18) - 0.02,
                        county = factor(county, levels = county_levels)) |>
                 ungroup() |>
                 arrange(county) |>
                 filter((row_number() / 10)  == floor(row_number() / 10)),
               aes(x = 0, xend = xend, y = county, yend = county),
               linewidth = 0.1, alpha = 0.3
               ) +
  geom_label(data = pop_u18_for_plot |> filter(year == 2020),
            aes(x = 0, y = county, label = county, fill = region),
            hjust = 0, size = 3, nudge_x = nudge, show.legend = TRUE, 
            label.size = 0, alpha = 0.3) +
  scale_x_continuous(labels = percent_format()) +
  scale_size_continuous(range = c(1, 15), labels = label_number(scale_cut = cut_short_scale())) +
  scale_color_viridis_d(end = 0.85) +
  scale_fill_manual(values =  region_color_levels,
                     labels = region_region_levels) +
  expand_limits(x = 0, y = length(county_levels) + 1) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.title = element_text(size = rel(2.0), face = "bold"),
        legend.position = c(0.35, 0.8)) + #c(0.3, 0.3)
  guides(fill = guide_legend(alpha = 0.7),
         size = "none"
  ) +
  labs(title = "Portion of NC county population under 18 (1970 - 2020)",
       subtitle = glue("Change each decennial census; ordered by percent of population under 18 in 2020", 
                       #"\nAverage {percent(avg_change_pct_u18)} reduction in portion of the population under 18",
                       "\nAverage portion under 18 dropped from {percent(change_in_pct_u18$ave_pop_pct_u18[[1]])}",
                         " to {percent(change_in_pct_u18$ave_pop_pct_u18[[2]])}",
                         " between {change_in_pct_u18$year[[1]]} and {change_in_pct_u18$year[[2]]} (not weighted by population)",
                       "\nCircles indicate relative population in 2020"),
       x = NULL,
       y = NULL,
       color = "Year",
       fill = "Region",
       caption = my_caption_nhgis)

```

<br>