33_lds_migration.r

library(tidyverse)
library(haven)
library(showtext)
library(magick)

font_add("Cairo_bold", "~/Library/Fonts/Cairo-Bold.ttf")
font_add("Cairo", "~/Library/Fonts/Cairo-Regular.ttf")
showtext::showtext_auto(enable = TRUE)
showtext::showtext_opts(dpi = 300)

data <- haven::read_sav("./data/GSS_spss/gss7224_r3.sav")

lds_base <- data |>
  filter(other == 64) |>
  mutate(
    reg16_num  = as.numeric(reg16_7222),
    region_num = as.numeric(region_7222),
    age_num    = as.numeric(age),
    educ_num   = as.numeric(educ)
  ) |>
  filter(
    reg16_num %in% 1:9,
    region_num %in% 1:9,
    !is.na(wtssps),
    !is.na(age_num)
  )

# Census region centroids and labels
region_coords <- tibble(
  reg   = 1:9,
  label = c("New\nEngland", "Mid\nAtlantic", "E. N.\nCentral", "W. N.\nCentral",
            "S.\nAtlantic", "E. S.\nCentral", "W. S.\nCentral", "Mountain", "Pacific"),
  lon   = c(-71.5, -76.0, -84.5, -96.0, -79.5, -86.5, -94.0, -111.5, -121.0),
  lat   = c( 44.5,  40.5,  43.0,  44.5,  32.5,  33.0,  31.0,   44.5,   44.5),
  hjust = c(  0.5,   0.5,   0.5,   0.5,   0.5,   0.5,   0.5,    0.5,    0.5),
  vjust = c( -0.7,  -0.7,  -0.7,  -0.7,   1.5,   1.5,   1.5,   -0.7,   -0.7)
)

# State → 9 Census division lookup
state_census <- tibble(
  state = c(
    "connecticut", "maine", "massachusetts", "new hampshire", "rhode island", "vermont",
    "new jersey", "new york", "pennsylvania",
    "illinois", "indiana", "michigan", "ohio", "wisconsin",
    "iowa", "kansas", "minnesota", "missouri", "nebraska", "north dakota", "south dakota",
    "delaware", "florida", "georgia", "maryland", "north carolina", "south carolina",
      "virginia", "west virginia",
    "alabama", "kentucky", "mississippi", "tennessee",
    "arkansas", "louisiana", "oklahoma", "texas",
    "arizona", "colorado", "idaho", "montana", "new mexico", "nevada", "utah", "wyoming",
    "california", "oregon", "washington"
  ),
  census_reg = c(
    rep(1, 6), rep(2, 3), rep(3, 5), rep(4, 7), rep(5, 8),
    rep(6, 4), rep(7, 4), rep(8, 8), rep(9, 3)
  )
)

region_pal <- c(
  "1" = "#CCEBC5", "2" = "#FED9A6", "3" = "#FFFFCC", "4" = "#FBB4AE",
  "5" = "#DECBE4", "6" = "#B3E2CD", "7" = "#FDCDAC", "8" = "#F4CAE4",
  "9" = "#E6F5C9"
)

us_map <- map_data("state") |>
  left_join(state_census, by = c("region" = "state"))

make_map <- function(lds, age_label) {
  flows <- lds |>
    group_by(from = reg16_num, to = region_num) |>
    summarise(weight = sum(wtssps, na.rm = TRUE), .groups = "drop") |>
    mutate(pct = weight / sum(weight)) |>
    filter(from != to, pct >= 0.005)

  flow_plot <- flows |>
    left_join(region_coords |> select(reg, from_lon = lon, from_lat = lat), by = c("from" = "reg")) |>
    left_join(region_coords |> select(reg, to_lon   = lon, to_lat   = lat), by = c("to"   = "reg")) |>
    mutate(to_mountain = to == 8)

  ggplot() +
    geom_polygon(
      data = us_map,
      aes(x = long, y = lat, group = group, fill = factor(census_reg)),
      color = "white", linewidth = 0.3
    ) +
    scale_fill_manual(values = region_pal, guide = "none") +
    geom_curve(
      data      = flow_plot,
      aes(x = from_lon, y = from_lat, xend = to_lon, yend = to_lat,
          linewidth = pct, alpha = pct, color = to_mountain),
      curvature = 0.4,
      arrow     = arrow(length = unit(0.12, "inches"), type = "closed")
    ) +
    scale_color_manual(values = c("FALSE" = "#4169e1", "TRUE" = "#e8543a"), guide = "none") +
    geom_text(
      data  = region_coords,
      aes(x = lon, y = lat, label = label, hjust = hjust, vjust = vjust),
      size       = 4.5, family = "Cairo", lineheight = 0.85, fontface = "bold"
    ) +
    scale_linewidth_continuous(range = c(0.6, 5), guide = "none") +
    scale_alpha_continuous(range = c(0.35, 0.95), guide = "none") +
    coord_fixed(ratio = 1.3, xlim = c(-125, -65), ylim = c(22, 50)) +
    labs(
      title    = paste0("LDS Migration — ", age_label),
      subtitle = "Arrow thickness = share of inter-regional moves  |  Orange = into Mountain region",
      caption  = "@mormon_metrics | Source: General Social Survey (GSS) 1972–2022"
    ) +
    theme(
      axis.ticks             = element_blank(),
      axis.text              = element_blank(),
      axis.title             = element_blank(),
      plot.background        = element_rect(fill = "grey95", color = "black", linewidth = 0.5),
      panel.background       = element_blank(),
      text                   = element_text(family = "Cairo", size = 10),
      plot.title             = element_text(size = 17, hjust = 0.5, family = "Cairo_bold"),
      plot.subtitle          = element_text(size = 12, hjust = 0.5, family = "Cairo_bold"),
      plot.title.position    = "plot",
      plot.subtitle.position = "plot",
      panel.grid             = element_blank(),
      plot.caption           = element_text(size = 9, family = "Cairo"),
      plot.margin            = margin(10, 10, 10, 10)
    )
}

