ct_selection_untitled-050EF50D.R

library(rio)
library(janitor)

pew <- import("E://data/atp124.sav") %>% clean_names() %>% as_tibble()


gg1 <- pew %>% 
  mutate(rel = uscomp_relig_w124) %>% 
  mutate(rel = frcode(rel == 1 ~ "More Religious", 
                      rel == 3 ~ "About the Same", 
                      rel == 2 ~ "Less Religious")) %>% 
  ct(rel, wt = weight_w124, show_na = FALSE) %>% 
  mutate(reltrad = "Entire Sample")


gg2 <- pew %>% 
  mutate(rel = uscomp_relig_w124) %>% 
  mutate(rel = frcode(rel == 1 ~ "More Religious", 
                      rel == 3 ~ "About the Same", 
                      rel == 2 ~ "Less Religious")) %>% 
  mutate(reltrad = f_religcat1) %>% 
  mutate(reltrad = frcode(reltrad == 1 & f_born == 1 ~ "Evangelical Protestant", 
                          reltrad == 1 & f_born == 2 ~ "Non-Evangelical Prot.", 
                          reltrad == 2 ~ "Catholic", 
                          reltrad == 3 ~ "No Religion", 
                          reltrad == 4 ~ "Something Else")) %>% 
  group_by(reltrad) %>% 
  ct(rel, wt = weight_w124, show_na = FALSE) %>% na.omit()

both <- bind_rows(gg1, gg2)


library(dplyr)
library(forcats)
library(ggplot2)
library(scales)

# --- Order facets: Entire Sample first, then by 'More Religious' (desc)
ord_levels <- both %>%
  filter(rel == "More Religious", reltrad != "Entire Sample") %>%
  arrange(desc(pct)) %>%
  pull(reltrad) %>%
  unique()

facet_levels <- c("Entire Sample", ord_levels)

both_ord <- both %>%
  mutate(reltrad = factor(reltrad, levels = facet_levels))

# --- Plot
both_ord %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(aes(x = 1, y = pct, fill = rel)) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ reltrad, ncol = 1, strip.position = "left") +
  # palette: green = More Religious, gray = About the Same, orange/red = Less Religious
  scale_fill_manual(values = c(
    "More Religious"  = "#1B9E77",
    "About the Same"  = "#9E9E9E",
    "Less Religious"  = "#D95F02"
  )) +
  theme_rb() +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = percent) +
  theme(strip.text.y.left = element_text(angle = 0, hjust = 1)) +
  guides(fill = guide_legend(reverse = TRUE, nrow = 1)) +
  theme(axis.title.y = element_blank(),
        axis.text.y  = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_blank(),
        plot.title = element_text(size = 16)) +
  geom_text(
    aes(label = ifelse(pct > .05, paste0(lab * 100, "%"), "")),
    position = position_stack(vjust = 0.5),
    size = 7, family = "font", color = "black"
  ) +
  labs(
    x = "", y = "",
    title = "How does the U.S. compare to other wealthy nations? The U.S. is…",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  )

save("us_relig_compare.png", wd = 9, ht = 4)




gg1 <- pew %>% 
  mutate(nat = nationality_relig_w124) %>% 
  mutate(nat = frcode(nat == 4 ~ "Not at all", 
                      nat == 3 ~ "Not too", 
                      nat == 2 ~ "Somewhat", 
                      nat == 1 ~ "Very")) %>% 
  ct(nat, wt = weight_w124, show_na = FALSE) %>% 
  mutate(reltrad = "Entire Sample")


gg2 <- pew %>% 
  mutate(nat = nationality_relig_w124) %>% 
  mutate(nat = frcode(nat == 4 ~ "Not at all", 
                      nat == 3 ~ "Not too", 
                      nat == 2 ~ "Somewhat", 
                      nat == 1 ~ "Very")) %>% 
  mutate(reltrad = f_religcat1) %>% 
  mutate(reltrad = frcode(reltrad == 1 & f_born == 1 ~ "Evangelical Protestant", 
                          reltrad == 1 & f_born == 2 ~ "Non-Evangelical Prot.", 
                          reltrad == 2 ~ "Catholic", 
                          reltrad == 3 ~ "No Religion", 
                          reltrad == 4 ~ "Something Else")) %>% 
  group_by(reltrad) %>% 
  ct(nat, wt = weight_w124, show_na = FALSE) %>% na.omit()

