2_tablesfigures.Rmd

---
title: "Wage Disparity Project"
date: "`r format(Sys.time(), '%d %B %Y')`"
output: 
  word_document:
    highlight: NULL
    reference_docx: word_styles3.docx
editor_options: 
  chunk_output_type: console
---

```{r, setup, include=FALSE}
knitr::opts_chunk$set(
  echo=FALSE, message=FALSE, warning=FALSE, comment=NA, 
  tab.cap.sep = ". ",
  fig.height=8, fig.width=8
)
options(max.print = 2000, tibble.print_max = 100)
```


```{r results='hide'}
#empty environment.
rm(list = ls())

# Libraries.
library(readr)
library(dplyr)
library(tidyr)
library(stringr)
library(janitor)
library(flextable)
library(officer)
library(survey)
library(ggplot2)
library(viridis)



# flextable custom theme
my_theme <- function(ft, pgwidth = 7){
  ft_out <- ft %>% 
    autofit() %>% 
    border_remove() %>% 
    hline_top(border = fp_border(color="black", width = 1), part = "header") %>% 
    hline_top(border = fp_border(color="black", width = 1)) %>% 
    hline_bottom(border = fp_border(color="black", width = 1)) %>% 
    bold(part = "header")
  ft_out <- width(ft_out, width = dim(ft_out)$widths*pgwidth /(flextable_dim(ft_out)$widths))
  return(ft_out)
}
set_flextable_defaults(
  line_spacing = 1.08,
  font.family = "Times New Roman",
  font.size = 12,
  padding = 3
)


# Folder path.
fpath <- "..."

# Data.
df0 <- readr::read_csv(file = file.path(fpath, "Data_final.csv"), show_col_types = F) %>% as.data.frame() %>% 
  # relevel the factors.
  mutate(EmplStat = factor(EmplStat, levels = c("FullTime","PartTime","Other","Unemployed"))) %>% 
  mutate(AGEG10LFS_T = factor(AGEG10LFS_T, levels = c("2534", "3544", "4554"))) %>% 
  mutate(RACETHN_4CAT = factor(RACETHN_4CAT, levels = c("White", "Black", "Hispanic", "Other"))) %>% 
  mutate(EDCAT3 = factor(EDCAT3, levels = c("Secondary","Non-tertiary","Tertiary"))) %>% 
  mutate(GENDER_R = factor(GENDER_R, levels = c("Female","Male"))) %>% 
  mutate(withkids = case_when(withkids==1 ~ "Yes",
                              withkids==0 ~ "No"))
```





```{r results='asis'}
# //TABLE 2a.
# Sample Size Dist by gender, edcat, and race.
df0 %>% 
  count(EDCAT3,GENDER_R,RACETHN_4CAT) %>%
  pivot_wider(names_from = "RACETHN_4CAT", values_from = "n") %>%
  adorn_totals(c("row","col")) %>%
  adorn_percentages("row") %>%
  adorn_pct_formatting(digits=0) %>%
  adorn_ns() %>%
  rename("Ed Attain" = "EDCAT3", "Gender"="GENDER_R") %>%
  flextable() %>% my_theme() %>% 
  align(j = c(-1,-2), align="right", part = "all") %>% 
  set_caption("TABLE 2a. Sample size distribution - All Respondents")
cat("\n \\newline \n")


# //TABLE 2b.
df0 %>% 
  filter(C_Q07==1) %>% # full-timers only
  count(EDCAT3,GENDER_R,RACETHN_4CAT) %>%
  pivot_wider(names_from = "RACETHN_4CAT", values_from = "n") %>%
  adorn_totals(c("row","col")) %>%
  adorn_percentages("row") %>%
  adorn_pct_formatting(digits=0) %>%
  adorn_ns() %>%
  rename("Ed Attain" = "EDCAT3", "Gender"="GENDER_R") %>%
  flextable() %>% my_theme() %>% 
  align(j = c(-1,-2), align="right", part = "all") %>% 
  set_caption("TABLE 2b. Sample size distribution - Full-time Respondents")
cat("\n\n\\pagebreak\n")
```