# Overall map — standalone image only
p_overall <- make_map(lds_base, "All Ages")
ggsave("./images/33_lds_migration_overall.png", p_overall, width = 2000, height = 1400, units = "px", dpi = 300)

# Age-group maps — GIF only
age_groups <- list(
  list(label = "Under 40", data = lds_base |> filter(age_num < 40)),
  list(label = "40 and Over",  data = lds_base |> filter(age_num >= 40))
)

png_files <- character(length(age_groups))

for (i in seq_along(age_groups)) {
  grp   <- age_groups[[i]]
  p     <- make_map(grp$data, grp$label)
  fname <- sprintf("./images/33_lds_migration_%d.png", i)
  ggsave(fname, p, width = 2000, height = 1400, units = "px", dpi = 300)
  png_files[i] <- fname
}

# Combine PNGs into a GIF — 500 centiseconds = 5 seconds per frame, loop forever
frames <- image_read(png_files)
gif    <- image_animate(frames, delay = 300, loop = 0)
image_write(gif, "./images/33_lds_migration.gif")

# Education-group maps — GIF only. educ is highest year of school; 16 = 4 years
# of college, the threshold for holding a bachelor's degree.
educ_base <- lds_base |> filter(!is.na(educ_num))

educ_groups <- list(
  list(label = "No Bachelor's Degree", data = educ_base |> filter(educ_num < 16)),
  list(label = "Bachelor's Degree or More", data = educ_base |> filter(educ_num >= 16))
)

educ_png_files <- character(length(educ_groups))

for (i in seq_along(educ_groups)) {
  grp   <- educ_groups[[i]]
  p     <- make_map(grp$data, grp$label)
  fname <- sprintf("./images/33_lds_migration_educ_%d.png", i)
  ggsave(fname, p, width = 2000, height = 1400, units = "px", dpi = 300)
  educ_png_files[i] <- fname
}

educ_frames <- image_read(educ_png_files)
educ_gif    <- image_animate(educ_frames, delay = 300, loop = 0)
image_write(educ_gif, "./images/33_lds_migration_educ.gif")

# Stay-vs-leave rates by religion ------------------------------------------
# For each religion, what share of members still live in the Census region
# they lived in at age 16 ("stayed") vs. moved to a different region ("left").
relig_base <- data |>
  mutate(
    reg16_num  = as.numeric(reg16_7222),
    region_num = as.numeric(region_7222),
    educ_num   = as.numeric(educ),
    denom_lbl  = as.character(as_factor(denom)),
    relig_grp  = case_when(
      other == 64                  ~ "Latter-day Saint",
      relig == 1 & denom_lbl == "baptist"      ~ "Baptist",
      relig == 1 & denom_lbl == "methodist"    ~ "Methodist",
      relig == 1 & denom_lbl == "lutheran"     ~ "Lutheran",
      relig == 1 & denom_lbl == "presbyterian" ~ "Presbyterian",
      relig == 1 & denom_lbl == "episcopal"    ~ "Episcopal",
      relig == 1                   ~ "Other Protestant",
      relig == 2                   ~ "Catholic",
      relig == 3                   ~ "Jewish",
      relig == 4                   ~ "No Religion",
      TRUE                         ~ NA_character_
    )
  ) |>
  filter(
    reg16_num %in% 1:9,
    region_num %in% 1:9,
    !is.na(wtssps)
  )

relig_rates <- relig_base |>
  mutate(stayed = reg16_num == region_num) |>
  filter(!is.na(relig_grp)) |>
  group_by(relig_grp) |>
  summarise(
    stay_rate = sum(wtssps[stayed], na.rm = TRUE) / sum(wtssps, na.rm = TRUE),
    .groups   = "drop"
  ) |>
  mutate(leave_rate = 1 - stay_rate)

overall_stay <- relig_base |>
  summarise(stay_rate = sum(wtssps[reg16_num == region_num], na.rm = TRUE) /
                        sum(wtssps, na.rm = TRUE)) |>
  mutate(relig_grp = "U.S. Overall", leave_rate = 1 - stay_rate)

relig_rates <- bind_rows(relig_rates, overall_stay)

relig_order <- relig_rates |> arrange(stay_rate) |> pull(relig_grp)

relig_plot_df <- relig_rates |>
  pivot_longer(c(stay_rate, leave_rate), names_to = "outcome", values_to = "rate") |>
  mutate(
    relig_grp = factor(relig_grp, levels = relig_order),
    outcome   = factor(outcome, levels = c("leave_rate", "stay_rate"),
                       labels = c("Left home region", "Stayed in home region"))
  )

mm_theme <- theme(
  axis.ticks              = element_blank(),
  axis.title              = element_blank(),
  axis.text.y             = element_text(size = 10, family = "Cairo"),
  axis.text.x             = element_blank(),
  plot.background         = element_rect(fill = "grey95", color = "black", linewidth = 0.5),
  panel.background        = element_blank(),
  text                    = element_text(family = "Cairo", size = 10),
  plot.title              = element_text(size = 17, hjust = 0.5, family = "Cairo_bold"),
  plot.subtitle           = element_text(size = 12, hjust = 0.5, family = "Cairo_bold"),
  plot.title.position     = "plot",
  plot.subtitle.position  = "plot",
  legend.background       = element_blank(),
  legend.box.background   = element_blank(),
  legend.position         = "top",
  legend.title            = element_blank(),
  legend.text             = element_text(size = 10, family = "Cairo"),
  panel.grid              = element_blank(),
  plot.caption            = element_text(size = 9, family = "Cairo")
)

p_relig <- ggplot(relig_plot_df, aes(x = rate, y = relig_grp, fill = outcome)) +
  geom_col(width = 0.7) +
  geom_text(
    aes(label = scales::percent(rate, accuracy = 1)),
    position = position_stack(vjust = 0.5),
    family = "Cairo_bold", size = 3.5, color = "white"
  ) +
  scale_fill_manual(values = c("Left home region" = "#e8543a",
                               "Stayed in home region" = "#4169e1")) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1), expand = c(0, 0)) +
  labs(
    title    = "Who Stays Close to Home?",
    subtitle = "Share still living in their age-16 Census region, by religion",
    x = NULL, y = NULL,
    caption  = "@mormon_metrics | Source: General Social Survey (GSS) 1972–2022"
  ) +
  mm_theme

ggsave("./images/33_lds_migration_religion.png", p_relig,
       width = 1500, height = 1500, units = "px", dpi = 300)

# College-educated rate by religion ----------------------------------------
# Share of each religion holding a bachelor's degree (educ >= 16).
educ_rates <- relig_base |>
  filter(!is.na(educ_num), !is.na(relig_grp)) |>
  group_by(relig_grp) |>
  summarise(
    college_rate = sum(wtssps[educ_num >= 16], na.rm = TRUE) /
                   sum(wtssps, na.rm = TRUE),
    .groups = "drop"
  )

overall_educ <- relig_base |>
  filter(!is.na(educ_num)) |>
  summarise(college_rate = sum(wtssps[educ_num >= 16], na.rm = TRUE) /
                           sum(wtssps, na.rm = TRUE)) |>
  mutate(relig_grp = "U.S. Overall")

educ_rates <- bind_rows(educ_rates, overall_educ) |>
  mutate(relig_grp = reorder(relig_grp, college_rate))

p_educ_relig <- ggplot(educ_rates, aes(x = college_rate, y = relig_grp)) +
  geom_col(width = 0.7, fill = "#4169e1") +
  geom_text(
    aes(label = scales::percent(college_rate, accuracy = 1)),
    hjust = -0.15, family = "Cairo_bold", size = 3.5
  ) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1),
                     expand = expansion(mult = c(0, 0.12))) +
  labs(
    title    = "Who Holds a College Degree?",
    subtitle = "Share with a bachelor's degree or more, by religion",
    x = NULL, y = NULL,
    caption  = "@mormon_metrics | Source: General Social Survey (GSS) 1972–2022"
  ) +
  mm_theme

