# If they are not yet installed, uncomment this line and run it.
# install.packages(c("tidyverse", "tidycensus", "censusapi",
#   "sf", "tigris", "Amelia", "viridis", "ggpubr", "ggtext"))

# Then load them!

# For data management
library(tidyverse)

# For querying census data
library(tidycensus)
library(censusapi)

# For spatial data management
library(sf)
library(tigris)

# For multiple imputation
library(Amelia)

# For data viz
library(viridis)
library(ggpubr)
library(ggtext)
library(knitr)

# Create a folder to host raw index indicators
dir.create("raw_data")

# Load api key to download data from US Census Bureau
mykey <- "YOUR_API_KEY_HERE"
census_api_key(mykey)

myunits <- expand_grid(
  state = c(state.abb, "DC"),
  year = 2011:2020) %>%
  mutate(id = 1:n())

myvars <- read_csv("indicators.csv")

# Get the many variables applicable
tract_vars <- myvars %>% 
  filter(survey == "acs5" & geography == "tract") %>%
  mutate(code = code %>%  str_extract(".*_[0-9]{3}")) %>%
  select(code, value)
  # Now cut of the 'E' at the end for 'Estimate';
  # tidycensus just wants the code

# Create a folder to hold our acs5 estimates
dir.create("raw_data/acs5")

# Then write a function to procure them
get_data = function(data, variables){
  result <- get_acs(
    # Inserting a distinct year
    year = data$year,
    # and a distinct state each time
    state = data$state,
    # but always grabbing the same set of tract variables from the acs5
    geography = "tract",
    variables = variables,
    survey = "acs5")
  
  result %>%
      # and each time, write it to file 
      # in the format of raw_data/acs5/YEAR_STATE.csv
      saveRDS(paste("raw_data/acs5/", data$year, "_", data$state, ".rds", sep = ""))
}

# Then run the function
myunits %>%
  # For each year-state pair
  split(.$id) %>%
  # Run the following get_acs() function many times,
  map(~get_data(., variables = tract_vars$code), .id = "id")

zip_vars <- myvars %>% 
  filter(survey == "zpb" & geography == "zipcode") %>%
  with(code)

expand_grid(year = as.character(2012:2016),
            var = zip_vars) %>% 
  mutate(id = 1:n()) %>%
  split(.$id) %>%
  map_dfr(~getCensus(
        name = "zbp",
        key = mykey,
        vintage = .$year,
        vars = c("ESTAB", "GEO_ID", "YEAR"),
        region = "zipcode:*",
        naics2012 = .$var))  %>%
  write_rds("raw_data/zip_data.rds")

# Get the population in 2016 for every zipcode tabulation area
data.frame(year = 2011:2020) %>%
  split(.$year) %>%
  map_dfr(~get_acs(
    year = .$year,
    geography = "zcta",
    variables = "B01003_001",
    survey = "acs5"), .id = "year") %>%
  # Get the appropriate zipcode
  mutate(geoid = str_sub(GEOID, start = -5, end = -1)) %>%
  # and simplify columns
  select(year, geoid, pop = estimate) %>%
  # and save to file
  write_rds("raw_data/zip_pop.rds")

read_rds("raw_data/zip_pop.rds") %>% 
  left_join(by = c("year", "geoid"),
            y = read_rds("raw_data/zip_data.rds") %>% 
              select(year = YEAR, geoid = zip_code, 
                     value = ESTAB, variable = NAICS2012) %>%
              mutate(variable = variable %>% dplyr::recode(
                "8131" = "religious_orgs",
                "8132" = "charitable_orgs",
                "8133" = "advocacy_orgs",
                "8134" = "civic_orgs",
                "8139" = "unions")) %>%
              pivot_wider(id_cols = c(year, geoid), names_from = variable, 
                          values_from = value)) %>%
  # Filter to years of data which are available from Zipcode Business Trends API
  filter(year %in% 2012:2016) %>%
  mutate_at(vars(charitable_orgs, advocacy_orgs, civic_orgs, unions, religious_orgs),
            list(~. / pop * 10000)) %>%
  # Now recode the variable numbers into names
  # Set any infinite values to NA too
  mutate_at(vars(charitable_orgs, advocacy_orgs, civic_orgs, unions, religious_orgs),
            list(~if_else(is.nan(.) | is.infinite(.), NA_real_, .))) %>%
  write_rds("raw_data/zip_rates.rds")

# For Zipcode to Tract conversion
get_crosswalk = function(from, to, year){
  # Build URL of HUD crosswalk file
  myurl <- paste("https://www.huduser.gov/portal/datasets/usps/", 
                 toupper(from), "_", toupper(to), "_03", year, ".xlsx", sep = "")
  # Create a temporary file to receive the download
  myfile <- tempfile()
  # Download the file from the URL to the temporary site
  download.file(url = myurl, destfile = myfile)
  # Import the file
  readxl::read_excel(myfile) %>% 
    select(from = 1, to = 2) %>% 
    mutate(year = year,
           type = paste(from, "->", to, sep = "")) %>%
    return()
}

data.frame(year = 2011:2020) %>%
  split(.$year) %>%
  # For each year, run the crosswalk function
  map_dfr(~get_crosswalk(from = "zip", to = "tract", year = .$year)) %>%
  # Save compressed version of file
  saveRDS("raw_data/crosswalk_zip_to_tract.rds")

remove(myfile)

read_rds("raw_data/crosswalk_zip_to_tract.rds") %>%
  # Let's join in the zipcode rates to our crosswalk
  left_join(by = c("from" = "geoid", "year"), 
            y = read_rds("raw_data/zip_rates.rds") %>%
              mutate(year = as.numeric(year))) %>%
  # Then for each year and tract,
  group_by(year, geoid = to) %>%
  # Calculate the average
  summarize_at(
    vars(charitable_orgs, advocacy_orgs, civic_orgs, unions, religious_orgs),
    list(~mean(., na.rm = TRUE))) %>%
  ungroup() %>%
  # Fix any NaN or Inf results to NA
   mutate_at(
    vars(charitable_orgs, advocacy_orgs, civic_orgs, unions, religious_orgs),
    list(~if_else(is.nan(.) | is.infinite(.), NA_real_, as.numeric(.)))) %>%
  saveRDS("raw_data/tract_rates.rds")

# Finally, to get ESRI's data on fraternal ties and political acts, 
# we just borrow it from the original index at the county level.
read_delim("raw_data/Kyne Aldrich Capturing Dataset 2020.tab") %>%
  select(geoid, fraternal, politicalacts) %>%
  mutate(geoid = geoid %>% str_sub(-5, nchar(geoid))) %>%
  write_rds("raw_data/county_esri.rds")

get_format = function(data){
  
  print(data$fileid)
  
  read_rds(data$fileid) %>%
    mutate(GEOID = as.character(GEOID),
           year = str_extract(data$fileid, pattern = "20[0-9]{2}")) %>%
    select(year, geoid = GEOID, variable, estimate) %>%
    pivot_wider(id_cols = c(geoid, year), 
                names_from = variable, values_from = estimate) %>%
    return()
}

data.frame(fileid = dir("raw_data/acs5", full.names = TRUE)) %>%
  split(.$fileid) %>%
  map_dfr(~get_format(.)) %>% 
  # Create a variable for joining in our 2012-2016 variables to the full range
  mutate(yearjoin = case_when(
    year < 2012 ~ 2012,
    year > 2016 ~ 2016,
    TRUE ~ as.numeric(year))) %>%
  # Join in zipcode-averaged rates
  left_join(by = c("geoid", "yearjoin" = "year"),
            y = read_rds("raw_data/tract_rates.rds") %>% 
              filter(year %in% 2012:2016)) %>%
  # Get rid of the now-unnecessary joining variable
  select(-yearjoin) %>%
  # Join in 2 county variables
  mutate(county = str_sub(geoid, 1,5)) %>%
  left_join(by = c("county" = "geoid"), y = read_rds("raw_data/county_esri.rds")) %>% 
  write_rds("raw_data/tract_data.rds")

remove(get_format)

