ct_selection_untitled-2D85D6D7.R

library(rio)
library(janitor)

pew <- import("E://data/rls24.sav")

pew <- pew %>% clean_names()

gg1 <- pew %>% 
   mutate(decade = frcode(
    birthdecade == 1 ~ "1940s",
    birthdecade == 2 ~ "1950s",
    birthdecade == 3 ~ "1960s",
    birthdecade == 4 ~ "1970s",
    birthdecade == 5 ~ "1980s",
    birthdecade == 6 ~ "1990s",
    birthdecade == 7 ~ "2000s"
  )) %>% 
  mutate(ed = frcode(educrec == 1 ~ "High School", 
                     educrec == 2 ~ "Some College")) %>% 
  mutate(rel = frcode(reltrad == 1100 | reltrad == 1200 | reltrad == 1300 ~ "Protestant", 
                      reltrad == 10000 ~ "Catholic", 
                      reltrad >= 20000 & reltrad <= 90002 ~ "Other World Religions", 
                      reltrad == 100000 ~ "Non-Religious")) %>% 
  group_by(decade) %>% 
  ct(rel, wt = weight, show_na = FALSE)
  


gg1 %>% 
  filter(decade != "NA") %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = 1, y = pct, fill = fct_rev(rel))) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ decade, ncol =1, strip.position = "left") +
  theme_rb() +
  scale_fill_manual(values = c(
    "Protestant" = "#1b9e77",        # Teal
    "Catholic" = "#d95f02",          # Burnt Orange
    "Other World Religions" = "#7570b3", # Purple
    "Non-Religious" = "#e7298a"      # Magenta
  )) +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = percent) +
  theme(strip.text.y.left = element_text(angle=0)) +
  guides(fill = guide_legend(reverse=T, nrow = 1)) +
  theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  theme(panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank()) +
  geom_text(aes(label = ifelse(pct >.05, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 8, family = "font", color = "black") +
  theme(plot.title = element_text(size = 16)) +
  theme(strip.text.y.left = element_text(angle = 0, hjust = 1)) +
  labs(x = "", y = "", title = "Religious Composition by Decade of Birth", caption = "@ryanburge | Data: Pew Religious Landscape Study, 2023-2024") 
save("pew_rls_decade_reltrad.png", wd = 9, ht = 5)  


gg1 <- pew %>% 
  mutate(decade = frcode(
    birthdecade == 1 ~ "1940s",
    birthdecade == 2 ~ "1950s",
    birthdecade == 3 ~ "1960s",
    birthdecade == 4 ~ "1970s",
    birthdecade == 5 ~ "1980s",
    birthdecade == 6 ~ "1990s",
    birthdecade == 7 ~ "2000s"
  )) %>% 
  mutate(cn = frcode(chrnat == 4 ~ "Strongly Oppose", 
                     chrnat == 3 ~ "Oppose", 
                     chrnat == 2 ~ "Favor", 
                     chrnat == 1 ~ "Strongly Favor")) %>% 
  group_by(decade) %>% 
  ct(cn, wt = weight, show_na = FALSE)


gg1 %>% 
  filter(decade != "NA") %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = 1, y = pct, fill = fct_rev(cn))) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ decade, ncol =1, strip.position = "left") +
  theme_rb() +
  scale_fill_manual(values = c(
    "Strongly Oppose" = "#4575b4",       # Blue
    "Oppose" = "#91bfdb",         # Light Blue
    "Favor" = "#fdae61",  # Orange
    "Strongly Favor" = "#d73027"     # Red
  )) +
  theme(legend.position = "bottom") +
  theme(legend.text = element_text(size = 18)) +
  scale_y_continuous(labels = percent) +
  theme(strip.text.y.left = element_text(angle=0)) +
  guides(fill = guide_legend(reverse=T, nrow = 1)) +
  theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  theme(panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank()) +
  geom_text(aes(label = ifelse(pct >.05, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 9, family = "font", color = "black") +
  geom_text(aes(label = ifelse(pct >.05, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 9, family = "font", color = "black") +
  theme(plot.title = element_text(size = 16)) +
  theme(strip.text.y.left = element_text(angle = 0, hjust = 1)) +
  labs(x = "", y = "", title = "The federal government declaring the U.S. a Christian nation: favor or oppose?", caption = "@ryanburge | Data: Pew Religious Landscape Study, 2023-2024") 
save("pew_rls_decade_xtn_nation.png", wd = 9, ht = 5)  



gg1 <- pew %>% 
  mutate(decade = frcode(
    birthdecade == 1 ~ "1940s",
    birthdecade == 2 ~ "1950s",
    birthdecade == 3 ~ "1960s",
    birthdecade == 4 ~ "1970s",
    birthdecade == 5 ~ "1980s",
    birthdecade == 6 ~ "1990s",
    birthdecade == 7 ~ "2000s"
  )) %>% 
  mutate(reltrad = frcode(reltrad == 1100 ~ "Evangelical", 
                          reltrad == 1200 ~ "Mainline", 
                          reltrad == 1300 ~ "Black Protestant", 
                          reltrad == 10000 ~ "Catholic", 
                          reltrad >= 20000 & reltrad <= 99999 ~ "Other World Religions", 
                          reltrad == 100000 ~ "Non-Religious")) %>% 
  mutate(chrnat = scpry1) %>% 
  mutate(cn = frcode(chrnat == 4 ~ "Strongly Oppose", 
                     chrnat == 3 ~ "Oppose", 
                     chrnat == 2 ~ "Favor", 
                     chrnat == 1 ~ "Strongly Favor")) %>% 
  group_by(reltrad) %>% 
  ct(cn, wt = weight, show_na = FALSE) %>% na.omit()


gg1 %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = 1, y = pct, fill = fct_rev(cn))) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ reltrad, ncol =1, strip.position = "left") +
  theme_rb() +
  scale_fill_manual(values = c(
    "Strongly Oppose" = "#4575b4",       # Blue
    "Oppose" = "#91bfdb",         # Light Blue
    "Favor" = "#fdae61",  # Orange
    "Strongly Favor" = "#d73027"     # Red
  )) +
  theme(legend.position = "bottom") +
  theme(legend.text = element_text(size = 18)) +
  scale_y_continuous(labels = percent) +
  theme(strip.text.y.left = element_text(angle=0)) +
  guides(fill = guide_legend(reverse=T, nrow = 1)) +
  theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  theme(panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank()) +
  geom_text(aes(label = ifelse(pct >.05, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 9, family = "font", color = "black") +
  geom_text(aes(label = ifelse(cn == "Strongly Oppose" & pct >.05, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 9, family = "font", color = "white") +
  geom_text(aes(label = ifelse(cn == "Strongly Favor" & pct >.05, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 9, family = "font", color = "white") +
  theme(plot.title = element_text(size = 16)) +
  theme(strip.text.y.left = element_text(angle = 0, hjust = 1)) +
  labs(x = "", y = "", title = "Allowing public school teachers to lead their classes in prayers that refer to Jesus: favor or oppose?", caption = "@ryanburge | Data: Pew Religious Landscape Study, 2023-2024") 
save("pew_rls_reltrad_pray_jesus.png", wd = 11, ht = 5)  



gg1 <- pew %>% 
  mutate(decade = frcode(
    birthdecade == 1 ~ "1940s",
    birthdecade == 2 ~ "1950s",
    birthdecade == 3 ~ "1960s",
    birthdecade == 4 ~ "1970s",
    birthdecade == 5 ~ "1980s",
    birthdecade == 6 ~ "1990s",
    birthdecade == 7 ~ "2000s"
  )) %>% 
  filter(reltrad <= 10000) %>% 
  mutate(opp = case_when(scpry1 == 3 | scpry1 == 4 ~ 1,
                         scpry1 == 1 | scpry1 == 2 ~ 0)) %>% 
  group_by(decade) %>% 
  mean_ci(opp, wt = weight, ci = .84) %>% na.omit()

lab_bar2 <- function(type, pos = 0, sz = 8, above = TRUE, threshold = 0.32) {
  type <- enquo(type)
  geom_text(
    aes(
      y = if (above) !!type + pos else pos,
      label = paste0(round(!!type, 2) * 100, '%'),
      color = ifelse(!!type > threshold, "white", "black")
    ),
    position = position_dodge(width = 0.9),
    size = sz,
    family = "font",
    show.legend = FALSE
  )
}

gg1 %>% 
  ggplot(aes(x = decade, y = mean, fill = mean)) +
  geom_col(color = 'black') +
  theme_rb() +
  y_pct() + 
  error_bar() +
  scale_fill_gradient(low = "#deebf7", high = "#08519c") +
  lab_bar2(type = mean, above = FALSE, pos = .036, sz = 9) +
  scale_color_identity() +    # Add scale_color_identity here once
  labs(x = "", y = "", title = "Share of Christians Who Oppose School Teachers\nLeading Prayers That Refer to Jesus", caption = "@ryanburge | Data: Pew Religious Landscape Study, 2023-2024")

save("jesus_prayer_christian.png", wd = 6)


gg1 <- pew %>% 
  mutate(decade = frcode(
    birthdecade == 1 ~ "1940s",
    birthdecade == 2 ~ "1950s",
    birthdecade == 3 ~ "1960s",
    birthdecade == 4 ~ "1970s",
    birthdecade == 5 ~ "1980s",
    birthdecade == 6 ~ "1990s",
    birthdecade == 7 ~ "2000s"
  )) %>% 
  mutate(reltrad = frcode(reltrad == 1100 ~ "Evangelical", 
                          reltrad == 1200 ~ "Mainline", 
                          reltrad == 1300 ~ "Black Protestant", 
                          reltrad == 10000 ~ "Catholic", 
                          reltrad >= 20000 & reltrad <= 99999 ~ "Other World Religions", 
                          reltrad == 100000 ~ "Non-Religious")) %>% 
  mutate(opp = case_when(scpry1 == 3 | scpry1 == 4 ~ 1,
                         scpry1 == 1 | scpry1 == 2 ~ 0)) %>% 
  group_by(decade, reltrad) %>% 
  mean_ci(opp, wt = weight, ci = .84) %>% na.omit()


gg1 %>% 
  ggplot(., aes(x = decade, y = mean, fill = reltrad)) +
  geom_col(color = 'black') +
  theme_rb() +
  y_pct() + 
  scale_fill_manual(values = c(
    "Evangelical" = "#4c78a8",         # Pew Blue
    "Mainline" = "#72b7b2",            # Soft Teal
    "Black Protestant" = "#f58518",    # Burnt Orange
    "Catholic" = "#e45756",            # Salmon Red
    "Other World Religions" = "#54a24b", # Muted Green
    "Non-Religious" = "#bab0ac"        # Gray
  )) +
  facet_wrap(~ reltrad) +
  error_bar() +
  theme(axis.text.x = element_text(size = 10)) +
  lab_bar_white(above = FALSE, type = mean, pos = .05, sz = 4.5) +
  labs(x = "Decade of Birth", y = "", title = "Share Who Oppose School Teachers Leading Prayers That Refer to Jesus", captoin = "@ryanburge | Data: Pew Religious Landscape Study, 2023-2024")
save("jesus_prayer_christian_reltrad.png")






gg1 <- pew %>% 
  mutate(decade = frcode(
    birthdecade == 1 ~ "1940s",
    birthdecade == 2 ~ "1950s",
    birthdecade == 3 ~ "1960s",
    birthdecade == 4 ~ "1970s",
    birthdecade == 5 ~ "1980s",
    birthdecade == 6 ~ "1990s",
    birthdecade == 7 ~ "2000s"
  )) %>% 
  mutate(reltrad = frcode(reltrad == 1100 ~ "Evangelical", 
                          reltrad == 1200 ~ "Mainline", 
                          reltrad == 1300 ~ "Black Protestant", 
                          reltrad == 10000 ~ "Catholic", 
                          reltrad >= 20000 & reltrad <= 99999 ~ "Other World Religions", 
                          reltrad == 100000 ~ "Non-Religious")) %>% 
  mutate(chrnat = scpry1) %>% 
  mutate(cn = frcode(chrnat == 4 ~ "Strongly Oppose", 
                     chrnat == 3 ~ "Oppose", 
                     chrnat == 2 ~ "Favor", 
                     chrnat == 1 ~ "Strongly Favor")) %>% 
  # group_by(decade) %>% 
  ct(cn, wt = weight, show_na = FALSE) %>% na.omit()


gg1 %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = 1, y = pct, fill = fct_rev(cn))) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ decade, ncol =1, strip.position = "left") +
  theme_rb() +
  scale_fill_manual(values = c(
    "Strongly Oppose" = "#4575b4",       # Blue
    "Oppose" = "#91bfdb",         # Light Blue
    "Favor" = "#fdae61",  # Orange
    "Strongly Favor" = "#d73027"     # Red
  )) +
  theme(legend.position = "bottom") +
  theme(legend.text = element_text(size = 18)) +
  scale_y_continuous(labels = percent) +
  theme(strip.text.y.left = element_text(angle=0)) +
  guides(fill = guide_legend(reverse=T, nrow = 1)) +
  theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  theme(panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank()) +
  geom_text(aes(label = ifelse(pct >.05, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 9, family = "font", color = "black") +
  geom_text(aes(label = ifelse(cn == "Strongly Oppose", paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 9, family = "font", color = "white") +
  geom_text(aes(label = ifelse(cn == "Strongly Favor", paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 9, family = "font", color = "white") +
  theme(plot.title = element_text(size = 16)) +
  theme(strip.text.y.left = element_text(angle = 0, hjust = 1)) +
  labs(x = "", y = "", title = "Allowing public school teachers to lead their classes in prayers that refer to Jesus: favor or oppose?", caption = "@ryanburge | Data: Pew Religious Landscape Study, 2023-2024") 
save("pew_rls_decades_pray_jesus.png", wd = 11, ht = 5)  



gg1 <- pew %>% 
  mutate(
    race = frcode(
      hisp == 2 & racecmb == 1 ~ "White",       # Non-Hispanic White
      hisp == 2 & racecmb == 2 ~ "Black",       # Non-Hispanic Black
      hisp == 1 ~ "Hispanic",                   # Hispanic, any race
      hisp == 2 & racecmb == 3 ~ "Asian",       # Non-Hispanic Asian
      hisp == 2 & racecmb == 4 ~ "Mixed Race",  # Non-Hispanic Mixed Race
      hisp == 2 & racecmb == 5 ~ "Other Race"   # Non-Hispanic Other Race
    )
  ) %>% 
  filter(reltrad <= 10000) %>% 
  mutate(opp = case_when(scpry1 == 3 | scpry1 == 4 ~ 1,
                         scpry1 == 1 | scpry1 == 2 ~ 0)) %>% 
  group_by(race) %>% 
  mean_ci(opp, wt = weight, ci = .84) %>% na.omit()


gg1 %>% 
  ggplot(., aes(x = race, y = mean, fill = race)) +
  geom_col(color = 'black') +
  theme_rb() +
  y_pct() + 
  error_bar() +
  scale_fill_manual(values = c(
    "White" = "#1f78b4",       # Bright Blue
    "Black" = "#e31a1c",       # Bold Red
    "Hispanic" = "#33a02c",    # Bright Green
    "Asian" = "#ff7f00",       # Orange
    "Mixed Race" = "#6a3d9a",  # Purple
    "Other Race" = "#b15928"   # Brown
  )) +
  lab_bar(above = FALSE, type = mean, pos = .025, sz = 11) +
  labs(x = "", y = "", title = "Share of Christians Who Oppose School Teachers\nLeading Prayers That Refer to Jesus", captoin = "@ryanburge | Data: Pew Religious Landscape Study, 2023-2024")
save("jesus_prayer_christian_race.png", wd = 6)


### REGRESSION #####



regg <- pew %>% 
  filter(reltrad <= 10000) %>% 
  mutate(opp = case_when(scpry1 == 3 | scpry1 == 4 ~ 1,
                         scpry1 == 1 | scpry1 == 2 ~ 0)) %>% 
  mutate(white = case_when(hisp == 2 & racecmb == 1 ~ 1, 
                           TRUE ~ 0)) %>% 
  mutate(male = case_when(gender == 1 ~ 1, TRUE ~ 0)) %>% 
  mutate(rep = case_when(party == 1 ~ 1, TRUE ~ 0)) %>% 
  mutate(educ = educrec) %>% 
  mutate(age = 8 - birthdecade) %>% 
  mutate(income = inc_sdt1) %>% 
  mutate(attend = 7 - attndperrls) %>% 
  select(opp, white, male, rep, educ, age, income, attend)
  

regg_clean <- pew %>% 
  filter(reltrad <= 10000) %>% 
  mutate(opp = case_when(scpry1 == 3 | scpry1 == 4 ~ 1,
                         scpry1 == 1 | scpry1 == 2 ~ 0)) %>% 
  mutate(white = case_when(hisp == 2 & racecmb == 1 ~ 1, 
                           TRUE ~ 0)) %>% 
  mutate(male = case_when(gender == 1 ~ 1, TRUE ~ 0)) %>% 
  mutate(rep = case_when(party == 1 ~ 1, TRUE ~ 0)) %>% 
  mutate(educ = educrec) %>% 
  mutate(age = 8 - birthdecade) %>% 
  mutate(income = inc_sdt1) %>% 
  mutate(attend = 7 - attndperrls) %>% 
  select(opp, white, male, rep, educ, age, income, attend) %>%
  filter(
    income != 99,      # Remove "don't know/refused" income responses
    educ != 99,        # Remove "don't know/refused" education responses  
    age != -91,        # Remove invalid age values
    attend >= 1,       # Remove invalid attendance values (should be 1-6)
    !is.na(opp)        # Remove missing values in the outcome variable
  )

# Check the cleaned data
nrow(regg_clean)  # See how many observations remain
summary(regg_clean)  # Check the range of values





# Load required libraries
library(tidyverse)
library(broom)
library(performance)
library(see)

# 1. Basic Logistic Regression
model1 <- glm(opp ~ white + male + rep + educ + age + income + attend, 
              data = regg_clean, 
              family = binomial)

# View results
summary(model1)
tidy(model1, conf.int = TRUE, exponentiate = TRUE)

# 2. Model with interactions (religion and politics often interact)
model2 <- glm(opp ~ white + male + rep + educ + age + income + attend + 
                rep*attend + white*rep, 
              data = regg_clean, 
              family = binomial)

# 3. Linear Probability Model (easier interpretation)
model3 <- lm(opp ~ white + male + rep + educ + age + income + attend, 
             data = regg_clean)

# Compare model performance
compare_performance(model1, model2, model3)




# Analyze the interaction effects from Model 2
library(tidyverse)
library(broom)

# Look at the interaction coefficients
tidy(model2, conf.int = TRUE, exponentiate = TRUE) %>%
  filter(str_detect(term, ":")) %>%
  mutate(
    interaction = case_when(
      term == "rep:attend" ~ "Republican × Religious Attendance",
      term == "white:rep" ~ "White × Republican"
    )
  ) %>%
  select(interaction, estimate, conf.low, conf.high, p.value)

# Create prediction scenarios to understand interactions
# 1. Republican effect by religious attendance level
rep_attend_data <- expand_grid(
  white = 1,     # White
  male = 0,      # Female  
  rep = c(0, 1), # Republican vs Non-Republican
  educ = 2,      # Some college
  age = 5,       # Middle age
  income = 4,    # Middle income
  attend = 1:6   # All attendance levels
)

rep_attend_data$pred_prob <- predict(model2, rep_attend_data, type = "response")

# Visualize Republican effect by attendance
rep_attend_data %>%
  mutate(
    Party = ifelse(rep == 1, "Republican", "Non-Republican"),
    Attendance = case_when(
      attend == 1 ~ "Never",
      attend == 2 ~ "Seldom", 
      attend == 3 ~ "Few times/year",
      attend == 4 ~ "Once/twice month",
      attend == 5 ~ "Once/week",
      attend == 6 ~ "More than weekly"
    )
  ) %>%
  ggplot(aes(x = factor(attend), y = pred_prob, fill = Party)) +
  geom_col(position = "dodge", color = "black") +
  scale_fill_manual(values = c("Republican" = "#d73027", "Non-Republican" = "#4575b4")) +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels = c("Never", "Seldom", "Few/year", "1-2/month", "Weekly", ">Weekly")) +
  labs(
    title = "Predicted Opposition to Jesus Prayers by Party and Attendance",
    subtitle = "Among White Christian Women (middle-aged, some college, middle income)",
    x = "Religious Attendance",
    y = "Predicted Probability of Opposition"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 2. White vs Non-White Republicans
race_party_data <- expand_grid(
  white = c(0, 1),  # White vs Non-White
  male = c(0, 1),   # Male vs Female
  rep = c(0, 1),    # Republican vs Non-Republican
  educ = 2,         # Some college
  age = 5,          # Middle age
  income = 4,       # Middle income
  attend = 3        # Moderate attendance
)

race_party_data$pred_prob <- predict(model2, race_party_data, type = "response")

# Visualize race-party interaction
race_party_data %>%
  mutate(
    Race = ifelse(white == 1, "White", "Non-White"),
    Gender = ifelse(male == 1, "Male", "Female"),
    Party = ifelse(rep == 1, "Republican", "Non-Republican")
  ) %>%
  unite("Profile", Race, Gender, sep = " ") %>%
  ggplot(aes(x = Profile, y = pred_prob, fill = Party)) +
  geom_col(position = "dodge", color = "black") +
  scale_fill_manual(values = c("Republican" = "#d73027", "Non-Republican" = "#4575b4")) +
  scale_y_continuous(labels = scales::percent) +
  labs(
    title = "Predicted Opposition by Race, Gender, and Party",
    subtitle = "Among Christians with moderate attendance",
    x = "",
    y = "Predicted Probability of Opposition"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 3. Show the main effects and interactions in a table
model2_summary <- tidy(model2, conf.int = TRUE, exponentiate = TRUE) %>%
  mutate(
    Variable = case_when(
      term == "(Intercept)" ~ "Intercept",
      term == "white" ~ "White (vs Non-White)",
      term == "male" ~ "Male (vs Female)",
      term == "rep" ~ "Republican (vs Non-Republican)", 
      term == "educ" ~ "Education Level",
      term == "age" ~ "Age Group",
      term == "income" ~ "Income Level",
      term == "attend" ~ "Religious Attendance",
      term == "rep:attend" ~ "Republican × Attendance",
      term == "white:rep" ~ "White × Republican",
      TRUE ~ term
    ),
    OR = round(estimate, 3),
    CI_Lower = round(conf.low, 3),
    CI_Upper = round(conf.high, 3),
    P_Value = round(p.value, 4),
    Significance = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**", 
      p.value < 0.05 ~ "*",
      p.value < 0.1 ~ ".",
      TRUE ~ ""
    )
  ) %>%
  select(Variable, OR, CI_Lower, CI_Upper, P_Value, Significance)

print("Model 2 Results (Odds Ratios):")
print(model2_summary)

cat("\n*** p<0.001, ** p<0.01, * p<0.05, . p<0.1")
cat("\nOdds Ratios > 1 = Increased odds of opposition")
cat("\nOdds Ratios < 1 = Decreased odds of opposition")





library(tidyverse)
library(scales)

# Create a clean comparison showing the white vs non-white difference
# Using average values for other predictors to isolate the race effect

comparison_data <- tibble(
  white = c(0, 1),
  male = 0.42,    # Average from your data (42% male)
  rep = 0.39,     # Average from your data (39% Republican) 
  educ = 2.5,     # Average education level
  age = 4.9,      # Average age group
  income = 4.5,   # Average income level
  attend = 3.6    # Average attendance level
)

# Get predicted probabilities
comparison_data$pred_prob <- predict(model2, comparison_data, type = "response")

# Also get confidence intervals using bootstrap or manual calculation
# For simplicity, we'll calculate standard errors manually
pred_logit <- predict(model2, comparison_data, type = "link", se.fit = TRUE)
comparison_data$pred_logit <- pred_logit$fit
comparison_data$se_logit <- pred_logit$se.fit

# Calculate confidence intervals on logit scale, then transform
comparison_data <- comparison_data %>%
  mutate(
    ci_lower_logit = pred_logit - 1.96 * se_logit,
    ci_upper_logit = pred_logit + 1.96 * se_logit,
    ci_lower = plogis(ci_lower_logit),
    ci_upper = plogis(ci_upper_logit),
    Race = ifelse(white == 1, "White Christians", "Non-White Christians"),
    difference = pred_prob[2] - pred_prob[1]  # Will calculate separately
  )

# Calculate the difference for annotation
white_prob <- comparison_data$pred_prob[comparison_data$white == 1]
nonwhite_prob <- comparison_data$pred_prob[comparison_data$white == 0]
percentage_point_diff <- (white_prob - nonwhite_prob) * 100
relative_diff <- (white_prob / nonwhite_prob - 1) * 100

# Create the main comparison plot
p1 <- comparison_data %>%
  ggplot(aes(x = Race, y = pred_prob, fill = Race)) +
  geom_col(color = "black", width = 0.6, alpha = 0.8) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), 
                width = 0.2, color = "black", size = 1) +
  scale_fill_manual(values = c(
    "White Christians" = "#d73027",      # Red
    "Non-White Christians" = "#4575b4"   # Blue
  )) +
  scale_y_continuous(labels = percent_format(accuracy = 1),
                     limits = c(0, max(comparison_data$ci_upper) * 1.1)) +
  labs(
    title = "Opposition to Jesus Prayers in Public Schools",
    subtitle = "Among Christians Only",
    y = "Predicted Probability of Opposition",
    x = "",
    caption = "Controlling for gender, party, education, age, income, and religious attendance\nError bars show 95% confidence intervals"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    plot.title = element_text(size = 16, face = "bold"),
    plot.subtitle = element_text(size = 12),
    axis.text.x = element_text(size = 12, face = "bold"),
    axis.text.y = element_text(size = 11),
    axis.title.y = element_text(size = 12),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  # Add percentage labels on bars
  geom_text(aes(label = paste0(round(pred_prob * 100, 1), "%")), 
            vjust = -0.5, size = 5, fontface = "bold") +
  # Add annotation about the difference
  annotate("text", x = 1.5, y = max(comparison_data$pred_prob) * 0.9,
           label = paste0("Difference: +", round(percentage_point_diff, 1), 
                          " percentage points\n(", round(relative_diff, 0), "% higher)"),
           size = 4, fontface = "bold", hjust = 0.5,
           color = "black")

print(p1)

# Also create a version broken down by party to show the interaction
party_race_data <- expand_grid(
  white = c(0, 1),
  male = 0.42,    
  rep = c(0, 1),   # Show both parties
  educ = 2.5,     
  age = 4.9,      
  income = 4.5,   
  attend = 3.6    
)

party_race_data$pred_prob <- predict(model2, party_race_data, type = "response")

party_race_data %>%
  mutate(
    Race = ifelse(white == 1, "White", "Non-White"),
    Party = ifelse(rep == 1, "Republican", "Non-Republican")
  ) %>%
  ggplot(aes(x = Race, y = pred_prob, fill = Party)) +
  geom_col(position = "dodge", color = "black", width = 0.7, alpha = 0.8) +
  scale_fill_manual(values = c(
    "Republican" = "#d73027",      # Red
    "Non-Republican" = "#4575b4"   # Blue
  )) +
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, .5)) +
  labs(
    title = "Opposition to Jesus Prayers: Race and Party Effects",
    subtitle = "Among Christians Only", 
    y = "Predicted Probability of Opposition",
    x = "",
    fill = "Party Affiliation"
  ) +
  theme_rb() +
  theme(
    plot.title = element_text(size = 16, face = "bold"),
    plot.subtitle = element_text(size = 12),
    axis.text.x = element_text(size = 12),
    legend.position = "bottom",
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  geom_text(aes(label = paste0(round(pred_prob * 100, 1), "%")), 
            position = position_dodge(width = 0.7),
            vjust = -0.5, size = 8, fontface = "bold", family = "font")
save("jesus_prayer_race_regression.png", wd = 6)

print(p2)

# Print the key statistics
cat("\n=== KEY STATISTICS ===")
cat(paste0("\nWhite Christians: ", round(white_prob * 100, 1), "% oppose"))
cat(paste0("\nNon-White Christians: ", round(nonwhite_prob * 100, 1), "% oppose"))
cat(paste0("\nDifference: +", round(percentage_point_diff, 1), " percentage points"))
cat(paste0("\nWhite Christians are ", round(relative_diff, 0), "% more likely to oppose"))
cat(paste0("\nOdds Ratio: ", round(1.94, 2), " (", round((1.94-1)*100, 0), "% higher odds)"))