# Make one clean dataset with all the outcome vars for the regression analysis
# Including matched DVs

dir.create('intermediate_data/regression_data')

# Bring in MSA to state xwalk
stxw = read_csv('raw_data/msa_to_state_xwalk/msa_to_states_xwalk.csv')

# Bring in MSA names
msanames = read.csv('raw_data/msas/msa_names_2020.csv')
msanames$msa = as.character(msanames$CBSA.Code)
msanames = select(msanames, -CBSA.Code, msaname = CBSA.Title)
msanames$shortname = as.character(msanames$shortname)
msanames$msaname = as.character(msanames$msaname)

# Bring in MSA-level data in long form

mdtaall = read_csv('intermediate_data/sumstats/sumstats_msa_all.csv') %>% select(-division, -region, msa) %>% filter(incdef == 'faminc_sqrt' | subpop == 'all' | (incdef == 'inctot' & subpop == 'emp_emp')) 

mdtar = read_csv('intermediate_data/sumstats/sumstats_msa_raceex.csv') %>% select(-division,-region, -racedef, msa)
mdtaf = bind_rows(mdtaall, mdtar) 
mdtaf$racedef = 'expansive'

mdtar1 = read_csv('intermediate_data/sumstats/sumstats_msa_racenr.csv') %>% select(-division,-region, -racedef,  msa)
mdtaf1 = bind_rows(mdtaall, mdtar1) 
mdtaf1$racedef = 'narrow'

mdtaf = bind_rows(mdtaf, mdtaf1)


# Make clean labels etc
mdtaf$incsamp = paste(mdtaf$subpop,mdtaf$incdef, sep = '_')
incsamps = unique(mdtaf$incsamp)

# Convert year to year income was measured = year before for Census, middle year for ACS
mdtaf$year = factor(mdtaf$year, levels = c(1970,1980,1990,2000,2009,2014,2019),
                    labels = c(1969,1979,1989,1999,2007,2012,2017)) %>% as.character() %>% as.numeric()


# Output dataset for predictor variables with just populations
# Note that we fill NAs created by pivot with 0, because total populatin is observed so it's just places that actually do have 0 pop of whatever race

preddata = mdtaf %>% select(year, racedef,incsamp, incdef,subpop,race,msa, pop ) %>% pivot_wider(id_cols = c(year, racedef, incsamp,incdef,subpop, msa),names_from = race,values_from = pop,names_prefix = 'pop_',values_fill = 0)

# Add on state for merging later
preddata = left_join(preddata, stxw, by = c('msa'))
stopifnot(!is.na(preddata$st))
stopifnot(!is.na(preddata$region))

# Add in msa names
preddata = left_join(preddata, msanames, by = c('msa'))
preddata = mutate(preddata, 
               shortname = coalesce(shortname,msaname),
               msaname = coalesce(msaname, msa))

write_csv(preddata, 'intermediate_data/regression_data/msa_pred_rawdata.csv')


# Make the ratios first, then append the matched sample, then make all the other variables 
wdtaf = mdtaf %>% filter(race %in% c('b','w'),incdef == 'faminc_sqrt' | subpop %in% c('all','gq_hh') | (incdef == 'inctot' & subpop == 'emp_emp')) %>% pivot_wider(
  id_cols = c(msa,incsamp,year,racedef),
  names_from = race, values_from = c(pop,mean,med,p05,p10,p25,p75,p90,p95,obs,obswt))


for(basevar in c('mean','med','p05','p10','p25','p75','p90','p95')){
  
  # Replace negatives with zeros
  wdtaf[wdtaf[,paste0(basevar,'_b')] < 0 & !is.na(wdtaf[,paste0(basevar,'_b')]),paste0(basevar,'_b')] = 0
  wdtaf[wdtaf[,paste0(basevar,'_w')] < 0 & !is.na(wdtaf[,paste0(basevar,'_w')]),paste0(basevar,'_w')] = 0
  
  wdtaf[,paste0(basevar,'_rat')] = wdtaf[,paste0(basevar,'_b')] / wdtaf[,paste0(basevar,'_w')]
  wdtaf[,paste0(basevar,'_ldiff')] = log((wdtaf[,paste0(basevar,'_b')] +1) / (wdtaf[,paste0(basevar,'_w')]+1))
  
  # If black population is 0 or white population is 0, say 1
  wdtaf[is.na(wdtaf[,'pop_b']) | is.na(wdtaf[,'pop_w']) , paste0(basevar,'_rat')] = 1
  wdtaf[is.na(wdtaf[,'pop_b']) | is.na(wdtaf[,'pop_w']) , paste0(basevar,'_ldiff')] = 0

  # If white var is 0 and black var is 0, say 1, 0 
  wdtaf[(!is.na(wdtaf$pop_b) & wdtaf[,paste0(basevar,'_b')] == 0) & (!is.na(wdtaf$pop_w) & wdtaf[,paste0(basevar,'_w')] == 0), paste0(basevar,'_rat')] = 1
  wdtaf[(!is.na(wdtaf$pop_b) & wdtaf[,paste0(basevar,'_b')] == 0) & (!is.na(wdtaf$pop_w) & wdtaf[,paste0(basevar,'_w')] == 0), paste0(basevar,'_ldiff')] = 0

  # If white is 0 and black > 0, say 10
  wdtaf[(!is.na(wdtaf$pop_b) & wdtaf[,paste0(basevar,'_b')] > 0) & (!is.na(wdtaf$pop_w) & wdtaf[,paste0(basevar,'_w')] == 0), paste0(basevar,'_rat')] = 10
  wdtaf[(!is.na(wdtaf$pop_b) & wdtaf[,paste0(basevar,'_b')] > 0) & (!is.na(wdtaf$pop_w) & wdtaf[,paste0(basevar,'_w')] == 0), paste0(basevar,'_ldiff')] = log(10)
  
  stopifnot(!is.na(wdtaf[,paste0(basevar,'_rat')]))
  stopifnot(!is.na(wdtaf[,paste0(basevar,'_ldiff')]))
  
}

