30_republican_disaffiliation.r

library(tidyverse)
library(haven)
library(car)
library(binom)
library(ggtext)
library(sjlabelled)
library(clipr)

# ── Load data ──────────────────────────────────────────────────────────────────

cumulative <- readRDS("./data/CES/cumulative_2006-2024.rds")

# Load 2025 individual year data and harmonize to match cumulative columns
# NOTE: verify variable names against data/CES/codebooks/CES25_Common.pdf
ces25_raw <- haven::read_dta("./data/CES/data/2025.dta")

# Inspect column names to confirm mapping before running
# names(ces25_raw)

ces25 <- ces25_raw |>
  transmute(
    year         = 2025,
    case_id      = caseid,
    weight       = commonweight,
    pid3         = pid3,       # 1=Dem, 2=Rep, 3=Ind, 4=Other, 5=Not sure
    pid7         = pid7,       # 1=Strong Dem ... 7=Strong Rep
    # pid3_leaner is not in the individual year file — derive from pid7:
    # 1-3 = Democrat+Lean, 4 = Independent, 5-7 = Republican+Lean
    pid3_leaner  = case_when(
      pid7 %in% 1:3 ~ 1L,
      pid7 == 4     ~ 3L,
      pid7 %in% 5:7 ~ 2L,
      TRUE          ~ NA_integer_
    ),
    religion     = as.numeric(religpew),    # 3 = Mormon/Latter-day Saint
    relig_church = as.numeric(pew_churatd),
    relig_imp    = as.numeric(pew_religimp),
    gender4      = as.numeric(gender4),
    prayer      = as.numeric(pew_prayer),
    # harmonize to binary gender to match cumulative's `gender` (1=Male, 2=Female)
    gender       = case_when(
      as.numeric(gender4) == 1 ~ 1L,
      as.numeric(gender4) == 2 ~ 2L,
      TRUE ~ NA_integer_
    ),
    race         = as.numeric(race),
    inputstate   = as.numeric(inputstate),
    educ         = as.numeric(educ),
    birthyr      = as.numeric(birthyr),
    marstat      = as.numeric(marstat),
    geography = case_when(
      inputstate == 49 ~ "Utah",
      !is.na(inputstate) ~ "Rest of U.S.",
      TRUE ~ NA_character_
  )
  )

# Cumulative stores some columns as sjlabelled factors — convert to numeric
# so bind_rows can combine them with the plain numeric columns from ces25
cumulative <- cumulative |>
  mutate(across(c(religion, relig_church, relig_imp, pid3, pid3_leaner, educ, gender, race, marstat),
                \(x) suppressWarnings(as.numeric(x)))) |>
  mutate(geography = case_when(
    state == "Utah" ~ "Utah",
    !is.na(state) ~ "Rest of U.S.",
    TRUE ~ NA_character_
  ))

cumulative <- cumulative |>
  mutate(id = paste0(year, "X", case_id)) |>
  left_join(read.csv("./data/CES/prayer.csv"), by = "id") |>
  mutate(year = year.x)

data <- bind_rows(cumulative, ces25)

# ── Recodes ────────────────────────────────────────────────────────────────────

data <- data |>
  mutate(
    lds = religion == 3,

    # Strict party ID (no leaners folded in)
    party = case_when(
      pid3 == 1 ~ "Democrat",
      pid3 == 2 ~ "Republican",
      pid3 == 3 ~ "Independent",
      TRUE ~ NA_character_
    ),
    party = factor(party, levels = c("Democrat", "Independent", "Republican")),

    # Party ID with leaners folded in (pid3_leaner: 1=Dem+Lean, 2=Rep+Lean, 3=True Ind)
    party_leaner = case_when(
      pid3_leaner == 1 ~ "Democrat+Lean",
      pid3_leaner == 2 ~ "Republican+Lean",
      pid3_leaner == 3 ~ "Independent",
      TRUE ~ NA_character_
    ),
    party_leaner = factor(party_leaner,
      levels = c("Democrat+Lean", "Independent", "Republican+Lean")),

    # Net Republican advantage using leaners (the "20-point shift" measure)
    net_rep = case_when(
      pid3_leaner == 2 ~  1,
      pid3_leaner == 3 ~  0,
      pid3_leaner == 1 ~ -1,
      TRUE ~ NA_real_
    ),

    # Church attendance
    relig_church_rc = car::recode(
      as.numeric(relig_church),
      "1:2 = 'Weekly or More'; 3:7 = 'Less than Weekly'"
    ),
    relig_church_rc = factor(relig_church_rc,
      levels = c("Weekly or More", "Less than Weekly")),

    # Prayer frequency
    prayer_rc = case_when(
      prayer == 1             ~ "Several Times a Day",
      prayer %in% 2:7         ~ "Less Than Several Times a Day",
      TRUE                    ~ NA_character_
    ),

    # Religious engagement clusters
    # 1. Devout Traditionalists: prayer several times/day + religion very important + weekly attendance
    # 2. Adaptive Believers: weekly OR monthly attenders not meeting all three devout criteria
    # 3. Cultural Mormons: attend less than monthly (few times/year, seldom, never)
    cluster = case_when(
      prayer_rc == "Several Times a Day" &
        as.numeric(relig_imp) == 1 &
        relig_church_rc == "Weekly or More"          ~ "Devout Traditionalist",
      as.numeric(relig_church) %in% 1:3              ~ "Adaptive Believer",
      as.numeric(relig_church) %in% 4:7              ~ "Cultural Mormon",
      TRUE                                           ~ NA_character_
    ),
    cluster = factor(cluster,
      levels = c("Devout Traditionalist", "Adaptive Believer", "Cultural Mormon")),

    # Birth cohort (for pseudo-panel analysis)
    birth_cohort = case_when(
      birthyr < 1950        ~ "Pre-1950",
      birthyr >= 1950 & birthyr < 1965 ~ "1950-1964",
      birthyr >= 1965 & birthyr < 1981 ~ "1965-1980",
      birthyr >= 1981 & birthyr < 1997 ~ "1981-1996",
      birthyr >= 1997       ~ "1997+",
      TRUE ~ NA_character_
    ),
    birth_cohort = factor(birth_cohort,
      levels = c("Pre-1950", "1950-1964", "1965-1980", "1981-1996", "1997+")),

    # Gender
    gender_rc = case_when(
      as.numeric(gender) == 1 ~ "Male",
      as.numeric(gender) == 2 ~ "Female",
      TRUE ~ NA_character_
    ),
    gender_rc = factor(gender_rc, levels = c("Male", "Female")),

    # Race
    race_rc = case_when(
      as.numeric(race) == 1 ~ "White",
      as.numeric(race) %in% 2:8 ~ "Non-White",
      TRUE ~ NA_character_
    ),
    race_rc = factor(race_rc, levels = c("White", "Non-White")),

    # Marital status (marstat: 1=Married, 2:6=Not married)
    married_rc = case_when(
      marstat == 1         ~ "Married",
      marstat %in% 2:6     ~ "Not Married",
      TRUE                 ~ NA_character_
    ),
    married_rc = factor(married_rc, levels = c("Married", "Not Married")),

    # Geography: Utah/Idaho ("Jello Belt") vs. rest of U.S.
    geography = factor(geography, levels = c("Utah", "Rest of U.S."))
  )

