Skip to content

Instantly share code, notes, and snippets.

@acbass49

acbass49/31_CES_2025_update.R

Last active April 13, 2026 02:10
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save acbass49/fbb45d652ee6267214592a45abeed522 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save acbass49/fbb45d652ee6267214592a45abeed522 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
31_CES_2025_update.R
library(tidyverse)
library(haven)
library(car)
library(binom)
library(ggtext)
library(sjlabelled)
library(clipr)
library(showtext)
font_add("Cairo_bold", "~/Library/Fonts/Cairo-Bold.ttf")
font_add("Cairo", "~/Library/Fonts/Cairo-Regular.ttf")
showtext::showtext_auto(enable = TRUE)
# ── Load data ──────────────────────────────────────────────
cumulative <- readRDS("./data/CES/cumulative_2006-2024.rds")
ces25_raw <- haven::read_dta("./data/CES/data/2025.dta")
ces25 <- ces25_raw |>
transmute(
year = 2025,
case_id = caseid,
weight = commonweight,
pid3 = pid3,
pid7 = pid7,
pid3_leaner = case_when(
pid7 %in% 1:3 ~ 1L,
pid7 == 4 ~ 3L,
pid7 %in% 5:7 ~ 2L,
TRUE ~ NA_integer_
),
religion = as.numeric(religpew),
relig_church = as.numeric(pew_churatd),
relig_imp = as.numeric(pew_religimp),
gender4 = as.numeric(gender4),
prayer = as.numeric(pew_prayer),
relig_protestant = as.numeric(religpew_protestant),
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: convert labelled columns to numeric
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)
# ── Religion recode ────────────────────────────────────────────────────────────
data$religion_rc <- dplyr::case_when(
data$religion == 1 & data$relig_protestant == 1 ~ "Baptist",
data$religion == 1 & data$relig_protestant == 2 ~ "Methodist",
data$religion == 1 & data$relig_protestant == 3 ~ "Non-Denominational",
data$religion == 1 & data$relig_protestant == 4 ~ "Lutheran",
data$religion == 1 & data$relig_protestant == 5 ~ "Presbyterian",
data$religion == 1 & data$relig_protestant == 6 ~ "Pentecostal",
data$religion == 1 & data$relig_protestant == 7 ~ "Episcopalian",
data$religion == 1 & data$relig_protestant %in% c(8:99) ~ "Other Protestant",
data$religion == 3 ~ "Latter-day Saint",
data$religion == 2 ~ "Catholic",
data$religion == 4 ~ "Orthodox Christian",
data$religion == 5 ~ "Jewish",
data$religion == 6 ~ "Muslim",
data$religion == 7 ~ "Buddhist",
data$religion == 8 ~ "Hindu",
data$religion == 9 ~ "Atheist",
data$religion == 10 ~ "Agnostic",
data$religion == 11 ~ "Nothing in particular",
data$religion == 12 ~ "Other religion"
)
# ── Recodes ────────────────────────────────────────────────────────────────────
data <- data |>
mutate(
lds = religion == 3,
pid3_leaner_rc = case_when(
pid3_leaner == 1 ~ "Democrat+Lean",
pid3_leaner == 2 ~ "Republican+Lean",
pid3_leaner == 3 ~ "Independent",
TRUE ~ NA_character_
),
pid3_leaner_rc = factor(pid3_leaner_rc,
levels = c("Democrat+Lean", "Independent", "Republican+Lean")),
net_rep = case_when(
pid3_leaner == 2 ~ 1,
pid3_leaner == 3 ~ 0,
pid3_leaner == 1 ~ -1,
TRUE ~ NA_real_
),
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")),
geography = factor(geography, levels = c("Utah", "Rest of U.S."))
)
lds <- data |> filter(lds, pid3_leaner %in% 1:3)
# ── Shared theme ───────────────────────────────────────────────────────────────
mm_theme <- theme(
axis.ticks = element_blank(),
axis.title = element_text(size = 30, family = "Cairo_bold"),
axis.text.y = element_text(size = 20),
axis.text.x = element_text(size = 20, margin = margin(t = -5)),
plot.background = element_rect(fill = "grey95", color = NA),
panel.background = element_blank(),
text = element_text(family = "Cairo", size = 20),
plot.title = element_text(size = 45, hjust = 0.5, family = "Cairo_bold"),
plot.subtitle = element_text(size = 30, hjust = 0.5, family = "Cairo_bold"),
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 = 20, family = "Cairo")
)
# ── Analysis 1: Party identification among LDS over time ──────────────────────
trend <- lds |>
filter(year >= 2008) |>
group_by(year, pid3_leaner_rc) |>
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 = pid3_leaner_rc, group = pid3_leaner_rc,
fill = pid3_leaner_rc)) +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
geom_line(linewidth = 1.2) +
geom_point(size = 2, stroke = 1.5) +
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(2008, 2025, 2)) +
labs(
title = "Latter-day Saint Party Identification, 2008–2025",
subtitle = "Leaners counted with their leaned party; shaded bands = 95% CI",
x = NULL, y = "Share", color = NULL,
caption = "@mormon_metrics | Source: Cooperative Election Study"
) +
mm_theme
ggsave("./images/31_CES_2025_update_1.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 2: Proportion of US identifying as LDS over time ─────────────────
lds_share <- data |>
filter(!is.na(religion), year >= 2008) |>
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(2008, 2025, 2)) +
scale_y_continuous(labels = scales::percent_format(accuracy = 0.1)) +
labs(
title = "LDS Share of U.S. Population, 2008–2025",
subtitle = "Share of CES respondents identifying as Latter-day Saint",
x = NULL, y = "Share identifying as LDS",
caption = "@mormon_metrics | Source: Cooperative Election Study"
) +
mm_theme
ggsave("./images/31_CES_2025_update_2.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 3: RD identification among LDS over time ─────────────────────────
rd_trend <- lds |>
filter(year >= 2008) |>
group_by(year) |>
summarise(
n_raw = n(),
n_wtd = sum(weight, na.rm = TRUE),
net = weighted.mean(net_rep, w = weight, na.rm = TRUE),
se = sd(net_rep, na.rm = TRUE) / sqrt(n_raw),
lo = net - 1.96 * se,
hi = pmin(net + 1.96 * se, 1)
)
plot <- rd_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") +
scale_x_continuous(breaks = seq(2008, 2025, 2)) +
scale_y_continuous(labels = scales::percent_format()) +
labs(
title = "LDS Net Republican Advantage, 2008–2025",
subtitle = "Positive = more Republican than Democrat; negative = more Democrat",
x = NULL, y = "Net Republican advantage",
caption = "@mormon_metrics | Source: Cooperative Election Study"
) +
mm_theme
ggsave("./images/31_CES_2025_update_3.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 4: Prayer frequency among LDS over time ──────────────────────────
# pew_prayer: 1 = Several times a day, 2 = Once a day, 3+ = less frequent
prayer_trend <- data |>
filter(lds == TRUE, !is.na(prayer), year >= 2008) |>
group_by(year) |>
summarise(
prop = weighted.mean(as.numeric(prayer) %in% 1, w = weight, na.rm = TRUE),
n = n(),
se = sqrt(prop * (1 - prop) / n),
lo = prop - 1.96 * se,
hi = pmin(prop + 1.96 * se, 1)
)
plot <- prayer_trend |>
ggplot(aes(x = year, y = prop)) +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = "#2166ac", alpha = 0.15) +
geom_line(linewidth = 1.2, color = "#2166ac") +
geom_point(size = 2.5, color = "#2166ac") +
scale_x_continuous(breaks = seq(2008, 2025, 2)) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
labs(
title = "LDS Who Pray Several Times A Day, 2008–2025",
subtitle = "Share reporting praying once a day or several times a day",
x = NULL, y = "Share",
caption = "@mormon_metrics | Source: Cooperative Election Study"
) +
mm_theme
ggsave("./images/31_CES_2025_update_4.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 5: Religion very important among LDS over time ───────────────────
# pew_religimp: 1 = Very important, 2 = Somewhat, 3 = Not too, 4 = Not at all
imp_trend <- data |>
filter(lds == TRUE, !is.na(relig_imp), year >= 2008) |>
group_by(year) |>
summarise(
prop = weighted.mean(as.numeric(relig_imp) == 1, w = weight, na.rm = TRUE),
n = n(),
se = sqrt(prop * (1 - prop) / n),
lo = prop - 1.96 * se,
hi = pmin(prop + 1.96 * se, 1)
)
plot <- imp_trend |>
ggplot(aes(x = year, y = prop)) +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = "#762a83", alpha = 0.15) +
geom_line(linewidth = 1.2, color = "#762a83") +
geom_point(size = 2.5, color = "#762a83") +
scale_x_continuous(breaks = seq(2008, 2025, 2)) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
labs(
title = "LDS Who Say Religion Is Very Important, 2008–2025",
subtitle = "Share saying religion is very important in their life",
x = NULL, y = "Share",
caption = "@mormon_metrics | Source: Cooperative Election Study"
) +
mm_theme
ggsave("./images/31_CES_2025_update_5.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 6: Weekly church attendance among LDS over time ──────────────────
# pew_churatd: 1 = More than once a week, 2 = Once a week, 3+ = less frequent
attend_trend <- data |>
filter(lds == TRUE, !is.na(relig_church), year >= 2008) |>
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 = prop - 1.96 * se,
hi = pmin(prop + 1.96 * se, 1)
)
plot <- attend_trend |>
ggplot(aes(x = year, y = prop)) +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = "#1a6faf", alpha = 0.15) +
geom_line(linewidth = 1.2, color = "#1a6faf") +
geom_point(size = 2.5, color = "#1a6faf") +
scale_x_continuous(breaks = seq(2008, 2025, 2)) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
labs(
title = "LDS Who Attend Church Weekly or More, 2008–2025",
subtitle = "Share reporting weekly or more-than-weekly attendance",
x = NULL, y = "Share",
caption = "@mormon_metrics | Source: Cooperative Election Study"
) +
mm_theme
ggsave("./images/31_CES_2025_update_6.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 7: Percent devout among LDS over time ────────────────────────────
# Devout = prays several times a day AND religion very important AND weekly+ attendance
devout_trend <- data |>
filter(lds == TRUE, year >= 2008) |>
filter(!is.na(prayer), !is.na(relig_imp), !is.na(relig_church)) |>
mutate(devout = as.numeric(prayer) == 1 &
as.numeric(relig_imp) == 1 &
as.numeric(relig_church) %in% 1:2) |>
group_by(year) |>
summarise(
prop = weighted.mean(devout, w = weight, na.rm = TRUE),
n = n(),
se = sqrt(prop * (1 - prop) / n),
lo = prop - 1.96 * se,
hi = pmin(prop + 1.96 * se, 1)
)
plot <- devout_trend |>
ggplot(aes(x = year, y = prop)) +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = "#b35806", alpha = 0.15) +
geom_line(linewidth = 1.2, color = "#b35806") +
geom_point(size = 2.5, color = "#b35806") +
scale_x_continuous(breaks = seq(2008, 2025, 2)) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
labs(
title = "Percent Devout Among Latter-day Saints, 2008–2025",
subtitle = "Devout = prays several times/day + religion very important + weekly attendance",
x = NULL, y = "Share",
caption = "@mormon_metrics | Source: Cooperative Election Study"
) +
mm_theme
ggsave("./images/31_CES_2025_update_7.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Gun control data: cumulative (2013–2024) + 2025 ───────────────────────────
gun_cumul <- readRDS("./data/CES/cumulative_2006-2024.rds") |>
mutate(id = paste0(year, "X", case_id)) |>
left_join(read.csv("./data/CES/gun_control.csv"), by = "id") |>
mutate(
year = year.x,
mormon = suppressWarnings(as.numeric(religion)) == 3,
pid3 = suppressWarnings(as.numeric(pid3))
)
# 2025 gun control recodes
# CC25_321a = Ban assault rifles (1=favor, 2=oppose)
# CC25_321b = Easier concealed carry (1=favor, 2=oppose)
# CC25_321c = Background checks (1=favor, 2=oppose)
ces25_gun <- ces25_raw |>
transmute(
year = 2025,
weight = commonweight,
mormon = as.numeric(religpew) == 3,
pid3 = as.numeric(pid3),
favor_background_checks = as.numeric(sjlabelled::as_numeric(CC25_321c)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(CC25_321a)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(CC25_321b)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit)
)
gun_data <- bind_rows(gun_cumul, ces25_gun)
gun_theme <- theme(
axis.ticks = element_blank(),
axis.title = element_text(size = 30, family = "Cairo_bold"),
axis.text.y = element_text(size = 20),
axis.text.x = element_text(size = 20, margin = margin(t = -5)),
plot.background = element_rect(fill = "grey95", color = NA),
panel.background = element_blank(),
text = element_text(family = "Cairo", size = 20),
plot.title = element_text(size = 45, hjust = 0.5, family = "Cairo_bold"),
plot.subtitle = element_text(size = 30, hjust = 0.5, family = "Cairo_bold"),
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(),
panel.grid = element_blank(),
legend.position = "bottom",
panel.grid.major.y = element_line(color = "lightgrey", linetype = "solid"),
plot.caption = element_text(size = 20, family = "Cairo")
)
gun_colors <- c(
"US Overall" = "grey75",
"Democrat" = "#9595ee",
"Independent" = "#ca8bf1",
"Republican" = "tomato",
"Mormon" = "goldenrod"
)
make_gun_plot <- function(var, title, y_limits) {
gun_data |>
group_by(year) |>
count(.data[[var]], wt = weight) |>
drop_na() |>
mutate(prop = n / sum(n), total_n = sum(n),
lower = binom.confint(n, total_n, method = "asymptotic")$lower,
upper = binom.confint(n, total_n, method = "asymptotic")$upper,
audience = "US Overall") |>
filter(.data[[var]] == 1) |>
select(-n, -total_n, -all_of(var)) |>
bind_rows(
gun_data |>
group_by(year, pid3) |>
count(.data[[var]], wt = weight) |>
drop_na() |>
mutate(prop = n / sum(n), total_n = sum(n),
lower = binom.confint(n, total_n, method = "asymptotic")$lower,
upper = binom.confint(n, total_n, method = "asymptotic")$upper,
audience = case_when(
as.numeric(pid3) == 1 ~ "Democrat",
as.numeric(pid3) == 2 ~ "Republican",
as.numeric(pid3) == 3 ~ "Independent",
TRUE ~ NA_character_
)) |>
filter(pid3 %in% 1:3 & .data[[var]] == 1) |>
ungroup() |>
select(-pid3, -n, -total_n, -all_of(var))
) |>
bind_rows(
gun_data |>
group_by(year, mormon) |>
count(.data[[var]], wt = weight) |>
drop_na() |>
mutate(prop = n / sum(n), total_n = sum(n),
lower = binom.confint(n, total_n, method = "asymptotic")$lower,
upper = binom.confint(n, total_n, method = "asymptotic")$upper,
audience = "Mormon") |>
filter(.data[[var]] == 1 & mormon == TRUE) |>
ungroup() |>
select(-mormon, -n, -total_n, -all_of(var))
) |>
mutate(audience = factor(audience,
levels = c("US Overall", "Democrat", "Independent", "Republican", "Mormon"))) |>
ggplot(aes(x = year, y = prop, group = audience, color = audience)) +
geom_line() +
geom_point(size = 3, stroke = 1.5) +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
geom_text(aes(label = scales::percent(prop, accuracy = 1)),
size = 7, family = "Cairo", fontface = "bold", color = "black") +
labs(
title = title,
subtitle = "By Political Party and Mormons",
x = "Year", y = "Proportion",
caption = "@mormon_metrics | Data: Cooperative Election Study (CES) 2013-25"
) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1),
limits = y_limits,
breaks = seq(y_limits[1], y_limits[2], by = 0.1)) +
scale_x_continuous(breaks = seq(2013, 2025, by = 2)) +
scale_color_manual(values = gun_colors) +
theme_minimal() +
gun_theme
}
# ── Analysis 8: Background checks ─────────────────────────────────────────────
plot <- make_gun_plot("favor_background_checks",
"Support for Background Checks on All Gun Sales",
c(0.6, 1))
ggsave("./images/31_CES_2025_update_8.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 9: Assault rifle ban ─────────────────────────────────────────────
plot <- make_gun_plot("favor_ban_assault_rifles",
"Support for Banning Assault Rifles",
c(0.3, 1))
ggsave("./images/31_CES_2025_update_9.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 10: Oppose easy concealed carry ───────────────────────────────────
plot <- make_gun_plot("oppose_easy_permit",
"Oppose Making it Easier To Get A Concealed Carry Permit",
c(0.3, 1))
ggsave("./images/31_CES_2025_update_10.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
# ── Analysis 11: Presidential approval over time ───────────────────────────────
# cumulative approval_pres: 1=Strongly Approve, 2=Approve, 3=Disapprove,
# 4=Strongly Disapprove, 5=Never Heard/Not Sure, 6=Neither
# 2025 CC25_312a: 1=Strongly Approve, 2=Approve, 3=Neither/Not Sure,
# 4=Disapprove, 5=Strongly Disapprove
appr_cumul <- readRDS("./data/CES/cumulative_2006-2024.rds") |>
mutate(
year = as.numeric(year),
mormon = suppressWarnings(as.numeric(religion)) == 3,
approve_pres = case_when(
suppressWarnings(as.numeric(approval_pres)) %in% 1:2 ~ 1,
suppressWarnings(as.numeric(approval_pres)) %in% 3:4 ~ 0,
TRUE ~ NA_real_
)
) |>
select(year, weight, mormon, approve_pres) |>
drop_na()
appr_25 <- ces25_raw |>
transmute(
year = 2025,
weight = commonweight,
mormon = as.numeric(religpew) == 3,
approve_pres = case_when(
as.numeric(CC25_312a) %in% 1:2 ~ 1,
as.numeric(CC25_312a) %in% 4:5 ~ 0,
TRUE ~ NA_real_
)
) |>
drop_na()
appr_data <- bind_rows(appr_cumul, appr_25)
appr_trend <- appr_data |>
filter(!is.na(approve_pres), !is.na(mormon), year >= 2008) |>
group_by(year, mormon) |>
summarise(
prop = weighted.mean(approve_pres, w = weight, na.rm = TRUE),
n = n(),
se = sqrt(prop * (1 - prop) / n),
lo = prop - 1.96 * se,
hi = pmin(prop + 1.96 * se, 1),
.groups = "drop"
) |>
mutate(audience = ifelse(mormon, "Latter-day Saint", "US Overall"))
plot <- appr_trend |>
ggplot(aes(x = year, y = prop, color = audience, fill = audience, group = audience)) +
annotate("rect", xmin = 2008, xmax = 2016, ymin = -Inf, ymax = Inf,
fill = "dodgerblue3", alpha = 0.05) +
annotate("rect", xmin = 2016, xmax = 2020, ymin = -Inf, ymax = Inf,
fill = "darkred", alpha = 0.05) +
annotate("rect", xmin = 2020, xmax = 2024, ymin = -Inf, ymax = Inf,
fill = "dodgerblue3", alpha = 0.05) +
annotate("rect", xmin = 2024, xmax = 2026, ymin = -Inf, ymax = Inf,
fill = "darkred", alpha = 0.05) +
annotate("text", x = 2012, y = 0.97, label = "Obama", size = 9, color = "grey40", family = "Cairo", fontface = "bold") +
annotate("text", x = 2018, y = 0.97, label = "Trump", size = 9, color = "grey40", family = "Cairo", fontface = "bold") +
annotate("text", x = 2022, y = 0.97, label = "Biden", size = 9, color = "grey40", family = "Cairo", fontface = "bold") +
annotate("text", x = 2025.3, y = 0.97, label = "Trump", size = 9, color = "grey40", family = "Cairo", fontface = "bold") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, color = NA) +
geom_line(linewidth = 1.2) +
geom_point(size = 2.5) +
scale_x_continuous(breaks = seq(2008, 2025, 2)) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
scale_color_manual(values = c("Latter-day Saint" = "goldenrod", "US Overall" = "grey50")) +
scale_fill_manual(values = c("Latter-day Saint" = "goldenrod", "US Overall" = "grey50"),
guide = "none") +
labs(
title = "Presidential Approval: LDS vs. U.S. Overall, 2008–2025",
subtitle = "Share approving or strongly approving; shaded bands = 95% CI",
x = NULL, y = "Approval rate", color = NULL,
caption = "@mormon_metrics | Source: Cooperative Election Study"
) +
mm_theme
ggsave("./images/31_CES_2025_update_11.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)
Raw
combining_ces_data_gun_control.r
library(tidyverse)
library(haven)
library(sjlabelled)
library(scales)
# The gun control battery I'm interested in starts in 2013
# I want to create 5 variables:
# 1. owngun (1 = someone in household owns a gun, 0 = No one in household owns a gun)
# 2. favor_background_checks (1 = favor, 0 = oppose) on all gun sales
# 3. favor_ban_assault_rifles (1 = favor, 0 = oppose)
# 4. oppose_easy_permit (1 = favor, 0 = oppose) # Make it easier to obtain a conceal and carry permit
# 5. bpsca (1 = favor, 0 = oppose) # Bipartisan Safer Communities Act
final <- list()
# 2013
data <- load('./data/CES/data/2013.RData')
ls()
data <- table
names(data)
final[['2013']] <- data |>
mutate(
id = paste0('2013X', caseid),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC13_320a)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC13_320d)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC13_320e)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2014
data <- load('./data/CES/data/2014.RData')
data <- x
names(data)
sjlabelled::get_labels(data$CC14_320a)
final[['2014']] <- data |>
mutate(
id = paste0('2014X', V101),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC14_320a)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 3:4, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC14_320d)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 3:4, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC14_320e)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 3:4, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2015
data <- load('./data/CES/data/2015.RData')
data <- x
names(data)
sjlabelled::as_numeric(data$CC15_320a)
final[['2015']] <- data |>
mutate(
id = paste0('2015X', V101),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC15_320a)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 3:4, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC15_320c)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 3:4, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC15_320d)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 3:4, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2016
data <- load('./data/CES/data/2016.RData')
data <- x
names(data)
final[['2016']] <- data |>
mutate(
id = paste0('2016X', V101),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC16_330a)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 3:4, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC16_330d)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 3:4, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC16_330e)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 3:4, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2017
data <- load('./data/CES/data/2017.RData')
ls()
data <- table
names(data)
final[['2017']] <- data |>
mutate(
id = paste0('2017X', V101),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC17_330a)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC17_330c)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC17_330d)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2018
data <- load('./data/CES/data/2018.RData')
data <- x
names(data)
final[['2018']] <- data |>
mutate(
id = paste0('2018X', caseid),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC18_320a)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC18_320c)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC18_320d)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2019
data <- load('./data/CES/data/2019.RData')
ls()
data <- table
names(data)
final[['2019']] <- data |>
mutate(
id = paste0('2019X', caseid),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC19_320a)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC19_320c)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC19_320d)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2020
data <- haven::read_dta('./data/CES/data/2020.dta')
names(data)
final[['2020']] <- data |>
mutate(
id = paste0('2020X', caseid),
owngun = dplyr::case_when(
sjlabelled::as_numeric(data$gunown) %in% 1:2 ~ "Gun in household",
sjlabelled::as_numeric(data$gunown) == 3 ~ "No gun in household",
sjlabelled::as_numeric(data$gunown) == 8 ~ "Unsure",
TRUE ~ NA_character_
),
owngun = factor(owngun, levels = c("Gun in household", "No gun in household", "Unsure")),
favor_background_checks = NA,
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC20_330b)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC20_330c)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2021
data <- load('./data/CES/data/2021.RData')
ls()
data <- table
names(data)
final[['2021']] <- data |>
mutate(
id = paste0('2021X', caseid),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC21_321f)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC21_321a)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC21_321b)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2022
data <- haven::read_dta('./data/CES/data/2022.dta')
names(data)
final[['2022']] <- data |>
mutate(
id = paste0('2022X', caseid),
owngun = dplyr::case_when(
sjlabelled::as_numeric(data$gunown) %in% 1:2 ~ "Gun in household",
sjlabelled::as_numeric(data$gunown) == 3 ~ "No gun in household",
sjlabelled::as_numeric(data$gunown) == 8 ~ "Unsure",
TRUE ~ NA_character_
),
owngun = factor(owngun, levels = c("Gun in household", "No gun in household", "Unsure")),
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC22_330e)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC22_330b)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC22_330c)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = NA
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2023
data <- haven::read_dta('./data/CES/data/2023.dta')
names(data)
final[['2023']] <- data |>
mutate(
id = paste0('2023X', caseid),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC23_321c)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC23_321a)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC23_321b)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = as.numeric(sjlabelled::as_numeric(data$CC23_340e)),
bpsca = ifelse(bpsca == 2, 0, bpsca),
bpsca = ifelse(bpsca %in% 8:9, NA, bpsca)
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# 2024
data <- haven::read_dta('./data/CES/data/2024.dta', encoding = 'UTF-8')
names(data)
final[['2024']] <- data |>
mutate(
id = paste0('2024X', caseid),
owngun = NA,
favor_background_checks = as.numeric(sjlabelled::as_numeric(data$CC24_321c)),
favor_background_checks = ifelse(favor_background_checks == 2, 0, favor_background_checks),
favor_background_checks = ifelse(favor_background_checks %in% 8:9, NA, favor_background_checks),
favor_ban_assault_rifles = as.numeric(sjlabelled::as_numeric(data$CC24_321a)),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles == 2, 0, favor_ban_assault_rifles),
favor_ban_assault_rifles = ifelse(favor_ban_assault_rifles %in% 8:9, NA, favor_ban_assault_rifles),
oppose_easy_permit = as.numeric(sjlabelled::as_numeric(data$CC24_321b)),
oppose_easy_permit = ifelse(oppose_easy_permit == 1, 0, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit == 2, 1, oppose_easy_permit),
oppose_easy_permit = ifelse(oppose_easy_permit %in% 8:9, NA, oppose_easy_permit),
bpsca = as.numeric(sjlabelled::as_numeric(data$CC24_321f)),
bpsca = ifelse(bpsca == 2, 0, bpsca),
bpsca = ifelse(bpsca %in% 8:9, NA, bpsca)
) |>
select(id, owngun, favor_background_checks, favor_ban_assault_rifles, oppose_easy_permit, bpsca)
# Combine all years into a single data frame
combined_data <- bind_rows(final, .id = 'year')
write.csv(combined_data, 'data/CES/gun_control.csv')
Raw
combining_ces_data_prayer.r
library(tidyverse)
library(haven)
library(sjlabelled)
library(scales)
# 2008
load('./data/CES/data/2008.RData')
# Check what objects were loaded
ls()
data <- x
final <- list()
#2008
final[['2008']] <- data |>
mutate(
id = paste0('2008X', V100),
prayer = as.numeric(sjlabelled::as_numeric(data$V218))
) |>
select(id, prayer)
# 2009
data <- haven::read_dta('./data/CES/data/2009.dta', encoding = 'UTF-8')
names(data)
data$v218
final[['2009']] <- data |>
mutate(
id = paste0('2009X', v100),
prayer = as.numeric(data$v218)
) |>
select(id, prayer)
# 2010
data <- haven::read_dta('./data/CES/data/2010.dta')
names(data)
final[['2010']] <- data |>
mutate(
id = paste0('2010X', V100),
prayer = as.numeric(data$V218)
) |>
select(id, prayer)
# 2011
data <- load('./data/CES/data/2011.RData')
data <- x
names(data)
final[['2011']] <- data |>
mutate(
id = paste0('2011X', V100),
prayer = as.numeric(sjlabelled::as_numeric(data$V218))
) |>
select(id, prayer)
# 2012
data <- load('./data/CES/data/2012.RData')
data <- x
names(data)
final[['2012']] <- data |>
mutate(
id = paste0('2012X', V101),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2013
data <- load('./data/CES/data/2013.RData')
ls()
data <- table
names(data)
final[['2013']] <- data |>
mutate(
id = paste0('2013X', caseid),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2014
data <- load('./data/CES/data/2014.RData')
data <- x
names(data)
final[['2014']] <- data |>
mutate(
id = paste0('2014X', V101),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2015
data <- load('./data/CES/data/2015.RData')
data <- x
names(data)
final[['2015']] <- data |>
mutate(
id = paste0('2015X', V101),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2016
data <- load('./data/CES/data/2016.RData')
data <- x
names(data)
final[['2016']] <- data |>
mutate(
id = paste0('2016X', V101),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2017
data <- load('./data/CES/data/2017.RData')
ls()
data <- table
names(data)
final[['2017']] <- data |>
mutate(
id = paste0('2017X', V101),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2018
data <- load('./data/CES/data/2018.RData')
data <- x
names(data)
final[['2018']] <- data |>
mutate(
id = paste0('2018X', caseid),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2019
data <- load('./data/CES/data/2019.RData')
ls()
data <- table
names(data)
final[['2019']] <- data |>
mutate(
id = paste0('2019X', caseid),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2020
data <- haven::read_dta('./data/CES/data/2020.dta')
names(data)
final[['2020']] <- data |>
mutate(
id = paste0('2020X', caseid),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2021
data <- load('./data/CES/data/2021.RData')
ls()
data <- table
names(data)
final[['2021']] <- data |>
mutate(
id = paste0('2021X', caseid),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2022
data <- haven::read_dta('./data/CES/data/2022.dta')
names(data)
final[['2022']] <- data |>
mutate(
id = paste0('2022X', caseid),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2023
data <- haven::read_dta('./data/CES/data/2023.dta')
names(data)
final[['2023']] <- data |>
mutate(
id = paste0('2023X', caseid),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# 2024
data <- haven::read_dta('./data/CES/data/2024.dta', encoding = 'UTF-8')
names(data)
final[['2024']] <- data |>
mutate(
id = paste0('2024X', caseid),
prayer = as.numeric(sjlabelled::as_numeric(data$pew_prayer))
) |>
select(id, prayer)
# Combine all years into a single data frame
combined_data <- bind_rows(final, .id = 'year')
write.csv(combined_data, 'data/CES/prayer.csv')
# combined_ces_data <- haven::read_dta('data/CES/cumulative_2006-2024.dta')
# #missing 1 case in 2023
# combined_ces_data |>
# mutate(id = paste0(year ,'X', case_id)) |>
# dplyr::inner_join(combined_data, by = 'id') |>
# count(year.x) |>
# summarize(sum(n))
# combined_ces_data |>
# filter(year %in% 2008:2024) |>
# count(year) |>
# summarize(sum(n))
@thaumkid

thaumkid commented Apr 11, 2026

Copy link
Copy Markdown

Straightforward enough to see that the cumulative data could be found here --> https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/II2DB6 and that the 2025.dta is probably from here --> https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/DFUGJR. But I can't immediately see what you're using for your 2025 prayer and gun control data -- "./data/CES/prayer.csv" and "./data/CES/gun_control.csv". Can you specify how you generated those files so I don't have to guess about transformations and formatting? Thanks

@acbass49

acbass49 commented Apr 11, 2026

Copy link
Copy Markdown
Author

@thaumkid thanks for engaging with this work. I updated the gist with the scripts used to generate the data files you referenced. I tried to upload the data files referenced directly, but they are too large. If you want to data files I used, feel free to send me your email and I can send you them directly. Let me know if theres anything else I can do to help. Cheers!

@acbass49

acbass49 commented Apr 11, 2026

Copy link
Copy Markdown
Author

And yes those are the links to where to download the two ces data files used (cumulative and the new 2025)

@thaumkid

thaumkid commented Apr 13, 2026

Copy link
Copy Markdown

Thanks for your help! I was able to replicate your code, and then I used linear interpolation from the 2000, 2010, and 2020 US Census population counts (and the wikipedia 2025 estimate) to convert your charts 2 and 4-7 to approximate absolute counts out of curiosity --

us_pop_2000 <- 281421906
us_pop_2010 <- 308745538
us_pop_2020 <- 331449281
us_est_pop_2025 <- 341784857
us_pop <- c(approx(c(us_pop_2000,us_pop_2010),n=11)$y,approx(c(us_pop_2010,us_pop_2020),n=11)$y[-1],approx(c(us_pop_2020,us_est_pop_2025),n=6)$y[-1])[-c(1:8)]

then multiplied us_pop by lds_share$share_lds and various geom_ribbon values to derive the following plot edits

plot 2:
31_CES_2025_update_2 1

plot 4:
31_CES_2025_update_4 1

plot 5:
31_CES_2025_update_5 1

plot 6:
31_CES_2025_update_6 1

plot 7:
31_CES_2025_update_7 1

@thaumkid

thaumkid commented Apr 13, 2026 •
edited
Loading

Copy link
Copy Markdown

the apparent numbers going up during election years was an interesting artifact of this approach. of course, mostly everything just follows the overall population numbers of the first chart, so not really very interesting, except for maybe the numbers being raised a bit due to the changing census values, and the overall devout number kind of stabilizing from 2019 onwards

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment