# Create one dataset with all the predictor variables for regressions

# Bring in output variables

mw = read_csv('intermediate_data/regression_data/msa_pred_rawdata.csv')
mw$st = str_pad(mw$st, 2,pad = '0')

## Make predictor variables
## Start with population shares calculated from outcome sample

# Baseline sample for income, population, etc.
controlvars = filter(mw, incsamp == 'all_faminc_sqrt')

# Percent of people with BA
mwcol = filter(mw, incsamp == 'edu_ba_faminc_sqrt') %>% select(year,racedef,msa,pop_all,pop_b,pop_o,pop_w)
controlvars = left_join(controlvars,mwcol, by = c('year','racedef','msa'),suffix = c('','_col_ba'))

# Immigrants = 1 - % born in US
mwimm = filter(mw, incsamp == 'cit_usborn_faminc_sqrt') %>% select(year,racedef,msa,pop_all,pop_b,pop_o,pop_w)
controlvars = left_join(controlvars,mwimm, by = c('year','racedef','msa'),suffix = c('','_cit_usborn'))

# Employment to population 
# mwemp = filter(mw, incsamp == 'emp_emp_inctot') %>% select(year,racedef,msa,pop_all,pop_b,pop_o,pop_w)
# controlvars = left_join(controlvars,mwemp, by = c('year','racedef','msa'),suffix = c('','_emp'))

# Group quarters
mwgc = filter(mw, incsamp == 'gq_hh_faminc_sqrt') %>% select(year,racedef,msa,pop_all,pop_b,pop_o,pop_w)
controlvars = left_join(controlvars,mwgc, by = c('year','racedef','msa'),suffix = c('','_gq_hh'))

# Institutionalized group quarters = one approach to 
# mwinst = filter(mw, incsamp == 'gq_inst_faminc_sqrt') %>% select(year,racedef,msa,pop_all,pop_b,pop_o,pop_w)
# controlvars = left_join(controlvars,mwinst, by = c('year','racedef','msa'),suffix = c('','_gq_inst'))

# Actually incarcerated group quarters - 1970 and 1980 only 
mwincarc = filter(mw, incsamp == 'gq_incarc_faminc_sqrt') %>% select(year,racedef,msa,pop_all,pop_b,pop_o,pop_w)
stopifnot(mwincarc$year %in% c(1969,1979))
controlvars = left_join(controlvars,mwincarc, by = c('year','racedef','msa'),suffix = c('','_gq_incarc'))


# Domestic migration
mwsamest = filter(mw, incsamp == 'mig_samest_faminc_sqrt') %>% select(year,racedef,msa,pop_all,pop_b,pop_o,pop_w)
mwdifst = filter(mw, incsamp == 'mig_difst_faminc_sqrt') %>% select(year,racedef,msa,pop_all,pop_b,pop_o,pop_w)
controlvars = left_join(controlvars,mwsamest, by = c('year','racedef','msa'),suffix = c('','_mig_samest'))
controlvars = left_join(controlvars,mwdifst, by = c('year','racedef','msa'),suffix = c('','_mig_difst'))

 
# Calculate percentages of population in each of these categories

for(r in c('b','o','w','all')) {
  # Percent of everyone
  controlvars[,paste('pct_',r,'_pop',sep = '')] = controlvars[,paste('pop_',r,sep = '')] / controlvars[,paste('pop_all',sep = '')]
  controlvars[,paste('pct_',r,'_gq',sep = '')] = 1 - controlvars[,paste('pop_',r,'_gq_hh',sep = '')] / controlvars[,paste('pop_',r,sep = '')]
  # controlvars[,paste('pct_',r,'_gq_inst',sep = '')] = controlvars[,paste('pop_',r,'_gq_inst',sep = '')] / controlvars[,paste('pop_',r,sep = '')]
  controlvars[,paste('pct_',r,'_col_ba',sep = '')] = controlvars[,paste('pop_',r,'_col_ba',sep = '')] / controlvars[,paste('pop_',r,sep = '')]
  controlvars[,paste('pct_',r,'_mig_samest',sep = '')] = controlvars[,paste('pop_',r,'_mig_samest',sep = '')] / controlvars[,paste('pop_',r,sep = '')]
  controlvars[,paste('pct_',r,'_mig_difst',sep = '')] = controlvars[,paste('pop_',r,'_mig_difst',sep = '')] / controlvars[,paste('pop_',r,sep = '')]
  controlvars[,paste('pct_',r,'_cit_imm',sep = '')] = 1 - (controlvars[,paste('pop_',r,'_cit_usborn',sep = '')] / controlvars[,paste('pop_',r,sep = '')])
  # controlvars[,paste('pct_',r,'_emp',sep = '')] = controlvars[,paste('pop_',r,'_emp',sep = '')] / controlvars[,paste('pop_',r,sep = '')]
  # controlvars[,paste('inc_trans_pct_',r,sep = '')] = 1 - (controlvars[,paste('mean_',r,'_notrans',sep = '')] / controlvars[,paste('mean_',r,sep = '')])
}


