# Calculate summary statistics by metro
# These are the inputs used for most of the remaining analysis

dir.create('intermediate_data/sumstats')



# Set up output data frames
nat_everything = NULL
nat_raceex_everything = NULL
nat_racenr_everything = NULL
nat_raceex_match_everything = NULL
msa_all_everything = NULL
msa_racenr_everything = NULL
msa_raceex_everything = NULL
msa_raceex_match_everything = NULL
msa_raceex_matchg_everything = NULL
region_all_everything = NULL
region_raceex_everything = NULL
region_racenr_everything = NULL

# Get Census region - South vs rest of country
regions = read.csv('raw_data/census_regions/state-geocodes-v2016.csv')
regions$south = factor(regions$Region, levels = c(1,2,3,4), labels = c('Rest of country','Rest of country','South','Rest of country'))
regions$region = factor(regions$Region, levels = c(1,2,3,4), labels = c('Northeast','Midwest','South','West'))
regions$division = factor(regions$Division, levels = 1:9, labels = c('New England','Middle Atlantic','East North Central','West North Central','South Atlantic','East South Central','West South Central','Mountain','Pacific'))
regions$st = str_pad(regions$State..FIPS.,width = 2,pad = '0')
regions = select(regions, st, south,region, division)



# Loop through years and calculate summary statistics 
for(yr in c(1970,1980,1990, 2000,2009,2014,2019)){ 
  
  print(yr)
  load(paste('intermediate_data/ipums/ipums_clean_',yr,'.Rdata',sep = ''))

  # Make age buckets for matching analysis
  dtaexp$agebucket10 = floor(dtaexp$age / 10) * 10 

  # Add on region 
  dtaexp = left_join(dtaexp, regions, by = 'st')
  
  # Income variables to loop through
  everyoneincvars = c('faminc_sqrt') # For most samples just use square root of family income
  empincvars = 'inctot' # For employed sample do individual income
  fullincvars = c('faminc_sqrt','faminc_norm','faminc_tot','hhinc_sqrt','hhinc_norm','hhinc_tot') # For baseline sample use more income definitions
  incvarsdisp = c('faminc_sqrt','spmuinc','spmuinc_sqrt','spmuinc_disp','spmuinc_disp_sqrt','spmuinc_posttx','spmuinc_posttx_sqrt')
  
  # Only do individual incomes for employed + adult
  # Only do adult-normalized for adult
  
  # Loop through subpops
  subpops = c('all','gq_hh','gq_incarc','gq_inst','edu_hs','edu_somecol','edu_ba','edu_over25','mig_samest','mig_difst','cit_usborn','emp_emp','emp_lf','adult') # mig difst col, afact1'mig_imm','emp_emp',,'adult','child','senior'
  
  for(spop in subpops){
    
    print(spop)
    
  # Assign income variables as configured above
  incvars = everyoneincvars
  if(spop %in% c('emp_emp','emp_lf')) incvars = empincvars
  if(spop == 'all') incvars = fullincvars
  if(spop == 'gq_hh' & yr %in% c(2014, 2019) ) incvars = incvarsdisp
   
  # Do the filter just once
  spopdta = dtaexp %>% filter(get(spop) == T) 
  
  if(nrow(spopdta) == 0) next()
  
  
  # Now get national stats for everyone
  nat_all = spopdta %>% summarise(
    across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
    across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
    across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
    across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop'),
    
    # Rank gap
    across(all_of(incvars), ~ {tryCatch({sum(pcwt * (race_white == T ) * (. <= wtd.quantile(. * (race_black_e == T), weights = pcwt * (race_black_e == T),probs = .5 ))) / sum(pcwt * (race_white == T) )}, error = function(e) NA)} , .names = '{.col}_regap'),
    across(all_of(incvars), ~ {tryCatch({sum(pcwt * (race_white == T ) * (. <= wtd.quantile(. * (race_black_n == T), weights = pcwt * (race_black_n == T),probs = .5 ))) / sum(pcwt * (race_white == T) )}, error = function(e) NA)} , .names = '{.col}_rngap'),
    
  ) %>% pivot_longer(everything())
  
  nat_all = mutate(nat_all, 
                   var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                   incdef = str_remove(name,('_.{1,5}$'))
  ) %>% pivot_wider(id_cols = c(incdef),names_from = var,values_from = value)
  
  nat_all$year = yr
  nat_all$msa = 'nat'
  nat_all$race = 'all'
  nat_all$subpop = spop
  
  
  nat_everything = rbind(nat_everything, nat_all)
  
  # Now get national stats by expansive race
  nat_raceex = spopdta %>% group_by(race = race_expansive) %>% summarise(
    across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
    across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
    across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
    across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop')
  ) %>% pivot_longer(-race)
  
  nat_raceex = mutate(nat_raceex, 
                      var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                      incdef = str_remove(name,('_.{1,5}$'))
  ) %>% pivot_wider(id_cols = c(race,incdef),names_from = var,values_from = value)
  
  
  
  nat_raceex$year = yr
  nat_raceex$msa = 'nat'
  nat_raceex$racedef = 'expansive'
  nat_raceex$subpop = spop
  
  nat_raceex_everything = rbind(nat_raceex_everything, nat_raceex)
  
  # Now get national stats by narrow race
  nat_racenr = spopdta %>% group_by(race = race_narrow) %>% summarise(
    across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
    across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
    across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
    across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop')
  ) %>% pivot_longer(-race)
  
  nat_racenr = mutate(nat_racenr, 
                      var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                      incdef = str_remove(name,('_.{1,5}$'))
  ) %>% pivot_wider(id_cols = c(race,incdef),names_from = var,values_from = value)
  
  nat_racenr$year = yr
  nat_racenr$msa = 'nat'
  nat_racenr$racedef = 'narrow'
  nat_racenr$subpop = spop
  
  nat_racenr_everything = rbind(nat_racenr_everything, nat_racenr)
  
  
  # Do by region as well
  # Now do by region
  region_all = spopdta %>% group_by(south) %>% summarise(
    across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
    across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
    across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
    across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop'),
    
    # Rank gap
    across(all_of(incvars), ~ {tryCatch({sum(pcwt * (race_white == T ) * (. <= wtd.quantile(. * (race_black_e == T), weights = pcwt * (race_black_e == T),probs = .5 ))) / sum(pcwt * (race_white == T) )}, error = function(e) NA)} , .names = '{.col}_regap'),
    across(all_of(incvars), ~ {tryCatch({sum(pcwt * (race_white == T ) * (. <= wtd.quantile(. * (race_black_n == T), weights = pcwt * (race_black_n == T),probs = .5 ))) / sum(pcwt * (race_white == T) )}, error = function(e) NA)} , .names = '{.col}_rngap'),
    
  )  %>% pivot_longer(cols = c(-south))
  
  region_all = mutate(region_all, 
                      var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                      incdef = str_remove(name,('_.{1,5}$'))
  ) %>% pivot_wider(id_cols = c(south,incdef),names_from = var,values_from = value)
  
  region_all$year = yr
  region_all$race = 'all'
  region_all$subpop = spop
  
  
  
  region_all_everything = rbind(region_all_everything, region_all)
  
  
  # Race expansive
  region_raceex = spopdta %>% group_by(south,race = race_expansive) %>% summarise(
    across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
    across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
    across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
    across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop')
  )  %>% pivot_longer(cols = c(-south,-race)) 
  
  region_raceex = mutate(region_raceex, 
                         var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                         incdef = str_remove(name,('_.{1,5}$'))
  ) %>% pivot_wider(id_cols = c(south,race,incdef),names_from = var,values_from = value)
  
  region_raceex$year = yr
  region_raceex$racedef = 'expansive'
  region_raceex$subpop = spop
  
  region_raceex_everything = rbind(region_raceex_everything, region_raceex)
  
  # Race narrow
  region_racenr = spopdta %>% group_by(south,race = race_narrow) %>% summarise(
    across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
    across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
    across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
    across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
    across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop')
  )  %>% pivot_longer(cols = c(-south,-race))
  
  region_racenr = mutate(region_racenr, 
                         var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                         incdef = str_remove(name,('_.{1,5}$'))
  ) %>% pivot_wider(id_cols = c(south,race,incdef),names_from = var,values_from = value)
  
  region_racenr$year = yr
  region_racenr$racedef = 'narrow'
  region_racenr$subpop = spop
  
  region_racenr_everything = rbind(region_racenr_everything, region_racenr)
  
  
      
  # Do by MSA

    msa_all = spopdta %>% group_by(msa) %>% summarise(
      across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
      across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
      across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
      across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop'),
      
      # Calculate rank gap = average rank of the median Black person within the white distribution
      across(all_of(incvars), ~ {tryCatch({sum(pcwt * (race_white == T ) * (. <= wtd.quantile(. * (race_black_e == T), weights = pcwt * (race_black_e == T),probs = .5 ))) / sum(pcwt * (race_white == T) )}, error = function(e) NA)} , .names = '{.col}_regap'),
      across(all_of(incvars), ~ {tryCatch({sum(pcwt * (race_white == T ) * (. <= wtd.quantile(. * (race_black_n == T), weights = pcwt * (race_black_n == T),probs = .5 ))) / sum(pcwt * (race_white == T) )}, error = function(e) NA)} , .names = '{.col}_rngap'),
      
      region = names(which.max(table(region))),
      division = names(which.max(table(division)))
    )  %>% pivot_longer(cols = c(-msa,-region,-division))

    msa_all = mutate(msa_all,
                     var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                     incdef = str_remove(name,('_.{1,5}$'))
    ) %>% pivot_wider(id_cols = c(msa,region,division,incdef),names_from = var,values_from = value)

    msa_all$year = yr
    msa_all$race = 'all'
    msa_all$subpop = spop

    

    # Combine into summary stats
    msa_all_everything = rbind(msa_all_everything, msa_all)


    # Race expansive
    msa_raceex = spopdta %>% group_by(msa,race = race_expansive) %>% summarise(
      across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
      across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
      across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
      across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop'),
      region = names(which.max(table(region))),
      division = names(which.max(table(division)))
    )  %>% pivot_longer(cols = c(-msa,-region,-division,-race))

    msa_raceex = mutate(msa_raceex,
                        var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                        incdef = str_remove(name,('_.{1,5}$'))
    ) %>% pivot_wider(id_cols = c(msa,region,division,race,incdef),names_from = var,values_from = value)

    msa_raceex$year = yr
    msa_raceex$racedef = 'expansive'
    msa_raceex$subpop = spop

    msa_raceex_everything = rbind(msa_raceex_everything, msa_raceex)

    # Race narrow
    msa_racenr = spopdta %>% group_by(msa,race = race_narrow) %>% summarise(
      across(all_of(incvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
      across(all_of(incvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
      across(all_of(incvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
    across(all_of(incvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
    across(all_of(incvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop'),
    region = names(which.max(table(region))),
    division = names(which.max(table(division)))
    )  %>% pivot_longer(cols = c(-msa,-region,-division,-race))

msa_racenr = mutate(msa_racenr,
                     var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                     incdef = str_remove(name,('_.{1,5}$'))
) %>% pivot_wider(id_cols = c(msa,region,division,race,incdef),names_from = var,values_from = value)

  msa_racenr$year = yr
  msa_racenr$racedef = 'narrow'
  msa_racenr$subpop = spop
  msa_racenr_everything = rbind(msa_racenr_everything, msa_racenr)

  # Matching based on sex and age
  # Only doing this for baseline
  if(spop == 'all'){
  msa_raceex_matchg = spopdta %>% filter(pcwt != 0) %>% group_by(msa,race = race_expansive,agebucket10,sex) %>% summarise(
    across(all_of(everyoneincvars), ~ wtd.mean(.,weights = pcwt), .names = '{.col}_mean'),
    across(all_of(everyoneincvars), ~ wtd.quantile(.,weights = pcwt, probs = .5), .names = '{.col}_med'),
    across(all_of(everyoneincvars), ~ wtd.quantile(.,weights = pcwt, probs = .05), .names = '{.col}_p05'),
    across(all_of(everyoneincvars), ~ wtd.quantile(.,weights = pcwt, probs = .1), .names = '{.col}_p10'),
    across(all_of(everyoneincvars), ~ wtd.quantile(.,weights = pcwt, probs = .25), .names = '{.col}_p25'),
    across(all_of(everyoneincvars), ~ wtd.quantile(.,weights = pcwt, probs = .75), .names = '{.col}_p75'),
    across(all_of(everyoneincvars), ~ wtd.quantile(.,weights = pcwt, probs = .9), .names = '{.col}_p90'),
    across(all_of(everyoneincvars), ~ wtd.quantile(.,weights = pcwt, probs = .95), .names = '{.col}_p95'),
    across(all_of(everyoneincvars), ~ sum(!is.na(.)), .names = '{.col}_obs'),
    across(all_of(everyoneincvars), ~ sum(afactor * !is.na(.)), .names = '{.col}_obswt'),
    across(all_of(everyoneincvars), ~ sum(pcwt * !is.na(.)), .names = '{.col}_pop'),
    region = names(which.max(table(region))),
    division = names(which.max(table(division)))
  )  %>% pivot_longer(cols = c(-msa,-region,-division,-race,-agebucket10,-sex))
  
  msa_raceex_matchg = mutate(msa_raceex_matchg,
                            var = str_extract(name,('_.{1,5}$')) %>% str_remove('_'),
                            incdef = str_remove(name,('_.{1,5}$'))
  ) %>% pivot_wider(id_cols = c(msa,region,division,race,sex,agebucket10,incdef),names_from = var,values_from = value)
  
  msa_raceex_matchg$year = yr
  msa_raceex_matchg$racedef = 'expansive'
  msa_raceex_matchg$subpop = spop

  msa_raceex_matchg_everything = rbind(msa_raceex_matchg_everything, msa_raceex_matchg)

}


  

  
  
  

  write_csv(nat_everything,'intermediate_data/sumstats/sumstats_nat_all.csv')
  write_csv(nat_raceex_everything,'intermediate_data/sumstats/sumstats_nat_raceex.csv')
  write_csv(nat_racenr_everything,'intermediate_data/sumstats/sumstats_nat_racenr.csv')
  write_csv(region_raceex_everything,'intermediate_data/sumstats/sumstats_region_raceex.csv')
  write_csv(region_racenr_everything,'intermediate_data/sumstats/sumstats_region_racenr.csv')
  write_csv(region_all_everything,'intermediate_data/sumstats/sumstats_region_all.csv')
  write_csv(msa_raceex_everything,'intermediate_data/sumstats/sumstats_msa_raceex.csv')
  write_csv(msa_raceex_matchg_everything,'intermediate_data/sumstats/sumstats_msa_raceex_matchg.csv')
  write_csv(msa_racenr_everything,'intermediate_data/sumstats/sumstats_msa_racenr.csv')
  write_csv(msa_all_everything,'intermediate_data/sumstats/sumstats_msa_all.csv')
  
  
  
  } # subpop
} # year

write_csv(nat_everything,'intermediate_data/sumstats/sumstats_nat_all.csv')
write_csv(nat_raceex_everything,'intermediate_data/sumstats/sumstats_nat_raceex.csv')
write_csv(nat_racenr_everything,'intermediate_data/sumstats/sumstats_nat_racenr.csv')

write_csv(region_raceex_everything,'intermediate_data/sumstats/sumstats_region_raceex.csv')
write_csv(region_racenr_everything,'intermediate_data/sumstats/sumstats_region_racenr.csv')
write_csv(region_all_everything,'intermediate_data/sumstats/sumstats_region_all.csv')

write_csv(msa_raceex_everything,'intermediate_data/sumstats/sumstats_msa_raceex.csv')
write_csv(msa_raceex_matchg_everything,'intermediate_data/sumstats/sumstats_msa_raceex_matchg.csv')
write_csv(msa_racenr_everything,'intermediate_data/sumstats/sumstats_msa_racenr.csv')
write_csv(msa_all_everything,'intermediate_data/sumstats/sumstats_msa_all.csv')

