ct_selection_untitled-3665FA19.R

fips <- data.frame(
  stringsAsFactors = FALSE,
  id = 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"),
  abb = 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")
) %>% 
  mutate(FIPS = as.numeric(FIPS))


aaa1 <- cces %>% 
  mutate(ath = case_when(religion == 9 ~ 1, 
                         TRUE ~ 0)) %>% 
  filter(year >= 2016) %>% 
  group_by(state) %>% 
  mean_ci(ath, wt = weight, ci = .84) %>% 
  mutate(type = "CES") %>% 
  rename(id = state)

aaa1 <- left_join(aaa1, fips)

aaa2 <- ns %>% 
  mutate(ath = case_when(religion == 10 ~ 1, 
                         TRUE ~ 0)) %>% 
  group_by(state) %>% 
  mean_ci(ath, wt = weight, ci = .84) %>% 
  mutate(type = "NS") %>% 
  rename(abb = state)

aaa2 <- left_join(aaa2, fips)


graph <- bind_rows(aaa1, aaa2)

graph %>% 
  filter(id != "NA") %>% 
  filter(id != "District of Columbia") %>% 
  ggplot(., aes(x = fct_rev(id), y = mean, color = type, group = type)) +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=0) +
  geom_point(fill = "white", stroke = 1, shape = 21) + 
  coord_flip() +
  scale_color_calc() +
  scale_y_continuous(labels = percent) +
  theme_rb(legend = TRUE)  +
  guides(color = guide_legend(reverse = TRUE)) +
  labs(x = "", y = "", title = "State Level Estimate of Atheists in Two Surveys", caption = "@ryanburge\nData: Cooperative Election Study, 2016-2022 + Nationscape 2019-2021")
save("atheist_compare_points.png", ht = 9, wd = 6)

library(geofacet)


geo <- graph %>% 
  filter(id != "NA") %>% 
  filter(id != "District of Columbia") 

geo %>% 
  mutate(lab = round(mean, 2)) %>% 
  ggplot(., aes(x = fct_rev(type), y = mean, fill= type)) +
  geom_col(color = "black") +
  facet_geo(~ id) +
  theme_rb() +
  scale_fill_calc() +
  scale_y_continuous(labels = percent) +
  theme(strip.text = element_text(size = 8)) +
  geom_text(aes(y = mean + .015, label = paste0(lab*100, '%')), position = position_dodge(width = .9), size = 5, family = "font") +
  labs(x = "", y =  "", title = "Share of the State That Identifies as Atheist", caption = "@ryanburge\nData: Cooperative Election Study, 2016-2022 + Nationscape 2019-2021")
save("geo_facet_atheist.png", ht = 10)



one <- aaa1 %>% 
  select(id, mean_ces = mean, n_ces = n)


two <- aaa2 %>% 
  select(id, mean_ns = mean, n_ns = n)

both <- left_join(one, two) %>% 
  mutate(pct_diff = mean_ces - mean_ns) %>% 
  mutate(mid_pct = (mean_ces + mean_ns)/2)

both <- both %>% 
  select(id, pct = mid_pct)

both$id <- tolower(both$id)


library(fiftystater)

data("fifty_states")

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

map <- mapp %>% 
  mutate(bins = frcode(pct >= .06 ~ "6%+",
                       pct < .06  & pct >= .04 ~ "4%-6%", 
                       pct < .04 ~ "<4%")) %>% 
  mutate(bins = fct_rev(bins)) %>% 
  filter(bins != "NA")


centers <- read_csv("D://state_centers_inset.csv") %>% rename(id = full)

centers$id <- tolower(centers$id)


centers <- left_join(map, centers, by = c("id")) %>% as_tibble()

centers <- centers %>% distinct(id, .keep_all = TRUE)

centers <- centers %>% 
  mutate(mean = round(pct, 3)) %>% 
  mutate(lab = mean*100) %>% 
  mutate(lab = paste0(lab, '%')) %>% 
  mutate(label = paste(state, lab, sep = ": "))

centers <- left_join(centers, mapp)


library(ggrepel)

mapp <- left_join(map, mapp)


mapp %>% 
  ggplot(., mapping = aes(x = long, y = lat, group = group, fill = bins)) +
  geom_polygon(color = "black", linewidth = 0.1) +
  coord_map(projection = "albers", lat0 = 41, lat1 = 44) +
  scale_fill_manual(values = c("#c91105", "#caf0f8",  "#03045e")) +
  ggrepel::geom_label_repel(data = centers,
                             aes(x = longitude, y = latitude, label = label, group = group), size = 2.5, fill = "white", 
                             seed = 1002, family = "font", force = .001, show.legend = FALSE) + 
  theme_rb(legend = TRUE) +
  theme(axis.line = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid = element_blank(),
        panel.spacing = unit(0, "lines"),
        plot.background = element_blank(),
        legend.position = "bottom",
        plot.title = element_text(size = 18)) +
  labs(title = "Best Estimate of Adult Population That are Atheists, State-Level",  caption = "@ryanburge\nData: Cooperative Election Study, 2016-2022 + Nationscape 2019-2021") 
