ct_selection_untitled-7D9EA7BA.R

graph1 <- gss %>% 
  filter(year != 1982) %>% 
  filter(year != 1987) %>% 
  filter(year <= 2018) %>% 
  mutate(nd = case_when(denom == 70 ~ 1, 
                        TRUE ~ 0)) %>% 
  group_by(year) %>% 
  mean_ci(nd, wt = wtssall) %>% 
  mutate(type = "Non-Denom.")

graph2 <- gss %>% 
  filter(year == 1982 | year == 1987) %>% 
  mutate(nd = case_when(denom == 70 ~ 1, 
                        TRUE ~ 0)) %>% 
  group_by(year) %>% 
  mean_ci(nd, wt = oversamp) %>% 
  mutate(type = "Non-Denom.")

graph3 <- gss %>% 
  filter(year == 2021) %>% 
  mutate(nd = case_when(denom == 70 ~ 1, 
                        TRUE ~ 0)) %>% 
  group_by(year) %>% 
  mean_ci(nd, wt = wtssnrps) %>% 
  mutate(type = "Non-Denom.")

all <- bind_rows(graph1, graph2, graph3)

all %>% 
  ggplot(., aes(x = year, y = mean)) +
  geom_line() +
  geom_point(stroke = 1, shape = 21, fill = "white") +
  theme_rb() +
  scale_y_continuous(labels = percent, limits = c(0, .13)) + 
  geom_smooth(se = FALSE, linetype = "twodash", color = "red", linewidth = 1) +
  labs(x = "Year", y = "", title = "Share of All Americans Who Identify as Non-Denominational", caption = "@ryanburge\nData: General Social Survey, 1972-2021")
save("non_denom_gss.png")


graph1 <- gss %>% 
  filter(year != 1982) %>% 
  filter(year != 1987) %>% 
  filter(year <= 2018) %>% 
  filter(relig == 1) %>% 
  mutate(den = frcode(denom == 14 ~ "SBC", 
                      denom == 22 ~ "UMC", 
                      denom == 70 ~ "ND")) %>% 
  group_by(year) %>% 
  ct(den, wt = wtssall) %>% 
  filter(den != "NA")

graph2 <- gss %>% 
  filter(year == 1982 | year == 1987) %>% 
  filter(relig == 1) %>% 
  mutate(den = frcode(denom == 14 ~ "SBC", 
                      denom == 22 ~ "UMC", 
                      denom == 70 ~ "ND")) %>% 
  group_by(year) %>% 
  ct(den, wt = oversamp) %>% 
  filter(den != "NA")

graph3 <- gss %>% 
  filter(year == 2021) %>% 
  filter(relig == 1) %>% 
  mutate(den = frcode(denom == 14 ~ "SBC", 
                      denom == 22 ~ "UMC", 
                      denom == 70 ~ "ND")) %>% 
  group_by(year) %>% 
  ct(den, wt = wtssnrps) %>% 
  filter(den != "NA")

all <- bind_rows(graph1, graph2, graph3)

library(ggtext)

all %>% 
  filter(year >= 1984) %>% 
  filter(year <= 2018) %>% 
  ggplot(., aes(x = year, y = pct, color = den, group = den)) +
  # geom_line() +
  # geom_point(stroke = 1, shape = 21, fill = "white", alpha = 1) +
  theme_rb() +
  scale_y_continuous(labels = percent, limits = c(0, .225)) + 
  geom_labelsmooth(aes(label = den),  method = "loess", formula = y ~ x, family = "font", linewidth = 1, text_smoothing = 30, size = 6, linewidth = 1, boxlinewidth = 0.3, hjust = .45) +
  labs(
    x = "Year",
    y = "",
    title = "<span style = 'font-size:14pt'>Share of Protestants Who are </span><span style = 'color:#0099B4FF; font-size:18pt'>Southern Baptists,</span><span style = 'color:#AD002AFF; font-size:18pt'>United Methodists,</span>or <span style = 'color:#42B540FF; font-size:18pt'>Non-Denominational</span>", 
    caption = "@ryanburge\nData: General Social Survey, 1984-2018")  +
  theme_rb() +
  scale_color_manual(values = c("#0099B4FF", "#AD002AFF", "#42B540FF")) +
  theme(
    text = element_text(family = "font"),
    # plot.title.position = "plot",
    plot.title = element_markdown(size = 14, lineheight = 1.2))
  # labs(x = "Year", y = "", title = "Share of Protestants Who are Southern Baptist, United Methodist, or Non-Denominational", caption = "@ryanburge\nData: General Social Survey, 1984-2018")
save("non_denom_gss_prot_comp.png")


library(readxl)
library(janitor)

cen20 <- read_excel("E:/data/rel_cen20.xlsx")

cen20 <- cen20 %>% 
  clean_names()


cen20 <- cen20 %>% 
  arrange(-adherents) %>% 
  filter(group_name != "NA") %>% 
  head(20) 

cen20$labs <- formatC(cen20$adherents, format = "d", big.mark = ",")

