# Correlation matrix of predictors

# Bring in outcomes and control variabes
outv = read_csv('intermediate_data/regression_data/msa_reg_outcomes.csv')
controlv = read_csv('intermediate_data/regression_data/msa_reg_predictors.csv') 

# Merge
fulldta = full_join(outv %>% select(-st, -region), controlv, by = c('year','racedef','msa'))
stopifnot(nrow(fulldta) == nrow(outv))
stopifnot(!is.na(fulldta$pop_all_full))


fulldta = fulldta %>% filter(racedef == 'expansive', incsamp == 'all_faminc_sqrt')
fulldta %>% filter(year == 2017) %>% count(msa)

# Log pop - make this here so that the name is formatted to be easy to extract
fulldta$logpop10 = log(fulldta$pop_all_full,base = 10)

# Make south variable
fulldta$south = fulldta$region == 'South'
fulldta$region = factor(fulldta$region, levels = c('Northeast','Midwest','South','West'))


regvars = c('logpop10','pct_b_pop','pct_o_pop','pct_all_cit_imm','pct_all_col_ba','pct_ba_lograt_bw',
'cemp_pct_manuf','cemp_pct_public','incar_comb_total','incar_comb_bwlrat','union_combined','hrrv','dissimilarity_index')



outcor = NULL
for(yr in unique(fulldta$year)){
wcormat = wtd.cor(fulldta %>% filter(year == yr) %>% select(regvars),weight = fulldta %>% filter(year == yr) %>% select(pop_b) %>% unlist())$correlation %>% data.frame()
write_csv(wcormat, paste0('output/regressions/cormat_pop_b_',yr,'.csv'))
  
wcormat$var1 = row.names(wcormat)
wcormat$year = yr
wcormat = wcormat %>% select(year,var1, everything())

outcor = bind_rows(outcor, wcormat)
}


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'
)

outcor = outcor %>% mutate(
  var1 = factor(var1, levels = fullvars, labels = cleanvars)
)

write_csv(outcor, 'output/regressions/cormat_allyears.csv')

outcor %>% filter(logpop10 > .7 | logpop10 < -.7)


