23 devout decline investigation

23_trend_investigation.r

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

# Load the dataset
data <- readRDS("./data/CES/cumulative_2006-2024.rds")

data$relig_church_rc <- car::recode(as.numeric(data$relig_church),"1:2 = 'Weekly or More'; 3 = 'Monthly'; 4:5 = 'A few times a year/Seldom'; 6 = 'Never'")
data$relig_church_rc <- factor(data$relig_church_rc, 
    levels = c("Weekly or More", "Monthly", "A few times a year/Seldom", "Never"))

data$relig_church_rc2 <- car::recode(as.numeric(data$relig_church),"1:2 = 'Weekly or More'; 3:4 = 'Monthly/A few times a year'; 5:6 = 'Seldom/Never'")
data$relig_church_rc2 <- factor(data$relig_church_rc2, 
    levels = c("Weekly or More", "Monthly/A few times a year", "Seldom/Never"))

data$year_4 <- car::recode(as.numeric(data$year), "2008:2012 = '2008-12'; 2013:2016 = '2013-16'; 2017:2020 = '2017-20'; 2021:2024 = '2021-24'")
data$year_4 <- factor(data$year_4, 
    levels = c("2008-12", "2013-16", "2017-20", "2021-24"))

data$year_2 <- car::recode(as.numeric(data$year), "2007:2008 = '2007-08'; 2009:2010 = '2009-10'; 2011:2012 = '2011-12'; 2013:2014 = '2013-14'; 2015:2016 = '2015-16'; 2017:2018 = '2017-18'; 2019:2020 = '2019-20'; 2021:2022 = '2021-22'; 2023:2024 = '2023-24'")
data$year_2 <- factor(data$year_2, 
    levels = c("2007-08", "2009-10", "2011-12", "2013-14", "2015-16", "2017-18", "2019-20", "2021-22", "2023-24"))

# prayer was not in the cumulative file, so I need to merge it in
# I went through each year and added the prayer variable manually.
data <- data |> 
    mutate(id = paste0(year, "X", case_id)) |> 
    left_join(read.csv("./data/CES/prayer.csv"), by = "id") |> 
    mutate(year = year.x)

data <- data |> 
    mutate(
        prayer_rc = case_when(
            prayer == 1 ~ "Several Times a Day",
            prayer %in% 2:7 ~ "Less Than Several Times a Day",
            TRUE ~ NA_character_
        ),
        prayer_rc2 = case_when(
            prayer %in% 1:2 ~ "Daily or More",
            prayer %in% 3:7 ~ "Less Than Several Times a Day",
            TRUE ~ NA_character_
        )
    )

data$devout <- dplyr::case_when(
    data$prayer_rc == "Several Times a Day" & 
    data$relig_imp == "Very Important" &
    data$relig_church_rc == "Weekly or More" ~ 1,
    TRUE ~ 0
)

data$Utah <- ifelse(data$state == "Utah", "Utah Resident", "Non-Utah Resident")
data$jello_belt <- ifelse(data$state %in% c("Utah", "Idaho", "Arizona"), 1, 0)

data$educ_rc <- dplyr::case_when(
    data$educ %in% 1:2 ~ "HS or Less",
    data$educ %in% 3:4 ~ "Some college",
    data$educ %in% 5 ~ "Bachelor",
    data$educ %in% 6 ~ "Graduate"
)

data$educ_rc <- factor(data$educ_rc, levels = c("HS or Less","Some college","Bachelor","Graduate"))

data$age_rc <- dplyr::case_when(
    data$age %in% 18:30 ~ "18-30",
    data$age %in% 31:45 ~ "31-45",
    data$age %in% 46:65 ~ "46-65",
    data$age %in% 66:110 ~ "66+"
)

data$age_rc <- factor(data$age_rc, levels = c("18-30","31-45","46-65","66+"))

data$race_rc <- ifelse(data$race == 1, "white", "non-white")
data$race_rc <- factor(data$race_rc, levels = c("white", "non-white"))

