# Make Figure 3, scatter plot of change in Black incomes vs national median against change in white incomes vs national median
# Also make maps of Black and white incomes vs national median

dir.create('output/maps/medinc_vs_nat')


# Download county shapes
cntys = counties(resolution = '20m', cb = T) %>% shift_geometry()

# Drop PR
cntys = filter(cntys, STATEFP != '72') #STATEFP != '02', STATEFP != '15',  # Still works here
cntys$fips = paste(cntys$STATEFP,cntys$COUNTYFP,sep = '')
cntys = rename(cntys, st = STATEFP)

# Get crosswalk from counties to MSAs
# Only including metros not micros
msas = read.csv('raw_data/msas/list1_2020.csv')
msas$msa = str_pad(msas$CBSA.Code,width = 5,pad = '0')
msas$st = str_pad(msas$FIPS.State.Code,width = 2, pad = '0')
msas$cnty = str_pad(msas$FIPS.County.Code,width = 3,pad = '0')
msas$fips = paste(msas$st,msas$cnty,sep = '')
msas$mettype = str_replace(msas$Metropolitan.Micropolitan.Statistical.Area,pattern = ' Statistical Area','')
msas = select(msas,fips,mettype,msa)


# Merge
cntys = left_join(cntys, msas, by = 'fips') 


# Assign leftover categories
cntys = cntys %>% mutate(msafull = case_when(mettype == 'Metropolitan' ~ msa,
                                          mettype == 'Micropolitan' ~ paste('nm',st,sep = '_'),
                                          is.na(mettype) ~ paste('nm',st,sep = '_'))) %>% select(-msa) %>% rename(msa = msafull)



# Bring in outcome variables
regout = read_csv('intermediate_data/regression_data/msa_reg_outcomes.csv')
outvars = unique(regout$incsamp)
rdefs = unique(regout$racedef)
regint = filter(regout, incsamp %in% c('all_faminc_sqrt'), racedef == 'expansive')
popcut = 10000


# Get national median
natmed = read_csv('intermediate_data/sumstats/sumstats_nat_raceex.csv') %>% select(-msa)
natmed$year = factor(natmed$year, levels = c(1970,1980,1990,2000,2009,2014,2019), labels = c(1969,1979,1989,1999,2007,2012,2017)) %>% as.character() %>% as.numeric()
natmed$incsamp = paste(natmed$subpop, natmed$incdef, sep = '_')
regint = left_join(regint, natmed %>% filter(race == 'b'), by = c('year','incsamp','racedef'))
regint = left_join(regint, natmed %>% filter(race == 'w'), by = c('year','incsamp','racedef'),suffix = c('_natb','_natw'))

stopifnot(!is.na(regint$med_natb))
stopifnot(!is.na(regint$med_natw))

regint$med_b_natpct = regint$med_b / regint$med_natb
regint$med_w_natpct = regint$med_w / regint$med_natw
regint[regint$med_w_natpct > 2, 'med_w_natpct'] = 2
regint[regint$med_b_natpct > 2, 'med_b_natpct'] = 2


# Get categorical variables
rcuts = c(-1,.7,.9,1.1,1.3,10000)
rlabs = c('Less than 70%','70-89%','90-109%','110%-129%','130% or greater')
regint$b_rat_med_cats = regint$med_b_natpct %>% cut(rcuts,labels = rlabs ,include.lowest = T)
regint$w_rat_med_cats = regint$med_w_natpct %>% cut(rcuts,labels = rlabs ,include.lowest = T)
regint[regint$pop_b < popcut, 'b_rat_med_cats'] = NA
regint[regint$pop_w < popcut, 'w_rat_med_cats'] = NA


# Dissolve
msashape = cntys %>% group_by(msa) %>% summarise()  # Dissolve takes a long time to do


rdef = 'expansive'
    
