# Run regressions

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

# Bring in outcomes and control variabes
outv = read_csv('intermediate_data/regression_data/msa_reg_outcomes.csv') %>% select(-st,-region)
controlv = read_csv('intermediate_data/regression_data/msa_reg_predictors.csv')

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

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

# Rename rgap to med_rgap
fulldta$med_rgap = fulldta$rgap

# Make flag for non CBSA
fulldta$state_resid = is.na(as.numeric(fulldta$msa))
fulldta %>% filter(state_resid == T) %>% count(msa) %>% print(n = 100)

# Going to loop through differnet income definitions and sub populations
incsamps = c('all_faminc_sqrt','gq_hh_faminc_sqrt','cit_usborn_faminc_sqrt','mig_samest_faminc_sqrt','mig_difst_faminc_sqrt','emp_emp_inctot',
             'all_hhinc_sqrt','all_faminc_norm','all_faminc_tot',
             'all_faminc_sqrt_matched_b_age_sex','gq_hh_spmuinc_sqrt','gq_hh_spmuinc_disp_sqrt' 
             )
yrs = unique(fulldta$year)
rdefs = unique(fulldta$racedef)

# Right now doing log income, log pop, 

# Models
models = c()
modnames = c()

i = 1



modnames[i] = 'Baseline'
models[i] = 'outvar ~ 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'

i = i + 1

# Alternate income definition - done in loop
# Alternate population subsets - done in loop
# Alternate sample specifications

modnames[i] = 'No population weights'
models[i] = 'outvar ~ 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'

i = i + 1

modnames[i] = 'Full population weights'
models[i] = 'outvar ~ 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'

i = i + 1

modnames[i] = 'Consistent MSA sample'
models[i] = 'outvar ~ 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'

i = i + 1

modnames[i] = 'No non-MSA observations'
models[i] = 'outvar ~ 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'

i = i + 1



# Alternate segregation measures

modnames[i] = 'Isolation index'
models[i] = 'outvar ~ 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 + iso_black'

i = i + 1

modnames[i] = 'Theil index'
models[i] = 'outvar ~ 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 + theil'

i = i + 1

modnames[i] = 'Exposure index'
models[i] = 'outvar ~ 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 + exposure_b_w'

i = i + 1


modnames[i] = 'Theil index (consistent tract boundaries)'
models[i] = 'outvar ~ 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 + theil_ltdb'

i = i + 1


modnames[i] = 'Dissimilarity index (consistent tract boundaries)'
models[i] = 'outvar ~ 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_ltdb'

i = i + 1

modnames[i] = 'Isolation index (consistent tract boundaries)'
models[i] = 'outvar ~ 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 + iso_black_ltdb'

i = i + 1

modnames[i] = 'Exposure index (consistent tract boundaries)'
models[i] = 'outvar ~ 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 + exposure_b_w_ltdb'

i = i + 1


# Alternate specifications - HRR

modnames[i] = 'HRR factor excluded'
models[i] = 'outvar ~ 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 + dissimilarity_index'

i = i + 1


modnames[i] = 'Baker HRR'
models[i] = 'outvar ~ 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 + hrrbaker:hrrvalid + dissimilarity_index'

i = i + 1


modnames[i] = '2-variable state-level HRR'
models[i] = 'outvar ~ 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 + hrr_st:hrrvalid + dissimilarity_index'

i = i + 1

modnames[i] = '4-variable state-level HRR'
models[i] = 'outvar ~ 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 + hrr4_st:hrrvalid + dissimilarity_index'

i = i + 1


modnames[i] = '4-component HRR'
models[i] = 'outvar ~ 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 + hrr4:hrrvalid + dissimilarity_index'

i = i + 1

modnames[i] = 'Slavery only'
models[i] = 'outvar ~ 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 + slave_pct:hrrvalid + dissimilarity_index'

i = i + 1

modnames[i] = 'Sharecropping only'
models[i] = 'outvar ~ 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 + cropfarms_pctblack:hrrvalid + dissimilarity_index'

i = i + 1

modnames[i] = 'Disfranchisement devices only'
models[i] = 'outvar ~ 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 + numdis:hrrvalid + dissimilarity_index'

