ct_selection_untitled-CEE13118.R

library(rio)
library(labelled)

fem <- import("E://data/fg_female.dta") %>% mutate(gender = "Female")
men <- import("E://data/fg_male.dta") %>% mutate(gender = "Male")

fg <- bind_rows(fem, men) %>% rename(weight = WGT2015_2017)



combined_df <- bind_rows(
  fg %>% 
    mutate(vv = pill) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Pill"),
  
  fg %>% 
    mutate(vv = condom) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Condom"),
  
  fg %>% 
    mutate(vv = vasect) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Vasectomy"),
  
  fg %>% 
    mutate(vv = depoprov) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Depo Shot"),
  
  fg %>% 
    mutate(vv = widrawal) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Withdrawal"),
  
  fg %>% 
    mutate(vv = rhythm) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Rhythm Method"),
  
  fg %>% 
    mutate(vv = patch) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Patch"),
  
  fg %>% 
    mutate(vv = ring) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Ring"),
  
  fg %>% 
    mutate(vv = mornpill) %>% 
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "Plan B"),
  
  fg %>% 
    mutate(vv = OTHRMETH01) %>% 
    mutate(vv = case_when(vv == 19 ~ 1, 
                          vv <= 98 ~ 0)) %>% 
    mean_ci(vv, wt = weight) %>% 
    mutate(type = "IUD")
)

#####

color_mapping <- c(
  "Pill" = "#28666e",
  "Condom" = "#5f2680",
  "Vasectomy" = "#B5B682",
  "Depo Shot" = "#7D3C98",
  "Withdrawal" = "#033f63",
  "Rhythm Method" = "#FEDC97",
  "Patch" = "#5E4A78",
  "Ring" = "#3E6B7A",
  "Plan B" = "#A4A37A",
  "IUD" = "#8C5A9E"
)


combined_df %>% 
  mutate(lab = round(mean, 2)) %>% 
  ggplot(., aes(x = reorder(type, mean), y = mean, fill = type)) +
  geom_col(color = "black") +
  coord_flip() +
  theme_rb()+
  y_pct() + 
  scale_fill_manual(values = color_mapping) +
  theme(plot.title = element_text(size = 18)) +
  geom_text(aes(y = mean + .045, label = ifelse(mean >= .095, paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 8, family = "font") +
  geom_text(aes(y = mean + .03, label = ifelse(mean < .094, paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 8, family = "font") +
  labs(x = "", y = "", title = "Have You Ever Used The Following Birth Control Methods", caption = '@ryanburge + @religiondata\nData: National Survey of Family Growth, 2015-2017')
save("fam_growth_all_bc.png")



# Define the reltrad labels
reltrad_labels <- c(
  "1" = "Evangelical",
  "2" = "Mainline",
  "3" = "Black Prot.",
  "4" = "Catholic",
  "5" = "Other religion",
  "6" = "No Religion",
  "8" = "Refused",
  "9" = "Don't know"
)

combined_df <- bind_rows(
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = pill) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Pill"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = condom) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Condom"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = vasect) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Vasectomy"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = depoprov) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Depo Shot"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = widrawal) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Withdrawal"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = rhythm) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Rhythm Method"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = patch) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Patch"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = ring) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Ring"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = mornpill) %>%
    mutate(vv = case_when(vv == 1 ~ 1, 
                          vv == 5 | vv == 8 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "Plan B"),
  
  fg %>%
    group_by(reltrad = dplyr::recode(reltrad, !!!reltrad_labels)) %>%
    mutate(vv = OTHRMETH01) %>%
    mutate(vv = case_when(vv == 19 ~ 1, 
                          vv <= 98 ~ 0)) %>%
    mean_ci(vv, wt = weight) %>%
    mutate(type = "IUD")
)

all <- combined_df %>% 
  filter(reltrad != "Don't know") %>% 
  filter(reltrad != "Refused")