both <- bind_rows(gg1, gg2)


library(dplyr)
library(forcats)
library(ggplot2)
library(scales)

# Ensure nat is ordered low→high so "Very" stacks on top
both_nat <- both %>%
  mutate(nat = factor(nat, levels = c("Not at all", "Not too", "Somewhat", "Very")))

# Order facets: Entire Sample first, then by 'Very' share (desc)
ord_levels_nat <- both_nat %>%
  filter(nat == "Very", reltrad != "Entire Sample") %>%
  arrange(desc(pct)) %>%
  pull(reltrad) %>%
  unique()

facet_levels_nat <- c("Entire Sample", ord_levels_nat)

plot_df <- both_nat %>%
  mutate(reltrad = factor(reltrad, levels = facet_levels_nat)) %>%
  mutate(lab = round(pct, 2))

# New color scheme (ordinal, color-blind friendly)
# darkest = "Very" → lightest = "Not at all"
fills_nat <- c(
  "Not at all" = "#E6F0FA",
  "Not too"    = "#B3CDE3",
  "Somewhat"   = "#5C88C4",
  "Very"       = "#27476E"
)

plot_df %>%
  ggplot(aes(x = 1, y = pct, fill = fct_rev(nat))) +
  geom_col(color = "black") +
  coord_flip() +
  facet_wrap(~ reltrad, ncol = 1, strip.position = "left") +
  scale_fill_manual(values = fills_nat, breaks = c("Very","Somewhat","Not too","Not at all")) +
  theme_rb() +
  theme(
    legend.position = "bottom",
    strip.text.y.left = element_text(angle = 0, hjust = 1),
    axis.title.y = element_blank(),
    axis.text.y  = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.y = element_blank(),
    plot.title = element_text(size = 16)
  ) +
  scale_y_continuous(labels = percent) +
  guides(fill = guide_legend(nrow = 1, reverse = TRUE)) +
  geom_text(
    aes(label = ifelse(pct > .05, paste0(lab * 100, "%"), "")),
    position = position_stack(vjust = 0.5),
    size = 7, family = "font", color = "black"
  ) +
  geom_text(
    aes(label = ifelse(pct > .05 & nat == "Very", paste0(lab * 100, "%"), "")),
    position = position_stack(vjust = 0.5),
    size = 7, family = "font", color = "white"
  ) +
  labs(
    x = "", y = "",
    title = "How important is being a Christian to be truly American?",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  )

save("truly_american_christian_importance.png", wd = 9, ht = 4)



# --- Build by AGE x RELTRAD ---

# 1) Entire-sample within each age group
gg1_age <- pew %>% 
  mutate(nat = nationality_relig_w124) %>% 
  mutate(nat = frcode(nat == 4 ~ "Not at all", 
                      nat == 3 ~ "Not too", 
                      nat == 2 ~ "Somewhat", 
                      nat == 1 ~ "Very")) %>% 
  mutate(agecat = frcode(f_agecat == 1 ~ "18–29",
                         f_agecat == 2 ~ "30–49",
                         f_agecat == 3 ~ "50–64",
                         f_agecat == 4 ~ "65+")) %>%
  group_by(agecat) %>%
  ct(nat, wt = weight_w124, show_na = FALSE) %>% 
  mutate(reltrad = "Entire Sample")