read_rds("raw_data/tract_data.rds") %>% 
  # Calculate the race fractionalization index.
  mutate(race_fractionalization = (1 - (
    (B02001_002 / B01003_001) ^ 2 + # square of % White
      (B02001_003 / B01003_001) ^ 2 + # square of % Black
      (B02001_004 / B01003_001) ^ 2 + # square of % American Indian/Alaskan Native 
      (B02001_005 / B01003_001) ^ 2 + # square of % Asian
      (B02001_006 / B01003_001) ^ 2 + # square of % native Hawaiian/Pacific Islander
      (B02001_007 / B01003_001) ^ 2 + # square of % some other race
      (B02001_008 / B01003_001) ^ 2) # square of % two or more races
  ),
  # Now calculate race similarity index, so that 1 = homophily and 0 = heterophily
  race_similarity = 1 - race_fractionalization) %>%
  # Calculate the ethnicity fractionalization index
  mutate(ethnicity_fractionalization = (1 - (
    (B03001_002 / B01003_001) ^ 2 + # square of % Not Hispanic/Latino
      (B03001_003 / B01003_001) ^ 2)  # square of % Hispanic/Latino
  ),
  ethnicity_similarity = 1 - ethnicity_fractionalization) %>% 
  #Calculate educational similarity
  # get absolute value of difference between college and high school educated persons
  # where a large gap indicates homophily and small gap indicates heterophily
   #"EDUCATIONAL ATTAINMENT FOR THE POPULATION 25 YEARS AND OVER"
  mutate(
    # Calculate % with less than high school education
    #  = 1 minus % with high school degree or more
    less_than_high_school = 1 - (
      (B15002_011 + #High School Graduate, GED, or alternative (Male),
         B15002_012 + #Some college, less than 1 year (Male),
         B15002_013 + #Some college, 1 or more years, no degree (Male),
         B15002_014 + #Associate's degree (Male),
         B15002_015 + #Bachelor's degree (Male),    
         B15002_016 + #Master's degree (Male),
         B15002_017 + #Professional school Degree (Male),
         B15002_018 + #Doctoral Degree (Male),
         # For women   
         B15002_019 + #Population 25 years and over (Female),
         B15002_028 + #High School Graduate, GED, or alternative (Female),
         B15002_029 + #Some college, less than 1 year (Female),
         B15002_030 + #Some college, 1 or more years, no degree (Female),
         B15002_031 + #Associate's degree (Female),
         B15002_032 + #Bachelor's degree (Female),    
         B15002_033 + #Master's degree (Female),
         B15002_034 + #Professional school Degree (Female),
         B15002_035 #Doctoral Degree (Female)
   ) / B15002_001), # total ages 24 or above
    # Calculate % with college education or more
    # sum of all types of education at or above college level, divided by total
    bachelors_or_more = (
      B15002_015 + #Bachelor's degree (Male),    
        B15002_016 + #Master's degree (Male),
        B15002_017 + #Professional school Degree (Male),
        B15002_018 + #Doctoral Degree (Male),
        B15002_032 + #Bachelor's degree (Female),    
        B15002_033 + #Master's degree (Female),
        B15002_034 + #Professional school Degree (Female),
        B15002_035 #Doctoral Degree (Female)
      ) / # doctoral degree
      B15002_001, # total ages 24 or above
    EE = abs(less_than_high_school - bachelors_or_more),
    # now subtract it from 1, where 1 now = heterophily and 0 = homophily
    education_equality = 1-EE) %>%
  # Bonding: Race/income equality
  # Get gini index
  # Gini coefficient (0 = perfect equality to 1 = perfect inequality)
  mutate(income_inequality = B19083_001,
         income_equality = 1 - income_inequality) %>%
  # Bonding: Employment equality 
  # Absolute difference between % of total employed and % of total unemployed labor force 
  mutate(
   percent_employed = (B23025_004 / B23025_003),
   percent_unemployed = (B23025_005 / B23025_003),
   employment_inequality = abs(percent_employed - percent_unemployed),
   # rescale to range from 0 to 1
   employment_inequality = scales::rescale(employment_inequality, to = c(0,1)),
  # Here, like with college,
  # 1 now = heterophily, while 0 = homophily.
  # So, let's rescale it so 1 = homophily
  employment_equality = 1 - employment_inequality) %>%
  # Get county mean for missing data
  group_by(county) %>%
  # If values are missing, impute them with the mean value from that variable
  mutate_at(vars(B19326_001, B19326_002, B19326_005),
            funs(if_else(is.na(.), mean(., na.rm = TRUE), .))) %>%
  ungroup() %>%
  #Calculate gender income similarity
  # Gender Income Similarity
  # Median income in the past 12 months (in 20XX inflation-adjusted dollars)
  # adjusted by Census Bureau to whatever year that census was (eg. 2010, 2011, 2012, etc.)
  # Gender Income Fractionalization (0 = homogeneity, 1 = heterogeneity)
  mutate(
    # Get percentage of total median income that is for men, squared
    MI = (B19326_002 / (B19326_002 + B19326_005)) ^ 2,
    # Get percentage of total median income that is for women, squared
    FI = (B19326_005 / (B19326_002 + B19326_005)) ^ 2,
    # How much space remains? Calculate fractionalization (low = homophily, high = heterophily)
    gender_income_fractionalization = (1 - (MI + FI)),
    # Now, to get a larger spanning variable, rescale.
    # smallest values become 0 (homophily) and 
    # largest values become 1 (heterophily)
    gender_income_fractionalization = (gender_income_fractionalization - min(gender_income_fractionalization, na.rm = TRUE)) / 
      # Divide by range
      (max(gender_income_fractionalization, na.rm = TRUE) - min(gender_income_fractionalization)),
    # Calculate gender income similarity, by subtracting from 1 so that
    # 1 = homophily and 0 = heterophily
    gender_income_similarity = 1 - gender_income_fractionalization
  ) %>%
  
  #Bonding: Language competency 
  # % of total population proficient English speakers
  mutate(lang = ((
      B16004_025 + # "Speak only English",
      B16004_027 + # "Speak Spanish & Speak English 'very well'",
      B16004_032 + # "Speak other Indo-European languages & Speak English 'very well'",
      B16004_037 + # "Speak Asian and Pacific Island languages & Speak English 'very well'",
      B16004_042 # "Speak other languages & Speak English 'very well'"
  ) / B16004_024) # "total 18 to 64 years old",
  ) %>%

  # Bonding: Communication capacity 
  # % of total households with a telephone
  mutate(telephone = ((B25043_003 + # "Owner occupied, with with telephone service available",
                       B25043_012 # "Renter occupied, with telephone service available",
                       ) / 
                        B25043_001) # total households
         ) %>%
  # Bonding: Non-elder population 
  # % of total population below 65 years of age
  mutate(
    # elder population
        elder = (
          # Sum of elder male residents (65+)
          B01001_020 + # Male 65-66
           B01001_021 + # Male 67-69
           B01001_022 + # Male 70-74
           B01001_023 + # Male 75-79
           B01001_024 + # Male 80-84
           B01001_025 + # Male 85+
           # Sum of elder female residents (65+)
           B01001_044 + # Female 65-66
           B01001_045 + # Female 67-69
           B01001_046 + # Female 70-74
           B01001_047 + # Female 75-79
           B01001_048)  # Female 85+
        # Divided by sum of male and female residents
        / (B01001_002 +   # Male
             B01001_026), # Female
        # non-elder population       
        non_elder = 1 - elder) %>%
  # % of Citizens who are of Voting Age
  mutate(voting_age = (
    # = Citizens of Voting Age / (Total Pop - Not Citizen)
      B05003_009 + # "Male, 18 years and over, Native (US Citizen)",
      B05003_011 + # "Male, 18 years and over, Foreign born, Naturalized U.S. citizen",
      B05003_020 + # "Female, 18 years and over, Native (US Citizen)",
      B05003_022   # "Female, 18 years and over, Foreign born, Naturalized U.S. citizen",
    ) / (
      B05001_001 - # "Total Population",
        B05001_006) # "Not a US Citizen"
    ) %>% 
  # % of government works at given level of government
  mutate(local_gov = B08128_006 / B08128_001, # Local govt works / total govt workers
         state_gov = B08128_007 / B08128_001, # State govt workers / total govt workers
         fed_gov = B08128_008 / B08128_001 # federal govt workers/ total govt workers
  ) %>%
  select(geoid, year, county, race_similarity, ethnicity_similarity, education_equality, 
         income_equality, employment_equality, gender_income_similarity, 
         lang, telephone, non_elder, 
         religious_orgs, charitable_orgs, advocacy_orgs, civic_orgs, unions, fraternal, 
         voting_age, local_gov, state_gov, fed_gov, politicalacts) %>%
  write_rds("raw_data/indicators.rds")