i = i + 1

modnames[i] = 'Southern Manifesto only'
models[i] = 'outvar ~ 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 + southman:hrrvalid + dissimilarity_index'

i = i + 1



# Other alternative models

# Add region dummy
modnames[i] = 'Region dummies'
models[i] = 'outvar ~ 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 + region + incar_comb_total + incar_comb_bwlrat +
  union_combined + hrrv + dissimilarity_index'

i = i + 1

# No HRR

modnames[i] = 'HRR excluded'
models[i] = 'outvar ~ 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 + dissimilarity_index'

i = i + 1

# No total incarceration rate 
modnames[i] = 'Incarceration inequality only'
models[i] = 'outvar ~ 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_bwlrat +
  union_combined + hrrv + dissimilarity_index'

i = i + 1

# Vera prison only (don't have to drop states)
modnames[i] = 'Imprisonment rate'
models[i] = 'outvar ~ 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 + total_prison_pop_rate + imprate_vera_lrat_bw +
  union_combined + hrrv + dissimilarity_index'

i = i + 1

modnames[i] = 'Incarceration difference'
models[i] = 'outvar ~ 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_bwdif +
  union_combined + hrrv + dissimilarity_index'

i = i + 1


# Drop segregation and incarceration for full sample in each year

modnames[i] = 'Segregation and incarceration excluded\n(full MSA sample in each year)'
models[i] = 'outvar ~ 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 + 
  union_combined + hrrv'

i = i + 1




# Run models 

allcoefsperiod = NULL
allr2s = data.frame()