lds <- data |> filter(lds, pid3_leaner %in% 1:3)

# ── Analysis 1: Year-by-year party ID trend ────────────────────────────────────
# The core "20-point shift" — look for a 2016 discontinuity

trend <- lds |>
  group_by(year, party_leaner) |>
  summarise(
    n_raw = n(),
    n_wtd = sum(weight, na.rm = TRUE),
    .groups = "drop"
  ) |>
  group_by(year) |>
  mutate(
    prop  = n_wtd / sum(n_wtd),
    se    = sqrt(prop * (1 - prop) / sum(n_raw)),
    lo    = prop - 1.96 * se,
    hi    = pmin(prop + 1.96 * se, 1)
  ) |>
  ungroup()

plot <- trend |>
  ggplot(aes(x = year, y = prop, color = party_leaner, group = party_leaner,
             fill = party_leaner)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2, stroke = 1.5) +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey60") +
  annotate("text", x = 2016.2, y = 0.7, label = "Trump\nnomination",
       hjust = 0, size = 3, color = "grey40", family = "Cairo") +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    limits = c(0, 1),
    breaks = seq(0, 1, by = 0.1)
  ) +
  scale_color_manual(values = c(
    "Democrat+Lean"   = "dodgerblue3",
    "Independent"     = "purple",
    "Republican+Lean" = "darkred"
  )) +
  scale_fill_manual(values = c(
    "Democrat+Lean"   = "dodgerblue3",
    "Independent"     = "purple",
    "Republican+Lean" = "darkred"
  ), guide = "none") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  labs(
    title = "Latter-day Saint Party Identification, 2006–2025",
    subtitle = "Leaners counted with their leaned party; shaded bands = 95% CI",
    x = NULL, y = "Share", color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_1.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 2: Net Republican advantage over time ─────────────────────────────
# Cleaner single-line view of the 20-point shift

net_trend <- lds |>
  group_by(year) |>
  summarise(
    net    = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n      = sum(!is.na(net_rep)),
    se     = sd(net_rep, na.rm = TRUE) / sqrt(n),
    lo     = net - 1.96 * se,
    hi     = pmin(net + 1.96 * se, 1)
  )

plot <- net_trend |>
  ggplot(aes(x = year, y = net)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "#d6604d", alpha = 0.15) +
  geom_line(linewidth = 1.3, color = "#d6604d") +
  geom_point(size = 2.5, color = "#d6604d") +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey50") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "LDS Net Republican Advantage, 2006–2025",
    subtitle = "Positive = more Republican than Democrat; negative = more Democrat",
    x = NULL, y = "Net Republican advantage",
    caption = "@mormon_metrics\nSource: Cooperative Election Study"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_2.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 3: Church attendance as moderator ─────────────────────────────────
# Does the shift happen among devout LDS or only among nominal members?

attend_trend <- lds |>
  filter(!is.na(relig_church_rc)) |>
  group_by(year, relig_church_rc) |>
  summarise(
    net = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n   = n(),
    se  = sd(net_rep, na.rm = TRUE) / sqrt(n),
    lo  = net - 1.96 * se,
    hi  = pmin(net + 1.96 * se, 1),
    .groups = "drop"
  )

plot <- attend_trend |>
  ggplot(aes(x = year, y = net, color = relig_church_rc, fill = relig_church_rc,
             group = relig_church_rc)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2) +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey50") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_color_manual(values = c("Weekly or More" = "#1a6faf", "Less than Weekly" = "#b2182b")) +
  scale_fill_manual(values = c("Weekly or More" = "#1a6faf", "Less than Weekly" = "#b2182b"),
                    guide = "none") +
  labs(
    title = "LDS Republican Disaffiliation by Church Attendance",
    subtitle = "If the shift is concentrated among low-attenders, it's largely a compositional story",
    x = NULL, y = "Net Republican advantage", color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_3.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 4: Cohort pseudo-panel ───────────────────────────────────────────
# Which birth cohorts are driving the shift?
# If it's generational replacement, older cohorts stay Republican,
# younger cohorts enter more Democratic.
# If cohorts themselves drift, that's within-cohort attitudinal change.

cohort_trend <- lds |>
  filter(!is.na(birth_cohort)) |>
  group_by(year, birth_cohort) |>
  summarise(
    net = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n   = n(),
    .groups = "drop"
  ) |>
  filter(n >= 30)   # drop cells with too few respondents

plot <- cohort_trend |>
  ggplot(aes(x = year, y = net, color = birth_cohort, group = birth_cohort)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_line(linewidth = 1.1) +
  geom_point(size = 2) +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey50") +
  scale_x_continuous(breaks = seq(2006, 2025, 4)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "LDS Republican Disaffiliation by Birth Cohort",
    subtitle = "Pseudo-panel: tracking the same birth cohorts across survey years",
    x = NULL, y = "Net Republican advantage", color = "Birth cohort",
    caption = "@mormon_metrics\nSource: Cooperative Election Study. Cells with n < 30 suppressed."
  ) +
  theme_minimal() +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_4.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 5: Demographic composition of LDS over time ──────────────────────
# Is the LDS pool in the CES changing? Compositional shift = disaffiliation story.

# Age composition
plot <- lds |>
  filter(!is.na(birth_cohort)) |>
  group_by(year, birth_cohort) |>
  summarise(n = sum(weight, na.rm = TRUE), .groups = "drop") |>
  group_by(year) |>
  mutate(prop = n / sum(n)) |>
  ggplot(aes(x = year, y = prop, fill = birth_cohort)) +
  geom_area(position = "stack") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Age Composition of LDS Respondents Over Time",
    subtitle = "Shifts in cohort composition would indicate disaffiliation or demographic change",
    x = NULL, y = "Share", fill = "Birth cohort",
    caption = "@mormon_metrics\nSource: Cooperative Election Study"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_5.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 6: Share of sample identifying as LDS ────────────────────────────
# If the LDS share of the CES drops over time, that suggests disaffiliation
# is reducing the pool — and selectively if politics of stayers differ from leavers.

lds_share <- data |>
  filter(!is.na(religion), pid3_leaner %in% 1:3) |>
  group_by(year) |>
  summarise(
    share_lds = weighted.mean(religion == 3, w = weight, na.rm = TRUE),
    n_total   = n(),
    se        = sqrt(share_lds * (1 - share_lds) / n_total),
    lo        = share_lds - 1.96 * se,
    hi        = pmin(share_lds + 1.96 * se, 1)
  )

plot <- lds_share |>
  ggplot(aes(x = year, y = share_lds)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "#4d9221", alpha = 0.15) +
  geom_line(linewidth = 1.2, color = "#4d9221") +
  geom_point(size = 2.5, color = "#4d9221") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 0.1)) +
  labs(
    title = "LDS Share of CES Sample Over Time",
    subtitle = "A declining share would be consistent with net disaffiliation from the Church",
    x = NULL, y = "Share identifying as LDS",
    caption = "Source: Cooperative Election Study"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_6.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 7: Party ID by gender ────────────────────────────────────────────

gender_trend <- lds |>
  filter(!is.na(gender_rc)) |>
  group_by(year, gender_rc) |>
  summarise(
    net = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n   = n(),
    se  = sd(net_rep, na.rm = TRUE) / sqrt(n),
    lo  = net - 1.96 * se,
    hi  = pmin(net + 1.96 * se, 1),
    .groups = "drop"
  )

plot <- gender_trend |>
  ggplot(aes(x = year, y = net, color = gender_rc, fill = gender_rc,
             group = gender_rc)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2, stroke = 1.5) +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey60") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_color_manual(values = c("Male" = "dodgerblue3", "Female" = "#f26c82")) +
  scale_fill_manual(values = c("Male" = "dodgerblue3", "Female" = "#f26c82"),
                    guide = "none") +
  labs(
    title = "LDS Republican Disaffiliation by Gender",
    subtitle = "Are men and women moving at the same rate?",
    x = NULL, y = "Net Republican advantage", color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study"
  ) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_7.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 8: Party ID by race ──────────────────────────────────────────────

race_trend <- lds |>
  filter(!is.na(race_rc)) |>
  group_by(year, race_rc) |>
  summarise(
    net = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n   = n(),
    se  = sd(net_rep, na.rm = TRUE) / sqrt(n),
    lo  = net - 1.96 * se,
    hi  = pmin(net + 1.96 * se, 1),
    .groups = "drop"
  ) |>
  filter(n >= 30)

plot <- race_trend |>
  ggplot(aes(x = year, y = net, color = race_rc, fill = race_rc,
             group = race_rc)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2, stroke = 1.5) +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey60") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_color_manual(values = c("White" = "#4d9221", "Non-White" = "#d6604d")) +
  scale_fill_manual(values = c("White" = "#4d9221", "Non-White" = "#d6604d"),
                    guide = "none") +
  labs(
    title = "LDS Republican Disaffiliation by Race",
    subtitle = "White vs. non-White Latter-day Saints",
    x = NULL, y = "Net Republican advantage", color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study. Cells with n < 30 suppressed."
  ) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_8.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 9: Party ID by geography ─────────────────────────────────────────

geo_trend <- lds |>
  filter(!is.na(geography)) |>
  group_by(year, geography) |>
  summarise(
    net = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n   = n(),
    se  = sd(net_rep, na.rm = TRUE) / sqrt(n),
    lo  = net - 1.96 * se,
    hi  = pmin(net + 1.96 * se, 1),
    .groups = "drop"
  ) |>
  filter(n >= 30)

plot <- geo_trend |>
  ggplot(aes(x = year, y = net, color = geography, fill = geography,
             group = geography)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2, stroke = 1.5) +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey60") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0,1)) +
  scale_color_manual(values = c("Utah" = "#1a6faf", "Rest of U.S." = "#b2182b")) +
  scale_fill_manual(values = c("Utah" = "#1a6faf", "Rest of U.S." = "#b2182b"),
                    guide = "none") +
  labs(
    title = "LDS Republican Disaffiliation by Geography",
    subtitle = "Utah vs. the rest of the U.S.",
    x = NULL, y = "Net Republican advantage", color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study. Cells with n < 30 suppressed."
  ) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text.y             = element_text(size = 10),
    axis.text.x             = element_text(size = 10, margin = margin(t = -5)),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.title            = element_blank(),
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_9.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 10: LDS church attendance over time ───────────────────────────────
# Among those who identify as LDS, are they becoming more nominal?
# A declining weekly attendance rate = internal disaffiliation even among identifiers.

attend_share <- data |>
  filter(lds, !is.na(relig_church)) |>
  group_by(year) |>
  summarise(
    prop = weighted.mean(as.numeric(relig_church) %in% 1:2, w = weight, na.rm = TRUE),
    n    = n(),
    se   = sqrt(prop * (1 - prop) / n),
    lo   = pmax(prop - 1.96 * se, 0),
    hi   = pmin(prop + 1.96 * se, 1)
  )

plot <- attend_share |>
  ggplot(aes(x = year, y = prop)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "#1a6faf", alpha = 0.15) +
  geom_line(linewidth = 1.3, color = "#1a6faf") +
  geom_point(size = 2.5, stroke = 1.5, color = "#1a6faf") +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey60") +
  annotate("text", x = 2016.2, y = 0.1, label = "Trump\nnomination",
           hjust = 0, size = 3, color = "grey40", family = "Cairo") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, 1), breaks = seq(0, 1, by = 0.1)) +
  labs(
    title = "LDS Weekly Church Attendance Over Time",
    subtitle = "Share of LDS identifiers attending church weekly or more; shaded band = 95% CI",
    x = NULL, y = "Share attending weekly or more",
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_10.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 11: LDS church attendance over time by political party ────────────
# Are Democrats/leaners becoming less active, or is it uniform across parties?

attend_party <- data |>
  filter(lds, !is.na(relig_church), !is.na(party_leaner)) |>
  group_by(year, party_leaner) |>
  summarise(
    prop = weighted.mean(as.numeric(relig_church) %in% 1:2, w = weight, na.rm = TRUE),
    n    = n(),
    se   = sqrt(prop * (1 - prop) / n),
    lo   = pmax(prop - 1.96 * se, 0),
    hi   = pmin(prop + 1.96 * se, 1),
    .groups = "drop"
  ) |>
  filter(n >= 30)

plot <- attend_party |>
  ggplot(aes(x = year, y = prop, color = party_leaner, fill = party_leaner,
             group = party_leaner)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2, stroke = 1.5) +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey60") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, 1), breaks = seq(0, 1, by = 0.1)) +
  scale_color_manual(values = c(
    "Democrat+Lean"   = "dodgerblue3",
    "Independent"     = "purple",
    "Republican+Lean" = "darkred"
  )) +
  scale_fill_manual(values = c(
    "Democrat+Lean"   = "dodgerblue3",
    "Independent"     = "purple",
    "Republican+Lean" = "darkred"
  ), guide = "none") +
  labs(
    title = "LDS Weekly Church Attendance by Political Party",
    subtitle = "Share attending weekly or more; shaded bands = 95% CI",
    x = NULL, y = "Share attending weekly or more", color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_11.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 12: Party ID by age group over time ──────────────────────────────
# Unlike the cohort chart, this tracks fixed age brackets over time.
# Shows whether young LDS are structurally more Democratic, or whether
# each age group has shifted.

age_trend <- lds |>
  filter(!is.na(birthyr)) |>
  mutate(
    age_at_survey = year - birthyr,
    age_group = case_when(
      age_at_survey >= 18 & age_at_survey <= 30 ~ "18-30",
      age_at_survey >= 31 & age_at_survey <= 45 ~ "31-45",
      age_at_survey >= 46 & age_at_survey <= 65 ~ "46-65",
      age_at_survey >= 66                        ~ "66+",
      TRUE ~ NA_character_
    ),
    age_group = factor(age_group, levels = c("18-30", "31-45", "46-65", "66+"))
  ) |>
  filter(!is.na(age_group)) |>
  group_by(year, age_group) |>
  summarise(
    net = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n   = n(),
    se  = sd(net_rep, na.rm = TRUE) / sqrt(n),
    lo  = net - 1.96 * se,
    hi  = pmin(net + 1.96 * se, 1),
    .groups = "drop"
  ) |>
  filter(n >= 30)

plot <- age_trend |>
  ggplot(aes(x = year, y = net, color = age_group, fill = age_group,
             group = age_group)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2, stroke = 1.5) +
  geom_vline(xintercept = 2016, linetype = "dashed", color = "grey60") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_color_manual(values = c(
    "18-30" = "#A8CDE9",
    "31-45" = "#74ADD1",
    "46-65" = "#4575B4",
    "66+"   = "#313695"
  )) +
  scale_fill_manual(values = c(
    "18-30" = "#A8CDE9",
    "31-45" = "#74ADD1",
    "46-65" = "#4575B4",
    "66+"   = "#313695"
  ), guide = "none") +
  labs(
    title = "LDS Party Identification by Age Group, 2006–2025",
    subtitle = "Net Republican advantage by age bracket; shaded bands = 95% CI",
    x = NULL, y = "Net Republican advantage", color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study. Cells with n < 30 suppressed."
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_12.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 13: Age group composition of LDS respondents over time ────────────
# Is the LDS sample getting older or younger? Shifts in composition help explain
# overall party ID trends — e.g., if younger (more Democratic) cohorts are growing.

age_comp <- lds |>
  filter(!is.na(birthyr)) |>
  mutate(
    age_at_survey = year - birthyr,
    age_group = case_when(
      age_at_survey >= 18 & age_at_survey <= 30 ~ "18-30",
      age_at_survey >= 31 & age_at_survey <= 45 ~ "31-45",
      age_at_survey >= 46 & age_at_survey <= 65 ~ "46-65",
      age_at_survey >= 66                        ~ "66+",
      TRUE ~ NA_character_
    ),
    age_group = factor(age_group, levels = c("18-30", "31-45", "46-65", "66+"))
  ) |>
  filter(!is.na(age_group)) |>
  group_by(year, age_group) |>
  summarise(n = sum(weight, na.rm = TRUE), .groups = "drop") |>
  group_by(year) |>
  mutate(prop = n / sum(n)) |>
  ungroup()

plot <- age_comp |>
  ggplot(aes(x = year, y = prop, fill = age_group)) +
  geom_area(position = "stack") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     breaks = seq(0, 1, by = 0.1)) +
  scale_fill_manual(values = c(
    "18-30" = "#A8CDE9",
    "31-45" = "#74ADD1",
    "46-65" = "#4575B4",
    "66+"   = "#313695"
  )) +
  labs(
    title = "Age Composition of LDS Respondents Over Time",
    subtitle = "Share of LDS sample in each age bracket by survey year",
    x = NULL, y = "Share", fill = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_13.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 14: Racial composition of LDS respondents over time ───────────────

race_comp <- lds |>
  filter(!is.na(race_rc)) |>
  group_by(year, race_rc) |>
  summarise(n = sum(weight, na.rm = TRUE), .groups = "drop") |>
  group_by(year) |>
  mutate(prop = n / sum(n)) |>
  ungroup()

plot <- race_comp |>
  ggplot(aes(x = year, y = prop, fill = race_rc)) +
  geom_area(position = "stack") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     breaks = seq(0, 1, by = 0.1)) +
  scale_fill_manual(values = c("White" = "#4d9221", "Non-White" = "#d6604d")) +
  labs(
    title = "Racial Composition of LDS Respondents Over Time",
    subtitle = "Share of LDS sample identifying as White vs. Non-White by survey year",
    x = NULL, y = "Share", fill = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_14.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Kitagawa Decomposition ─────────────────────────────────────────────────────
# Decomposes the change in net Republican advantage into:
#   - Composition effect: how much changed because group SIZES changed
#   - Rate effect: how much changed because within-group ATTITUDES changed

kitagawa <- function(data, group_var, early_years = 2006:2010, late_years = 2020:2025) {
  
  group_name <- deparse(substitute(group_var))  # get string name upfront
  
  early <- data |>
    filter(year %in% early_years, !is.na({{ group_var }})) |>
    group_by({{ group_var }}) |>
    summarise(
      share_early = sum(weight, na.rm = TRUE),
      rate_early  = weighted.mean(net_rep, w = weight, na.rm = TRUE),
      .groups = "drop"
    ) |>
    mutate(share_early = share_early / sum(share_early))
  
  late <- data |>
    filter(year %in% late_years, !is.na({{ group_var }})) |>
    group_by({{ group_var }}) |>
    summarise(
      share_late = sum(weight, na.rm = TRUE),
      rate_late  = weighted.mean(net_rep, w = weight, na.rm = TRUE),
      .groups = "drop"
    ) |>
    mutate(share_late = share_late / sum(share_late))
  
  decomp <- early |>
    left_join(late, by = group_name) |>  # plain string works here
    mutate(
      composition = (share_late - share_early) * ((rate_early + rate_late) / 2),
      rate        = ((share_early + share_late) / 2) * (rate_late - rate_early)
    )
  
  list(
    detail = decomp,
    summary = decomp |>
      summarise(
        composition_effect = sum(composition),
        rate_effect        = sum(rate),
        total_explained    = sum(composition) + sum(rate)
      )
  )
}

# Run for each dimension
decomp_attendance <- kitagawa(lds, relig_church_rc)
decomp_race       <- kitagawa(lds, race_rc)
decomp_gender     <- kitagawa(lds, gender_rc)
decomp_age        <- kitagawa(lds, birth_cohort)

# Print summaries
cat("=== ATTENDANCE ===\n"); print(decomp_attendance$summary)
cat("=== RACE ===\n");       print(decomp_race$summary)
cat("=== GENDER ===\n");     print(decomp_gender$summary)
cat("=== COHORT ===\n");     print(decomp_age$summary)

# ── Plot the decomposition results ────────────────────────────────────────────

decomp_results <- bind_rows(
  decomp_attendance$summary |> mutate(dimension = "Church Attendance"),
  decomp_race$summary       |> mutate(dimension = "Race"),
  decomp_gender$summary     |> mutate(dimension = "Gender"),
  decomp_age$summary        |> mutate(dimension = "Age Cohort")
) |>
  pivot_longer(
    cols = c(composition_effect, rate_effect),
    names_to = "effect_type",
    values_to = "effect"
  ) |>
  mutate(
    effect_type = case_when(
      effect_type == "composition_effect" ~ "Composition Effect\n(who's in the pool)",
      effect_type == "rate_effect"        ~ "Rate Effect\n(attitude change)"
    ),
    dimension = factor(dimension, 
                       levels = c("Church Attendance", "Race", "Age Cohort", "Gender"))
  )

# Also compute total observed change for reference line
total_change <- lds |>
  filter(year %in% c(2006:2010)) |>
  summarise(net = weighted.mean(net_rep, w = weight, na.rm = TRUE)) |>
  pull(net) -
  lds |>
  filter(year %in% c(2020:2025)) |>
  summarise(net = weighted.mean(net_rep, w = weight, na.rm = TRUE)) |>
  pull(net)

plot <- decomp_results |>
  ggplot(aes(x = dimension, y = effect, fill = effect_type)) +
  geom_col(position = "dodge") +
  geom_hline(yintercept = 0, color = "grey40") +
  geom_hline(yintercept = -total_change, linetype = "dashed", 
             color = "grey40", linewidth = 0.8) +
  annotate("text", x = 0.5, y = -total_change + 0.01, 
           label = "Total observed\nshift", hjust = 0, 
           size = 3, color = "grey40", family = "Cairo") +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c(
    "Composition Effect\n(who's in the pool)" = "#d6604d",
    "Rate Effect\n(attitude change)"           = "#4575b4"
  )) +
  labs(
    title    = "Decomposing the LDS Republican Shift, 2006–2025",
    subtitle = "How much of the shift is who's in the pool vs. attitude change within groups?",
    x = NULL, y = "Contribution to net Republican shift", fill = NULL,
    caption  = "@mormon_metrics\nSource: Cooperative Election Study"
  ) +
  theme(
    axis.ticks              = element_blank(),
    axis.title              = element_text(size = 14),
    axis.text               = element_text(size = 10),
    plot.background         = element_rect(fill = "grey95", color = NA),
    panel.background        = element_blank(),
    text                    = element_text(face = "bold", family = "Cairo"),
    plot.title              = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle           = element_text(size = 11, hjust = 0.5),
    plot.title.position     = "plot",
    plot.subtitle.position  = "plot",
    legend.background       = element_blank(),
    legend.box.background   = element_blank(),
    legend.key              = element_blank(),
    legend.position         = "top",
    panel.grid              = element_blank(),
    panel.grid.major.y      = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption            = element_text(size = 8, family = "Cairo")
  )

ggsave("./images/30_republican_disaffiliation_decomp_15.png", plot,
       width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 16: Cluster composition over time ────────────────────────────────

cluster_comp <- lds |>
  filter(!is.na(cluster)) |>
  group_by(year, cluster) |>
  summarise(n = sum(weight, na.rm = TRUE), .groups = "drop") |>
  group_by(year) |>
  mutate(prop = n / sum(n)) |>
  ungroup()

plot <- cluster_comp |>
  ggplot(aes(x = year, y = prop, fill = cluster)) +
  geom_area(position = "stack") +
  scale_x_continuous(breaks = seq(2006, 2025, 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     breaks = seq(0, 1, by = 0.1)) +
  scale_fill_manual(values = c(
    "Devout Traditionalist" = "#619CFF",
    "Adaptive Believer"     = "#F8766D",
    "Cultural Mormon"       = "#00BA38"
  )) +
  labs(
    title = "LDS Religious Engagement Clusters Over Time",
    subtitle = "Share of LDS respondents in each cluster by survey year",
    x = NULL, y = "Share", fill = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_16.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 17: Net Republican advantage by cluster over time ─────────────────

cluster_trend <- lds |>
  filter(!is.na(cluster)) |>
  mutate(
    period = case_when(
      year %in% 2008:2012 ~ "2008-12",
      year %in% 2013:2016 ~ "2013-16",
      year %in% 2017:2020 ~ "2017-20",
      year %in% 2021:2025 ~ "2021-25",
      TRUE ~ NA_character_
    ),
    period = factor(period, levels = c("2008-12", "2013-16", "2017-20", "2021-25"))
  ) |>
  filter(!is.na(period)) |>
  group_by(period, cluster) |>
  summarise(
    net = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n   = n(),
    se  = sd(net_rep, na.rm = TRUE) / sqrt(n),
    lo  = net - 1.96 * se,
    hi  = pmin(net + 1.96 * se, 1),
    .groups = "drop"
  ) |>
  filter(n >= 30)

plot <- cluster_trend |>
  ggplot(aes(x = period, y = net, color = cluster, group = cluster)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.2, linewidth = 0.8) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3, stroke = 1.5) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_color_manual(values = c(
    "Devout Traditionalist" = "#619CFF",
    "Adaptive Believer"     = "#F8766D",
    "Cultural Mormon"       = "#00BA38"
  )) +
  labs(
    title = "LDS Republican Disaffiliation by Religious Cluster",
    subtitle = "Net Republican advantage by five-year period; error bars = 95% CI",
    x = NULL, y = "Net Republican advantage", color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_17.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 18: Devout Traditionalist share vs. Republican ID over time ───────

period_levels <- c("2008-12", "2013-16", "2017-20", "2021-25")

add_period <- function(df) {
  df |> mutate(
    period = case_when(
      year %in% 2008:2012 ~ "2008-12",
      year %in% 2013:2016 ~ "2013-16",
      year %in% 2017:2020 ~ "2017-20",
      year %in% 2021:2025 ~ "2021-25",
      TRUE ~ NA_character_
    ),
    period = factor(period, levels = period_levels)
  ) |> filter(!is.na(period))
}

year2_levels <- c("2006-07", "2008-09", "2010-11", "2012-13", "2014-15",
                  "2016-17", "2018-19", "2020-21", "2022-23", "2024-25")

dt_trend <- lds |>
  filter(!is.na(cluster)) |>
  add_period() |>
  group_by(period) |>
  summarise(
    value = weighted.mean(cluster == "Devout Traditionalist", w = weight, na.rm = TRUE),
    n     = n(),
    se    = sqrt(value * (1 - value) / n),
    lo    = pmax(value - 1.96 * se, 0),
    hi    = pmin(value + 1.96 * se, 1),
    metric = "Devout Traditionalist Share",
    .groups = "drop"
  )

rep_trend <- lds |>
  add_period() |>
  group_by(period) |>
  summarise(
    value = weighted.mean(net_rep, w = weight, na.rm = TRUE),
    n     = n(),
    se    = sd(net_rep, na.rm = TRUE) / sqrt(n),
    lo    = value - 1.96 * se,
    hi    = pmin(value + 1.96 * se, 1),
    metric = "Net Republican Advantage",
    .groups = "drop"
  )

combo <- bind_rows(dt_trend, rep_trend) |>
  mutate(metric = factor(metric,
    levels = c("Devout Traditionalist Share", "Net Republican Advantage")))

# Compute R-squared using two-year groupings
add_year2 <- function(df) {
  df |> mutate(year2 = case_when(
    year %in% 2006:2007 ~ "2006-07",
    year %in% 2008:2009 ~ "2008-09",
    year %in% 2010:2011 ~ "2010-11",
    year %in% 2012:2013 ~ "2012-13",
    year %in% 2014:2015 ~ "2014-15",
    year %in% 2016:2017 ~ "2016-17",
    year %in% 2018:2019 ~ "2018-19",
    year %in% 2020:2021 ~ "2020-21",
    year %in% 2022:2023 ~ "2022-23",
    year %in% 2024:2025 ~ "2024-25",
    TRUE ~ NA_character_
  ))
}

dt_y2 <- lds |>
  filter(!is.na(cluster)) |>
  add_year2() |>
  filter(!is.na(year2)) |>
  group_by(year2) |>
  summarise(dt_share = weighted.mean(cluster == "Devout Traditionalist",
                                     w = weight, na.rm = TRUE), .groups = "drop")

rep_y2 <- lds |>
  add_year2() |>
  filter(!is.na(year2)) |>
  group_by(year2) |>
  summarise(net_rep_avg = weighted.mean(net_rep, w = weight, na.rm = TRUE),
            .groups = "drop")

r_sq <- dt_y2 |>
  inner_join(rep_y2, by = "year2") |>
  summarise(r2 = cor(dt_share, net_rep_avg)^2) |>
  pull(r2)

plot <- combo |>
  ggplot(aes(x = period, y = value, color = metric, group = metric)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.2, linewidth = 0.8) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3, stroke = 1.5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     breaks = seq(0, 1, by = 0.1)) +
  scale_color_manual(values = c(
    "Devout Traditionalist Share" = "#619CFF",
    "Net Republican Advantage"    = "darkred"
  )) +
  labs(
    title = "Devout Traditionalists and Republican ID Among LDS",
    subtitle = "Share of LDS who are Devout Traditionalists vs. net Republican advantage",
    x = NULL, y = NULL, color = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_18.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Analysis 19: Cluster breakdown by age group ────────────────────────────────

age_cluster <- lds |>
  filter(!is.na(cluster), !is.na(birthyr)) |>
  mutate(
    age_at_survey = year - birthyr,
    age_group = case_when(
      age_at_survey >= 18 & age_at_survey <= 30 ~ "18-30",
      age_at_survey >= 31 & age_at_survey <= 45 ~ "31-45",
      age_at_survey >= 46 & age_at_survey <= 65 ~ "46-65",
      age_at_survey >= 66                        ~ "66+",
      TRUE ~ NA_character_
    ),
    age_group = factor(age_group, levels = c("18-30", "31-45", "46-65", "66+"))
  ) |>
  filter(!is.na(age_group)) |>
  group_by(age_group, cluster) |>
  summarise(n = sum(weight, na.rm = TRUE), .groups = "drop") |>
  group_by(age_group) |>
  mutate(prop = n / sum(n)) |>
  ungroup()

plot <- age_cluster |>
  ggplot(aes(x = age_group, y = prop, fill = cluster)) +
  geom_col(position = "stack") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     breaks = seq(0, 1, by = 0.1)) +
  scale_fill_manual(values = c(
    "Devout Traditionalist" = "#619CFF",
    "Adaptive Believer"     = "#F8766D",
    "Cultural Mormon"       = "#00BA38"
  )) +
  labs(
    title = "LDS Religious Cluster Breakdown by Age Group",
    subtitle = "Share of each age group falling into each religious engagement cluster",
    x = NULL, y = "Share", fill = NULL,
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 10, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.title           = element_blank(),
    legend.background      = element_blank(),
    legend.box.background  = element_blank(),
    legend.key             = element_blank(),
    legend.position        = "top",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_19.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# ── Chart 20: % Devout Traditionalist over time by education ───────────────────

educ_levels <- c("No College", "Some College / 2-Year", "4-Year Degree", "Postgraduate")

plot <- lds |>
  mutate(
    educ_rc = case_when(
      educ %in% 1:3 ~ "No College",
      educ == 4     ~ "Some College / 2-Year",
      educ == 5     ~ "4-Year Degree",
      educ == 6     ~ "Postgraduate",
      TRUE          ~ NA_character_
    ),
    educ_rc = factor(educ_rc, levels = educ_levels)
  ) |>
  add_period() |>
  group_by(educ_rc, period) |>
  count(cluster, wt = weight) |>
  drop_na() |>
  mutate(
    prop     = n / sum(n),
    total_n  = sum(n),
    lower    = binom.confint(x = n, n = total_n, method = "asymptotic")$lower,
    upper    = pmin(binom.confint(x = n, n = total_n, method = "asymptotic")$upper, 1)
  ) |>
  filter(cluster == "Devout Traditionalist") |>
  ggplot(aes(x = period, y = prop, color = educ_rc, group = educ_rc)) +
  geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  geom_text(aes(label = scales::percent(prop, accuracy = 1)),
            vjust = -2, size = 3, family = "Cairo", fontface = "bold", color = "black") +
  facet_grid(~ educ_rc) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, 0.8), breaks = seq(0, 0.8, by = 0.1)) +
  labs(
    title    = "Devout Traditionalists Over Time by Education",
    subtitle = "% of LDS identifying as Devout Traditionalist\n(pray several times daily, religion very important, attend church weekly+)",
    x = NULL, y = "Proportion",
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 7, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.position        = "none",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    strip.background       = element_blank(),
    strip.text             = element_text(size = 9, face = "bold", family = "Cairo"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_20.png", plot, width = 2400, height = 1500, units = "px", dpi = 300)

# ── Chart 21: % Devout Traditionalist over time by gender × marital status ─────

gender_marital_levels <- c("Married Male", "Married Female", "Single Male", "Single Female")

plot <- lds |>
  mutate(
    gender_marital = case_when(
      gender_rc == "Male"   & married_rc == "Married"     ~ "Married Male",
      gender_rc == "Female" & married_rc == "Married"     ~ "Married Female",
      gender_rc == "Male"   & married_rc == "Not Married" ~ "Single Male",
      gender_rc == "Female" & married_rc == "Not Married" ~ "Single Female",
      TRUE ~ NA_character_
    ),
    gender_marital = factor(gender_marital, levels = gender_marital_levels)
  ) |>
  add_period() |>
  group_by(gender_marital, period) |>
  count(cluster, wt = weight) |>
  drop_na() |>
  mutate(
    prop    = n / sum(n),
    total_n = sum(n),
    lower   = binom.confint(x = n, n = total_n, method = "asymptotic")$lower,
    upper   = pmin(binom.confint(x = n, n = total_n, method = "asymptotic")$upper, 1)
  ) |>
  filter(cluster == "Devout Traditionalist") |>
  ggplot(aes(x = period, y = prop, color = gender_marital, group = gender_marital)) +
  geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  geom_text(aes(label = scales::percent(prop, accuracy = 1)),
            vjust = -2, size = 3, family = "Cairo", fontface = "bold", color = "black") +
  facet_grid(~ gender_marital) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, 0.8), breaks = seq(0, 0.8, by = 0.1)) +
  labs(
    title    = "Devout Traditionalists Over Time by Gender and Marital Status",
    subtitle = "% of LDS identifying as Devout Traditionalist\n(pray several times daily, religion very important, attend church weekly+)",
    x = NULL, y = "Proportion",
    caption = "@mormon_metrics\nSource: Cooperative Election Study cumulative file + 2025 individual year data"
  ) +
  theme(
    axis.ticks             = element_blank(),
    axis.title             = element_text(size = 14),
    axis.text.y            = element_text(size = 10),
    axis.text.x            = element_text(size = 7, margin = margin(t = -5)),
    plot.background        = element_rect(fill = "grey95", color = NA),
    panel.background       = element_blank(),
    text                   = element_text(face = "bold", family = "Cairo"),
    plot.title             = element_text(size = 15, face = "bold", hjust = 0.5),
    plot.subtitle          = element_text(size = 11, hjust = 0.5),
    plot.title.position    = "plot",
    plot.subtitle.position = "plot",
    legend.position        = "none",
    panel.grid             = element_blank(),
    panel.grid.major.y     = element_line(color = "lightgrey", linetype = "solid"),
    strip.background       = element_blank(),
    strip.text             = element_text(size = 9, face = "bold", family = "Cairo"),
    plot.caption           = element_text(size = 8, family = "Cairo")
  )
ggsave("./images/30_republican_disaffiliation_21.png", plot, width = 2400, height = 1500, units = "px", dpi = 300)