# Several variables are actually rates, so to normalize them, we log-transform them.
# First, in the event of zeros, we're going to replace these with a small-but-realistic value,
# since you can't log-transform the value zero.
# We're going to identify two points - the smallest non-zero value and half that value.
# that lower value will be our replacement for zero.
replace_zero = function(x){
  xhat <- x[x != 0]
  upper <- sort(unique(x))[2]
  lower <- upper / 2
  return(lower)
  
}


# Grab a list of all census tracts available
read_rds("raw_data/indicators.rds") %>%
 # Let's impute the county-mean for our rates of organizations, 
  # to try to reduce missing data in a spatially conscious way, 
  # since these are spatially aggregated variables from the zipcode level anyways. 
  group_by(year, county) %>%
  mutate_at(vars(religious_orgs, charitable_orgs, advocacy_orgs, civic_orgs),
            funs(if_else(
              condition = is.na(.),
              true = mean(., na.rm = TRUE),
              false = .))) %>%
  ungroup() %>%
  # recode any Infinite responses as NA
  mutate_at(vars(race_similarity:politicalacts),
            list(~if_else(is.infinite(.) | is.nan(.), NA_real_, .))) %>%
  # Next, a few of these variables are actually rates, so they are naturally right-skewed.
  # In order to be nice and comparable to other indicators, 
  # we're going to log-transform them to make them normally distributed again.
  # To do so, we need to impute 0s with a very small but realistic value, as described above.
  mutate_at(vars(advocacy_orgs, civic_orgs,religious_orgs, unions, 
                 charitable_orgs, fraternal, local_gov, state_gov, fed_gov),
            list(~case_when(. == 0 ~  replace_zero(.),   TRUE ~ .)))  %>%
  # Now take these variables and log-transform them
  mutate_at(vars(charitable_orgs,advocacy_orgs, civic_orgs,religious_orgs, unions, 
                 fraternal, local_gov, state_gov, fed_gov),
            list(~log(.))) %>% 
  # Cap all values above or below the 95% most common values,
  # so as not to unduly bias our estimates of the majority of communities
  mutate_at(vars(race_similarity:politicalacts),
            list(~case_when(
              . < quantile(., probs = 0.025, na.rm = TRUE) ~ quantile(., probs = 0.025, na.rm = TRUE),
              . > quantile(., probs = 0.975, na.rm = TRUE) ~ quantile(., probs = 0.975, na.rm = TRUE),
              TRUE ~ .))) %>%
  # Now rescale ALL variables (with their new polarities) from a min of 0 to a max of 1
  mutate_at(vars(race_similarity:politicalacts),
            list(~scales::rescale(.))) %>%
  saveRDS("raw_data/indicators_rescaled.rds")