cen20 %>% 
  ggplot(., aes(x = reorder(group_name, adherents), y = adherents, fill = group_name)) +
  geom_col(color = "black") +
  coord_flip() +
  geom_text(aes(y = adherents + 6100000, label = labs), position = position_dodge(width = .9), size = 5, family = "font") +
  theme_rb() + 
  scale_fill_manual(values=met.brewer("Hiroshige", 20)) +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6), limits = c(0, 71000000)) +
  labs(x = "", y = "", title = "The 20 Largest Religious Traditions in the United States", caption = "@ryanburge\nData: 2020 Religion Census")
save("rel_cen_20_largest.png")




cen20 <- read_excel("E:/data/rel_cen20.xlsx", sheet = 3)

cen20 <- cen20 %>% 
  clean_names()


nd <- cen20 %>% 
  filter(group_code == "500") %>% 
  mutate(fips = as.numeric(fips))

nd <- nd %>% 
  mutate(bins = adherents_as_percent_of_total_adherents) %>% 
  mutate(bins = frcode(bins >= .20 ~ "20%+",
                       bins >= .15 & bins <= .199 ~ "15%-20%",
                       bins >= .10 & bins <= .1499 ~ "10%-15%",
                       bins >= .05 & bins <= .0999 ~ "5%-10%",
                       bins <= .0499 ~ "<5%")) %>%  
  filter(bins != "NA")

library(urbnmapr)
library(urbnthemes)


territories_counties <- get_urbn_map(map = "territories_counties")

territories_counties <-  territories_counties %>% 
  mutate(fips =  as.numeric(county_fips))

joined <- left_join(territories_counties, nd) 

library(paletteer) 
paletteer_d("ggsci::lanonc_lancet")



joined %>% 
  filter(bins != 'NA') %>% 
  ggplot(aes(long, lat, group = group, fill = bins)) +
  geom_polygon(color = "black", linewidth = 0.05) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  labs(title  = "      Percent of All Adherents Who Are Non-Denominational", caption = "@ryanburge\nData: 2020 Religious Census", fill = "") +
  theme(legend.title = element_text(), legend.key.width = unit(.5, "in")) +
  urbnthemes::theme_urbn_map() +
  scale_fill_paletteer_d("PNWColors::Sunset2") +
  theme(plot.title = element_text(size=20, family= "font", face= "bold")) +
  theme(plot.caption = element_text(size = 14, family="font")) +
  theme(text=element_text(size=18, family="font")) 
save("relcen2020_map_nd_new.png")






all <- read_csv("E://data/county_pop22.csv")

age <- all %>% 
  filter(year == 4) %>% 
  mutate(fips = paste0(state, county)) %>% 
  select(fips, median_age_tot, pop = popestimate, stname)

library(readxl)
library(janitor)

cen <- read_excel("E://data/rel_cen20.xlsx", sheet = 3) %>% 
  clean_names()

rel <- cen %>% 
  filter(group_code == "500") %>% 
  select(fips, pct = adherents_as_percent_of_total_adherents)

graph <- left_join(rel, age)


fips <- data.frame(
  stringsAsFactors = FALSE,
  stname = c("Alabama","Alaska","Arizona",
             "Arkansas","California","Colorado","Connecticut",
             "Delaware","Florida","Georgia","Hawaii","Idaho",
             "Illinois","Indiana","Iowa","Kansas","Kentucky","Louisiana",
             "Maine","Maryland","Massachusetts","Michigan",
             "Minnesota","Mississippi","Missouri","Montana","Nebraska",
             "Nevada","New Hampshire","New Jersey","New Mexico",
             "New York","North Carolina","North Dakota","Ohio",
             "Oklahoma","Oregon","Pennsylvania","Rhode Island",
             "South Carolina","South Dakota","Tennessee","Texas","Utah",
             "Vermont","Virginia","Washington","West Virginia",
             "Wisconsin","Wyoming","American Samoa","Guam",
             "Northern Mariana Islands","Puerto Rico","Virgin Islands"),
  Postal.Code = c("AL","AK","AZ","AR","CA",
                  "CO","CT","DE","FL","GA","HI","ID","IL","IN","IA",
                  "KS","KY","LA","ME","MD","MA","MI","MN","MS",
                  "MO","MT","NE","NV","NH","NJ","NM","NY","NC","ND",
                  "OH","OK","OR","PA","RI","SC","SD","TN","TX",
                  "UT","VT","VA","WA","WV","WI","WY","AS","GU","MP",
                  "PR","VI"),
  FIPS = c("01","02","04","05","06",
           "08","09","10","12","13","15","16","17","18","19",
           "20","21","22","23","24","25","26","27","28",
           "29","30","31","32","33","34","35","36","37","38",
           "39","40","41","42","44","45","46","47","48",
           "49","50","51","53","54","55","56","60","66","69",
           "72","78")
)


graph1 <- left_join(graph, fips) %>% 
  mutate(FIPS = as.numeric(FIPS))