for(incs  in incsamps){
  
  dir.create(paste('output/regressions/fulloutput/',incs,sep = ''))
  
      # For baseline sample and income definition, also run extended models and percentiles, and expanded race definition
        
      totmods = c(1)
      rdefs = 'expansive'
      
        if(incs %in% c('all_faminc_sqrt')){
          rdefs = c('expansive','narrow')
        }
      if(incs %in% c('all_faminc_sqrt')){
        totmods = 1:length(models)
      }

      for(rdef in rdefs){
        
        if(rdef == 'expansive' & incs == 'all_faminc_sqrt'){
          msalistdta = filter(fulldta, incsamp == incs,  racedef == rdef,pop_b > 0)
          
          segmsas = filter(msalistdta, !is.na(dissimilarity_index), year == 1969) %>% select(msa) %>% unlist()
          incarmsas = filter(msalistdta, !is.na(incar_comb_bwlrat), year == 2017) %>% select(msa) %>% unlist()
          msas_consist = intersect(segmsas,incarmsas)
        }
        
        for(modnum in totmods){

          mod = models[modnum] %>% formula()
          modname = modnames[modnum]
          
          outquants = 'med'
          if(modname == 'Baseline' & incs == 'all_faminc_sqrt') outquants = c('med','mean','p05','p10','p25','p75','p90','p95') #,'mean','p05','p10','p25','p75','p90','p95')
          
          for(outq in outquants){
            
            # Model median white, median black, and difference
            outvs = c('w','b','ldiff')
            
            if(modname == 'Baseline' & incs == 'all_faminc_sqrt' & outq == 'med') outvs = c('w','b','ldiff','rgap')

            for(outthing in outvs){ 
          
                outv = paste(outq,outthing, sep = '_')

                outregs = list()
                outregsr = list()
                outregss = list()
                position = 1
                outcoefs = NULL
                vifs = NULL
            
                # For vera incarceration by race, limit to >1990
                # For dispinc, limit to 2019
                regyrs = yrs
                if(modname %in% c('Imprisonment rate')) regyrs = c(1989,1999,2007,2012,2017)
                if(incs %in% c('gq_hh_spmuinc','gq_hh_spmuinc_disp','gq_hh_spmuinc_sqrt','gq_hh_spmuinc_disp_sqrt','gq_hh_spmuinc_posttx','gq_hh_spmuinc_posttx_sqrt')) regyrs = c(2012,2017)
                
                for(yr in regyrs){
                  
                  rdta = filter(fulldta, incsamp == incs, racedef == rdef,year == yr, pop_b > 0)
                  rdta = rdta %>% rename(outvar = !!sym(outv))
                  
                  # Log the medians
                  if(outv %in% c('med_w','med_b')) rdta$outvar = log(rdta$outvar + 1)

                  # Run regressions
                  lm1 = lm_robust(mod,  
                           data = rdta, weights = pop_b )
                  # Get normalized coefficients
                  lm2 = lm(mod, data = rdta, weights = pop_b ) %>% lm.beta()
                  
                  if(modname == 'Consistent MSA sample') {
                    lm1 = lm_robust(mod, data = rdta %>% filter(msa %in% msas_consist), weights = pop_b )
                    lm2 = lm(mod, data = rdta %>% filter(msa %in% msas_consist), weights = pop_b ) %>% lm.beta()
                  }
                  if(modname == 'No non-MSA observations'){
                    lm1 = lm_robust(mod, data = rdta %>% filter(state_resid == F), weights = pop_b )
                    lm2 = lm(mod, data = rdta %>% filter(state_resid == F), weights = pop_b ) %>% lm.beta()
                    
                  }
                  if(modname == 'Full population weights'){
                    lm1 = lm_robust(mod, data = rdta, weights = pop_all_full )
                    lm2 = lm(mod, data = rdta, weights = pop_all_full ) %>% lm.beta()
                  }
                  if(modname == 'No population weights'){
                    lm1 = lm_robust(mod, data = rdta )
                    lm2 = lm(mod, data = rdta) %>% lm.beta()
                  }
                  outregs[[position]] = lm1
                  position = position + 1
                  
                  # Plan: get everything from robust, but separately estimate normal and get standard coef from there
                  coefgraphs = cbind(names(lm1$coefficients),summary(lm1)$coefficients) %>% data.frame()
                  scoef = summary(lm2)$coefficients[,'Standardized'] 
                  coefgraphs = bind_cols(coefgraphs, scoef)
                  names(coefgraphs) = c('var','coef','se','t','pvalue','cil','cih','df','standardcoef')
                  coefgraphs$coef = as.numeric(coefgraphs$coef)
                  coefgraphs$standardcoef = as.numeric(coefgraphs$standardcoef)
                  coefgraphs$pvalue = as.numeric(coefgraphs$pvalue)
                  coefgraphs$pcat = coefgraphs$pvalue < .05
                  coefgraphs$mod = modname
                  coefgraphs$year = yr
                  coefgraphs$out = outv
                  coefgraphs$racedef = rdef
                  coefgraphs$outvar = outv
                  coefgraphs$incsamp = incs
                  
                  outcoefs = rbind(outcoefs, coefgraphs)
                  allcoefsperiod = bind_rows(allcoefsperiod,coefgraphs)
                  
                  cr2 = summary(lm1)$r.squared
                  cr2a = model.frame(lm1) %>% nrow() 
                  outrd = c(rdef,incs,outv,modname,yr,cr2,cr2a)
                  allr2s = rbind(allr2s, outrd)
                  names(allr2s) = c('racedef','incsamp','outvar','mod','year','r2','nobs')
                  
                  # For baseline only get VIFs
                  if(modname == 'Baseline' & outv == 'med_ldiff' & rdef == 'expansive' & incs == 'all_faminc_sqrt'){
                    vifbase = vif(lm(mod, data = rdta, weights = pop_b ))
                    vifrow = cbind(names(vifbase), vifbase) %>% as.data.frame() 
                    names(vifrow) = c(paste0('var',yr),paste0('vif',yr))
                    vifs = bind_cols(vifs, vifrow)
                    
                    write_csv(vifs, 'output/regressions/vifs_baseline.csv')
                  }
                  
                  
                  } # year
                  
  write_csv(outcoefs,paste('output/regressions/fulloutput/',incs,'/coefs_r_',incs,'_',outv,'_',rdef,'_',modname,'.csv',sep = ''))
  write_csv(allcoefsperiod,paste('output/regressions/all_coefs.csv',sep = ''))
  write_csv(allr2s,paste('output/regressions/all_r2s.csv',sep = ''))
  
} # model
        } #outthing
          
    
    
} # outq
    } # rdef
  } # incsamp