save("atheist_state_map_estimate.png")



one <- aaa1 %>% 
  select(id, mean_ces = mean, n_ces = n)


two <- aaa2 %>% 
  select(id, mean_ns = mean, n_ns = n)

both <- left_join(one, two) %>% 
  mutate(pct_diff = mean_ces - mean_ns) %>% 
  mutate(mid_pct = (mean_ces + mean_ns)/2) %>% 
  select(id, pct = mid_pct, n = n_ns) %>% 
  na.omit()

age <- read_csv("E://data/state_age.csv")

graph <- left_join(both, age)

graph <- 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() %>% as_tibble()

graph <- left_join(graph, cc)

graph %>% 
  ggplot(., aes(x = age, y = pct)) +
  geom_smooth(se = FALSE, method = lm, linetype = "twodash", color = "black", linewidth = .5) +
  geom_point(data = graph, aes(size  = n, color = region), stroke = 1.25, shape = 21, fill = "white") +
  geom_text_repel(data = graph, aes(x = age, y = pct, label = abb), family = "font") +
  theme_rb() +
  scale_color_manual(values = c("#ff595e", "#8ac926", '#1982c4', "#6a4c93")) +
  scale_y_continuous(labels = percent) +
  labs(x = "Median Age", y = "Share Who are Atheists", title = "The Relationship Between Age and Atheists %, State Level", caption = "@ryanburge\nData: Cooperative Election Study, 2016-2022 + Nationscape 2019-2021 + US Census Bureau 2019")
save("ath_age_state.png")


lm(pct ~ age, data = graph) ## Five years older = 1% more atheist 

library(googlesheets4)

trump <- read_sheet("https://docs.google.com/spreadsheets/d/1B0ImV7kqEh6NjM1AaEbz822agJcPQUaNfTR7okwUuns/edit?usp=sharing")

graph <- left_join(graph, trump)

graph <- graph %>% 
  mutate(tbins = frcode(trump  <  .40 ~ "a", 
                        trump >= .40 & trump < .45 ~ "b", 
                        trump >= .45 & trump < .50 ~ "c", 
                        trump >= .50 & trump < .55 ~ "d", 
                        trump >= .55 & trump < .60 ~ "e", 
                        trump >= .60 ~ "f"))

graph %>% 
  ggplot(., aes(x = trump, y = pct)) +
  geom_smooth(se = FALSE, method = lm, linetype = "twodash", color = "black", linewidth = .5) +
  geom_point(data = graph, aes(size  = n, color = tbins), stroke = 1.25, shape = 21, fill = "white") +
  geom_text_repel(data = graph, aes(x = trump, y = pct, label = abb), family = "font") +
  theme_rb() +
  scale_color_manual(values = c("#0000ff", "#3300cc", "#660099", "#990066", "#cc0033", "#ff0000")) +
  scale_y_continuous(labels = percent) +
  scale_x_continuous(labels = percent) +
  labs(x = "Trump's Vote Share in 2020", y = "Share Who are Atheists", title = "The Relationship Between Age and Atheists %, State Level", caption = "@ryanburge\nData: Cooperative Election Study, 2016-2022 + Nationscape 2019-2021 + American Presidency Project")
save("ath_trump_state.png")

lm(pct ~ trump, data = graph) ## 10 pt increase for Trump = .9% lower atheist rate 

white <- read_sheet("https://docs.google.com/spreadsheets/d/18Yps7aLor4VUS3i-ERfIStGuH7iKvUVtbUgJTUKFnzQ/edit?usp=sharing")

graph <- left_join(graph, white)

reg <- lm(pct ~ trump + age + white + bach + income + male, data = graph)

library(jtools)



coef_names <- c("Trump %" = "trump",
                "Education" = "bach", 
                "Age" = "age", 
                "Income" = "income", 
                "White" = "white", 
                "Male" = "male")

out <- plot_summs(reg, scale = TRUE, coefs = coef_names)

out + 
  theme_rb() + 
  labs(x = "", y = "", title = "What Factors Predict More Atheists at the State Level?", caption = "@ryanburge") +
  add_text(x = .0075, y = 5.5, word = "More\nAtheists", sz = 6) +
  add_text(x =-.013, y = 1.5, word = "Fewer\nAtheists", sz = 6) 
save("reg_atheist_state.png", ht = 4)