ggsave("./images/33_lds_migration_college.png", p_educ_relig,
       width = 1500, height = 1500, units = "px", dpi = 300)

# Where do Pacific leavers go? --------------------------------------------
# Among people who lived in the Pacific region (reg16 == 9) at age 16 and
# later moved to a different Census region, what share landed in each
# destination region? Broken down by U.S. overall vs. Mormons only.
region_names <- c(
  "1" = "New England", "2" = "Mid Atlantic", "3" = "E. N. Central",
  "4" = "W. N. Central", "5" = "S. Atlantic", "6" = "E. S. Central",
  "7" = "W. S. Central", "8" = "Mountain", "9" = "Pacific"
)

pacific_flow <- function(df, label) {
  df |>
    filter(reg16_num == 9, region_num != 9) |>
    group_by(to = region_num) |>
    summarise(weight = sum(wtssps, na.rm = TRUE), .groups = "drop") |>
    mutate(
      group       = label,
      destination = region_names[as.character(to)],
      pct         = weight / sum(weight)
    ) |>
    arrange(desc(pct))
}

pacific_us  <- pacific_flow(relig_base, "U.S. Overall")
pacific_lds <- pacific_flow(lds_base,   "Latter-day Saint")

pacific_destinations <- bind_rows(pacific_us, pacific_lds) |>
  select(group, destination, pct) |>
  mutate(pct = scales::percent(pct, accuracy = 0.1))

cat("\n=== Where Pacific leavers go ===\n")
cat("\n-- U.S. Overall --\n")
print(pacific_destinations |> filter(group == "U.S. Overall") |> select(-group), n = Inf)
cat("\n-- Latter-day Saint --\n")
print(pacific_destinations |> filter(group == "Latter-day Saint") |> select(-group), n = Inf)

# Side-by-side table PNG of Pacific-leaver destinations
pacific_table <- bind_rows(pacific_us, pacific_lds) |>
  select(group, destination, pct) |>
  pivot_wider(names_from = group, values_from = pct, values_fill = 0) |>
  arrange(desc(`U.S. Overall`))

dest_order <- pacific_table$destination
n_rows     <- nrow(pacific_table)

table_long <- pacific_table |>
  mutate(destination = factor(destination, levels = rev(dest_order))) |>
  pivot_longer(c(`U.S. Overall`, `Latter-day Saint`),
               names_to = "group", values_to = "pct") |>
  mutate(
    group = factor(group, levels = c("U.S. Overall", "Latter-day Saint")),
    label = scales::percent(pct, accuracy = 0.1)
  )

p_pacific_table <- ggplot(table_long, aes(x = group, y = destination)) +
  geom_tile(aes(fill = pct), color = "white", linewidth = 1) +
  geom_text(aes(label = label), family = "Cairo_bold", size = 4.2,
            color = ifelse(table_long$pct > 0.35, "white", "grey20")) +
  scale_fill_gradient(low = "#e8eefc", high = "#1f3f8f", guide = "none") +
  scale_x_discrete(position = "top", expand = expansion(mult = c(0.05, 0.05))) +
  scale_y_discrete(expand = expansion(mult = c(0.02, 0.02))) +
  labs(
    title    = "Where Do Pacific Leavers Go?",
    subtitle = "Destination share of Pacific leavers",
    x = NULL, y = NULL,
    caption  = "@mormon_metrics | Source: General Social Survey (GSS) 1972–2022"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_blank(),
    axis.text.y            = element_text(size = 11, family = "Cairo_bold",
                                          color = "grey20", hjust = 1),
    axis.text.x.top        = element_text(size = 12, family = "Cairo_bold", color = "grey20"),
    plot.background        = element_rect(fill = "grey95", color = "black", linewidth = 0.5),
    panel.background       = element_blank(),
    text                   = element_text(family = "Cairo", size = 10),
    plot.title             = element_text(size = 17, hjust = 0.5, family = "Cairo_bold"),
    plot.subtitle          = element_text(size = 12, hjust = 0.5, family = "Cairo_bold"),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    panel.grid             = element_blank(),
    plot.caption           = element_text(size = 9, family = "Cairo"),
    plot.margin            = margin(15, 20, 10, 20)
  )

ggsave("./images/33_lds_migration_pacific_table.png", p_pacific_table,
       width = 1500, height = 1500, units = "px", dpi = 300)

# gist: https://gist.github.com/acbass49/bb60eb6cf872b1f1125c6000d86f5452