16 Cluster Confirmation

16_cluster_confirmation.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_
        )
    )

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

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

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

#creating charts...
# decline of devout Mormons...
mormon_data |> 
    group_by(year_4) |> 
    count(devout) |> 
    drop_na() |> 
    summarise(total = sum(n)) |> 
    pull(total) |> sum()

plot <- mormon_data |> 
    group_by(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, group = 1)) +
    geom_line(color = "darkorange") +
    geom_point(
        size = 3,
        stroke = 1.5,
        color = "darkorange",
    ) +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, color = "#FF8C00") +
    geom_text(aes(label = scales::percent(prop, accuracy = 1)), vjust = -1, size = 4, family = "Cairo", fontface = "bold", color = "black") +
    labs(
        title = "Proportion of US Devout Latter-day Saints Declining",
        subtitle = "Devout defined as praying several times a day,\nreligion very important,\nand attending church weekly or more",
        x = "Year Group",
        y = "Proportion",
        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)) +
    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"),  # 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.title = element_blank(),
        legend.background = element_blank(),
        legend.box.background = 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")
    )

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

# I included a few other ways of defining devout Mormons to see if the trend holds
# It does across the two other definitions I tested below.

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

mormon_data |> 
    group_by(year_4) |> 
    count(devout_2, 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_2 == 1)

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

mormon_data |> 
    group_by(year_4) |> 
    count(devout_3, 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_3 == 1)


# Rise of cultural Mormons...
mormon_data |> 
    group_by(year_4) |> 
    count(relig_imp, 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
    )

plot <- mormon_data |> 
    group_by(year_4) |> 
    count(relig_church_rc2, 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
    ) |> 
    ggplot(aes(x = year_4, y = prop, color = relig_church_rc2, group = relig_church_rc2)) +
    geom_line() +
    geom_point() +
    geom_text(aes(label = scales::percent(prop, accuracy = 1)), vjust = -1, size = 3, family = "Cairo", fontface = "bold", color = "black") +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
    labs(
        title = "Proportion of US Latter-day Saints by Church Attendance",
        subtitle = "Church attendance recoded to three groups:\nWeekly or More, Monthly/A few times a year, Seldom/Never",
        x = "Year Group",
        y = "Proportion",
        color = "Church Attendance",
        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)) +
    theme(
        axis.ticks = element_blank(),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 10),
        axis.text.x = element_text(size = 8, margin = margin(t = -5)),
        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.title = element_text(size = 12, face = "bold"),
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        legend.position = "top",
        panel.grid.major.y = element_line(color = "lightgrey", linetype = 'solid'),
        plot.caption = element_text(size = 8,family = "Cairo")
    ) +
    facet_wrap(~relig_church_rc2)
ggsave("./images/16_confirmation_2.png", plot, width = 2000, height = 1500, units = "px", dpi = 300)

# Rise of the adaptive Mormons...

mormon_data |> 
    filter(relig_church_rc == "Weekly or More") |>
    group_by(year_4) |> 
    count(pid3) |> 
    drop_na() |> 
    pull(n) |> sum()
#6302 weekly attenders

plot <- mormon_data |> 
    filter(relig_church_rc == "Weekly or More") |>
    group_by(year_4) |> 
    count(pid3, wt = weight) |> 
    drop_na() |>
    filter(!pid3 %in% 4:5) |>
    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,
        pid3 = sjlabelled::as_label(pid3),
        pid3 = factor(pid3, levels = c("Democrat", "Independent", "Republican"))
    ) |> 
    ggplot(aes(x = year_4, y = prop, color = pid3, group = pid3)) +
    geom_line() +
    geom_point() +
    geom_text(aes(label = scales::percent(prop, accuracy = 1)), vjust = -1.5, size = 3, family = "Cairo", fontface = "bold", color = "black") +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
    labs(
        title = "Proportion of US Latter-day Saints by Party ID",
        subtitle = "Among Weekly or More Church Attending Latter-day Saints",
        x = "Year Group",
        y = "Proportion",
        color = "Party ID",
        caption = "@mormon_metrics\nData: Cooperative Election Study (CES) 2008-24. N=6,302 LDS Weekly Attenders"
    ) +
    theme_minimal() +
    scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0,1) , breaks = seq(0, 1, by = 0.1)) +
    theme(
        axis.ticks = element_blank(),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 10),
        axis.text.x = element_text(size = 8, margin = margin(t = -5)),
        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.title = element_text(size = 12, face = "bold"),
        legend.background = element_blank(),
        legend.box.background = 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_wrap(~pid3) + 
    scale_color_manual(values = c("blue", "grey40", "red"))
ggsave("./images/16_confirmation_3.png", plot, width = 2000, height = 1500, units = "px", dpi = 300)

mormon_data |> 
    filter(relig_church_rc == "Weekly or More") |>
    group_by(year) |> 
    count(abortion_always_allow, 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(abortion_always_allow == 1) |>
    print(n=50)

plot <- mormon_data |> 
    filter(relig_church_rc == "Weekly or More") |>
    group_by(year_2) |> 
    count(abortion_always_allow, 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(abortion_always_allow == 1) |> 
    ggplot(aes(x = year_2, y = prop, group = 1)) +
    geom_line(color = "darkorange") +
    geom_point(
        size = 3,
        stroke = 1.5,
        color = "darkorange",
    ) +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, color = "#FF8C00") +
    geom_text(aes(label = scales::percent(prop, accuracy = 1)), vjust = -2, size = 4, family = "Cairo", fontface = "bold", color = "black") +
    labs(
        title = "Share of US LDS Who Say 'Abortion Should Always be Allowed'",
        subtitle = "Among Weekly or More Church Attending Latter-day Saints",
        x = "Year Group",
        y = "Proportion",
        caption = "@mormon_metrics\nData: Cooperative Election Study (CES) 2008-24. N=6,302 LDS Weekly Attenders"
    ) +
    theme_minimal() +
    scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0,0.4) , breaks = seq(0, 0.4, by = 0.05)) +
    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"),  # 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.title = element_text(size = 12, face = "bold"),
        legend.background = element_blank(),
        legend.box.background = 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")
    )
ggsave("./images/16_confirmation_4.png", plot, width = 2000, height = 1500, units = "px", dpi = 300)
# link: https://gist.github.com/acbass49/ae50bccdc1a376b89bda0d829559c62f

mormon_data |> 
    filter(relig_church_rc == "Weekly or More") |>
    group_by(year_2) |> 
    count(abortion_always_allow, 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(abortion_always_allow == 1) |>
    print(n=50)