# Take matched results and compute weighted average (by black and by white pop) for each metro

dta = read_csv('intermediate_data/sumstats/sumstats_msa_raceex_matchg.csv')

# Get local ratios for black and white, drop unneeded specifications for speed
dta = filter(dta, incdef %in% c('faminc_sqrt','faminc_tot'), subpop == 'all')
dta = dta %>% filter(race %in% c('b','w')) 


# This sometimes exhausts vector memory limit on laptop
wdta = dta %>% pivot_wider(
  id_cols = c(msa,region,division,agebucket10,sex,incdef,year,racedef,subpop),
  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 0
  wdta[(!is.na(wdta$pop_b) & wdta[,paste0(basevar,'_b')] < 0), paste0(basevar,'_b')] = 0
  wdta[(!is.na(wdta$pop_w) & wdta[,paste0(basevar,'_w')] < 0), paste0(basevar,'_w')] = 0
  
  # Compute ratios
  wdta[,paste0(basevar,'_rat')] = wdta[,paste0(basevar,'_b')] / wdta[,paste0(basevar,'_w')]
  wdta[,paste0(basevar,'_ldiff')] = log((wdta[,paste0(basevar,'_b')] +1) / (wdta[,paste0(basevar,'_w')]+1))
  
  # If black population is 0 or white population is 0, say ratio = 1, logdif = 0
  wdta[is.na(wdta[,'pop_b']) | is.na(wdta[,'pop_w']) , paste0(basevar,'_rat')] = 1
  wdta[is.na(wdta[,'pop_b']) | is.na(wdta[,'pop_w']) , paste0(basevar,'_ldiff')] = 0
  #  # | wdta[,paste0(basevar,'_b')] == 0)  | wdta[,paste0(basevar,'_w')] == 0 
  
  # If white var is 0 and black var is 0, say 1, 0 
  wdta[(!is.na(wdta$pop_b) & wdta[,paste0(basevar,'_b')] == 0) & (!is.na(wdta$pop_w) & wdta[,paste0(basevar,'_w')] == 0), paste0(basevar,'_rat')] = 1
  wdta[(!is.na(wdta$pop_b) & wdta[,paste0(basevar,'_b')] == 0) & (!is.na(wdta$pop_w) & wdta[,paste0(basevar,'_w')] == 0), paste0(basevar,'_ldiff')] = 0
  
  # If white is 0 nad black > 0, say 10
  wdta[(!is.na(wdta$pop_b) & wdta[,paste0(basevar,'_b')] > 0) & (!is.na(wdta$pop_w) & wdta[,paste0(basevar,'_w')] == 0), paste0(basevar,'_rat')] = 10
  wdta[(!is.na(wdta$pop_b) & wdta[,paste0(basevar,'_b')] > 0) & (!is.na(wdta$pop_w) & wdta[,paste0(basevar,'_w')] == 0), paste0(basevar,'_ldiff')] = log(10)
  
  stopifnot(!is.na(wdta[,paste0(basevar,'_rat')]))
  stopifnot(!is.na(wdta[,paste0(basevar,'_ldiff')]))
  # wdta[wdta[,paste0(basevar,'_b')] > 0 & wdta[,paste0(basevar,'_w')] == 0, paste0(basevar,'_rat')] = 10
  
}


# Now aggregate based on black pop
# Missing populations are 0

wdta[is.na(wdta$pop_b),'pop_b'] = 0
wdta[is.na(wdta$pop_w),'pop_w'] = 0


aggbgdta = wdta %>% group_by(msa,incdef,year,racedef,subpop) %>% 
  summarise(#across(ends_with('_rat'),~ wtd.mean(.x,weights = pop_b)),
    across(matches('^(mean|med|p05|p10|p25|p75|p90|p95)'),~ wtd.mean(.x,weights = pop_b)),
    across(matches('^(pop|obs|obswt)'), ~ sum(.x,na.rm = T))
  )

aggwgdta = wdta %>% group_by(msa,incdef,year,racedef,subpop) %>% 
  summarise(#across(ends_with('_rat'),~ wtd.mean(.x,weights = pop_w)),
    across(matches('^(mean|med|p05|p10|p25|p75|p90|p95)'),~ wtd.mean(.x,weights = pop_w)),
    across(matches('^(pop|obs|obswt)'), ~ sum(.x,na.rm = T))
  )    

aggbgdta$matchby = 'age_sex'
aggwgdta$matchby = 'age_sex'


write_csv(aggbgdta,'intermediate_data/sumstats/msa_matched_bpop.csv')
write_csv(aggwgdta,'intermediate_data/sumstats/msa_matched_wpop.csv')