```{r results='asis'}
# //TABLE 3a.
# Sample Size Dist by gender, edcat, and OccupCat.
df0 %>% 
  count(EDCAT3,GENDER_R,OccupCat) %>%
  pivot_wider(names_from = "OccupCat", values_from = "n") %>%
  adorn_totals(c("row","col")) %>%
  adorn_percentages("row") %>%
  adorn_pct_formatting(digits=0) %>%
  adorn_ns() %>%
  rename("Ed Attain" = "EDCAT3", "Gender"="GENDER_R") %>%
  flextable() %>% my_theme() %>% 
  align(j = c(-1,-2), align="right", part = "all") %>% 
  set_caption("TABLE 3a. Sample size distribution - All Respondents")
cat("\n \\newline \n")


# //TABLE 3b.
df0 %>% 
  filter(C_Q07==1) %>% # full-timers only
  count(EDCAT3,GENDER_R,OccupCat) %>%
  pivot_wider(names_from = "OccupCat", values_from = "n") %>%
  adorn_totals(c("row","col")) %>%
  adorn_percentages("row") %>%
  adorn_pct_formatting(digits=0) %>%
  adorn_ns() %>%
  rename("Ed Attain" = "EDCAT3", "Gender"="GENDER_R") %>%
  flextable() %>% my_theme() %>% 
  align(j = c(-1,-2), align="right", part = "all") %>% 
  set_caption("TABLE 3b. Sample size distribution - Full-time Respondents")
cat("\n\n\\pagebreak\n")
```




```{r results='asis'}
# //TABLE 4a.
# Weighted Descriptives by subgroups - All respondents
dsgnsubgps_allRs <- svydesign(ids = ~1, strata = NULL, weights = ~SPFWT0, 
                              data = df0 %>% 
                                mutate(subgrps = paste(EDCAT3,RACETHN_4CAT,GENDER_R, sep=":")))
# weighted means
tt_means <- svybys(~EARNMTHBONUSUS_C, by=~subgrps, dsgnsubgps_allRs, svymean, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  rename("mean.weighted"="EARNMTHBONUSUS_C") %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","RACETHN_4CAT","GENDER_R"), remove=T) 
# weighted sds
tt_sds <- svybys(~EARNMTHBONUSUS_C, by=~subgrps, dsgnsubgps_allRs, svyvar, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  rename("var.weighted"="EARNMTHBONUSUS_C") %>% 
  mutate(sd.weighted = sqrt(var.weighted)) %>% select(-var.weighted) %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","RACETHN_4CAT","GENDER_R"), remove=T)
# simple stats
tt_simplestats <- df0 %>% 
  group_by(EDCAT3,RACETHN_4CAT,GENDER_R) %>% 
  summarise_at(c("EARNMTHBONUSUS_C"), list(~ min(., na.rm = T), ~ median(., na.rm = TRUE), ~ max(., na.rm = T), ~ n())) %>% 
  ungroup() %>% as.data.frame() 

# merge all
tt_allstats <- tt_simplestats %>%
  full_join(tt_means, by = c("EDCAT3","RACETHN_4CAT","GENDER_R")) %>%
  full_join(tt_sds, by = c("EDCAT3","RACETHN_4CAT","GENDER_R"))

# table form
tt_allstats %>% 
  relocate(c(min:n), .after = "sd.weighted") %>% 
  mutate(across(4:8, round, -2)) %>%
  mutate(across(4:8, ~./100)) %>%
  mutate(across(c(mean.weighted:max), format, nsmall = 0)) %>%
  adorn_totals(c("row")) %>%
  rename("Gender"="GENDER_R","Race/Ethnicity"="RACETHN_4CAT", "Ed Attain" = "EDCAT3",
         "Mean (weighted)"="mean.weighted", "SD (weighted)"="sd.weighted",
         "Min"="min", "Median"="median", "Max"="max", "N"="n") %>%
  flextable() %>% my_theme() %>% 
  align(j = c(-1,-2,-3), align="center", part = "all") %>% 
  set_caption("TABLE 4a. Descriptive statistics of top-coded monthly income with bonus (EARNMTHBONUSUS_C) by subgroups - All Respondents")
cat("\n\n\\pagebreak\n")


# save final table to use in bar charts
tt_allstats_copy <- tt_allstats
```