# Get variables
vars_bonding <- c("race_similarity", "ethnicity_similarity", "education_equality", 
         "income_equality", "employment_equality", "gender_income_similarity", 
         "lang", "telephone", "non_elder") 
# Get variables
vars_bridging <- c("civic_orgs", "advocacy_orgs",  "charitable_orgs",
                      "religious_orgs", "unions", "fraternal")
# Get variables
vars_linking <- c("voting_age", "local_gov", "state_gov", "fed_gov", "politicalacts")


library(tidyverse)
library(shadowtext)

obslong <- read_rds("raw_data/indicators_rescaled.rds") %>%
  sample_n(size = 10000, replace = FALSE) %>%
  pivot_longer(cols = c(race_similarity:politicalacts), 
               names_to = "variable", values_to = "value") %>%
  mutate(group = case_when(
    variable %in% vars_bonding ~ "bonding",
    variable %in% vars_bridging ~ "bridging",
    variable %in% vars_linking ~ "linking")) %>%
  mutate(variable = variable %>% dplyr::recode_factor(
    # Bonding
    "race_similarity" = "Race\nSimilarity",
    "ethnicity_similarity" = "Ethnicity\nSimilarity",
    "education_equality" = "Educational\nEquality",
    "income_equality" = "Race/Income\nEquality",
    "employment_equality" = "Employment\nEquality",
    "gender_income_similarity" = "Gender Income\nSimilarity",
    "lang" = "Language\nCompetency",
    "telephone" = "Communication\nCapacity",
    "non_elder" = "Non-Elder\nPopulation",
    # Bridging
    "religious_orgs" = "Religious\nOrgs", 
    "civic_orgs" = "Civic\nOrgs", 
    "advocacy_orgs" = "Social\nAdvocacy Orgs",
    "charitable_orgs" = "Charitable\nOrgs", 
    "fraternal" = "Fraternal\nOrders", 
    "unions" = "Unions",
    # Linking
    "voting_age" = "Political\nLinkage", 
    "local_gov" = "Local Govt.\nEmployees",
    "state_gov" = "State Govt.\nEmployees", 
    "fed_gov" = "Fed. Govt.\nEmployees",
    "politicalacts" = "Political\nActivities") %>%
      factor(levels = levels(.)))


# Get variables of interest

mylabels <- obslong %>%
  group_by(variable, group) %>%
  summarize(y = 0.5,
            x = min(value, na.rm = TRUE)) 



g1 <- ggplot(data = obslong) +
  geom_density(mapping = aes(x = value, y= ..scaled.., fill = group, group = variable)) +
  geom_shadowtext(data = mylabels, mapping = aes(x = x, y = y, sample = NULL, 
                                                 color = group, label = variable),
                  bg.color = "white", color = "black", bg.r = 0.2, hjust = 0, size = 3, fontface = "bold") +
  facet_wrap(~variable, scales = "free") + 
  theme_classic(base_size = 14) +
  theme(panel.border = element_blank()) +
  theme(legend.position = "bottom", strip.text.x = element_blank(),
        plot.caption = element_text(hjust=  0.5),
        plot.subtitle = element_text(hjust = 0.5)) +
  guides(fill = "none", color = "none") +
  scale_fill_manual(values = c("#648FFF",  "#DC267F", "#FFB000")) +
  scale_color_manual(values = c("#648FFF",  "#DC267F", "#FFB000")) +
  labs(x = "Observed Values",
       y = "Frequency (%)",
       subtitle = "Distributions of Indicators",
       caption = "Based on a random sample of 10,000 census-tract years, due to large sample size.") +
  scale_x_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1),
                     labels = c("0", ".25", ".5", ".75", "1")) +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1),
                     labels = c("0", ".25", ".5", ".75", "1"))

ggsave(g1, filename = "images/densityplot.png", dpi = 500, width = 8, height = 6)

include_graphics("images/densityplot.png")

remove(g1, obslong, mylabels, mymissing)

# Get variables
vars_bonding <- c("race_similarity", "ethnicity_similarity", "education_equality", 
         "income_equality", "employment_equality", "gender_income_similarity", 
         "lang", "telephone", "non_elder") 
# Get variables
vars_bridging <- c("religious_orgs", "charitable_orgs", "advocacy_orgs", 
                   "civic_orgs", "unions", "fraternal")

# Get variables
vars_linking <- c("voting_age", "local_gov", "state_gov", "fed_gov", "politicalacts")


# Import indicators
total <- read_rds("raw_data/indicators_rescaled.rds")

# Calculate social capital index using mean
data.frame(
  year = total$year,
  geoid = total$geoid,
  bonding = matrixStats::rowMeans2(x = total[, vars_bonding] %>% as.matrix()),
  bridging = matrixStats::rowMeans2(x = total[, vars_bridging] %>% as.matrix()),
  linking = matrixStats::rowMeans2(x = total[, vars_linking] %>% as.matrix())) %>%
  mutate(social_capital = matrixStats::rowMeans2(.[, c("bonding", "bridging", "linking")] %>% as.matrix())) %>%
  select(year, geoid, social_capital, bonding, bridging, linking) %>%
  saveRDS("raw_data/index_tract_mean.rds")

# Calculate social capital index using medians
data.frame(
  year = total$year,
  geoid = total$geoid,
  bonding = matrixStats::rowMedians(x = total[, vars_bonding] %>% as.matrix()),
  bridging = matrixStats::rowMedians(x = total[, vars_bridging] %>% as.matrix()),
  linking = matrixStats::rowMedians(x = total[, vars_linking] %>% as.matrix())) %>%
  mutate(social_capital = matrixStats::rowMedians(.[, c("bonding", "bridging", "linking")] %>% as.matrix())) %>%
  select(year, geoid, social_capital, bonding, bridging, linking) %>%
  saveRDS("raw_data/index_tract_median.rds")

remove(total)

bind_rows(
  read_rds("raw_data/index_tract_mean.rds") %>% mutate(type = "mean"),
  read_rds("raw_data/index_tract_median.rds") %>% mutate(type = "median")) %>%
  pivot_longer(cols = c(social_capital:linking), names_to = "variable", values_to = "value") %>%
  pivot_wider(id_cols = c(geoid, year, variable), names_from = type, values_from = value) %>%
  group_by(variable) %>%
  summarize(r = cor(mean, median, use = "pairwise.complete.obs")) %>%
  write_rds("raw_data/correlation.rds")