# 2) Age x Reltrad breakdown
gg2_age <- pew %>% 
  mutate(nat = nationality_relig_w124) %>% 
  mutate(nat = frcode(nat == 4 ~ "Not at all", 
                      nat == 3 ~ "Not too", 
                      nat == 2 ~ "Somewhat", 
                      nat == 1 ~ "Very")) %>% 
  mutate(agecat = frcode(f_agecat == 1 ~ "18–29",
                         f_agecat == 2 ~ "30–49",
                         f_agecat == 3 ~ "50–64",
                         f_agecat == 4 ~ "65+")) %>%
  mutate(reltrad = f_religcat1) %>% 
  mutate(reltrad = frcode(reltrad == 1 & f_born == 1 ~ "Evangelical Protestant", 
                          reltrad == 1 & f_born == 2 ~ "Non-Evangelical Prot.", 
                          reltrad == 2 ~ "Catholic", 
                          reltrad == 3 ~ "No Religion", 
                          reltrad == 4 ~ "Something Else")) %>% 
  group_by(agecat, reltrad) %>% 
  ct(nat, wt = weight_w124, show_na = FALSE) %>% 
  na.omit()

both_age <- bind_rows(gg1_age, gg2_age)

# --- Plot prep (order, palette, labels) ---
library(dplyr)
library(forcats)
library(ggplot2)
library(scales)

# Ensure stack order: light→dark with "Very" on top in legend
both_age <- both_age %>%
  mutate(nat = factor(nat, levels = c("Not at all", "Not too", "Somewhat", "Very")))

# Order reltrad columns: Entire Sample first, then by overall share "Very" (desc)
ord_cols <- both_age %>%
  filter(nat == "Very", reltrad != "Entire Sample") %>%
  arrange(desc(pct)) %>%
  pull(reltrad) %>%
  unique()

reltrad_levels <- c("Entire Sample", ord_cols)

plot_df <- both_age %>%
  mutate(reltrad = factor(reltrad, levels = reltrad_levels)) %>%
  mutate(lab = round(pct, 2))

fills_nat <- c(
  "Not at all" = "#E6F0FA",
  "Not too"    = "#B3CDE3",
  "Somewhat"   = "#5C88C4",
  "Very"       = "#27476E"
)

plot_df <- plot_df %>%
  mutate(reltrad = factor(reltrad,
                          levels = c("Entire Sample",
                                     "Evangelical Protestant",
                                     "Non-Evangelical Prot.",
                                     "Catholic",
                                     "No Religion")))


# --- Plot: grid of Age (rows) × Reltrad (cols) ---
plot_df %>%
  filter(agecat != "NA") %>% 
  filter(reltrad != "Something Else") %>% 
  ggplot(aes(x = fct_rev(agecat), y = pct, fill = fct_rev(nat))) +
  geom_col(color = "black") +
  coord_flip() +
  facet_wrap(~ reltrad, ncol = 1) +
  scale_fill_manual(values = fills_nat,
                    breaks = c("Very","Somewhat","Not too","Not at all")) +
  theme_rb() +
  theme(
    legend.position = "bottom",
    strip.placement = "outside",
    strip.text.y.left = element_text(angle = 0, hjust = 1),
    axis.title.y = element_blank(),
    # axis.text.y  = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.y = element_blank(),
    plot.title = element_text(size = 22),
    plot.subtitle = element_text(size = 12), 
    legend.text = element_text(size = 24)
  ) +
  scale_y_continuous(labels = percent) +
  guides(fill = guide_legend(nrow = 1, reverse = TRUE)) +
  # Black labels on all large slices
  geom_text(aes(label = ifelse(pct > .05, paste0(lab * 100, "%"), "")),
            position = position_stack(vjust = 0.5),
            size = 6, family = "font", color = "black") +
  # White label overlay for 'Very' to ensure contrast
  geom_text(aes(label = ifelse(pct > .05 & nat == "Very", paste0(lab * 100, "%"), "")),
            position = position_stack(vjust = 0.5),
            size = 6, family = "font", color = "white") +
  labs(
    x = "", y = "",
    title = "How important is being a Christian to be truly American?",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  )

save("truly_american_christian_importance_by_age_reltrad.png", wd = 8, ht = 12.5)





pew2 <- pew %>% 
  mutate(
    born_imp = frcode(nationality_born_w124 %in% c(1,2) ~ 1, 
                      nationality_born_w124 %in% c(3,4) ~ 0),
    lang_imp = frcode(nationality_lang_w124 %in% c(1,2) ~ 1, 
                      nationality_lang_w124 %in% c(3,4) ~ 0),
    christ_imp = frcode(nationality_relig_w124 %in% c(1,2) ~ 1, 
                        nationality_relig_w124 %in% c(3,4) ~ 0)
  )