```{r results='asis'}
# //TABLE 4b.
# Weighted Descriptives by subgroups - Full-time respondents
dsgnsubgps_FTRs <- svydesign(ids = ~1, strata = NULL, weights = ~SPFWT0, 
                             data = df0 %>% 
                               filter(C_Q07==1) %>% # full-timers only
                               mutate(subgrps = paste(EDCAT3,RACETHN_4CAT,GENDER_R, sep=":")))

# weighted means
tt_means <- svybys(~EARNMTHBONUSUS_C, by=~subgrps, dsgnsubgps_FTRs, svymean, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  rename("mean.weighted"="EARNMTHBONUSUS_C") %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","RACETHN_4CAT","GENDER_R"), remove=T) 
# weighted sds
tt_sds <- svybys(~EARNMTHBONUSUS_C, by=~subgrps, dsgnsubgps_FTRs, svyvar, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  rename("var.weighted"="EARNMTHBONUSUS_C") %>% 
  mutate(sd.weighted = sqrt(var.weighted)) %>% select(-var.weighted) %>%  
  separate(subgrps, sep=":", into = c("EDCAT3","RACETHN_4CAT","GENDER_R"), remove=T)
# simple stats
tt_simplestats <- df0 %>% filter(C_Q07==1) %>% # full-timers only
  group_by(EDCAT3,RACETHN_4CAT,GENDER_R) %>% 
  summarise_at(c("EARNMTHBONUSUS_C"), list(~ min(., na.rm = T), ~ median(., na.rm = TRUE), ~ max(., na.rm = T), ~ n())) %>% 
  ungroup() %>% as.data.frame() 

# merge all
tt_allstats <- tt_simplestats %>%
  full_join(tt_means, by = c("EDCAT3","RACETHN_4CAT","GENDER_R")) %>%
  full_join(tt_sds, by = c("EDCAT3","RACETHN_4CAT","GENDER_R"))

# table form
tt_allstats %>% 
  relocate(c(min:n), .after = "sd.weighted") %>% 
  mutate(across(4:8, round, -2)) %>%
  mutate(across(4:8, ~./100)) %>%
  mutate(across(c(mean.weighted:max), format, nsmall = 0)) %>%
  adorn_totals(c("row")) %>%
  rename("Gender"="GENDER_R","Race/Ethnicity"="RACETHN_4CAT", "Ed Attain" = "EDCAT3",
         "Mean (weighted)"="mean.weighted", "SD (weighted)"="sd.weighted",
         "Min"="min", "Median"="median", "Max"="max", "N"="n") %>%
  flextable() %>% my_theme() %>% 
  align(j = c(-1,-2,-3), align="center", part = "all") %>% 
  set_caption("TABLE 4b. Descriptive statistics of top-coded monthly income with bonus (EARNMTHBONUSUS_C) by subgroups - Full-time Respondents")
cat("\n\n\\pagebreak\n")
```




```{r results='asis', fig.height=6, fig.width=6}
# //FIGURE 1.
# Boxplots.
df0 %>% 
  ggplot(aes(x=RACETHN_4CAT, y=round(EARNMTHBONUSUS_C, -2)/100, fill = RACETHN_4CAT)) + 
  geom_boxplot() +
  scale_fill_viridis(discrete = T, option = "D") +
  theme_light() + xlab("\nRace/Ethnicity") + ylab("Montly income with bonus\n") + theme(legend.position = "none") +
  theme(plot.title = element_text(size=10),
        plot.subtitle = element_text(size=8))
cat("\n\n\\pagebreak\n")
```