data$ageXgender <- dplyr::case_when(
    data$age %in% 18:30 & data$gender == 1 ~ "M: 18-30",
    data$age %in% 31:45 & data$gender == 1 ~ "M: 31-45",
    data$age %in% 46:65 & data$gender == 1 ~ "M: 46-65",
    data$age %in% 66:110 & data$gender == 1 ~ "M: 66+",
    data$age %in% 18:30 & data$gender == 2 ~ "F: 18-30",
    data$age %in% 31:45 & data$gender == 2 ~ "F: 31-45",
    data$age %in% 46:65 & data$gender == 2 ~ "F: 46-65",
    data$age %in% 66:110 & data$gender == 2 ~ "F: 66+"
)

data$ageXgender <- factor(data$ageXgender, levels = c("M: 18-30","M: 31-45","M: 46-65","M: 66+","F: 18-30","F: 31-45","F: 46-65","F: 66+"))

data$generation <- case_when(
    data$birthyr < 1946 ~ "Greatest/Silent",
    data$birthyr >= 1946 & data$birthyr < 1965 ~ "Boomers",
    data$birthyr >= 1965 & data$birthyr < 1981 ~ "GenX",
    data$birthyr >= 1981 & data$birthyr < 1997 ~ "Mil.",
    data$birthyr >= 1997 ~ "GenZ"
)

data$generation <- factor(
    data$generation, 
    levels = 
        c("Greatest/Silent", "Boomers", "GenX", "Mil.", "GenZ")
    )

data$cohort <- dplyr::case_when(
    data$birthyr %in% 1940:1949 ~ "Born in 1940s",
    data$birthyr %in% 1950:1959 ~ "Born in 1950s",
    data$birthyr %in% 1960:1969 ~ "Born in 1960s",
    data$birthyr %in% 1970:1979 ~ "Born in 1970s",
    data$birthyr %in% 1980:1989 ~ "Born in 1980s",
    data$birthyr %in% 1990:1999 ~ "Born in 1990s"
)

data$cohort <- factor(data$cohort, levels = c(
    "Born in 1940s",
    "Born in 1950s",
    "Born in 1960s",
    "Born in 1970s",
    "Born in 1980s",
    "Born in 1990s"
))

data$age_buckets <- dplyr::case_when(
    data$age %in% 18:24 ~ "18-24",
    data$age %in% 25:30 ~ "25-30",
    data$age %in% 31:35 ~ "31-35",
    data$age %in% 36:40 ~ "36-40",
    data$age %in% 41:45 ~ "41-45",
    data$age %in% 46:50 ~ "46-50",
    data$age %in% 51:55 ~ "51-55",
    data$age %in% 56:60 ~ "56-60",
    data$age %in% 61:65 ~ "61-65",
    data$age %in% 66:70 ~ "66-70",
    data$age %in% 66:70 ~ "71-75"
)

data$age_buckets <- factor(data$age_buckets, levels = c(
    "18-24",
    "25-30",
    "31-35",
    "36-40",
    "41-45",
    "46-50",
    "51-55",
    "56-60",
    "61-65",
    "66-70",
    "71-75"
))

data$Republican <- ifelse(data$pid3 == 2, "Republican", "Other")

data$pid3_rc <- ifelse(data$pid3 %in% 1, "Democrat", ifelse(data$pid3 %in% 2, "Republican", "Independent"))
data$pid3_rc <- factor(data$pid3_rc, levels = c("Republican", "Independent", "Democrat"))

mormon_data <- data |> filter(data$religion == 3)

# gender, age, education, race, party, ageXgender

mormon_data |> 
    group_by(gender,year_4) |> 
    count(devout, 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 = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1)

mormon_data |> 
    group_by(educ_rc,year_4) |> 
    count(devout, 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 = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1)

mormon_data |> 
    group_by(age_rc,year_4) |> 
    count(devout, 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 = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1)

mormon_data |> 
    group_by(race_rc,year_4) |> 
    count(devout, 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 = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1)

mormon_data |> 
    group_by(pid3,year_4) |> 
    count(devout, wt = weight) |> 
    drop_na() |>
    filter(pid3 %in% 1:3) |> 
    mutate(
        prop = n / sum(n),
        total_n = sum(n),
        lower = binom.confint(x = n, n = total_n, method = "asymptotic")$lower,
        upper = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1)

