Skip to content

Instantly share code, notes, and snippets.

@ryanburge
Created March 4, 2024 01:08
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ryanburge/190eb8bfb5665554b072ea5d8b88cb95 to your computer and use it in GitHub Desktop.
Save ryanburge/190eb8bfb5665554b072ea5d8b88cb95 to your computer and use it in GitHub Desktop.
library(rio)
library(janitor)
hp <- read_csv("E://data/hps.csv") %>% clean_names()
gg <- hp %>%
mutate(msa = est_msa) %>%
mutate(msa = frcode(msa == 35620 ~ "NYC",
msa == 31080 ~ "Charlotte",
msa == 16980 ~ "Chicago",
msa == 19100 ~ "Dallas",
msa == 26420 ~ "Houston",
msa == 47900 ~ "DC",
msa == 33100 ~ "Miami",
msa == 37980 ~ "Philadelphia",
msa == 12060 ~ "Atlanta",
msa == 38060 ~ "Phoenix",
msa == 14460 ~ "Boston",
msa == 41860 ~ "San Francisco",
msa == 40140 ~ "Riverside",
msa == 19820 ~ "Detroit",
msa == 42660 ~ "Seattle")) %>%
mutate(att = support3) %>%
mutate(att = frcode(att == 1 ~ "Less than Once a Year",
att == 2 ~ "1-3/Yr.",
att == 3 ~ "4-11/Yr.",
att == 4 ~ "12+/Yr.")) %>%
group_by(msa) %>%
ct(att, wt = pweight, show_na = FALSE)
lvl <- gg %>%
filter(att == "Less than Once a Year") %>%
select(msa, sort = pct)
gg <- left_join(gg, lvl)
gg$msa <- fct_reorder(gg$msa, gg$sort, .desc = TRUE)
gg %>%
filter(msa != "NA") %>%
mutate(lab = round(pct, 2)) %>%
ggplot(., aes(x = 1, y = pct, fill = fct_rev(att))) +
geom_col(color = "black") +
coord_flip() +
facet_wrap(~ msa, ncol =1, strip.position = "left") +
theme_rb() +
scale_fill_manual(values = c(moma.colors("ustwo", 4))) +
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 >.03, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 4, family = "font", color = "black") +
geom_text(aes(label = ifelse(att == "Less than Once a Year", paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 4, family = "font", color = "white") +
# geom_text(aes(label = ifelse(att == "12+/Yr.", paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 4, 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 = "Religious Attendance in the 15 Largest Metro Areas", caption = "@ryanburge\nData: Household Pulse Survey, 2024")
save("att_msa_pulse.png", wd = 9, ht = 6)
gg <- hp %>%
mutate(msa = est_msa) %>%
mutate(msa = frcode(msa == 35620 ~ "NYC",
msa == 31080 ~ "Charlotte",
msa == 16980 ~ "Chicago",
msa == 19100 ~ "Dallas",
msa == 26420 ~ "Houston",
msa == 47900 ~ "DC",
msa == 33100 ~ "Miami",
msa == 37980 ~ "Philadelphia",
msa == 12060 ~ "Atlanta",
msa == 38060 ~ "Phoenix",
msa == 14460 ~ "Boston",
msa == 41860 ~ "San Francisco",
msa == 40140 ~ "Riverside",
msa == 19820 ~ "Detroit",
msa == 42660 ~ "Seattle")) %>%
mutate(educ = eeduc) %>%
mutate(educ = frcode(educ == 1 | educ == 2 | educ == 3 ~ "HS\nor\nLess",
educ == 4 | educ == 5 ~ "Some\nColl.",
educ == 6 ~ "4 Yr.",
educ == 7 ~ "Grad\nSchool")) %>%
mutate(att = support3) %>%
mutate(att = case_when(att == 1 ~ 1,
att <= 4 ~ 0)) %>%
group_by(msa, educ) %>%
mean_ci(att, wt = pweight, ci = .84) %>% filter(msa != "NA")
gg %>%
mutate(lab = round(mean, 2)) %>%
ggplot(., aes(x = educ, y = mean, fill = msa, group = msa)) +
geom_col(color = "black") +
facet_wrap(~ msa) +
theme_rb() +
error_bar() +
y_pct() +
scale_fill_manual(values = c(moma.colors("Warhol", 15))) +
lab_bar(top = FALSE, pos = .065, sz = 4, type = lab) +
geom_text(aes(y = .065, label = ifelse(msa == "NYC" | msa == "DC" | msa == "Boston", paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font", color = "white") +
geom_text(aes(y = .065, label = ifelse(msa == "San Francisco" | msa == "Riverside" | msa == "Seattle", paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font", color = "white") +
labs(x = "Highest Level of Education", y = "", title = "Share Never Attending Religious Services", caption = "@ryanburge\nData: Household Pulse Survey, 2024")
save("msa_ed_never_att.png", wd = 6, ht = 8)
gg <- hp %>%
filter(income > 0) %>%
mutate(att = support3) %>%
mutate(att = case_when(att == 1 ~ 1,
att <= 4 ~ 0)) %>%
mutate(msa = est_msa) %>%
mutate(msa = frcode(msa == 35620 ~ "NYC",
msa == 31080 ~ "Charlotte",
msa == 16980 ~ "Chicago",
msa == 19100 ~ "Dallas",
msa == 26420 ~ "Houston",
msa == 47900 ~ "DC",
msa == 33100 ~ "Miami",
msa == 37980 ~ "Philadelphia",
msa == 12060 ~ "Atlanta",
msa == 38060 ~ "Phoenix",
msa == 14460 ~ "Boston",
msa == 41860 ~ "San Francisco",
msa == 40140 ~ "Riverside",
msa == 19820 ~ "Detroit",
msa == 42660 ~ "Seattle")) %>%
group_by(msa, income) %>%
mean_ci(att, wt = pweight, ci = .84) %>% filter(msa != "NA")
gg %>%
mutate(lab = round(mean, 2)) %>%
ggplot(., aes(x = income, y = mean, color = msa, group = msa)) +
error_bar() +
geom_point(stroke = .75, shape = 21, fill = "white") +
geom_smooth(se = FALSE, method = lm, aes(color = msa)) +
facet_wrap(~ msa) +
theme_rb() +
y_pct() +
scale_x_continuous(breaks = c(2,4,6,8), labels = c("$25K", "$50K", "$100K", "$200K")) +
scale_color_manual(values = c(moma.colors("Warhol", 15))) +
theme(axis.text.x = element_text(size = 8)) +
# lab_bar(top = FALSE, pos = .065, sz = 4, type = lab) +
# geom_text(aes(y = .065, label = ifelse(msa == "NYC" | msa == "DC" | msa == "Boston", paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font", color = "white") +
# geom_text(aes(y = .065, label = ifelse(msa == "San Francisco" | msa == "Riverside" | msa == "Seattle", paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font", color = "white") +
labs(x = "Household Income", y = "", title = "Share Never Attending Religious Services", caption = "@ryanburge\nData: Household Pulse Survey, 2024")
save("msa_income_never_att.png", wd = 6, ht = 8)
gg <- hp %>%
mutate(hisp = rhispanic) %>%
mutate(hisp = case_when(hisp == 2 ~ 1,
hisp == 1 ~ 0)) %>%
mutate(race = rrace) %>%
mutate(race = frcode(race == 1 & hisp == 0 ~ "White",
race == 2 ~ "Black",
hisp == 1 ~ "Hisp.",
race == 3 ~ "Asian")) %>%
mutate(att = support3) %>%
mutate(att = case_when(att == 1 ~ 1,
att <= 4 ~ 0)) %>%
mutate(msa = est_msa) %>%
mutate(msa = frcode(msa == 35620 ~ "NYC",
msa == 31080 ~ "Charlotte",
msa == 16980 ~ "Chicago",
msa == 19100 ~ "Dallas",
msa == 26420 ~ "Houston",
msa == 47900 ~ "DC",
msa == 33100 ~ "Miami",
msa == 37980 ~ "Philadelphia",
msa == 12060 ~ "Atlanta",
msa == 38060 ~ "Phoenix",
msa == 14460 ~ "Boston",
msa == 41860 ~ "San Francisco",
msa == 40140 ~ "Riverside",
msa == 19820 ~ "Detroit",
msa == 42660 ~ "Seattle")) %>%
group_by(msa, race) %>%
mean_ci(att, wt = pweight, ci = .84) %>%
filter(msa != "NA") %>%
filter(race != "NA")
gg %>%
mutate(lab = round(mean, 2)) %>%
ggplot(., aes(x = race, y = mean, fill = msa, group = msa)) +
geom_col(color = "black") +
facet_wrap(~ msa) +
theme_rb() +
error_bar() +
y_pct() +
scale_fill_manual(values = c(moma.colors("Warhol", 15))) +
lab_bar(top = FALSE, pos = .065, sz = 4, type = lab) +
theme(axis.text.x = element_text(size = 8)) +
geom_text(aes(y = .065, label = ifelse(msa == "NYC" | msa == "DC" | msa == "Boston", paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font", color = "white") +
geom_text(aes(y = .065, label = ifelse(msa == "San Francisco" | msa == "Riverside" | msa == "Seattle", paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font", color = "white") +
labs(x = "Race of Respondent", y = "", title = "Share Never Attending Religious Services", caption = "@ryanburge\nData: Household Pulse Survey, 2024")
save("msa_race_never_att.png", wd = 6, ht = 8)
gg <- hp %>%
mutate(birthyr = tbirth_year) %>%
mutate(age = 2024 - birthyr) %>%
mutate(age2 = frcode(age <= 35 ~ "18-35",
age >= 36 & age <= 49 ~ "36-49",
age >= 50 & age <= 64 ~ "50-64",
age >= 65 ~ "65+")) %>%
mutate(att = support3) %>%
mutate(att = case_when(att == 1 ~ 1,
att <= 4 ~ 0)) %>%
mutate(msa = est_msa) %>%
mutate(msa = frcode(msa == 35620 ~ "NYC",
msa == 31080 ~ "Charlotte",
msa == 16980 ~ "Chicago",
msa == 19100 ~ "Dallas",
msa == 26420 ~ "Houston",
msa == 47900 ~ "DC",
msa == 33100 ~ "Miami",
msa == 37980 ~ "Philadelphia",
msa == 12060 ~ "Atlanta",
msa == 38060 ~ "Phoenix",
msa == 14460 ~ "Boston",
msa == 41860 ~ "San Francisco",
msa == 40140 ~ "Riverside",
msa == 19820 ~ "Detroit",
msa == 42660 ~ "Seattle")) %>%
group_by(msa, age2) %>%
mean_ci(att, wt = pweight, ci = .84) %>%
filter(msa != "NA") %>%
filter(age2 != "NA")
gg %>%
mutate(lab = round(mean, 2)) %>%
ggplot(., aes(x = age2, y = mean, fill = msa, group = msa)) +
geom_col(color = "black") +
facet_wrap(~ msa) +
theme_rb() +
error_bar() +
y_pct() +
scale_fill_manual(values = c(moma.colors("Warhol", 15))) +
lab_bar(top = FALSE, pos = .065, sz = 4, type = lab) +
theme(axis.text.x = element_text(size = 8)) +
geom_text(aes(y = .065, label = ifelse(msa == "NYC" | msa == "DC" | msa == "Boston", paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font", color = "white") +
geom_text(aes(y = .065, label = ifelse(msa == "San Francisco" | msa == "Riverside" | msa == "Seattle", paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font", color = "white") +
labs(x = "Age of Respondent", y = "", title = "Share Never Attending Religious Services", caption = "@ryanburge\nData: Household Pulse Survey, 2024")
save("msa_age2_never_att.png", wd = 6, ht = 8)
regg <- hp %>%
filter(income > 0) %>%
mutate(att = support3) %>%
mutate(att = case_when(att == 1 ~ 1,
att <= 4 ~ 0)) %>%
mutate(msa = est_msa) %>%
mutate(msa = frcode(msa == 35620 ~ "NYC",
msa == 31080 ~ "Charlotte",
msa == 16980 ~ "Chicago",
msa == 19100 ~ "Dallas",
msa == 26420 ~ "Houston",
msa == 47900 ~ "DC",
msa == 33100 ~ "Miami",
msa == 37980 ~ "Philadelphia",
msa == 12060 ~ "Atlanta",
msa == 38060 ~ "Phoenix",
msa == 14460 ~ "Boston",
msa == 41860 ~ "San Francisco",
msa == 40140 ~ "Riverside",
msa == 19820 ~ "Detroit",
msa == 42660 ~ "Seattle")) %>%
mutate(birthyr = tbirth_year) %>%
mutate(age = 2024 - birthyr) %>%
mutate(hisp = rhispanic) %>%
mutate(hisp = case_when(hisp == 2 ~ 1,
hisp == 1 ~ 0)) %>%
mutate(race = rrace) %>%
mutate(white = case_when(race == 1 & hisp == 0 ~ 1,
TRUE ~ 0)) %>%
mutate(educ = eeduc) %>%
select(att, msa, age, white, income, educ)
list_of_msa_datasets <- split(regg, regg$msa)
names(list_of_msa_datasets) <- gsub(" ", "_", names(list_of_msa_datasets))
for(msa in names(list_of_msa_datasets)) {
assign(msa, list_of_msa_datasets[[msa]])
}
library(jtools)
coef_names <- c("Education" = "educ",
"Age" = "age",
"Income" = "income",
"White" = "white")
library(patchwork)
plots_list <- list()
# Loop through each MSA dataset
for (msa_name in names(list_of_msa_datasets)) {
# Fit the model
model <- glm(att ~ age + white + income + educ, family = "binomial", data = list_of_msa_datasets[[msa_name]])
# Create the plot with plot_summs
plot <- plot_summs(model, scale = TRUE, robust = "HC3", coefs = coef_names) +
theme_rb() +
labs(x = "", y = "", title = msa_name) # Use MSA name as the title
# Add the plot to the list
plots_list[[msa_name]] <- plot
}
# Combine all plots using patchwork
combined_plot <- reduce(plots_list, `+`)
# To arrange the plots in a more structured manner, you can specify the layout
# For example, to arrange in a 4x4 grid (adjust based on the number of MSAs)
combined_plot + plot_layout(ncol = 4, nrow = 4)
# Display the combined plot
print(combined_plot)
combined_plot <- combined_plot +
plot_annotation(title = "What Factors Increase the Likelihood of Never Attending?",
theme = theme(plot.title = element_text(hjust = 0.5)))
ggsave("E://graphs23/patched.png", combined_plot, width = 10, height = 10)
regg <- hp %>%
filter(income > 0) %>%
mutate(att = support3) %>%
mutate(att = case_when(att == 1 ~ 1,
att <= 4 ~ 0)) %>%
mutate(msa = est_msa) %>%
mutate(msa = frcode(msa == 35620 ~ "NYC",
msa == 31080 ~ "Charlotte",
msa == 16980 ~ "Chicago",
msa == 19100 ~ "Dallas",
msa == 26420 ~ "Houston",
msa == 47900 ~ "DC",
msa == 33100 ~ "Miami",
msa == 37980 ~ "Philadelphia",
msa == 12060 ~ "Atlanta",
msa == 38060 ~ "Phoenix",
msa == 14460 ~ "Boston",
msa == 41860 ~ "San Francisco",
msa == 40140 ~ "Riverside",
msa == 19820 ~ "Detroit",
msa == 42660 ~ "Seattle")) %>%
mutate(birthyr = tbirth_year) %>%
mutate(age = 2024 - birthyr) %>%
mutate(hisp = rhispanic) %>%
mutate(hisp = case_when(hisp == 2 ~ 1,
hisp == 1 ~ 0)) %>%
mutate(race = rrace) %>%
mutate(race = frcode(race == 1 & hisp == 0 ~ "White",
race == 2 ~ "Black",
hisp == 1 ~ "Hisp.",
race == 3 ~ "Asian")) %>%
mutate(educ = eeduc) %>%
select(att, msa, age, race, income, educ) %>%
filter(age <= 75) %>%
filter(msa != "NA")
int <- glm(att ~ race*age + income + educ, family = 'binomial', data = regg)
library(interactions)
graph <- interact_plot(int, pred = age, modx = race, interval = TRUE, int.width = .76)
graph +
theme_rb(legend = TRUE) +
y_pct() +
labs(x = "Age", y = "Estimate of Share Never Attending", title = "The Interaction of Race and Age on Never Attending",
subtitle = "In Major Metro Areas - Controls for Education and Income",
caption = "@ryanburge\nData: Household Pulse Survey, 2024")
ggsave("E://graphs23/never_attend_pulse_interact.png", bg = "white", width = 7, height = 6)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment