# Compute segregation with LTDB data


# Bring in MSAs
msas = read.csv('raw_data/msas/list1_2020.csv')
msas$msa = str_pad(msas$CBSA.Code,width = 5,pad = '0')
msas$mettype = str_replace(msas$Metropolitan.Micropolitan.Statistical.Area,pattern = ' Statistical Area','')

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 = filter(msas, mettype != '', st != '72')
msas = select(msas, st,cnty,mettype, msa)


# Load LTDB

d70 = read_csv('raw_data/ltdb_std_all_fullcount/LTDB_Std_1970_fullcount.csv') %>% mutate(
  pop_o = POP70 - WHITE70 - BLACK70) %>% select(
  fips = TRTID10, pop_w = WHITE70, pop_b = BLACK70,pop_o
) %>% mutate(fips = str_pad(fips, 11,pad = '0'), st = substr(fips,1,2), cnty = substr(fips, 3,5), year = 1969)

d80 = read_csv('raw_data/ltdb_std_all_fullcount/LTDB_Std_1980_fullcount.csv') %>% mutate(
  pop_o = POP80 - NHWHT80 - NHBLK80) %>% select(
  fips = TRTID10, pop_w = NHWHT80, pop_b = NHBLK80, pop_o
) %>% mutate(fips = str_pad(fips, 11,pad = '0'), st = substr(fips,1,2), cnty = substr(fips, 3,5), year = 1979)

d90 = read_csv('raw_data/ltdb_std_all_fullcount/LTDB_Std_1990_fullcount.csv') %>% mutate(
  pop_o = POP90 - NHWHT90 - NHBLK90) %>% select(
  fips = TRTID10, pop_w = NHWHT90, pop_b = NHBLK90, pop_o
) %>% mutate(fips = str_pad(fips, 11,pad = '0'), st = substr(fips,1,2), cnty = substr(fips, 3,5), year = 1989)

d00 = read_csv('raw_data/ltdb_std_all_fullcount/LTDB_Std_2000_fullcount.csv') %>% mutate(
  pop_o = POP00 - NHWHT00 - NHBLK00) %>% select(
  fips = TRACTID10, pop_w = NHWHT00, pop_b = NHBLK00, pop_o
) %>% mutate(fips = str_pad(fips, 11,pad = '0'), st = substr(fips,1,2), cnty = substr(fips, 3,5), year = 1999)

d10 = read_csv('raw_data/ltdb_std_all_fullcount/LTDB_Std_2010_fullcount.csv') %>% mutate(
  pop_o = pop10 - nhwht10 - nhblk10) %>% select(
  fips = tractid, pop_w = nhwht10, pop_b = nhblk10, pop_o
) %>% mutate(fips = str_pad(fips, 11,pad = '0'), st = substr(fips,1,2), cnty = substr(fips, 3,5), year = 2009)

d20 = read_csv('raw_data/ltdb_std_all_fullcount/ltdb_std_2020_fullcount.csv') %>% mutate(
  pop_o = pop20 - nhwt20 - nhblk20) %>% select(
  fips = TRTID2010, pop_w = nhwt20, pop_b = nhblk20, pop_o
) %>% mutate(fips = str_pad(fips, 11,pad = '0'), st = substr(fips,1,2), cnty = substr(fips, 3,5), year = 2019)

ltdb = bind_rows(d70,d80,d90,d00,d10,d20)

# Add on MSA
ltdb = left_join(ltdb,msas, by = c('st','cnty'))

ltdb = ltdb %>% 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)



  
  ltdblong = ltdb %>% select(year,msa,fips,pop_b,pop_w) %>% pivot_longer(
    cols = c(pop_w, pop_b),
    names_to = "race",
    values_to = "population"
  ) 
  
  
  ltdblongtot = ltdb %>% select(year,msa,fips,pop_b,pop_w,pop_o) %>% pivot_longer(
    cols = c(pop_w, pop_b,pop_o),
    names_to = "race",
    values_to = "population"
  ) 
  


filter(ltdblong, is.na(population)) # These are actual missings in the data, mostly 2000 Census
filter(ltdblong, is.na(population)) %>% count(year) 