mormon_data |> 
    group_by(ageXgender,year_4) |> 
    count(devout, 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 = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1) |> 
    print(n=50)

mormon_data |> 
    group_by(generation,year_4) |> 
    count(devout, 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 = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1) |> 
    print(n=50)


summary(lm(devout ~ age*year + educ*year + gender*year + race_rc*year + Republican*year + year, data = mormon_data, weights = mormon_data$weight))

mormon_data |> 
    group_by(cohort, age_buckets) |> 
    summarise("devoutness" = mean(devout), n = n()) |> 
    drop_na() |> 
    filter(n>75) |> 
    ggplot(aes(x = age_buckets, y = devoutness, color = cohort, group = cohort)) +
    geom_smooth(method = "loess", se = FALSE)


plot <- mormon_data |> 
    group_by(age_rc,year_4) |> 
    count(devout, 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 = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1) |> 
    ggplot(aes(x = year_4, y = prop, color = age_rc, group = age_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") +
    labs(
        title = "Decline Of 'Devout' Seems Driven By Age",
        subtitle = "Devout defined as praying several times a day,\nreligion very important,\nand attending church weekly or more",
        x = "Year Group",
        y = "Proportion of Devout",
        caption = "@mormon_metrics\nData: Cooperative Election Study (CES) 2008-24. N=9,189 LDS"
    ) +
    theme_minimal() +
    scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0,0.8) , breaks = seq(0, 0.8, by = 0.1)) +
    scale_color_manual(values = c("18-30" = "#A8CDE9", "31-45" = "#74ADD1", "46-65"="#4575B4", "66+"="#313695")) +
    theme(
        axis.ticks = element_blank(),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 10),
        axis.text.x = element_text(size = 6, margin = margin(t = -5), angle = 45),
        plot.background = element_rect(fill = "grey95"),  # Change entire plot background color here
        panel.background = element_blank(),
        text = element_text(face = "bold", family = "Cairo"),
        plot.title = element_text(size = 15, face = "bold"),
        plot.subtitle = element_text(size = 11),
        plot.title.position = "plot",
        plot.subtitle.position = "plot",
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.title = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        legend.position = "none",
        panel.grid.major.y = element_line(color = "lightgrey", linetype = 'solid'),
        plot.caption = element_text(size = 8,family = "Cairo")
    ) +
    facet_grid(~age_rc)
ggsave("./images/23_trend_investigation_1.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

plot <- mormon_data |> 
    group_by(pid3_rc,year_4) |> 
    count(devout, 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 = binom.confint(x = n, n = total_n, method = "asymptotic")$upper
    ) |> 
    filter(devout == 1) |> 
    ggplot(aes(x = year_4, y = prop, color = pid3_rc, group = pid3_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") +
    labs(
        title = "Republican 'Devoutness' Decreasing Slightly Faster",
        subtitle = "Devout defined as praying several times a day,\nreligion very important,\nand attending church weekly or more",
        x = "Year Group",
        y = "Proportion of Devout",
        caption = "@mormon_metrics\nData: Cooperative Election Study (CES) 2008-24. N=9,189 LDS"
    ) +
    theme_minimal() +
    scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0,0.8) , breaks = seq(0, 0.8, by = 0.1)) +
    scale_color_manual(values = c("Republican" = "darkred", "Independent" = "purple", "Democrat"="dodgerblue3")) +
    theme(
        axis.ticks = element_blank(),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 10),
        axis.text.x = element_text(size = 6, margin = margin(t = -5), angle = 45),
        plot.background = element_rect(fill = "grey95"),  # Change entire plot background color here
        panel.background = element_blank(),
        text = element_text(face = "bold", family = "Cairo"),
        plot.title = element_text(size = 15, face = "bold"),
        plot.subtitle = element_text(size = 11),
        plot.title.position = "plot",
        plot.subtitle.position = "plot",
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.title = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        legend.position = "none",
        panel.grid.major.y = element_line(color = "lightgrey", linetype = 'solid'),
        plot.caption = element_text(size = 8,family = "Cairo")
    ) +
    facet_grid(~pid3_rc)
ggsave("./images/23_trend_investigation_2.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)