pew2 <- pew2 %>% 
  mutate(native_eng_christ = ifelse(born_imp == 1 & lang_imp == 1 & christ_imp == 1, 1, 0))


pew2 <- pew2 %>%
  mutate(reltrad = f_religcat1) %>% 
  mutate(reltrad = frcode(reltrad == 1 & f_born == 1 ~ "Evangelical Protestant", 
                          reltrad == 1 & f_born == 2 ~ "Non-Evangelical Prot.", 
                          reltrad == 2 ~ "Catholic", 
                          reltrad == 3 ~ "No Religion", 
                          reltrad == 4 ~ "Something Else"))



results <- pew2 %>% 
  group_by(reltrad) %>% 
  mean_ci(native_eng_christ, wt = weight_w124) %>% na.omit()

results %>% 
  ggplot(., aes(x = reorder(reltrad, mean), y = mean, fill = reltrad)) +
  geom_col(color = 'black') +
  coord_flip() +
  theme_rb() +
  y_pct() + 
  scale_fill_manual(values = c(
    "Evangelical Protestant"   = "#1B9E77",  # teal green
    "Non-Evangelical Prot."    = "#7570B3",  # purple
    "Catholic"                 = "#E7298A",  # magenta
    "No Religion"              = "#D95F02",  # orange
    "Something Else"           = "#66A61E"   # olive green
  )) +
  lab_bar(above = FALSE, pos = .04, sz = 10, type = mean) +
  error_bar() +
  labs(x = "", y = "", title = "The Share Saying That Being a Native-Born, English Speaking,\nChristian is Important to Being an American", 
       caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124")
save("native_born_american_english.png", ht = 4)




gg1 <- pew %>% 
  mutate(country = countries_visited_w124) %>% 
  mutate(country = frcode(country == 1 ~ "None", 
                          country == 2 ~ "1", 
                          country == 3 ~ "2", 
                          country == 4 ~ "3-4", 
                          country == 5 ~ "5-9", 
                          country == 6 ~ "10+")) %>% 
  ct(country, wt = weight_w124, show_na = FALSE) %>% 
  mutate(reltrad = "Entire Sample")


gg2 <- pew %>% 
  mutate(country = countries_visited_w124) %>% 
  mutate(country = frcode(country == 1 ~ "None", 
                          country == 2 ~ "1", 
                          country == 3 ~ "2", 
                          country == 4 ~ "3-4", 
                          country == 5 ~ "5-9", 
                          country == 6 ~ "10+")) %>% 
  mutate(reltrad = f_religcat1) %>% 
  mutate(reltrad = frcode(reltrad == 1 & f_born == 1 ~ "Evangelical Protestant", 
                          reltrad == 1 & f_born == 2 ~ "Non-Evangelical Prot.", 
                          reltrad == 2 ~ "Catholic", 
                          reltrad == 3 ~ "No Religion", 
                          reltrad == 4 ~ "Something Else")) %>% 
  group_by(reltrad) %>% 
  ct(country, wt = weight_w124, show_na = FALSE) %>% na.omit()


both <- bind_rows(gg1, gg2)


library(dplyr)
library(forcats)
library(ggplot2)
library(scales)

# --- Order facets: Entire Sample first, then by 'More Religious' (desc)
ord_levels <- both %>%
  filter(country == "None", reltrad != "Entire Sample") %>%
  arrange(desc(pct)) %>%
  pull(reltrad) %>%
  unique()

facet_levels <- c("Entire Sample", ord_levels)

both_ord <- both %>%
  mutate(reltrad = factor(reltrad, levels = facet_levels))

# --- Plot
both_ord %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(aes(x = 1, y = pct, fill = fct_rev(country))) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ reltrad, ncol = 1, strip.position = "left") +
  # palette: green = More Religious, gray = About the Same, orange/red = Less Religious
  theme_rb() +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = percent) +
  scale_fill_manual(values = c(
    "None"  = "#f7f4f9",
    "1"     = "#e7e1ef",
    "2"     = "#d4b9da",
    "3-4"   = "#c994c7",
    "5-9"   = "#df65b0",
    "10+"   = "#ce1256"
  )) +
  theme(strip.text.y.left = element_text(angle = 0, hjust = 1)) +
  guides(fill = guide_legend(reverse = TRUE, nrow = 1)) +
  theme(axis.title.y = element_blank(),
        axis.text.y  = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_blank(),
        legend.text = element_text(size = 22), 
        plot.title = element_text(size = 16)) +
  geom_text(
    aes(label = ifelse(pct > .05, paste0(lab * 100, "%"), "")),
    position = position_stack(vjust = 0.5),
    size = 7, family = "font", color = "black"
  ) +
  labs(
    x = "", y = "",
    title = "How many countries have you traveled to outside the United States? ",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  )