read_rds("raw_data/correlation.rds")

read_rds("raw_data/index_tract_mean.rds") %>% 
  filter(year == 2011) %>%
  select(social_capital:linking) %>%
  summary()

remove(mycor)

rm(list = ls())

# We're going to work straight from the imputed data
library(tidyverse)
library(Amelia)

# Import and restructure imputations as tibble
read_rds("mi/mi.rds")$imputations %>% 
  map_dfr(~as_tibble(.), .id = "imp") %>%
  write_rds("mi/mi_imp.rds")

# Load imputed cells
mymi <- read_rds("mi/mi_imp.rds")

################
# Bonding
################

# Get variables
vars_bonding <- c("race_similarity", "ethnicity_similarity", "education_equality", 
         "income_equality", "employment_equality", "gender_income_similarity", 
         "lang", "telephone", "non_elder") 
# Get stats
mybonding <- tibble(
  mean = mymi[,vars_bonding] %>% as.matrix() %>% matrixStats::rowMeans2(., na.rm = TRUE),
  sd = mymi[,vars_bonding] %>% as.matrix() %>% matrixStats::rowSds(., na.rm = TRUE))



################
# Bridging
################

# Get variables
vars_bridging <- c("religious_orgs", "charitable_orgs", "advocacy_orgs", 
                   "civic_orgs", "unions", "fraternal")

# Get stats
mybridging <- tibble(
  mean = mymi[,vars_bridging] %>% as.matrix() %>% matrixStats::rowMeans2(., na.rm = TRUE),
  sd = mymi[,vars_bridging] %>% as.matrix() %>% matrixStats::rowSds(., na.rm = TRUE))


################
# Linking
################

# Get variables
vars_linking <- c("voting_age", "local_gov", "state_gov", "fed_gov", "politicalacts")

# Get stats
mylinking <- tibble(
  mean = mymi[,vars_linking] %>% as.matrix() %>% matrixStats::rowMeans2(., na.rm = TRUE),
  sd = mymi[,vars_linking] %>% as.matrix() %>% matrixStats::rowSds(., na.rm = TRUE))


# Bind together the results
total <- bind_rows(
  mybonding %>%
    mutate(imp = mymi$imp,
           geoid = mymi$geoid,
           year = mymi$year,
           type = "bonding"),
  mybridging %>%
    mutate(imp = mymi$imp,
           geoid = mymi$geoid,
           year = mymi$year,
           type = "bridging"),
  mylinking %>%
    mutate(imp = mymi$imp,
           geoid = mymi$geoid,
           year = mymi$year,
           type = "linking")) %>%
  # Pivot so that imputations are in columns, not rows
  pivot_wider(id_cols = c(type, geoid, year), names_from = imp, values_from = c(mean, sd))

# Remove the excess data
remove(mymi, mybonding, mybridging, mylinking,
       vars_bonding, vars_bridging, vars_linking)

dir.create("index")

total %>%
  # Now meld the results for the five imputations
  with(
    # Now, we want to generate single bonding, bridging, and linking indices,
    # combining them and incorporating error within and between imputations 
    # by merging them with Rubin's Rules
    Amelia::mi.meld(
      q = select(., mean_imp1:mean_imp5) %>%
        as.matrix(),
      se = select(., sd_imp1:sd_imp5) %>%
        as.matrix(),
      byrow = FALSE) ) %>%
 # Now report them in a final tibble
  with(
    tibble(
      # Add in unique identifiers
      type = total$type,
      geoid = total$geoid,
      year = total$year,
      # And the vectors of melded results
      score = .$q.mi %>% as.vector(),
      se = .$se.mi %>% as.vector()
    )
  )  %>%
  # Now pivot so each index is a column
  pivot_wider(
    id_cols = c(year, geoid),
    names_from = type,
    values_from = c(score, se)) %>%
  rename(bonding = score_bonding,
         bridging = score_bridging,
         linking = score_linking,
         bonding_se = se_bonding,
         bridging_se = se_bridging,
         linking_se = se_linking) %>%
  # Get social capital index
  mutate(social_capital = (bonding + bridging + linking)/3 ) %>%
  select(year, geoid, social_capital, bonding, bridging, linking,
         bonding_se, bridging_se, linking_se) %>%
  # Save to file
  saveRDS("index/index_tract_V3_04_10_2022.rds")