# Difference in college 
controlvars$pct_ba_dif_bw = controlvars$pct_w_col_ba - controlvars$pct_b_col_ba
controlvars[!is.na(controlvars$pct_b_col_ba) & controlvars$pct_b_col_ba == 0, 'pct_b_col_ba'] = 0.0003616
controlvars$pct_ba_lograt_bw = log(controlvars$pct_w_col_ba / (controlvars$pct_b_col_ba) )


# # Get statewide institutionalization
# mwinst_st = filter(mw, incsamp == 'gq_inst_faminc_sqrt') %>% select(year,racedef,st,pop_all,pop_b,pop_o,pop_w) %>% group_by(year,racedef,st) %>% summarise(
#   across(starts_with('pop_'), sum)
# )
mwtot_st = filter(mw, incsamp == 'all_faminc_sqrt') %>% select(year,racedef,st,pop_all,pop_b,pop_o,pop_w) %>% group_by(year,racedef,st) %>% summarise(
  across(starts_with('pop_'), sum)
)
#   
# inst_st = left_join(mwtot_st, mwinst_st, by = c('year','racedef','st'), suffix = c('_tot','_inst'))
# inst_st = inst_st %>% mutate(
#   instst_all = pop_all_inst / pop_all_tot,
#   instst_b = pop_b_inst / pop_b_tot,
#   instst_o = pop_o_inst / pop_o_tot,
#   instst_w = pop_w_inst / pop_w_tot,
# ) %>% ungroup()
#   
# # This doesn't look exatly right - not much total change over time, rate today lower than in 1970 in california
# # Might be getting switch from mental institutions to prison, which census groups
# # Confirmed by https://www.journals.uchicago.edu/doi/full/10.1086/667773 that total institutionalization rate was comparable in 1960s and 70s to today but mostly mental institutions rather than jails 
# # ggplot(inst_st %>% filter(msadef == 'msa3', racedef == 'expansive'), aes(x = year, y = instst_all, alpha = pop_all_tot, group = st)) + geom_line()
# 
# inst_st %>% select(year,racedef,st,starts_with('instst_'))
# 
# inst_st$log_instst_bw = log(inst_st$instst_b / inst_st$instst_w)
# inst_st[inst_st$log_instst_bw == -Inf,'log_instst_bw'] = log(.1) # If the true BW ratio is 0 (no black incarceration, some white), treat as 10% (lowest otherwise is 0.42). Happens in several states in 1969, in Vermont as late as 1989
# controlvars = left_join(controlvars,inst_st, by = c('year','racedef','st'),suffix = c('','_gq_inst')) 


# Get state incarceration rate for 1969 and 1979
mwincar_st = filter(mw, incsamp == 'gq_incarc_faminc_sqrt') %>% select(year,racedef,st,pop_all,pop_b,pop_o,pop_w) %>% group_by(year,racedef,st) %>% summarise(
  across(starts_with('pop_'), sum)
)

cenincar_st = left_join(mwtot_st %>% filter(year %in% c(1969,1979)), mwincar_st, by = c('year','racedef','st'), suffix = c('_tot','_incar'))
cenincar_st = cenincar_st %>% mutate(
  cenincarst_all = pop_all_incar / pop_all_tot * 100000,
  cenincarst_b = pop_b_incar / pop_b_tot * 100000,
  cenincarst_o = pop_o_incar / pop_o_tot * 100000,
  cenincarst_w = pop_w_incar / pop_w_tot * 100000,
) %>% ungroup() # Do rates per 100000 for comparison with Vera

cenincar_st %>% select(year,racedef,st,starts_with('cenincarst_'))

cenincar_st$incar_cen_lrat_bw = log(cenincar_st$cenincarst_b / cenincar_st$cenincarst_w)
filter(cenincar_st, incar_cen_lrat_bw != -Inf) %>% select(incar_cen_lrat_bw) %>% exp() %>% summary() # Minimum is 0.2785
cenincar_st[cenincar_st$incar_cen_lrat_bw == -Inf,'incar_cen_lrat_bw'] = log(.1) # If the true BW ratio is 0 (no black incarceration, some white), treat as 10% (lowest otherwise is 0.2785)
controlvars = left_join(controlvars,cenincar_st, by = c('year','racedef','st'),suffix = c('','_gq_cenincar')) 