save("us_relig_countries_visited.png", wd = 9, ht = 4)



# 1) Composite: native-born + English-speaking + Christian (Very/Somewhat = important)
pew_nec <- pew %>% 
  mutate(
    born_imp   = frcode(nationality_born_w124  %in% c(1,2) ~ 1, nationality_born_w124  %in% c(3,4) ~ 0),
    lang_imp   = frcode(nationality_lang_w124  %in% c(1,2) ~ 1, nationality_lang_w124  %in% c(3,4) ~ 0),
    christ_imp = frcode(nationality_relig_w124 %in% c(1,2) ~ 1, nationality_relig_w124 %in% c(3,4) ~ 0),
    native_eng_christ = as.integer(born_imp == 1 & lang_imp == 1 & christ_imp == 1)
  )

# 2) Countries visited categories (match your earlier bins)
country_levels <- c("None","1","2","3-4","5-9","10+")
pew_nec <- pew_nec %>% 
  mutate(country = frcode(
    countries_visited_w124 == 1 ~ "None",
    countries_visited_w124 == 2 ~ "1",
    countries_visited_w124 == 3 ~ "2",
    countries_visited_w124 == 4 ~ "3-4",
    countries_visited_w124 == 5 ~ "5-9",
    countries_visited_w124 == 6 ~ "10+"
  )) %>%
  mutate(country = factor(country, levels = country_levels))

# 3) Weighted share by countries visited
nec_by_country <- pew_nec %>%
  group_by(country) %>%
  mean_ci(native_eng_christ, wt = weight_w124) %>%
  filter(!is.na(country))

# 4) Plot (sequential blue palette; light = fewer countries, dark = more)
fills_country <- c(
  "None" = "#f7fbff", "1" = "#deebf7", "2" = "#c6dbef",
  "3-4"  = "#9ecae1", "5-9" = "#6baed6", "10+" = "#2171b5"
)

nec_by_country %>%
  ggplot(aes(x = country, y = mean, fill = country)) +
  geom_col(color = "black") +
  error_bar() +
  lab_bar(above = FALSE, pos = 0.03, sz = 9, type = mean) +
  y_pct() +
  scale_fill_manual(values = fills_country) +
  coord_flip() +
  theme_rb(legend = FALSE) +
  theme(plot.title = element_text(size = 16)) +
  labs(
    x = "", y = "",
    title = "Propensity to Say It’s Important to Be a Native-Born, English-Speaking Christian",
    subtitle = "By number of countries visited outside the U.S.",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  )

save("nec_by_countries_visited.png", wd = 7.5, ht = 4.5)



gg1 <- pew %>% 
  filter(want_travel_w124 != "NA") %>% 
  mutate(reltrad = f_religcat1) %>% 
  mutate(reltrad = frcode(reltrad == 1 & f_born == 1 ~ "Evangelical Protestant", 
                          reltrad == 1 & f_born == 2 ~ "Non-Evangelical Prot.", 
                          reltrad == 2 ~ "Catholic", 
                          reltrad == 3 ~ "No Religion", 
                          reltrad == 4 ~ "Something Else")) %>% 
  mutate(want = want_travel_w124) %>% 
  mutate(want = case_when(want == 1 ~ 1, 
                          want == 2 ~ 0)) %>% 
  group_by(reltrad) %>% 
  mean_ci(want, wt = weight_w124, ci = .84) %>% filter(reltrad != "NA")