```{r results='asis'}
# //FIGURE 2.
# Across ed attainment levels; within gender categories
tt_allstats_copy %>% select(GENDER_R,RACETHN_4CAT,EDCAT3,mean.weighted) %>% 
  ggplot(aes(x=factor(RACETHN_4CAT, levels = c("White", "Black", "Hispanic", "Other")), 
             y = round(mean.weighted,-2)/100, 
             fill=factor(EDCAT3, levels = c("Secondary","Non-tertiary","Tertiary")))) +
  geom_bar(position = "dodge", stat="identity") +
  scale_fill_viridis(discrete = T, option = "C") +
  facet_wrap(.~factor(GENDER_R), ncol=1) +
  geom_text(aes(label = format(round(mean.weighted,-2)/100,nsmall = 0)),
            vjust = -0.5, size = 3.5, position = position_dodge(width = .9)) +
  ylim(0,83) +
  theme_light() +
  theme(legend.position="top",
        strip.text = element_text(colour = 'black'), 
        plot.title = element_text(size=12)) +
  labs(y="Weighted means of top-coded monthly income with bonus\n", x=element_blank(),
       fill = "Ed Attainment")



# //FIGURE 3.
# Within ed attainment levels; across race categories 
tt_allstats_copy %>% select(GENDER_R,RACETHN_4CAT,EDCAT3,mean.weighted) %>% 
  ggplot(aes(x=GENDER_R, 
             y = round(mean.weighted,-2)/100, 
             fill=factor(RACETHN_4CAT, levels = c("White", "Black", "Hispanic", "Other")))) +
  geom_bar(position = "dodge", stat="identity") +
  scale_fill_viridis(discrete = T, option = "D") +
  facet_wrap(.~factor(EDCAT3, levels = c("Secondary","Non-tertiary","Tertiary")), ncol=1) +
  geom_text(aes(label = format(round(mean.weighted,-2)/100,nsmall = 0)),
            vjust = -0.5, size = 3.5, position = position_dodge(width = .9)) +
  ylim(0,83) +
  theme_light() +
  theme(legend.position="top",
        strip.text = element_text(colour = 'black'), 
        plot.title = element_text(size=12)) +
  labs(y="Weighted means of top-coded monthly income with bonus\n", x=element_blank(),
       fill = "Race/ethnicity")



# //FIGURE 4.
# Within ed attainment levels; across gender categories
tt_allstats_copy %>% select(GENDER_R,RACETHN_4CAT,EDCAT3,mean.weighted) %>% 
  ggplot(aes(x=factor(RACETHN_4CAT, levels = c("White", "Black", "Hispanic", "Other")), 
             y = round(mean.weighted,-2)/100, 
             fill=GENDER_R)) +
  geom_bar(position = "dodge", stat="identity") +
  scale_fill_viridis(discrete = T, option = "E") +
  facet_wrap(.~factor(EDCAT3, levels = c("Secondary","Non-tertiary","Tertiary")), ncol=1) +
  geom_text(aes(label = format(round(mean.weighted,-2)/100,nsmall = 0)),
            vjust = -0.5, size = 3.5, position = position_dodge(width = .9)) +
  ylim(0,83) +
  theme_light() +
  theme(legend.position="top",
        strip.text = element_text(colour = 'black'), 
        plot.title = element_text(size=12)) +
  labs(y="Weighted means of top-coded monthly income with bonus\n", x=element_blank(),
       fill = "Gender")
cat("\n\n\\pagebreak\n")
```






