Map to compare percentage of married population and percent of people who gave birth in past 12 mo to US average

birth_marriage_counties.R

library(tidycensus)
library(tidyverse)
library(sf)
library(tmap)

birth <- get_acs(geography = "county", year = 2023, table = "B13002", geometry = TRUE,
                 summary_var = "B13002_001")
marriage <- get_acs(geography = "county", year = 2023, table = "B12001",
                    summary_var = "B12001_001")

birth <- birth |> 
  filter(variable %in% c("B13002_002")) |> 
  pivot_wider(id_cols = c("GEOID", "NAME", "summary_est", "geometry"), names_from = variable, values_from = estimate) |> 
  mutate(birth_12_mo = B13002_002/summary_est) |> 
  mutate(avg_birth = mean(birth_12_mo, na.rm = TRUE)) |> 
  mutate(above_avg_birth = if_else(birth_12_mo > avg_birth, TRUE, FALSE))

marriage <- marriage |> 
  filter(variable %in% c(
    "B12001_004",
    "B12001_013"
  )) |> 
  group_by(GEOID) |> 
  reframe(GEOID, NAME, married = sum(estimate), summary_est) |> 
  unique() |> 
  mutate(pct_married = married / summary_est,
         avg_married = mean(pct_married, na.rm = TRUE)) |> 
  mutate(above_avg_married = if_else(pct_married > avg_married, TRUE, FALSE))

birth |> left_join(marriage, join_by(GEOID, NAME)) |> 
  mutate(above_avg_birth_married = if_else(above_avg_birth == TRUE & above_avg_married == TRUE, TRUE, FALSE)) |> 
  filter(!str_detect(NAME, "Hawaii|Alaska|Puerto")) |> 
  tm_shape() +
  tm_polygons(fill = "above_avg_birth_married",lwd = .5, col = "white",
              fill.legend = tm_legend(title = "Above average births and marriage")) +
  tm_title("US counties with birth and marriage rates, compared to national average") +
  tm_credits("Data: 2019-2023 American Community Survey 5-year estimates. Tables B13002, B12001")

Great map! I am unsure on the wording so I will replace/try the & (AND) with | (OR) (L31).

Edit: Now I think I am stupid, kids need to be born under marriage, in which century I think I am.