gg1 %>%
  ggplot(aes(x = reorder(reltrad, mean), y = mean, fill = reltrad)) +
  geom_col(color = "black") +
  error_bar() +
  lab_bar(above = FALSE, pos = 0.05, sz = 9, type = mean) +
  y_pct() +
  scale_fill_manual(values = c(
    "Evangelical Protestant"   = "#1B9E77",  # teal green
    "Non-Evangelical Prot."    = "#7570B3",  # purple
    "Catholic"                 = "#E7298A",  # magenta
    "No Religion"              = "#D95F02",  # orange
    "Something Else"           = "#66A61E"   # olive green
  )) +
  coord_flip() +
  theme_rb(legend = FALSE) +
  theme(plot.title = element_text(size = 16)) +
  labs(
    x = "", y = "",
    title = "If you had the opportunity, would you like to travel outside the United States? ",
    subtitle = "Among People Who Have Never Traveled Internationally",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  )

save("want_to_travel_reltrad.png", wd = 8.5, ht = 4)


pew %>% ct(f_inc_sdt1, wt = weight_w124)


plot_df <- pew %>%
  mutate(
    inc4 = frcode(
      f_inc_sdt1 == 1           ~ "<$30k",
      f_inc_sdt1 %in% 2:4       ~ "$30–59k",
      f_inc_sdt1 %in% 5:8       ~ "$60–99k",
      f_inc_sdt1 == 9           ~ "$100k+")) %>% 
  mutate(country = countries_visited_w124) %>% 
  mutate(country = frcode(country == 1 ~ "None", 
                          country == 2 ~ "1", 
                          country == 3 ~ "2", 
                          country == 4 ~ "3-4", 
                          country == 5 ~ "5-9", 
                          country == 6 ~ "10+")) %>% 
  group_by(inc4) %>% 
  ct(country, wt = weight_w124, show_na = FALSE)



fills_country <- c(
  "None" = "#f7fbff", "1" = "#deebf7", "2" = "#c6dbef",
  "3-4"  = "#9ecae1", "5-9" = "#6baed6", "10+" = "#2171b5"
)

plot_df %>%
  filter(inc4 != "NA") %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(aes(x = fct_rev(inc4), y = pct, fill = fct_rev(country))) +
  geom_col(color = "black") +
  coord_flip() +
  scale_fill_manual(values = fills_country,
                    breaks = c("10+","5-9","3-4","2","1","None")) +
  scale_y_continuous(labels = percent) +
  theme_rb() +
  theme(
    legend.position = "bottom",
    plot.title = element_text(size = 16)
  ) +
  guides(fill = guide_legend(nrow = 1, reverse = TRUE, title = NULL)) +
  # inside-bar labels (only if slice > 5%)
  geom_text(
    aes(label = ifelse(pct > .05, paste0(lab * 100, "%"), "")),
    position = position_stack(vjust = 0.5),
    size = 8, family = "font", color = "black"
  ) +
  labs(
    x = "", y = "",
    title = "Countries Visited Outside the U.S., by Family Income",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  )

save("countries_by_income_stacked.png", wd = 8.5, ht = 4)



pew_ct <- pew %>% 
  mutate(
    country = frcode(
      countries_visited_w124 == 1 ~ "None",
      countries_visited_w124 == 2 ~ "1",
      countries_visited_w124 == 3 ~ "2",
      countries_visited_w124 == 4 ~ "3-4",
      countries_visited_w124 == 5 ~ "5-9",
      countries_visited_w124 == 6 ~ "10+"
    ),
    reltrad = frcode(
      f_religcat1 == 1 & f_born == 1 ~ "Evangelical Protestant",
      f_religcat1 == 1 & f_born == 2 ~ "Non-Evangelical Prot.",
      f_religcat1 == 2               ~ "Catholic",
      f_religcat1 == 3               ~ "No Religion",
      f_religcat1 == 4               ~ "Something Else"
    ),
    inc4 = frcode(
      f_inc_sdt1 == 1           ~ "<$30k",
      f_inc_sdt1 %in% 2:4       ~ "$30–59k",
      f_inc_sdt1 %in% 5:8       ~ "$60–99k",
      f_inc_sdt1 == 9           ~ "$100k+"
    )
  ) %>%
  group_by(reltrad, inc4) %>%
  ct(country, wt = weight_w124, show_na = FALSE) %>%
  na.omit()


