13 2 hour church policy change?

13_2hour_policy_change.r

library(tidyverse)
library(tidysynth)
library(sjlabelled)

data  <- readRDS('../client_work/jana_reiss/cumulative_2006-2024.rds')

data$religion_rc <- dplyr::case_when(

    # ---- Jehovah's Witnesses ----
    as_label(data$relig_protestant) == "Jehovah's Witness" ~ "Jehovah's Witnesses",

    # ---- Historically Black Protestant ----
    as_label(data$religion) == "Protestant" &
        data$race == 2 ~ "Historically Black Protestant",  # race == 2 is typically Black in CES

    # ---- Evangelical Protestant ----
    as_label(data$religion) == "Protestant" &
        as_label(data$relig_bornagain) == "Yes" ~ "Evangelical Protestant",

    # ---- Mainline Protestant ----
    as_label(data$religion) == "Protestant" &
        as_label(data$relig_bornagain) != "Yes" ~ "Mainline Protestant",

    # ---- Catholic ----
    as_label(data$religion) == "Roman Catholic" ~ "Catholic",

    # ---- LDS ----
    as_label(data$religion) == "Mormon" ~ "LDS",

    # ---- Orthodox ----
    as_label(data$religion) == "Eastern or Greek Orthodox" ~ "Orthodox",

    # ---- Other Christian ----
    as_label(data$religion) == "Protestant" ~ "Other Christian",

    # ---- Jewish ----
    as_label(data$religion) == "Jewish" ~ "Jewish",

    # ---- Muslim ----
    as_label(data$religion) == "Muslim" ~ "Muslim",

    # ---- Buddhist ----
    as_label(data$religion) == "Buddhist" ~ "Buddhist",

    # ---- Hindu ----
    as_label(data$religion) == "Hindu" ~ "Hindu",

    # ---- Atheist ----
    as_label(data$religion) == "Atheist" ~ "Atheist",

    # ---- Agnostic ----
    as_label(data$religion) == "Agnostic" ~ "Agnostic",

    # ---- Nothing in particular ----
    as_label(data$religion) == "Nothing in Particular" ~ "Nothing in particular",

    # ---- All other cases ----
    TRUE ~ "Other"
    )

data$relig_church_rc <- car::recode(as.numeric(data$relig_church),
    "1:2=1; 3:6=0; 7:8=NA"
)

data$educ_rc <- as_label(data$educ)
data$educ_rc <- factor(data$educ_rc, 
    levels = c("No HS", "High School Graduate", "Some College", "2-Year", "4-Year", "Post-Grad"))

data$ba_plus <- case_when(
    data$educ_rc %in% c("4-Year", "Post-Grad") ~ 1,
    data$educ_rc %in% c("2-Year", "No HS", "High School Graduate", "Some College") ~ 0,
    TRUE ~ NA_real_
)

data$pct_white <- case_when(
    data$race == 1 ~ 1,
    data$race %in% 2:8 ~ 0,
    TRUE ~ NA_real_
)

data$pct_married <- case_when(
    data$marstat == 1 ~ 1, # Married
    data$marstat %in% 2:6 ~ 0, # Not married
    TRUE ~ NA_real_
)

data$pct_rep <- case_when(
    data$pid3 == 2 ~ 1, # Republican
    data$pid3 %in% c(1,3,4,5) ~ 0, # Not Republican
    TRUE ~ NA_real_
)

data$pct_has_kids <- case_when(
    data$has_child == "Yes" ~ 1, # Has kids
    data$has_child %in% c("No", "Don't know", "Refused") ~ 0, # No kids
    TRUE ~ NA_real_
)

data$pct_female <- ifelse(data$sex == 2, 1, 0)

data$relig_imp_rc <- dplyr::case_when(
    data$relig_imp == "Very Important" ~ 1,
    data$relig_imp == "Somewhat Important" ~ 0,
    data$relig_imp == "Not Too Important" ~ 0,
    data$relig_imp == "Not at All Important" ~ 0,
    TRUE ~ NA_real_
)

# OK I just want to test the aggregated data
df_church <- data |> 
    filter(religion %in% c(3)) |>
    group_by(year) |>
    count(relig_church_rc, wt = weight) |> 
    drop_na() |> 
    mutate(
        prop = n / sum(n)
    ) |>
    filter(relig_church_rc == 1) |> 
    mutate(
        post_2019 = ifelse(year %in% 2019:2024, 1, 0)
    )