```{r results='asis'}
# //TABLE 6a.
# Cog Skill distr within ed attain
dsgnsubgps <- svydesign(ids = ~1, strata = NULL, weights = ~SPFWT0, 
                        data = df0 %>% 
                          group_by(EDCAT3) %>% 
                          mutate(grp = paste0("Cog_Q",ntile(STLITNUM1, 3))) %>% 
                          ungroup() %>% 
                          mutate(subgrps = paste(EDCAT3,grp, sep=":")))

tt_lit1_means <- svybys(~(PVLIT1), by=~subgrps, dsgnsubgps, svymean, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","PR"), remove=T) %>% 
  mutate(PR = str_replace(PR, "Cog", "PVLIT1")) %>% 
  pivot_wider(names_from = PR, values_from = PVLIT1)

tt_num1_means <- svybys(~(PVNUM1), by=~subgrps, dsgnsubgps, svymean, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","PR"), remove=T) %>% 
  mutate(PR = str_replace(PR, "Cog", "PVNUM1")) %>% 
  pivot_wider(names_from = PR, values_from = PVNUM1)

tt_stlitnum1_means <- svybys(~(STLITNUM1), by=~subgrps, dsgnsubgps, svymean, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","PR"), remove=T) %>% 
  mutate(PR = str_replace(PR, "Cog", "STLITNUM1")) %>% 
  pivot_wider(names_from = PR, values_from = STLITNUM1)


tt_simplestats %>% select(EDCAT3) %>% distinct() %>% 
  full_join(tt_lit1_means, by = c("EDCAT3")) %>%
  full_join(tt_num1_means, by = c("EDCAT3")) %>%
  full_join(tt_stlitnum1_means, by = c("EDCAT3")) %>% 
  mutate(across(c(-1), round, 1)) %>%
  flextable() %>% my_theme() %>% 
  align(j = c(-1), align="center", part = "all") %>% 
  set_caption("TABLE 6a. Average weighted means of skills by terciles and educational attainment - All Respondents")
cat("\n \\newline \n")
```


```{r results='asis'}
# //TABLE 6b.
# Cog Skill distr within ed attain - Fulltimers
dsgnsubgps <- svydesign(ids = ~1, strata = NULL, weights = ~SPFWT0, 
                        data = df0 %>% 
                          filter(C_Q07==1) %>% # full-timers only
                          group_by(EDCAT3) %>% 
                          mutate(grp = paste0("Cog_Q",ntile(STLITNUM1, 3))) %>% 
                          ungroup() %>% 
                          mutate(subgrps = paste(EDCAT3,grp, sep=":")))

tt_lit1_means <- svybys(~(PVLIT1), by=~subgrps, dsgnsubgps, svymean, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","PR"), remove=T) %>% 
  mutate(PR = str_replace(PR, "Cog", "PVLIT1")) %>% 
  pivot_wider(names_from = PR, values_from = PVLIT1)

tt_num1_means <- svybys(~(PVNUM1), by=~subgrps, dsgnsubgps, svymean, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","PR"), remove=T) %>% 
  mutate(PR = str_replace(PR, "Cog", "PVNUM1")) %>% 
  pivot_wider(names_from = PR, values_from = PVNUM1)

tt_stlitnum1_means <- svybys(~(STLITNUM1), by=~subgrps, dsgnsubgps, svymean, na.rm=T) %>% 
  as.data.frame() %>% tibble::rownames_to_column() %>% 
  select(-rowname, -se) %>% 
  separate(subgrps, sep=":", into = c("EDCAT3","PR"), remove=T) %>% 
  mutate(PR = str_replace(PR, "Cog", "STLITNUM1")) %>% 
  pivot_wider(names_from = PR, values_from = STLITNUM1)


tt_simplestats %>% select(EDCAT3) %>% distinct() %>% 
  full_join(tt_lit1_means, by = c("EDCAT3")) %>%
  full_join(tt_num1_means, by = c("EDCAT3")) %>%
  full_join(tt_stlitnum1_means, by = c("EDCAT3")) %>% 
  mutate(across(c(-1), round, 1)) %>%
  flextable() %>% my_theme() %>% 
  align(j = c(-1), align="center", part = "all") %>% 
  set_caption("TABLE 6b. Average weighted means of skills by terciles and educational attainment - Full-time Respondents")
cat("\n\n\\pagebreak\n")
```