# Replace missings with 0
ltdblong = ltdblong %>% mutate(population = replace_na(population,0))
ltdblongtot = ltdblongtot %>% mutate(population = replace_na(population,0))

dissimilarity_results = ltdblong %>%
  # Group by CBSA
  group_by(msa,year) %>%
  # Calculate multiple segregation indices for each CBSA
  summarise(
    # Dissimilarity Index (evenness)
    dissimilarity_index = if(n_distinct(fips) <= 1) {
      0  # Return 0 if only one tract (no segregation possible)
    } else {
      dissimilarity(
        data = cur_data(),
        group = "race",
        unit = "fips",
        weight = "population"
      )$est  # Extract the estimate from the dissimilarity object
    },
    
    
    exposure_b_w = if(n_distinct(fips) <= 1 ) {
      #   # For single tract, isolation equals the proportion Black in that tract
      sum(population[race == "pop_w"], na.rm = TRUE) / sum(population, na.rm = TRUE)
    } else if(sum(population * (race == 'pop_b')) == 0) {
      0
    } else {
      exposure(data = cur_data(), 
               group = "race",
               unit = "fips",
               weight = "population") %>% 
        filter(of == 'pop_b',to == 'pop_w') %>% select(exposure) %>% unlist()
    },
    
    theil = if(n_distinct(fips) <= 1 ) {
      0
    } else if(sum(population * (race == 'pop_b')) == 0) {
      0
    } else {
      mutual_total(data = cur_data(), 
                   group = "race",
                   unit = "fips",
                   weight = "population") %>% 
        filter(stat == 'H') %>% select(est) %>% unlist()
    },
    
    total_population = sum(population, na.rm = TRUE),
    total_black_pop = sum(population[race == "pop_b"], na.rm = TRUE),
    total_white_pop = sum(population[race == "pop_w"], na.rm = TRUE),
    num_tracts = n_distinct(fips),
  )  %>% ungroup()


iso_results = ltdblongtot %>%
  # Group by CBSA
  group_by(msa,year) %>%
  # Calculate multiple segregation indices for each CBSA
  summarise(
    iso_black =  if(n_distinct(fips) <= 1 ) {
      #   # For single tract, isolation equals the proportion Black in that tract
      sum(population[race == "pop_b"], na.rm = TRUE) / sum(population, na.rm = TRUE)
    } else if(sum(population * (race == 'pop_b')) == 0) {
      0
    } else {
      
      isolation(
        data = cur_data(),
        group = "race",
        unit = "fips",
        weight = "population"
      ) %>% filter(race == 'pop_b') %>% select(isolation) %>% unlist()
    }
  )

seg_results = full_join(dissimilarity_results, iso_results, by = c('msa','year'))
stopifnot(nrow(seg_results) == nrow(dissimilarity_results))


write_csv(seg_results, 'intermediate_data/predictors/msa_segregation_ltdb.csv')

# Check difference
yby = read_csv('intermediate_data/predictors/msa_segregation_nhgis.csv')
ltdbs = read_csv('intermediate_data/predictors/msa_segregation_ltdb.csv')

yby = yby %>% mutate(yearm = case_when(
  year %in% c(2007,2012) ~ 2009,
  year == 2017 ~ 2019,
  .default = year
))

test = left_join(yby,ltdbs, by = c('msa', 'yearm' = 'year'), suffix = c('','_lt'))

cor(test$dissimilarity_index,test$dissimilarity_index_lt, use = 'co') #0.8162797
cor(test$iso_black,test$iso_black_lt,use = 'co') #0.9906358
cor(test$exposure_b_w,test$exposure_b_w_lt,use = 'co') #0.9744922
cor(test$theil,test$theil_lt,use = 'co') #0.9565022


wtd.cor(test$dissimilarity_index,test$dissimilarity_index_lt,test$total_population) #0.9140499
wtd.cor(test$dissimilarity_index,test$dissimilarity_index_lt,test$total_black_pop) #0.9449086

# High correlation for all but dissimilarity, that is higher when weighting by population
# Differences seem to be places with low black populations

test$dif = test$dissimilarity_index - test$dissimilarity_index_lt

test %>% filter(dif > .5) %>% count(year) # All 1969 and 1979
test %>% filter(dif > .15) %>% count(year) # All years 