all %>% 
  mutate(lab = round(mean, 2)) %>% 
  ggplot(., aes(x = reorder(type, mean), y = mean, fill = type)) +
  geom_col(color = "black") +
  coord_flip() +
  theme_rb()+
  facet_wrap(~ reltrad) +
  y_pct() + 
  scale_fill_manual(values = color_mapping) +
  theme(plot.title = element_text(size = 18)) +
  geom_text(aes(y = mean + .07, label = ifelse(mean >= 0, paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 4, family = "font") +
  theme(strip.text = element_text(size = 20))+
  labs(x = "", y = "", title = "Have You Ever Used The Following Birth Control Methods", caption = '@ryanburge + @religiondata\nData: National Survey of Family Growth, 2015-2017')
save("fam_growth_all_bc_reltrad.png")


graph <- fg %>% 
  filter(oppyearnum <= 20) %>% 
  group_by(ager, gender) %>% 
  mean_ci(oppyearnum, wt = weight)


graph %>% 
  ggplot(., aes(x = ager, y = mean, color = gender, group = gender)) +
  geom_point(stroke = 1, shape = 21) + 
  geom_labelsmooth(aes(label = gender),  method = "loess", formula = y ~ x, family = "font", linewidth = 1, text_smoothing = 30, size = 6, linewidth = 1, boxlinewidth = 0.3, hjust = .75) +
  scale_color_manual(values = c("#B5B682", "#033f63")) +
  theme_rb() + 
  labs(x = "", y = "", title = "Number of opposite-sex partners in last 12 months", caption = '@ryanburge + @religiondata\nData: National Survey of Family Growth, 2015-2017')
save("num_partners_gender_age.png")



reltrad_labels <- c(
  "1" = "Evangelical",
  "2" = "Mainline",
  "3" = "Black Prot.",
  "4" = "Catholic",
  "5" = "Other religion",
  "6" = "No Religion",
  "8" = "Refused",
  "9" = "Don't know"
)

gg <- fg %>% 
  mutate(marstat = frcode(marstat == 1 ~ "Married", 
                          marstat == 2 ~ "Co-habitating", 
                          marstat == 6 ~ "Never been married", 
                          marstat == 4 | marstat == 5 ~ "Divorced/Separated"),
         reltrad = factor(dplyr::recode(reltrad, !!!reltrad_labels),
                          levels = c("Evangelical", "Mainline", "Black Prot.", "Catholic", "Other religion", "No Religion"))) %>% 
  group_by(reltrad) %>%
  ct(marstat, wt = weight, show_na = FALSE) %>% 
  filter(reltrad != "Refused") %>% 
  filter(reltrad != "Don't know")


color_mapping <- c(
  "Evangelical" = "#28666e",
  "Mainline" = "#5f2680",
  "Black Prot." = "#B5B682",
  "Catholic" = "#7D3C98",
  "Other religion" = "#033f63",
  "No Religion" = "#FEDC97"
)

gg %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = reltrad, y = pct, fill = reltrad)) +
  geom_col(color = "black") +
  facet_wrap(~ marstat) +
  theme_rb() +
  scale_fill_manual(values = color_mapping) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(labels = percent) +
  theme(strip.text = element_text(size = 20)) +
  geom_text(aes(y = pct + .04, label = paste0(lab*100, '%')), position = position_dodge(width = .9), size = 8, family = "font") +
  labs(x = "", y = "", title = "What is your current marital or cohabiting status?", caption = '@ryanburge + @religiondata\nData: National Survey of Family Growth, 2015-2017')
save("marstat_reltrad_arda.png", ht = 8)

gg <- fg %>% 
  filter(oppyearnum <= 20) %>% 
  filter(marstat == 6) %>% 
  filter(ager >= 25) %>% 
  mutate(reltrad = factor(dplyr::recode(reltrad, !!!reltrad_labels),
                          levels = c("Evangelical", "Mainline", "Black Prot.", "Catholic", "Other religion", "No Religion"))) %>% 
  mutate(partners = oppyearnum) %>% 
  mutate(partners = frcode(partners == 0 ~ "Zero", 
                           partners == 1 ~ "One", 
                           partners == 2 ~ "Two", 
                           partners == 3 ~ "Three or More")) %>% 
  group_by(reltrad) %>% 
  ct(partners, wt = weight, show_na = FALSE) %>% filter(reltrad != 'NA')