# Fill in missing populations with 0s
wdtaf[is.na(wdtaf[,'pop_b']) ,'pop_b'] = 0
wdtaf[is.na(wdtaf[,'pop_w']) ,'pop_w'] = 0




# Add on rank gap 
mrank = mdtaf %>% filter(race == 'all', racedef == 'expansive') %>% select(msa,incdef,year, subpop,regap, rngap) %>% 
  pivot_longer(cols = c(regap,rngap), values_to = 'rgap', names_to = 'racedef')
mrank = mrank %>% mutate(
  racedef = factor(racedef, levels = c('regap','rngap'), labels = c('expansive','narrow') %>% as.character()),
  incsamp = paste(subpop,incdef, sep = '_')
) %>% select(-subpop, -incdef)


wdtaf = left_join(wdtaf, mrank, by = c('msa','incsamp','year','racedef'))

stopifnot(!is.na(wdtaf$rgap) | wdtaf$pop_b == 0 | wdtaf$pop_w == 0 | str_detect(wdtaf$incsamp, 'spmuinc') )


### Bring in matched results
matched_b = read_csv('intermediate_data/sumstats/msa_matched_bpop.csv')
matched_b$incsamp = paste(matched_b$subpop,matched_b$incdef,'matched_b',matched_b$matchby,sep = '_')

matched_b = filter(matched_b, incdef %in% c('faminc_sqrt','faminc_tot'), subpop == 'all') %>% select(-matchby,-incdef,-subpop)

matched_w = read_csv('intermediate_data/sumstats/msa_matched_wpop.csv')
matched_w$incsamp = paste(matched_w$subpop,matched_w$incdef,'matched_w',matched_w$matchby,sep = '_')
matched_w = filter(matched_w, incdef %in% c('faminc_sqrt','faminc_tot'), subpop == 'all') %>% select(-matchby,-incdef,-subpop)

matched_b$year = factor(matched_b$year, levels = c(1970,1980,1990,2000,2009,2014,2019),
                   labels = c(1969,1979,1989,1999,2007,2012,2017)) %>% as.character() %>% as.numeric()

matched_w$year = factor(matched_w$year, levels = c(1970,1980,1990,2000,2009,2014,2019),
                   labels = c(1969,1979,1989,1999,2007,2012,2017)) %>% as.character() %>% as.numeric()


wdtaf = bind_rows(wdtaf, matched_b)
wdtaf = bind_rows(wdtaf, matched_w)

stopifnot(!is.na(wdtaf$rgap) | wdtaf$pop_b == 0 | wdtaf$incsamp %in% c('all_faminc_sqrt_matched_b_age_sex','all_faminc_sqrt_matched_w_age_sex','all_faminc_tot_matched_b_age_sex','all_faminc_tot_matched_w_age_sex') | wdtaf$pop_w == 0 | str_detect(wdtaf$incsamp, 'spmuinc') )


# Add in MSA names
wdtaf = left_join(wdtaf, msanames, by = 'msa')
wdtaf = mutate(wdtaf, 
            shortname = coalesce(shortname,msa),
            msaname = coalesce(msaname, msa))


# Add on state for merging later
wdtaf = left_join(wdtaf, stxw, by = c('msa'))
stopifnot(!is.na(wdtaf$st))
stopifnot(!is.na(wdtaf$region))


wdtaf = select(wdtaf, year, incsamp, msa, msaname, shortname, region, st, everything())


# Save output
write_csv(wdtaf, 'intermediate_data/regression_data/msa_reg_outcomes.csv')