# Loop through years
  for(yr in c(1969,1979,1989,1999,2007,2012,2017)){

    # Get data
    dtaout = filter(regint, year == yr,  racedef == rdef)
    mdta = left_join(msashape, dtaout, by = 'msa')
    
    
    #Make maps of income categories
    # Black
    pdf(paste('output/maps/medinc_vs_nat/map_medb_cat_',yr,'_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
    gplt = ggplot(mdta, aes(fill = b_rat_med_cats)) +
      geom_sf(size = .1, color = 'black') + 
      scale_fill_brewer(name = 'Black median family income as fraction of national Black median',
                        palette = 'YlGnBu',
                        labels = c(levels(mdta$b_rat_med_cats),paste('Black population < 10,000')),
                        guide = guide_legend(title.position = 'top',title.hjust = .5)) +
      theme(legend.position = 'bottom')+
      theme(
        panel.background = element_blank(),
        panel.grid = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        plot.title = element_text(hjust = .5)
      ) 
    print(gplt) 
    dev.off()
    
    
    # White
    pdf(paste('output/maps/medinc_vs_nat/map_medw_cat_',yr,'_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
    gplt = ggplot(mdta, aes(fill = w_rat_med_cats)) +
      geom_sf(size = .1, color = 'black') + 
      scale_fill_brewer(name = 'White median family income as fraction of national white median',
                        palette = 'YlGnBu',
                        labels = c(levels(mdta$w_rat_med_cats),paste('White population < 10,000')),
                        guide = guide_legend(title.position = 'top',title.hjust = .5)) +
      theme(legend.position = 'bottom')+
      theme(
        panel.background = element_blank(),
        panel.grid = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        plot.title = element_text(hjust = .5)
      ) 
    print(gplt) 
    dev.off()
    
  }



# Now do change
mch = regint %>% filter(year %in% c(1969,2017)) %>% select(year,msa,racedef,pop_b,pop_w,med_b_natpct,med_w_natpct,med_rat,med_b, med_w,med_natb,med_natw) %>%
  pivot_wider(id_cols = c(msa,racedef),names_from = year,values_from = c(pop_b,pop_w,med_b_natpct,med_w_natpct,med_rat,med_b,med_w,med_natb,med_natw))

mch$change_b = mch$med_b_natpct_2017 - mch$med_b_natpct_1969
mch$change_w = mch$med_w_natpct_2017 - mch$med_w_natpct_1969

mch$changep_b = mch$med_b_natpct_2017 / mch$med_b_natpct_1969
mch$changep_w = mch$med_w_natpct_2017 / mch$med_w_natpct_1969
mch$change_rat = mch$med_rat_2017 - mch$med_rat_1969

# For MSAs where med_b_1969 is 0, divide by 100 to avoid divide-by-zero; otherwise use the actual value
mch$grow_b = (mch$med_b_2017 / ifelse(!is.na(mch$med_b_1969) & mch$med_b_1969 == 0, 100, mch$med_b_1969)) - 1
mch$grow_w = (mch$med_w_2017 / mch$med_w_1969) - 1

mch$grownat_b = (mch$med_natb_2017 / mch$med_natb_1969) - 1
mch$grownat_w = (mch$med_natw_2017 / mch$med_natw_1969) - 1

mch$growp_b = mch$grow_b / mch$grownat_b
mch$growp_w = mch$grow_w / mch$grownat_w


# Place next to abluquerque is missing - replace with ABQ which is MSA 10740
filter(mch, msa == '10740')
mch[mch$msa == '24380', 'changep_w'] = 0.948
mch[mch$msa == '24380', 'changep_b'] = 1.13

mch[mch$msa == '24380', 'change_w'] = -0.0530
mch[mch$msa == '24380', 'change_b'] = 0.138
mch[mch$msa == '24380', 'growp_w'] = -0.942
mch[mch$msa == '24380', 'growp_b'] = 1.10
mch[mch$msa == '24380', 'change_rat'] = 0.187
mch[mch$msa == '24380', 'pop_b_1969'] = 6355

popcut = 10000
chcuts = c(-100000, -.15,-.05,.05,.15,10000)
chlabs = c('Decline over 15 pp','Decline 5-15 pp','Minimal change','Increase 5-15 pp','Increase over 15 pp')
mch$changecat_b = mch$change_b %>% cut(chcuts,labels = chlabs ,include.lowest = T)
mch$changecat_w = mch$change_w %>% cut(chcuts,labels = chlabs ,include.lowest = T)
mch[mch$pop_b_1969 < popcut, 'changecat_b'] = NA
mch[mch$pop_w_1969 < popcut, 'changecat_w'] = NA


mdta = left_join(msashape, mch %>% filter(racedef == rdef), by = 'msa')

# Black
pdf(paste('output/maps/medinc_vs_nat/map_paper_change_b_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
gplt = ggplot(mdta, aes(fill = changecat_b)) +
  geom_sf(size = .1, color = 'black') + coord_sf(lims_method = 'geometry_bbox',default_crs = NULL) +
  # facet_wrap(.~race_long,ncol = 2) +
  scale_fill_brewer(palette = 'RdYlBu',
                    name = paste('1969-2017 change in Black median family income vs nation'),
                    labels = c(levels(mdta$changecat_b),paste('Black population < 10,000 in 1969')),
                    guide = guide_legend(title.position = 'top',title.hjust = .5)) +
  # scale_fill_brewer(name = '',palette = 'YlGnBu', labels = c(levels(mdta$bw_rat_med_cats_popcut),paste('Black population < 10,000'))) +
  theme(legend.position = 'bottom')+
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    plot.title = element_text(hjust = .5)
  )
print(gplt)
dev.off()

# White
pdf(paste('output/maps/medinc_vs_nat/map_paper_change_w_',rdef,'_',popcut,'.pdf',sep = ''),height = 5)
gplt = ggplot(mdta, aes(fill = changecat_w)) +
  geom_sf(size = .1, color = 'black') + coord_sf(lims_method = 'geometry_bbox',default_crs = NULL) +
  # facet_wrap(.~race_long,ncol = 2) +
  scale_fill_brewer(palette = 'RdYlBu',
                    name = paste('1969-2017 change in white median family income vs nation'),
                    labels = c(levels(mdta$changecat_w),paste('White population < 10,000 in 1969')),
                    guide = guide_legend(title.position = 'top',title.hjust = .5)) +
  # scale_fill_brewer(name = '',palette = 'YlGnBu', labels = c(levels(mdta$bw_rat_med_cats_popcut),paste('Black population < 10,000'))) +
  theme(legend.position = 'bottom')+
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    plot.title = element_text(hjust = .5)
  )
print(gplt)
dev.off()

## Now make scatterplots for Figure 3 of text
# 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$st = str_extract(msanames$msaname,pattern = ',.*') %>% str_replace(pattern = ', ',replacement = '')
msanames$shortname = as.character(msanames$shortname)
msanames$msaname = as.character(msanames$msaname)
msanames$shortclean = paste(msanames$shortname, msanames$st, sep = ', ')

mch = left_join(mch, msanames, by = 'msa')

mch$better = cut(mch$change_rat, breaks = c(-10000,-.05,.05,10000),labels = c('Decrease over 5%\n(more inequality)','Minimal change','Increase over 5%\n(less inequality)')) 
# mch$better = mch$change_rat >= 0 #cut(mch$change_rat, breaks = c(-10000,-.05,.05,10000),labels = c('Change < -5%','Minimal change','Change > 5%')) 

mchnames = filter(mch, pop_b_2017 > 200000 | shortname == 'Austin')


pdf('output/maps/fig3_sc_grow.pdf',height = 6,width = 9)
gplt = ggplot(mch, aes(x = growp_b, y = growp_w, size = pop_b_2017,col=better)) + 
  geom_point(alpha = .4) + 
  scale_size_continuous(range = c(0,7),labels = comma,name = 'Black population, 2017') + 
  scale_y_continuous(label = percent,limits = c(.75,1.5),breaks = c(.8,1,1.2,1.4),
                     name = '1969-2017 growth in local white median income as fraction\nof growth in national white median income') + #,breaks = c(.7,.8,.9,1,1.1,1.25,1.4)) + 
  scale_color_manual(values = c('#d7191c','#e6e600','#2c7bb6'),name = '1969-2017 change in local\nBlack-white median income ratio')+
  geom_hline(yintercept = 1,lty = 3,alpha = .4) +
  geom_vline(xintercept = 1,lty = 3,alpha = .4) +
  geom_text_repel(data = mchnames, aes(label = shortclean),size = 1.9,color = 'black',seed = 11,min.segment.length = .05,segment.size = .2, box.padding = .2, max.overlaps = 14) + 
  scale_x_continuous(label = percent,limits = c(.4,1.7), breaks = c(.6,1,1.4),name = '1969-2017 growth in local Black median income as fraction\nof growth in national Black median income') #,breaks = c(.5,.7,1,1.4)) 

print(gplt)
dev.off()


pdf('output/maps/fig3_sc_grow_bw.pdf',height = 6,width = 9)
gplt = ggplot(mch, aes(x = growp_b, y = growp_w, size = pop_b_2017,col=better)) + 
  geom_point(alpha = .4) + 
  scale_size_continuous(range = c(0,7),labels = comma,name = 'Black population, 2017') + 
  scale_y_continuous(label = percent,limits = c(.75,1.5),breaks = c(.8,1,1.2,1.4),
                     name = '1969-2017 growth in local white median income as fraction\nof growth in national white median income') + #,breaks = c(.7,.8,.9,1,1.1,1.25,1.4)) + 
  scale_color_manual(values = c('grey20','grey80','grey50'),name = '1969-2017 change in local\nBlack-white median income ratio')+
  geom_hline(yintercept = 1,lty = 3,alpha = .4) +
  geom_vline(xintercept = 1,lty = 3,alpha = .4) +
  geom_text_repel(data = mchnames, aes(label = shortclean),size = 1.9,color = 'black',seed = 11,min.segment.length = .05,segment.size = .2, box.padding = .2, max.overlaps = 14) + 
  scale_x_continuous(label = percent,limits = c(.4,1.7), breaks = c(.6,1,1.4),name = '1969-2017 growth in local Black median income as fraction\nof growth in national Black median income') #,breaks = c(.5,.7,1,1.4)) 

print(gplt)
dev.off()



write_csv(mch,'output/trends/mchange.csv')
mch$bbetter = mch$changep_b > 1
mch$wbetter = mch$changep_w > 1
mch$bothbetter = mch$bbetter == T & mch$wbetter == T
mch %>% ungroup() %>% group_by(better,wbetter) %>% summarise(pb = sum(pop_b_2017))
mch %>% ungroup() %>% group_by(better,bbetter) %>% summarise(pb = sum(pop_b_2017))
mch %>% ungroup() %>% group_by(better,bothbetter) %>% summarise(pb = sum(pop_b_2017))