gg %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = reltrad, y = pct, fill = reltrad)) +
  geom_col(color = "black") +
  facet_wrap(~ partners) +
  theme_rb() +
  scale_fill_manual(values = color_mapping) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(labels = percent) +
  theme(strip.text = element_text(size = 20)) +
  geom_text(aes(y = pct + .04, label = paste0(lab*100, '%')), position = position_dodge(width = .9), size = 8, family = "font") +
  labs(x = "", y = "", title = "Number of opposite-sex partners in last 12 months", 
       subtitle = "Among Never Been Married Respondents, 25-50 Years Old", 
       caption = '@ryanburge + @religiondata\nData: National Survey of Family Growth, 2015-2017')
save("num_partners_unmarried_reltrad.png", ht = 8)
  


gg <- fg %>% 
  mutate(reltrad = factor(dplyr::recode(reltrad, !!!reltrad_labels),
                          levels = c("Evangelical", "Mainline", "Black Prot.", "Catholic", "Other religion", "No Religion"))) %>% 
  mutate(intact = case_when(INTACT18 == 1 ~ 1, 
                            INTACT18 == 2 ~ 0)) %>% 
  group_by(reltrad) %>% 
  mean_ci(intact, wt = weight) %>% filter(reltrad != 'NA')


gg %>% 
  mutate(lab = round(mean, 2)) %>% 
  ggplot(., aes(x = reorder(reltrad, mean), y = mean, fill = reltrad)) +
  geom_col(color = "black") +
  coord_flip() +
  theme_rb()+
  y_pct() + 
  scale_fill_manual(values = moma.colors("Panton", 6)) +
  theme(plot.title = element_text(size = 18)) +
  geom_text(aes(y = mean + .045, label = ifelse(mean >= .095, paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 11, family = "font") +
  geom_text(aes(y = mean + .03, label = ifelse(mean < .094, paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 8, family = "font") +
  labs(x = "", y = "", title = "Share who lived in an intact family from birth to age 18 ", caption = '@ryanburge\nData: National Survey of Family Growth, 2015-2017')
save("fam_growth_intact_fam_new.png", ht = 5)


gg <- fg %>% 
  mutate(rel = frcode(reltrad == 1  | reltrad == 2 | reltrad == 3 ~ "Protestant", 
                      reltrad == 4 ~ "Catholic", 
                      reltrad == 5 ~ "Other Religion", 
                      reltrad == 6 ~ "No Religion")) %>% 
  mutate(preg = case_when(numpregs == 0 ~ 1, 
                          numpregs <= 20 ~ 0)) %>% 
  group_by(ager, rel) %>% 
  mean_ci(preg, wt = weight) %>% 
  na.omit()


color_mapping <- c(
  "Evangelical" = "#28666e",
  "Mainline" = "#5f2680",
  "Black Prot." = "#B5B682",
  "Catholic" = "#7D3C98",
  "Other Religion" = "#033f63",
  "No Religion" = "#FEDC97"
)

gg %>% 
  ggplot(., aes(x = ager, y= mean, color = rel, group = rel)) +
  geom_point(stroke = .25, shape = 21, alpha = .5) +
  geom_smooth(se = FALSE) +
  theme_rb(legend = TRUE) +
  scale_color_manual(values = c("#033f63", "#B5B682", "#5f2680", "#FEDC97")) +
  scale_y_continuous(labels = percent, limits = c(0, 1)) +
  labs(x = "", y = "", title = "How Many Times Have You Been Pregnant in Your Life? - Share Saying Zero", caption = '@ryanburge + @religiondata\nData: National Survey of Family Growth, 2015-2017')
save("never_preg.png")



gg <- fg %>% 
  filter(ager >= 35) %>% 
  mutate(rel = frcode(reltrad == 1  | reltrad == 2 | reltrad == 3 ~ "Protestant", 
                      reltrad == 4 ~ "Catholic", 
                      reltrad == 5 ~ "Other Religion", 
                      reltrad == 6 ~ "No Religion")) %>% 
  mutate(num = frcode(numbabes == 0 ~ "0", 
                      numbabes == 1 ~ "1", 
                      numbabes == 2 ~ "2", 
                      numbabes == 3 ~ "3", 
                      numbabes >= 4 & numbabes <= 14 ~ "4+")) %>% 
  group_by(rel) %>% 
  ct(num, wt = weight, show_na = FALSE)


gg %>% 
  filter(rel != "NA") %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = fct_rev(num), y = pct, fill = rel)) +
  geom_col(color = "black") +
  coord_flip() +
  theme_rb()+
  facet_wrap(~ rel) +
  scale_y_continuous(labels = percent, limits = c(0, .43)) + 
  scale_fill_manual(values = color_mapping) +
  theme(plot.title = element_text(size = 18)) +
  geom_text(aes(y = pct + .035, label = ifelse(pct >= 0, paste0(lab*100, '%'), "")), position = position_dodge(width = .9), size = 8, family = "font") +
  theme(strip.text = element_text(size = 20))+
  labs(x = "", y = "", title = "Number of Babies Born Alive, Women Aged 35-50", caption = '@ryanburge + @religiondata\nData: National Survey of Family Growth, 2015-2017')
save("number_babies_rel.png")


regg <- fg %>% 
  filter(ager >= 35) %>%
  mutate(marstat = case_when(marstat == 6 ~ 1,
                             marstat <= 5 ~ 0)) %>% 
  mutate(educ = higrade) %>% 
  mutate(income = totincr) %>% 
  mutate(white = case_when(race == 1 ~ 1, 
                          TRUE ~ 0)) %>% 
  mutate(rel = frcode(reltrad == 1  ~ "Evangelical", 
                      reltrad == 2 ~ "Mainline", 
                      reltrad == 3 ~ "Black Prot.", 
                      reltrad == 4 ~ "Catholic", 
                      reltrad == 5 ~ "Other Religion", 
                      reltrad == 6 ~ "No Religion")) %>% 
  mutate(number = numbabes) %>% 
  select(income, educ, white, rel, number, marstat) %>% as_tibble()
  
regg <- regg %>% 
  mutate(rel = relevel(rel, ref = "Evangelical"))

library(jtools)

model <- lm(number ~ white + educ + income + rel + marstat, data = regg)
summary(model)




regg <- fg %>% 
  mutate(age = ager) %>% 
  mutate(marstat = case_when(marstat == 6 ~ 1,
                             marstat <= 5 ~ 0)) %>% 
  mutate(educ = higrade) %>% 
  mutate(income = totincr) %>% 
  mutate(white = case_when(race == 1 ~ 1, 
                           TRUE ~ 0)) %>% 
  mutate(rel2 = frcode(reltrad <= 5 ~ "Religious", 
                       reltrad == 6 ~ "No Religion")) %>% 
  mutate(number = numbabes) %>% 
  select(income, educ, white, rel2, number, marstat, age) %>% as_tibble()



library(jtools)

model <- lm(number ~ age*rel2 + white + educ + income + marstat, data = regg)
summary(model)

library(interactions)

gg <- interact_plot(model, pred = age, modx = rel2, int.width = .76, interval = TRUE) 

gg + 
  theme_rb() +
  scale_color_manual(values = c("#B5B682", "#033f63")) +
  scale_fill_manual(values = c("#B5B682", "#033f63")) +
  add_text(x = 35, y = 1.70, word = "Religious", sz = 7) +
  add_text(x = 35, y = 1.0, word = "Non-Religious", sz = 7) +
  labs(x = "Age", y = "Predicted Number of Children", 
       title = "Predicted Number of Babies Born Alive for Women by Age and Religion", 
       subtitle = "Controls for Race, Education, Income, and Marital Status", 
       caption = '@ryanburge + @religiondata\nData: National Survey of Family Growth, 2015-2017')
save("reg_numkids.png")