## Merge on variables from other sources

# Incarceration from vera
incar_vera = read_csv('raw_data/incarceration/incarceration_trends_state.csv')
incar_vera %>% filter(year == 1969) %>% select(state_name,total_incarceration_rate)

# Use 1970 (not missing much) instead of 1969 (missing lots)
incar_vera = incar_vera %>% filter(year >= 1970)
incar_vera[incar_vera$year == 1970,'year'] = 1969
incar_vera_total = incar_vera %>% select(st = state_fips, year,
                                   total_incarceration_rate)

# Get race-specific rates for 1990-onward (including assigning 1990 value to 1989)
incar_vera_race = incar_vera %>% filter(year >= 1990) %>% select(st = state_fips, year, black_prison_pop_rate, white_prison_pop_rate, total_prison_pop_rate, total_prison_jur_pop_rate, black_jail_pop_rate, white_jail_pop_rate,total_jail_pop_rate)
incar_vera_race[incar_vera_race$year == 1990,'year'] = 1989
incar_vera_race = incar_vera_race %>% mutate(
  black_incar_pop_rate = black_prison_pop_rate + black_jail_pop_rate,
  white_incar_pop_rate = white_prison_pop_rate + white_jail_pop_rate,
  total_incar_pop_rate = total_prison_pop_rate + total_jail_pop_rate,
  
)

incar_vera_race$imprate_vera_lrat_bw = log(incar_vera_race$black_prison_pop_rate / incar_vera_race$white_prison_pop_rate)
incar_vera_race$incrate_vera_lrat_bw = log(incar_vera_race$black_incar_pop_rate / incar_vera_race$white_incar_pop_rate)

controlvars = left_join(controlvars, incar_vera_total, by = c('st','year'))
controlvars = left_join(controlvars, incar_vera_race, by = c('st','year'))

# Make combined incarceration rate variable = Census incarceration data for 1969 and 1979, Vera imprisonment data for remaining years
controlvars = controlvars %>% mutate(
  incar_comb_total = case_when(
    year %in% c(1969,1979) ~ cenincarst_all,
    year %in% c(1989,1999,2007, 2012, 2017) ~ total_incar_pop_rate
  ),
  incar_comb_bwlrat = case_when(
    year %in% c(1969,1979) ~ incar_cen_lrat_bw,
    year %in% c(1989,1999,2007, 2012, 2017) ~ incrate_vera_lrat_bw
  ) ,
  incar_comb_b = case_when(
    year %in% c(1969,1979) ~ cenincarst_b,
    year %in% c(1989,1999,2007, 2012, 2017) ~ black_incar_pop_rate
  ),
  incar_comb_w = case_when(
    year %in% c(1969,1979) ~ cenincarst_w,
    year %in% c(1989,1999,2007, 2012, 2017) ~ white_incar_pop_rate
  ),
  imp_comb_total = case_when(
    year %in% c(1969,1979) ~ cenincarst_all,
    year %in% c(1989,1999,2007, 2012, 2017) ~ total_prison_pop_rate
  ),
  imp_comb_bwlrat = case_when(
    year %in% c(1969,1979) ~ incar_cen_lrat_bw,
    year %in% c(1989,1999,2007, 2012, 2017) ~ imprate_vera_lrat_bw
  ) ,
  imp_comb_b = case_when(
    year %in% c(1969,1979) ~ cenincarst_b,
    year %in% c(1989,1999,2007, 2012, 2017) ~ black_prison_pop_rate
  ),
  imp_comb_w = case_when(
    year %in% c(1969,1979) ~ cenincarst_w,
    year %in% c(1989,1999,2007, 2012, 2017) ~ white_prison_pop_rate
  ),
  incar_comb_bwdif = incar_comb_b - incar_comb_w
)


# Employment by industry
cen_emp = read_csv('intermediate_data/predictors/msa_census_emp_all.csv')
controlvars = left_join(controlvars, cen_emp, by = c('msa','year'))
stopifnot(!is.na(controlvars$cemp_pct_public))
stopifnot(!is.na(controlvars$cemp_tot))


# Segregation
seg = read_csv('intermediate_data/predictors/msa_segregation_nhgis.csv') %>% select(year,msa,dissimilarity_index,iso_black,exposure_b_w,theil)
controlvars = left_join(controlvars,seg, by = c('year','msa'))
filter(controlvars, is.na(dissimilarity_index)) %>% count(year) # Lots of places missing segregation in 1969 because not tracted yet