```{r results='asis'}
# //TABLE 7a.
# skills premium: gain in income for greater skills within subgroups
dsgnsubgps_allRs <- svydesign(ids = ~1, strata = NULL, weights = ~SPFWT0, 
                              data = df0 %>% 
                                group_by(EDCAT3,RACETHN_4CAT,GENDER_R) %>% 
                                mutate(grp = paste0("Cog_Q",ntile(STLITNUM1, 3))) %>% 
                                ungroup() %>% 
                                mutate(subgrps = paste(EDCAT3,RACETHN_4CAT,GENDER_R,grp, sep=":")))

tt_simplestats %>% select(EDCAT3,RACETHN_4CAT,GENDER_R) %>% 
  full_join(svybys(~EARNMTHBONUSUS_C, by=~subgrps, dsgnsubgps_allRs, svymean, na.rm=T) %>% as.data.frame() %>% 
              tibble::rownames_to_column() %>% 
              select(-rowname, -se) %>% 
              separate(subgrps, sep=":", into = c("EDCAT3","RACETHN_4CAT","GENDER_R","PR"), remove=T) %>% 
              pivot_wider(names_from = PR, values_from = EARNMTHBONUSUS_C), 
            by = c("EDCAT3","RACETHN_4CAT","GENDER_R")) %>% 
  mutate(across(4:6, round, -2)) %>%
  mutate(across(4:6, ~./100)) %>%
  mutate(across(4:6, format, nsmall = 0)) %>% 
  rename("Gender"="GENDER_R", "Race/Ethnicity"="RACETHN_4CAT", "Ed Attain"="EDCAT3",
         "Low Proficiency"="Cog_Q1", "Medium Proficiency"="Cog_Q2", "High Proficiency"="Cog_Q3") %>% 
  flextable() %>% my_theme() %>% 
  align(j = c(-1,-2,-3), align="right", part = "all") %>% 
  set_caption("TABLE 7a. Average weighted means of top-coded monthly income with bonus (EARNMTHBONUSUS_C) by terciles of cognitive skills (STLITNUM1) and subgroups - All Respondents")
cat("\n\n\\pagebreak\n")


# //TABLE 7b.
dsgnsubgps_FTRs <- svydesign(ids = ~1, strata = NULL, weights = ~SPFWT0, 
                             data = df0 %>% 
                               filter(C_Q07==1) %>% # full-timers only
                               group_by(EDCAT3,RACETHN_4CAT,GENDER_R) %>% 
                               mutate(grp = paste0("Cog_Q",ntile(STLITNUM1, 3))) %>% 
                               ungroup() %>% 
                               mutate(subgrps = paste(EDCAT3,RACETHN_4CAT,GENDER_R,grp, sep=":")))

tt_simplestats %>% select(EDCAT3,RACETHN_4CAT,GENDER_R) %>% 
  full_join(svybys(~EARNMTHBONUSUS_C, by=~subgrps, dsgnsubgps_FTRs, svymean, na.rm=T) %>% as.data.frame() %>%
              tibble::rownames_to_column() %>% 
              select(-rowname, -se) %>% 
              separate(subgrps, sep=":", into = c("EDCAT3","RACETHN_4CAT","GENDER_R","PR"), remove=T) %>% 
              pivot_wider(names_from = PR, values_from = EARNMTHBONUSUS_C), 
            by = c("EDCAT3","RACETHN_4CAT","GENDER_R")) %>% 
  mutate(across(4:6, round, -2)) %>%
  mutate(across(4:6, ~./100)) %>%
  mutate(across(4:6, format, nsmall = 0)) %>% 
  rename("Gender"="GENDER_R", "Race/Ethnicity"="RACETHN_4CAT", "Ed Attain"="EDCAT3",
         "Low Proficiency"="Cog_Q1", "Medium Proficiency"="Cog_Q2", "High Proficiency"="Cog_Q3") %>% 
  flextable() %>% my_theme() %>% 
  align(j = c(-1,-2,-3), align="right", part = "all") %>% 
  set_caption("TABLE 7b. Average weighted means of top-coded monthly income with bonus (EARNMTHBONUSUS_C) by terciles of cognitive skills (STLITNUM1) and subgroups - Full-time Respondents")
```