pew_ct %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = fct_rev(reltrad), y = pct, fill = fct_rev(country))) +
  geom_col(color = "black") +
  coord_flip() +
  facet_wrap(~ inc4, ncol = 2) +
  scale_y_continuous(labels = percent) +
  scale_fill_manual(values = fills_country) +
  theme_rb() +
  theme(legend.position = "bottom") +
  guides(fill = guide_legend(reverse = TRUE, nrow = 1)) +
  geom_text(
    aes(label = ifelse(pct > .07, paste0(lab * 100, "%"), "")),
    position = position_stack(vjust = 0.5),
    size = 7, family = "font", color = "black"
  ) +
  labs(x = "", y = "", 
       title = "Countries Visited by Religious Tradition, Stratified by Family Income",
       caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124")
save("income_reltrad_countries_visited.png", wd = 12, ht = 7)



gg1 <- pew %>% 
  mutate(country = countries_visited_w124) %>% 
  mutate(country = frcode(country == 1 ~ "None", 
                          country == 2 ~ "1", 
                          country == 3 ~ "2", 
                          country == 4 ~ "3-4", 
                          country == 5 ~ "5-9", 
                          country == 6 ~ "10+")) %>% 
  cces_attend(f_attend) %>% 
  group_by(att) %>% 
  ct(country, wt = weight_w124, show_na = FALSE)




# --- Plot
gg1 %>% 
  filter(att != 'NA') %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(aes(x = 1, y = pct, fill = fct_rev(country))) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ att, ncol = 1, strip.position = "left") +
  # palette: green = More Religious, gray = About the Same, orange/red = Less Religious
  theme_rb() +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = percent) +
  scale_fill_manual(values = c(
    "None"  = "#f7f4f9",
    "1"     = "#e7e1ef",
    "2"     = "#d4b9da",
    "3-4"   = "#c994c7",
    "5-9"   = "#df65b0",
    "10+"   = "#ce1256"
  )) +
  theme(strip.text.y.left = element_text(angle = 0, hjust = 1)) +
  guides(fill = guide_legend(reverse = TRUE, nrow = 1)) +
  theme(axis.title.y = element_blank(),
        axis.text.y  = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_blank(),
        legend.text = element_text(size = 22), 
        plot.title = element_text(size = 16)) +
  geom_text(
    aes(label = ifelse(pct > .05, paste0(lab * 100, "%"), "")),
    position = position_stack(vjust = 0.5),
    size = 7, family = "font", color = "black"
  ) +
  labs(
    x = "", y = "",
    title = "How many countries have you traveled to outside the United States? ",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  )

save("us_relig_countries_visited_attendance.png", wd = 9, ht = 4)




pew_reg <- pew %>%
  mutate(
    never_travel = ifelse(countries_visited_w124 == 1, 1, 0),  # 1 = None, 0 = any travel
    # Religion
    reltrad = frcode(
      f_religcat1 == 1 & f_born == 1 ~ "Evangelical Protestant",
      f_religcat1 == 1 & f_born == 2 ~ "Non-Evangelical Prot.",
      f_religcat1 == 2               ~ "Catholic",
      f_religcat1 == 3               ~ "No Religion",
      f_religcat1 == 4               ~ "Something Else"
    ),
    # Income (collapsed as we set up earlier)
    inc4 = frcode(
      f_inc_sdt1 == 1           ~ "<$30k",
      f_inc_sdt1 %in% 2:4       ~ "$30–59k",
      f_inc_sdt1 %in% 5:8       ~ "$60–99k",
      f_inc_sdt1 == 9           ~ "$100k+"
    ),
    # Gender
    gender = frcode(f_gender == 1 ~ "Male",
                    f_gender == 2 ~ "Female"),
    # Education (using the Pew coding, F_EDUCCAT2 = standard bins)
    educ = frcode(f_educcat2 == 1 ~ "HS or less",
                  f_educcat2 == 2 ~ "Some college",
                  f_educcat2 == 3 ~ "BA",
                  f_educcat2 == 4 ~ "Post-grad"),
    # Race (simple 3-way; check codebook for exact coding)
    race = frcode(f_racecmb == 1 ~ "White",
                  f_racecmb == 2 ~ "Black",
                  f_racecmb == 3 ~ "Hispanic",
                  f_racecmb %in% 4:6 ~ "Other")
  ) %>%
  filter(!is.na(never_travel), !is.na(reltrad), !is.na(inc4),
         !is.na(gender), !is.na(educ), !is.na(race), !is.na(f_agecat))



