# Make regression coefficient plots for paper

dir.create('output/regressions/coef_plots_paper')

# Main coef plots: 
coefs = read_csv('output/regressions/all_coefs.csv')

# Make clean names
fullvars = c("(Intercept)","logpop10",
              "pct_b_pop","pct_o_pop",  "pct_all_cit_imm",
             "union_combined","cemp_pct_manuf","cemp_pct_public",
              "dissimilarity_index",'pct_ba_lograt_bw',
              'incar_comb_bwlrat','hrrv',
              "regionMidwest","regionSouth","regionWest", 
              "pct_all_col_ba",
             'incar_comb_total','incar_comb_bwdif','total_prison_pop_rate','imprate_vera_lrat_bw',
              'dissimilarity_index_ltdb','iso_black','iso_black_ltdb','exposure_b_w','exposure_b_w_ltdb','theil','theil_ltdb',
              'slave_pct:hrrvalid','cropfarms_pctblack:hrrvalid','hrrbaker:hrrvalid','hrr4:hrrvalid','numdis:hrrvalid','southman:hrrvalid','hrr_st:hrrvalid','hrr4_st:hrrvalid'
             )

cleanvars = c('Intercept','MSA population (log10)', 
              '% Black','% neither Black nor white','% foreign-born',
              'Unionization rate','% Manufacturing','% Public administration',
              'BW dissimilarity index','Logged BA ratio',
              'Logged incarc. ratio','HRR factor',
              'Midwest','South','West',
              "% with bachelor's degree",
              'State incarceration rate','BW difference in incarc.','State imprisionment rate','Logged imp. ratio',
              'BW dissimilarity index','Black isolation index','Black isolation index','BW exposure index','BW exposure index','BW Theil index','BW Theil index',
              '% Enslaved, 1860','% Black\nsharecroppers, 1930','HRR factor','HRR factor','Number of\ndisfranchisement devices','Support for\nSouthern Manifesto','HRR factor','HRR factor'
)

coefs$clean_names = factor(coefs$var, levels = (fullvars), 
                           labels = (cleanvars))


# Create bounds on confidence intervals
coefs$smin = coefs$coef - coefs$se * 1.96
coefs$smax = coefs$coef + coefs$se * 1.96

# Also get min and max for standardized
coefs$stmin = coefs$standardcoef - coefs$standardcoef / coefs$t * 1.96 
coefs$stmax = coefs$standardcoef + coefs$standardcoef / coefs$t * 1.96 


## F4 - % Black, with and without HRR
incs = 'all_faminc_sqrt'
outvs = c('med_ldiff','med_b','med_w')
rdef = 'expansive'
modints = c('Baseline','HRR factor excluded')

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              racedef == rdef,
              var == 'pct_b_pop',
              incsamp == incs)


pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$outvarclean = factor(pdta$outvar, levels = c('med_ldiff','med_b','med_w'), labels = c('Median income ratio','Black median income','White median income'))

pdf(paste0('output/regressions/coef_plots_paper/f4_pct_black.pdf'),height = 4,width = 6)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(mod ~ outvarclean) +
  geom_vline(alpha = .5,xintercept = 0)+
  scale_x_continuous(name = 'Standardized coefficient') +
  scale_y_discrete(name = '')

print(gplt)
dev.off()


# Figure 5 - Industry and union

varsint = c('union_combined','cemp_pct_manuf','cemp_pct_public')
modints = c('Baseline')

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              racedef == rdef,
              var %in% varsint,
              incsamp == incs)


pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$outvarclean = factor(pdta$outvar, levels = c('med_ldiff','med_b','med_w'), labels = c('Median income ratio','Black median income','White median income'))

pdf(paste0('output/regressions/coef_plots_paper/f5_formal_labor.pdf'),height = 5.5,width = 6)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ outvarclean) +
  geom_vline(alpha = .5,xintercept = 0)+
  scale_x_continuous(name = 'Standardized coefficient') +
  scale_y_discrete(name = '')