df_imp <- data |> 
    filter(religion %in% c(3)) |>
    group_by(year) |>
    count(relig_imp_rc, wt = weight) |> 
    drop_na() |> 
    mutate(
        prop = n / sum(n)
    ) |>
    filter(relig_imp_rc == 1) |> 
    mutate(
        post_2019 = ifelse(year %in% 2019:2024, 1, 0)
    )

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_
        )
    )

df_prayer <- data |> 
    filter(religion %in% c(3)) |>
    group_by(year) |>
    count(prayer_rc, wt = weight) |> 
    drop_na() |> 
    mutate(
        prop = n / sum(n)
    ) |>
    filter(prayer_rc == "Several Times a Day") |> 
    mutate(
        post_2019 = ifelse(year %in% 2019:2024, 1, 0)
    )

lm(prop ~ year + post_2019, data = df_prayer)

# Add predictions to dataframes
model_church <- lm(prop ~ year + post_2019, data = df_church)
df_church$predicted <- predict(model_church)

model_imp <- lm(prop ~ year + post_2019, data = df_imp)
df_imp$predicted <- predict(model_imp)

model_prayer <- lm(prop ~ year + post_2019, data = df_prayer)
df_prayer$predicted <- predict(model_prayer)

plot <- df_church |> 
    mutate(
        after_2_hour_change = factor(ifelse(year >= 2019, "Post-2019", "Pre-2019"), levels = c("Pre-2019", "Post-2019"))
    ) |> 
    ggplot(aes(x = year, y = predicted)) +
    geom_line(aes(color = after_2_hour_change)) +
    geom_point(aes(y = prop), shape = 1, size = 2) +
    scale_color_manual(values = c("Pre-2019" = "steelblue", "Post-2019" = "darkorange")) +
    labs(
        title = "Mormon Only: Church Attendance by Year",
        subtitle = "Percent Answering 'Weekly or More'",
        x = "Year",
        y = "Proportion",
        caption = "@mormon_metrics\nSource: Cooperative Election Study (CES) 2008-24"
    ) +
    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_x_continuous(breaks = seq(2008, 2024, by = 2)) +
    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 = 16, face = "bold"),
        legend.title = element_blank(),
        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")
    )

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

plot <- df_imp |> 
    mutate(
        after_2_hour_change = factor(ifelse(year >= 2019, "Post-2019", "Pre-2019"), levels = c("Pre-2019", "Post-2019"))
    ) |> 
    ggplot(aes(x = year, y = predicted)) +
    geom_line(aes(color = after_2_hour_change)) +
    geom_point(aes(y = prop), shape = 1, size = 2) +
    scale_color_manual(values = c("Pre-2019" = "steelblue", "Post-2019" = "darkorange")) +
    labs(
        title = "Mormon Only: Religious Importance by Year",
        subtitle = "Percent Answering 'Very Important'",
        x = "Year",
        y = "Proportion",
        caption = "@mormon_metrics\nSource: Cooperative Election Study (CES) 2008-24"
    ) +
    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_x_continuous(breaks = seq(2008, 2024, by = 2)) +
    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 = 16, face = "bold"),
        legend.title = element_blank(),
        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")
    )

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

plot <- df_prayer |> 
    mutate(
        after_2_hour_change = factor(ifelse(year >= 2019, "Post-2019", "Pre-2019"), levels = c("Pre-2019", "Post-2019"))
    ) |> 
    ggplot(aes(x = year, y = predicted)) +
    geom_line(aes(color = after_2_hour_change)) +
    geom_point(aes(y = prop), shape = 1, size = 2) +
    scale_color_manual(values = c("Pre-2019" = "steelblue", "Post-2019" = "darkorange")) +
    labs(
        title = "Mormon Only: Prayer Frequency by Year",
        subtitle = "Percent Answering 'Several Times a Day'",
        x = "Year",
        y = "Proportion",
        caption = "@mormon_metrics\nSource: Cooperative Election Study (CES) 2008-24"
    ) +
    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_x_continuous(breaks = seq(2008, 2024, by = 2)) +
    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 = 16, face = "bold"),
        legend.title = element_blank(),
        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")
    )

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