# Segregation from LTDB consistent tracts
seg_ltdb = read_csv('intermediate_data/predictors/msa_segregation_ltdb.csv') %>% select(year,msa,dissimilarity_index,iso_black,exposure_b_w,theil)

# Map years
controlvars = controlvars %>% mutate(yearm = case_when(
    year %in% c(2007,2012) ~ 2009,
    year == 2017 ~ 2019,
    .default = year
  ))
  
controlvars = left_join(controlvars,seg_ltdb, by = c('yearm' = 'year','msa'),suffix = c('','_ltdb')) %>% select(-yearm)
filter(controlvars, is.na(dissimilarity_index_ltdb)) %>% count(year) #1969 and 1979 missing some

# Unionization rate

# Statewide rate for 1969 and 1979
unions = read_csv('intermediate_data/predictors/union_statewide.csv')
controlvars = left_join(controlvars, unions, by = c('year','st' = 'fips'))
stopifnot(!is.na(controlvars$unionrate))
controlvars = rename(controlvars, union_statewide = unionrate)

# CPS for 1989 on
unions_50one_met = read_csv('intermediate_data/predictors/union_50oneyr_cbsa.csv')
unions_50one_met$METFIPS = as.character(unions_50one_met$METFIPS)
unions_50one_st = read_csv('intermediate_data/predictors/union_50oneyr_st.csv')

controlvars = left_join(controlvars, unions_50one_met %>% select(YEAR,METFIPS,union_cps_50one_met = union_pct), by = c('year' = "YEAR",'msa' = 'METFIPS'))
controlvars = left_join(controlvars, unions_50one_st %>% select(YEAR,st,union_cps_50one_st = union_pct), by = c('year' = "YEAR",'st'))

# Combine state and metro
controlvars = mutate(controlvars, union_cps_50one = coalesce(union_cps_50one_met,union_cps_50one_st))
controlvars %>% filter(is.na(union_cps_50one)) %>% count(year) 
stopifnot(!is.na(controlvars$union_cps_50one) | controlvars$year < 1989 ) # Only missing are first two decades before union was in CPS

controlvars %>% group_by(year) %>% summarise(popmi = sum(pop_all * is.na(union_cps_50one_met)), pop = sum(pop_all),pctmiss = popmi / pop) # Have metro level for about 50% of population

controlvars = controlvars %>% select(-union_cps_50one_met,-union_cps_50one_st) # Remove extra variables

# Make combined
controlvars = mutate(controlvars, 
                     union_combined = coalesce(union_cps_50one,union_statewide))

stopifnot(!is.na(controlvars$union_combined))


# Racial oppression
racialoppr = read_csv('intermediate_data/predictors/racialoppr_msa_wfactor.csv') %>% select(-st,-region)
controlvars = left_join(controlvars, racialoppr, by = c('msa'))

# Make flag for HRR observed
controlvars$hrrvalid = (controlvars$region == "South" | controlvars$st %in% c('29','54')) * 1
controlvars$hrrinvalid = 1- controlvars$hrrvalid
controlvars$hrrv = controlvars$hrr * controlvars$hrrvalid

# No MSAs are missing from racial oppression data
not_in_ro = filter(controlvars, is.na(pop_tot_1860  & is.na(cropfarms_pctblack))) %>% select(year,msa,msaname,shortname)
stopifnot(nrow(not_in_ro) == 0)

# State HRR factor
racialoppr_st = read_csv('intermediate_data/predictors/racialoppr_st_wfactor.csv') %>% select(st, slave_pct_st = slave_pct, cropfarms_pctblack_st = cropfarms_pctblack, hrr_st, hrr4_st)
controlvars = left_join(controlvars,racialoppr_st,by = 'st')

stopifnot(!is.na(controlvars$hrr_st) | controlvars$st %in% c('02','15')) # Alaska and Hawaii are missing, nothing else. Those should just get the minimum
controlvars = mutate(controlvars, 
                     hrr_st = case_when(!(st %in% c('02','15')) ~ hrr_st,
                                        st %in% c('02','15') ~ -0.8331151),
                     hrr4_st =  case_when(!(st %in% c('02','15')) ~ hrr4_st,
                                          st %in% c('02','15') ~ -1.103453))





# Clean up and export 
controlvars = select(controlvars,-msaname,-shortname, -incsamp) 
controlvars = rename(controlvars, pop_b_full = pop_b, pop_w_full = pop_w,pop_o_full = pop_o, pop_all_full = pop_all) 

# Make pop10 variable
controlvars$logpop10 = log(controlvars$pop_all_full,base = 10)

write_csv(controlvars, 'intermediate_data/regression_data/msa_reg_predictors.csv')