remove(total)

# Delete this hulking dataset of imputations, since the original is stored in mi.rds
unlink("mi/mi_imp.rds")

rm(list = ls())

# For Zipcode to Tract conversion
get_crosswalk = function(from, to, year){
  # Build URL of HUD crosswalk file
  myurl <- paste("https://www.huduser.gov/portal/datasets/usps/", 
                 toupper(from), "_", toupper(to), "_03", year, ".xlsx", sep = "")
  # Create a temporary file to receive the download
  myfile <- tempfile()
  # Download the file from the URL to the temporary site
  download.file(url = myurl, destfile = myfile)
  # Import the file
  readxl::read_excel(myfile) %>% 
    select(from = 1, to = 2) %>% 
    mutate(year = year,
           type = paste(from, "->", to, sep = "")) %>%
    return()
}

data.frame(year = 2011:2020) %>%
  split(.$year) %>%
  # For each year, run the crosswalk function
  map_dfr(~get_crosswalk(from = "tract", to = "zip", year = .$year)) %>%
  # Save compressed version of file
  saveRDS("raw_data/crosswalk_tract_to_zip.rds")

read_rds("raw_data/crosswalk_tract_to_zip.rds") %>% 
  select(from, to, year) %>%
  # Extract the state from each census tract
  mutate(state_code = str_sub(from, 1,2)) %>%
  # Join in the state name for each tract
  left_join(by = "state_code", y = tigris::fips_codes %>%
              select(state, state_code) %>%
              distinct()) %>%
  # Join in the census tract index scores
  left_join(by = c("from" = "geoid", "year"),
            y = read_rds("index/index_tract_V3_04_10_2022.rds")) %>%
  # Get the average score for each zipcode per state-year
  group_by(state, geoid = to, year) %>%
  summarize(social_capital = mean(social_capital, na.rm = TRUE),
            bonding = mean(bonding, na.rm = TRUE),
            bridging = mean(bridging, na.rm = TRUE),
            linking = mean(linking, na.rm = TRUE)) %>%
  ungroup() %>%
  # Some zipcodes could not be estimated, due to their overlapping census tracts 
  # having no having no census data available for measuring indices. We cut them here.
  filter(!is.na(social_capital)) %>%
  write_rds("index/index_zipcode_V3_04_10_2022.rds")

read_rds("raw_data/crosswalk_tract_to_csub.rds") %>% 
  select(from, to) %>%
  # Extract the state from each census tract
  mutate(state_code = str_sub(from, 1,2)) %>%
  # Join in the state name for each tract
  left_join(by = "state_code", y = tigris::fips_codes %>%
              select(state, state_code) %>%
              distinct()) %>%
  # Join in the census tract index scores
  left_join(by = c("from" = "geoid"),
            y = read_rds("index/index_tract_V3_04_10_2022.rds")) %>%
  # Get average score per census county subdivision, per year, per state
  group_by(state, geoid = to, year) %>%
  summarize(social_capital = mean(social_capital, na.rm = TRUE),
            bonding = mean(bonding, na.rm = TRUE),
            bridging = mean(bridging, na.rm = TRUE),
            linking = mean(linking, na.rm = TRUE)) %>%
  ungroup() %>%
  # 10 subdivisions in rural NY end up without a corresponding year or values; 
  # They were not estimatable from current data. we cut these.
  filter(!is.na(year)) %>%
  write_rds("index/index_csub_V3_04_10_2022.rds")

library(tidyverse)

# Export census tracts
read_rds("index/index_tract_V3_04_10_2022.rds") %>% 
  mutate(state_code = str_sub(geoid, 1,2)) %>%
  left_join(by = "state_code", y = tigris::fips_codes %>%
              select(state, state_code) %>%
              distinct()) %>%
  select(state, year, geoid, social_capital:linking_se) %>%
  write_csv("index/index_tracts_V3_04_10_2022.csv")

# Export zipcodes
read_rds("index/index_zipcode_V3_04_10_2022.rds") %>% 
  write_csv("index/index_zipcode_V3_04_10_2022.csv")

# Export county subdivisions
read_rds("index/index_csub_V3_04_10_2022.rds") %>% 
  write_csv("index/index_csub_V3_04_10_2022.csv")