cc <- cces20 %>% 
  select(FIPS = inputstate, region) %>% 
  mutate(region = frcode(region == 1 ~ "Northeast", 
                         region == 2 ~ "Midwest", 
                         region == 3 ~ "South", 
                         region == 4 ~ "West")) %>% 
  distinct() %>% 
  mutate(FIPS = as.numeric(FIPS)) %>% 
  as_tibble()

graph2 <- left_join(graph1, cc) %>% filter(region != "NA")



graph2 %>% 
  ggplot(., aes(x = median_age_tot, y = pct)) +
  geom_point(data = graph2, aes(size = pop, color = region), stroke = .5, shape = 21) + 
  geom_smooth(se = FALSE, method = lm, color = "black", linetype = "twodash", linewidth = .5) +
  guides(size = FALSE) +
  y_pct() + 
  theme_rb() +
  scale_color_tableau() +
  theme(legend.position = c(.85, .55)) +
  labs(x = "Median Age", y = "Share of All Adherents Who Are Non-Denominational", title = "Are Non-Denominationals More Prominent in Younger Counties?", caption = "@ryanburge\nData: Religion Census 2020 + US Census Data")
save("med_age_scatter_nd.png")



cen20 <- read_excel("E:/data/rel_cen20.xlsx", sheet = 2)

cen20 <- cen20 %>% 
  clean_names()

state <- cen20 %>% 
  group_by(state_name) %>% 
  slice_max(order_by = adherents, n = 1) %>% 
  select(state_name, group_name)



library(fiftystater)

data("fifty_states")

state$id <- tolower(state$state_name)


mapp <- left_join(fifty_states, state, by = c("id")) %>% as_tibble() 

mapp <- mapp %>% 
  mutate(group_name = str_replace(group_name, "Church of Jesus Christ of Latter-day Saints", "LDS")) %>% 
  mutate(group_name = str_replace(group_name, "Non-denominational Christian Churches", "Non-Denominational")) 


mapp %>% 
  ggplot(., aes(map_id = id)) +
  geom_map(aes(fill = group_name), map = fifty_states, color = "black") + 
  expand_limits(x = fifty_states$long, y = fifty_states$lat) +
  coord_map() +
  scale_x_continuous(breaks = NULL) + 
  scale_y_continuous(breaks = NULL) +
  labs(x = "", y = "") +
  fifty_states_inset_boxes() +
  scale_fill_manual(values = c("Catholic Church" =  "#00468BFF", 
                               "LDS" = "#ED0000FF", 
                               "Non-Denominational" = "#42B540FF", 
                               "Southern Baptist Convention" = "#0099B4FF")) +
  theme_rb(legend = TRUE) +
  theme(legend.position = "bottom", 
        legend.text = element_text(size = 16), 
        panel.background = element_blank()) +
  labs(title = "The Denomination with the Most Adherents in Each State", caption = "@ryanburge\nData: Religion Census 2020") 
save("fifty_state_maps_denom2020.png")



top <- cen20 %>% 
  group_by(state_name) %>% 
  top_n(adherents, n = 2) %>% 
  ungroup(state_name)

graph <- top %>% 
  group_by(state_name) %>% 
  top_n(-1) %>% 
  select(state_name, group_name)


graph$id <- tolower(graph$state_name)


library(fiftystater)

data("fifty_states")

mapp <- left_join(fifty_states, graph, by = c("id")) %>% as_tibble() %>% filter(group_name != "NA")

mapp <- mapp %>% 
  mutate(group_name = str_replace(group_name, "Church of Jesus Christ of Latter-day Saints", "LDS")) %>% 
  mutate(group_name = str_replace(group_name, "Non-denominational Christian Churches", "Non-Denominational")) %>% 
  mutate(group_name = str_replace(group_name, "Evangelical Lutheran Church in America", "ELCA")) %>% 
  mutate(group_name = str_replace(group_name, "Lutheran Church--Missouri Synod", "Lutheran - MO Synod"))



mapp %>% 
  ggplot(., aes(map_id = id)) +
  geom_map(aes(fill = group_name), map = fifty_states, color = "black") + 
  expand_limits(x = fifty_states$long, y = fifty_states$lat) +
  coord_map() +
  scale_x_continuous(breaks = NULL) + 
  scale_y_continuous(breaks = NULL) +
  labs(x = "", y = "") +
  theme(legend.position = "bottom", 
        panel.background = element_blank()) +
  fifty_states_inset_boxes() +
  theme_rb(legend = TRUE) +
  theme(legend.text = element_text(size = 10)) +
  scale_fill_manual(values = c("Catholic Church" =  "#00468BFF", 
                               "LDS" = "#ED0000FF", 
                               "Non-Denominational" = "#42B540FF", 
                               "Southern Baptist Convention" = "#0099B4FF", 
                               "Muslim Estimate" = "#925E9FFF", 
                               "ELCA" = "#FDAF91FF", 
                               "United Methodist Church" = "#AD002AFF", 
                               "Lutheran - MO Synod" = "#ADB6B6FF")) +
  labs(title = "Religious Group with the Second Most Adherents by State", caption = "@ryanburge\nData: Religion Census, 2020") 
save("second_most_pop2020.png")