model1 <- glm(
  never_travel ~ reltrad + f_agecat + inc4 + race + educ + gender,
  family = binomial(),
  weights = weight_w124,
  data = pew_reg
)

summary(model1)


library(broom)
tidy(model1, exponentiate = TRUE, conf.int = TRUE) %>%
  select(term, estimate, conf.low, conf.high, p.value)

library(ggeffects)
library(dplyr)
library(ggplot2)
library(scales)

# Get marginal effects for reltrad and convert to data frame
preds <- ggeffect(model1, terms = "reltrad") %>% as.data.frame()

# Ensure desired order & a reltrad column to map fills
preds <- preds %>%
  mutate(
    reltrad = factor(x, levels = c(
      "Evangelical Protestant",
      "Non-Evangelical Prot.",
      "Catholic",
      "No Religion",
      "Something Else"
    ))
  )

ggplot(preds, aes(x = reorder(reltrad, predicted), y = predicted, fill = reltrad)) +
  geom_col(color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  geom_text(aes(y = .03, label = scales::percent(predicted, accuracy = 1)),
            size = 10, family = "font") +
  coord_flip(ylim = c(0, max(preds$conf.high) * 1.1)) +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c(
    "Evangelical Protestant"   = "#1B9E77",
    "Non-Evangelical Prot."    = "#7570B3",
    "Catholic"                 = "#E7298A",
    "No Religion"              = "#D95F02",
    "Something Else"           = "#66A61E"
  )) +
  labs(
    x = "", y = "Predicted probability",
    title = "Predicted Probability of Never Traveling Outside the U.S.",
    subtitle = "By religious tradition, holding age, income, race, education, and gender constant",
    caption = "@ryanburge + @religiondata | Data: Pew American Trends Panel, Wave 124"
  ) +
  theme_rb(legend = FALSE)


save("never_travel_predprob_reltrad.png", wd = 8.5, ht = 4)

pew %>% 
 filter(f_religcat1 == 2) %>% 
  group_by(f_hisp) %>% 
  mutate(country = countries_visited_w124) %>% 
  mutate(country = frcode(country == 1 ~ "None", 
                          country == 2 ~ "1", 
                          country == 3 ~ "2", 
                          country == 4 ~ "3-4", 
                          country == 5 ~ "5-9", 
                          country == 6 ~ "10+")) %>% 
  ct(country, wt = weight_w124, show_na = FALSE)



plot_df <- pew %>%
  mutate(
    inc4 = frcode(
      f_inc_sdt1 == 1           ~ "<$30k",
      f_inc_sdt1 %in% 2:4       ~ "$30–59k",
      f_inc_sdt1 %in% 5:8       ~ "$60–99k",
      f_inc_sdt1 == 9           ~ "$100k+")) %>% 
  mutate(country = countries_visited_w124) %>% 
  mutate(country = frcode(country == 1 ~ "None", 
                          country == 2 ~ "1", 
                          country == 3 ~ "2", 
                          country == 4 ~ "3-4", 
                          country == 5 ~ "5-9", 
                          country == 6 ~ "10+")) %>% 
  group_by(inc4) %>% 
  ct(country, wt = weight_w124, show_na = FALSE)