print(gplt)
dev.off()


## Figure 6 - Segregation and college and incarceration and HRR
varsint = c('dissimilarity_index','pct_ba_lograt_bw','incar_comb_bwlrat','hrrv')
modints = c('Baseline')
pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              racedef == rdef,
              var %in% varsint,
              incsamp == incs)


pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$outvarclean = factor(pdta$outvar, levels = c('med_ldiff','med_b','med_w'), labels = c('Median income ratio','Black median income','White median income'))


pdf(paste0('output/regressions/coef_plots_paper/f6_newmech.pdf'),height = 7.5,width = 6)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ outvarclean) +
  geom_vline(alpha = .5,xintercept = 0)+
  scale_x_continuous(name = 'Standardized coefficient') +
  scale_y_discrete(name = '')

print(gplt)
dev.off()



## Appendix figures 

## 1 - Alternative income definitions

modints = c('Baseline')
incs = c('all_faminc_sqrt','all_hhinc_sqrt','all_faminc_norm','all_faminc_tot','gq_hh_spmuinc_sqrt','gq_hh_spmuinc_disp_sqrt')
rdef = c('expansive') 
outvs = c('med_ldiff','med_rgap')
varsint = c('pct_b_pop','union_combined','cemp_pct_manuf','cemp_pct_public','dissimilarity_index','pct_ba_lograt_bw','hrrv','incar_comb_bwlrat')

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              var %in% varsint,
              racedef %in% rdef,
              incsamp %in% incs
)


pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$cleansamp = factor(pdta$incsamp, 
                        levels = incs ,
                        labels = c('Baseline','Household income','Family income\nnormalized by\nnumber of members','Family income\n(not normalized)','SPM sample, baseline\nincome definition','SPM sample,\ndisposable income')) #,'Narrow definition\nof "Black," US-born\ncitizens only')) 'Matched on age\nand sex','Matched on age,\nfamily size, and sex'
pdta = mutate(pdta,scen = case_when(
  incsamp != 'all_faminc_sqrt' ~ as.character(cleansamp),
  incsamp == 'all_faminc_sqrt' & outvar == 'med_ldiff' ~ 'Baseline',
  incsamp == 'all_faminc_sqrt' & outvar == 'med_rgap' ~ 'Income rank gap',

)
)

pdta$scen = factor(pdta$scen, 
                   levels = c('Baseline','Narrow definition\nof "Black"',#'Narrow definition\nof "Black," US-born\ncitizens only',
                              'Household income','Family income\nnormalized by\nnumber of members','Family income\n(not normalized)','SPM sample, baseline\nincome definition','SPM sample,\ndisposable income','Income rank gap'))

pdf(paste0('output/regressions/coef_plots_paper/figa_incdefs_robust.pdf'),height = 12,width = 10)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ scen) +
  geom_vline(alpha = .5,xintercept = 0)+
  theme(strip.text = element_text(size = 7)) +
  scale_x_continuous(name = 'Standardized coefficient', breaks = c(-.4,0,.4)) +
  scale_y_discrete(name = '')

print(gplt)
dev.off()


## 2 - Alternative population samples

modints = c('Baseline') #'White population weights','No population weights',
incs = c('all_faminc_sqrt','gq_hh_faminc_sqrt','emp_emp_inctot','cit_usborn_faminc_sqrt','mig_samest_faminc_sqrt','mig_difst_faminc_sqrt')
rdef = c('expansive','narrow')
outvs = 'med_ldiff'
varsint = c('pct_b_pop','union_combined','cemp_pct_manuf','cemp_pct_public','dissimilarity_index','pct_ba_lograt_bw','hrrv','incar_comb_bwlrat')

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              var %in% varsint,
              racedef %in% rdef,
              incsamp %in% incs
)


pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$cleansamp = factor(pdta$incsamp, 
                        levels = c('all_faminc_sqrt','gq_hh_faminc_sqrt','cit_usborn_faminc_sqrt','mig_samest_faminc_sqrt','mig_difst_faminc_sqrt','emp_emp_inctot'),
                        labels = c('Baseline','Group quarters\nexcluded','US-born citizens only','Residents of\nbirth state only','Residents of\nnon-birth states only','Employed persons only\n(individual income)'))
pdta = mutate(pdta,scen = case_when(
  incsamp != 'all_faminc_sqrt' ~ as.character(cleansamp),
  incsamp == 'all_faminc_sqrt' & mod == 'Baseline' & racedef == 'expansive' ~ 'Baseline',
  incsamp == 'all_faminc_sqrt' & mod == 'Baseline' & racedef == 'narrow' ~ 'Narrow definition\nof "Black"',
  incsamp == 'all_faminc_sqrt' & mod == 'Consistent MSA sample' ~ 'Consistent\nMSA sample'
)
)

pdta$scen = factor(pdta$scen, 
                   levels = c('Baseline','Group quarters\nexcluded','US-born citizens only','Residents of\nbirth state only','Residents of\nnon-birth states only',
                              'Consistent\nMSA sample','Employed persons only\n(individual income)','Narrow definition\nof "Black"')) #,'White population\nweights','No population weights'

pdf(paste0('output/regressions/coef_plots_paper/figa_popsamp_robust.pdf'),height = 12,width = 10)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ scen) +
  geom_vline(alpha = .5,xintercept = 0)+
  theme(strip.text = element_text(size = 7)) +
  scale_x_continuous(name = 'Standardized coefficient', breaks = c(-.4,0,.4)) +
  scale_y_discrete(name = '')

print(gplt)
dev.off()


## 3 Alternative samples and weights

modints = c('Baseline','Consistent MSA sample','No non-MSA observations','No population weights','Full population weights')
incs = c('all_faminc_sqrt','all_faminc_sqrt_matched_b_age_sex') 
rdef = c('expansive') 
outvs = 'med_ldiff'
varsint = c('pct_b_pop','union_combined','cemp_pct_manuf','cemp_pct_public','dissimilarity_index','pct_ba_lograt_bw','hrrv','incar_comb_bwlrat')

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              var %in% varsint,
              racedef %in% rdef,
              incsamp %in% incs
)

pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$cleansamp = factor(pdta$incsamp,
                        levels = incs ,
                        labels = c('Baseline','Matched on age\nand sex')) #,'Narrow definition\nof "Black," US-born\ncitizens only')) 'Matched on age\nand sex','Matched on age,\nfamily size, and sex'
pdta = mutate(pdta,scen = case_when(
  incsamp != 'all_faminc_sqrt' ~ as.character(cleansamp),
  incsamp == 'all_faminc_sqrt' ~ mod 
)
)

pdta$scen = factor(pdta$scen, 
                   levels = c('Baseline','Consistent MSA sample','No non-MSA observations','No population weights','Full population weights','Matched on age\nand sex'))
                   
pdf(paste0('output/regressions/coef_plots_paper/figa_msasamp_robust.pdf'),height = 12,width = 9)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ scen) +
  geom_vline(alpha = .5,xintercept = 0)+
  theme(strip.text = element_text(size = 7)) +
  scale_x_continuous(name = 'Standardized coefficient', breaks = c(-.4,0,.4)) +
  scale_y_discrete(name = '')

print(gplt)
dev.off()

## 4 - alternative models
modints = c('Baseline','Region dummies','HRR excluded','Segregation and incarceration excluded\n(full MSA sample in each year)','Incarceration inequality only','Incarceration difference','Imprisonment rate')
incs = c('all_faminc_sqrt')
rdef = c('expansive')
outvs = 'med_ldiff'
varsint = c('pct_b_pop','union_combined','cemp_pct_manuf','cemp_pct_public','dissimilarity_index','pct_ba_lograt_bw','hrrv','incar_comb_bwlrat','incar_comb_bwdif','imprate_vera_lrat_bw',"regionMidwest","regionSouth","regionWest")

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              var %in% varsint,
              racedef %in% rdef,
              incsamp %in% incs
)

pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))

pdta$mod = factor(pdta$mod, 
                  levels = c('Baseline','Region dummies','HRR excluded','Segregation and incarceration excluded\n(full MSA sample in each year)','Incarceration inequality only','Incarceration difference','Imprisonment rate'),
                  labels = c('Baseline','Region dummies','HRR excluded','Segregation and incarceration\nexcluded (full MSA sample\nin each year)','Incarceration inequality only','Incarceration difference','Imprisonment rate')
)


pdf(paste0('output/regressions/coef_plots_paper/figa_mods_robust.pdf'),height = 16,width = 12.5)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ mod) +
  geom_vline(alpha = .5,xintercept = 0)+
  theme(strip.text = element_text(size = 7)) +
  scale_x_continuous(name = 'Standardized coefficient', breaks = c(-.4,0,.4)) +
  scale_y_discrete(name = '')

print(gplt)
dev.off()


## 5 - Exploration of segregation
modints = c('Baseline','Isolation index','Exposure index','Theil index','Dissimilarity index (consistent tract boundaries)','Isolation index (consistent tract boundaries)','Theil index (consistent tract boundaries)')
mdefint = 'msa3'
incs = c('all_faminc_sqrt')
rdef = c('expansive')
outvs = 'med_ldiff'
varsint = c('pct_b_pop','union_combined','cemp_pct_manuf','cemp_pct_public','dissimilarity_index','pct_ba_lograt_bw','incar_comb_bwlrat','hrrv','dissimilarity_index_ltdb','exposure_b_w','iso_black','iso_black_ltdb','theil','theil_ltdb')

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              var %in% varsint,
              racedef %in% rdef,
              incsamp %in% incs
)

pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$mod = factor(pdta$mod, 
                  levels = modints,
                  labels = c('Baseline','Isolation index','Exposure index','Theil index','Dissimilarity index\n(consistent tract boundaries)','Isolation index\n(consistent tract boundaries)','Theil index\n(consistent tract boundaries)')
)

pdf(paste0('output/regressions/coef_plots_paper/figa_seg_robust.pdf'),height = 14,width = 11)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ mod) +
  geom_vline(alpha = .5,xintercept = 0)+
  theme(strip.text = element_text(size = 7)) +
  scale_x_continuous(name = 'Standardized coefficient', breaks = c(-.4,0,.4)) +
  scale_y_discrete(name = '')

print(gplt)
dev.off()


## 6 - Exploration of HRR
modints = c('Baseline','4-component HRR','4-variable state-level HRR','Baker HRR','2-variable state-level HRR')
incs = c('all_faminc_sqrt')
rdef = c('expansive')
outvs = 'med_ldiff'
varsint = c('pct_b_pop','union_combined','cemp_pct_manuf','cemp_pct_public','dissimilarity_index','pct_ba_lograt_bw','incar_comb_bwlrat','hrrv','hrrbaker:hrrvalid','hrr:hrrvalid','hrrbaker:hrrvalid','hrr_st:hrrvalid','hrr4:hrrvalid','hrr4_st:hrrvalid','slave_pct:hrrvalid','cropfarms_pctblack:hrrvalid','numdis:hrrvalid','southman:hrrvalid')

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              var %in% varsint,
              racedef %in% rdef,
              incsamp %in% incs
)

pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$mod = factor(pdta$mod, levels = c('Baseline','4-component HRR','2-variable state-level HRR','4-variable state-level HRR','Baker HRR','Slavery only','Sharecropping only','Disfranchisement devices only','Southern Manifesto only'))

pdf(paste0('output/regressions/coef_plots_paper/figa_hrr_robust.pdf'),height = 12,width = 9)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ mod) +
  geom_vline(alpha = .5,xintercept = 0)+
  theme(strip.text = element_text(size = 7)) +
  scale_x_continuous(name = 'Standardized coefficient', breaks = c(-.4,0,.4)) +
  scale_y_discrete(name = '')

print(gplt)
dev.off()


## 6 - Exploration of HRR components
modints = c('Baseline','Slavery only','Sharecropping only','Disfranchisement devices only','Southern Manifesto only')
incs = c('all_faminc_sqrt')
rdef = c('expansive')
outvs = 'med_ldiff'
varsint = c('pct_b_pop','union_combined','cemp_pct_manuf','cemp_pct_public','dissimilarity_index','pct_ba_lograt_bw','incar_comb_bwlrat','hrr:hrrvalid','hrrbaker:hrrvalid','hrr:hrrvalid','hrrbaker:hrrvalid','hrr_st:hrrvalid','hrr4:hrrvalid','hrr4_st:hrrvalid','slave_pct:hrrvalid','cropfarms_pctblack:hrrvalid','numdis:hrrvalid','southman:hrrvalid')

pdta = filter(coefs, clean_names != "Intercept",
              mod %in% modints,
              outvar %in% outvs,
              var %in% varsint,
              racedef %in% rdef,
              incsamp %in% incs
)

pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$mod = factor(pdta$mod, levels = c('Baseline','4-component HRR','2-variable state-level HRR','4-variable state-level HRR','Baker HRR','Slavery only','Sharecropping only','Disfranchisement devices only','Southern Manifesto only'))

pdf(paste0('output/regressions/coef_plots_paper/figa_hrr_comp.pdf'),height = 14,width = 10.5)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ mod) +
  geom_vline(alpha = .5,xintercept = 0)+
  theme(strip.text = element_text(size = 7)) +
  scale_x_continuous(name = 'Standardized coefficient', breaks = c(-.4,0,.4)) +
  scale_y_discrete(name = '')

print(gplt)
dev.off()



## Appendix 4 - alternative percentiles
modint = 'Baseline'
incs = 'all_faminc_sqrt'
rdef = 'expansive'

outvs = c('med_ldiff','mean_ldiff','p10_ldiff','p25_ldiff','p75_ldiff','p90_ldiff')
varsint = c('pct_b_pop','union_combined','cemp_pct_manuf','cemp_pct_public','dissimilarity_index','pct_ba_lograt_bw','incar_comb_bwlrat','hrrv')

pdta = filter(coefs, clean_names != "Intercept",
              mod == modint,
              outvar %in% outvs,
              var %in% varsint,
              racedef == rdef,
              incsamp %in% c('all_faminc_sqrt')
)



pdta$year = factor(pdta$year, levels = c(2017,2012,2007,1999,1989,1979,1969))
pdta$outvar = factor(pdta$outvar, levels = c('p05_ldiff','p10_ldiff','p25_ldiff','med_ldiff','p75_ldiff','p90_ldiff','p95_ldiff','mean_ldiff'),
                     labels = c('5th pctile','10th pctile','25th pctile','Median (baseline)','75th pctile','90th pctile','95th pctile','Mean'))

pdf(paste0('output/regressions/coef_plots_paper/fig_a4_pctiles_robust.pdf'),height = 14,width = 10)
gplt = ggplot(pdta, aes(x = standardcoef, y = year)) +
  geom_point() +
  geom_segment(aes(x = stmin, xend = stmax,yend = year)) +
  facet_grid(clean_names ~ outvar) +
  geom_vline(alpha = .5,xintercept = 0)+
  scale_x_continuous(name = 'Standardized coefficient', breaks = c(-.4,0,.4)) +
  scale_y_discrete(name = '')

print(gplt)
dev.off()

