NIH_impacts.Rmd

---
title: "NIH_impacts"
output: html_document
date: "2025-05-05"
---


# Setup

```{r Packages}

if (!require('rstudioapi')) BiocManager::install("rstudioapi");library(rstudioapi)

# change working directory to folder with script
setwd(dirname(rstudioapi::getSourceEditorContext()$path)) 

library(tidyverse)
library(readxl)
library(sf)
library(magrittr)
library(data.table)

options(scipen=999)
options(datatable.use.logical64 = TRUE)
```

```{r Loading OSF data}
#2024 NIH data
NIH <- read_csv("https://osf.io/4trj7/download")
```

```{r Data Wrangling - for IDCS}

#format column names to match
names(NIH) <- make.names(names(NIH))

#select relevant variables
NIH_clean <- select(NIH, "Application.ID", "Organization.Name", "Project.Title", "Total.Cost", "Direct.Cost.IC", "Indirect.Cost.IC", "Organization.State", "Organization.City", "Latitude", "Longitude", "CFDA.Code", "Project.Number")

#rename columns
colnames(NIH_clean) <- c("ID", "org_name", "title", "total", "direct", "indirect", 
                         "state", "city", "lat", "lon", "CFDA","proj.num")

#exclude rows with NA for CFDA (these indicate sub-grants like "cores" within a center grant)
NIH_clean <- filter(NIH_clean, !is.na(NIH_clean$CFDA))

#replace NA values with zeros and sum directs and indirects
NIH_clean$direct[is.na(NIH_clean$direct)] <- 0
NIH_clean$indirect[is.na(NIH_clean$indirect)] <- 0
NIH_clean$grant_funds <- NIH_clean$direct + NIH_clean$indirect

#fix an issue in the raw data where some grants have total values (direct) but no breakdown by direct/indirect
for (i in 1:nrow(NIH_clean)){
  if (NIH_clean$grant_funds[i] == 0){
    NIH_clean$grant_funds[i] <- NIH_clean$total[i]
    NIH_clean$direct[i] <- NIH_clean$total[i]
  }
}

#deletes first character of project number
NIH_clean$proj.num <- str_sub(NIH_clean$proj.num, 2)

```

```{r Get FIPS Codes}
#merge in FIPS codes (geolocated separately)
NIH_tract_fips <- read.csv("data/tract_dictionary.csv") %>%
            mutate(tractid = str_pad(tractid, 11, pad = "0")) %>%
            mutate(FIPS = substr(tractid, 0, 5)) %>%
            filter(lat!=0 & !is.na(FIPS)) %>%
            select(tractid, FIPS, lat, lon) %>%
            distinct()

NIH_geoid <- read.csv("data/geoid_dictionary_July4.csv") %>%
            select(-GEOID) %>%
            mutate(GEOID = str_pad(New119Fips, 4, pad = "0")) %>%
            filter(lat!=0 & !is.na(GEOID)) %>%
            select(GEOID, NAMELSAD, lat, lon) %>%
            distinct()

NIH_clean_fips <- left_join(NIH_clean, NIH_tract_fips) %>%
                    left_join(NIH_geoid)

NIH_clean_fips <- filter(NIH_clean_fips, !is.na(FIPS))

NIH_clean_fips %>%
  filter(FIPS == "06059")
```

```{r terminated grants}
term_tab_orig <- read_csv("https://osf.io/txmjw/download") %>%
  mutate(award_remaining = as.numeric(str_replace(total_estimated_remaining, "\\$", ""))) %>%
  mutate(org_name = toupper(org_name)) %>%
  mutate(award_remaining_noself = ifelse(cancellation_source == "Self reported", 0 , award_remaining)) %>%
  # exclude reinstated grants
  filter(is.na(reinstated_est_date) & is.na(unfrozen_date)) %>%
    mutate(award_remaining = ifelse(award_remaining <0 , 0 , award_remaining),
           award_remaining_noself = ifelse(award_remaining_noself <0 , 0 , award_remaining_noself)) 

filter(term_tab, org_name== "UNIVERSITY OF CALIROFNIA-IRVINE")

term_tab <- term_tab_orig

# combine subinstitutions
repeated_dict <- read.csv("data/repeated_orgs.csv") 
# patch bc university of colorado schools are all called the same thing
unique(term_tab$org_congdist[term_tab$org_name == "UNIVERSITY OF COLORADO"])
term_tab$org_name[term_tab$org_name=="UNIVERSITY OF COLORADO" & term_tab$org_city != "Boulder"] <- "UNIVERSITY OF COLORADO - COLORADO SPRINGS"
term_tab$org_name[term_tab$org_name=="HARVARD UNIVERSITY D/B/A HARVARD SCHOOL OF PUBLIC HEALTH"] <- "HARVARD SCHOOL OF PUBLIC HEALTH"

# use FY2024 grants to geolocate terminated grants (by matching to org name)
loc_tab <- NIH_clean_fips %>%
  mutate(FIPS = str_pad(FIPS, 5, pad = "0")) %>%
  mutate(org_name = toupper(org_name)) %>%
  mutate(lat = round(as.numeric(lat), 5),
         lon = round(as.numeric(lon), 5)) %>%
  select(org_name, lat, lon, FIPS, tractid) %>%
  distinct() %>%
  filter(!is.na(org_name))

loc_tab <-read.csv("data/geoid_dictionary_July4.csv") %>%
            select(-GEOID) %>%
            mutate(GEOID = str_pad(New119Fips, 4, pad = "0"),
                   lat = round(lat, 5),
                   lon = round(lon, 5)) %>%
            select(lat, lon, GEOID) %>%
            distinct() %>%
            left_join(loc_tab) %>%
            select(org_name, lat, lon, GEOID, FIPS, tractid) %>%
            filter(!is.na(org_name))

loc_tab$org_name[loc_tab$org_name=="UNIVERSITY OF COLORADO" & loc_tab$GEOID == "0805"] <- "UNIVERSITY OF COLORADO - COLORADO SPRINGS"
loc_tab$org_name[loc_tab$org_name=="UNIVERSITY OF COLORADO" & loc_tab$GEOID == "0802"] <- "UNIVERSITY OF COLORADO - BOULDER"

# find termination grants that aren't in the location table
# (org didn't receive funds in FY2024)
canada_prov <- c("AB", "BC", "MB", "NB", "NL", "NS", "ON", "PE", "QC", "SK", "NT", "NU", "YT")

term_tab %>%
  filter(!org_name %in% loc_tab$org_name) %>%
  # remove grants from Canada, territories
  filter(! is.na(org_state) & ! org_state %in% canada_prov & org_state != "PR" & org_state != "VI") %>%
  select(org_name) %>% unique() -> missing_loc

static_loss <- read.csv("output/static_loss_county.csv") %>%
mutate(FIPS = str_pad(FIPS, 5, pad = "0"))
filter(static_loss, terminated_econ_loss>0) %>% select(FIPS) %>% distinct() %>% nrow()

NIH_prior <- read_csv("https://osf.io/ex89k/download")
names(NIH_prior) <- make.names(names(NIH_prior))

# find the lat/lon for these orgs based on prior year funding
missing_coords <- NIH_prior %>%
  mutate(org_name = toupper(Organization.Name)) %>%
  select(org_name, Longitude, Latitude) %>%
  rename(lon = Longitude,
         lat = Latitude) %>%
  distinct() %>%
  filter(!is.na(org_name)) %>%
  filter(org_name %in% missing_loc$org_name)

# if there are more organizations added with missing coords, will need to relocate
# using API call: https://geo.fcc.gov/api/census/block/find?latitude=lat&longitude=long&showall=true&format=json

missing_fips <- data.frame(org_name =  c("ECOHEALTH ALLIANCE, INC.",
                                         #"FUNCTION PROMOTING THERAPIES, LLC",
                                         "BENTLEY UNIVERSITY",
                                         "RESILIENT GAMES STUDIO, LLC",
                                         "HEBREW HOME FOR THE AGED AT RIVERDALE",
                                         "NEW YORK UNIVERSITY D/B/A NYU LONG ISLAND SCHOOL OF MEDICINE",
                                         "FAMILIES AGAINST NARCOTICS, INC. D/B/A FACE ADDICTION NOW"),
                                         #"PRIOGEN CORP"),
                           FIPS = c("36061", "25017", "17031", "36005", "36059", "26099"),
                           GEOID = c("3612", "2505", "1701", "3615", "3603", "2610"),
                           tractid = c("36061010300","25017369000","17031410600","36005030900", "36059303600", "26099242500")
                           ) %>%
                right_join(missing_coords)

missing_coords$lat[missing_coords$org_name=="NEW YORK UNIVERSITY D/B/A NYU LONG ISLAND SCHOOL OF MEDICINE"]<-40.741556066276
missing_coords$lon[missing_coords$org_name=="NEW YORK UNIVERSITY D/B/A NYU LONG ISLAND SCHOOL OF MEDICINE"]<--73.640354915859
missing_coords$lat[missing_coords$org_name=="FAMILIES AGAINST NARCOTICS, INC. D/B/A FACE ADDICTION NOW"]<-42.6191401
missing_coords$lon[missing_coords$org_name=="FAMILIES AGAINST NARCOTICS, INC. D/B/A FACE ADDICTION NOW"]<--82.9550412


loc_tab %<>% rbind(missing_fips)

# Append FIPS code to terminated grants
term_tab %<>%
  distinct() %>%
  left_join(loc_tab) %>%
  filter(!is.na(FIPS)) %>%
  mutate(state_FIPS = substr(FIPS, 1, 2))

term_tab <- repeated_dict %>%
  right_join(term_tab) %>%
  mutate(org_name = ifelse(fix_name=="" | is.na(fix_name), org_name, fix_name))

# read in state FIPS to name dictionary
state_dict <- read.csv("data/state_and_county_fips_master.csv") %>%
  mutate(state_FIPS = str_pad(fips, 5, pad = "0")) %>%
  filter(endsWith(state_FIPS, "000")) %>%
  mutate(state_FIPS = substr(state_FIPS, 1, 2)) %>%
  mutate(state = str_to_title(tolower(name))) %>%
  select(-name)

# manually fix issue with capitalization of "of" 
state_dict$state <- ifelse(state_dict$state == "District Of Columbia", "District of Columbia" , state_dict$state)

#add state name and abbreviation
term_tab <- read.csv("data/state_abbrev.csv") %>%
  rename(state_code = Abbreviation,
         state = State) %>%
  left_join(state_dict) %>%
  right_join(term_tab) 

# summary statistics at the point (organization) level
term_points <- term_tab %>%
  group_by(lat, lon, org_name) %>%
  summarize(terminated_loss = sum(award_remaining, na.rm=TRUE),
            terminated_num = n(),
            terminated_loss_noself = sum(award_remaining_noself, na.rm=TRUE),
            terminated_num_noself = sum(award_remaining_noself!=0 & !is.na(award_remaining_noself))) %>%
  mutate(org_name = str_to_title(org_name) ) %>%
  filter(terminated_loss != 0)

# if new organizations were added, set the number of new orgs below and then uncomment and run 
# through line 217 to identify new repeats

# some blanket fixes on org name
replacements <- c("Of" = "of", "At" = "at", "And" = "and",
                  "Univ " = "University ", "Ctr" = "Center", 
                  "Col " = "College ", "Hlth" = "Health", "Scis" = "Sciences",
                  "Sci " = "Science ", "Llc" = "LLC",
                  "The" = "the", "For" = "for", "Res " = "Research ", 
                  "Hosp " = "Hospital ", "Med " = "Medical ")

term_points$org_name <- str_replace_all(term_points$org_name, replacements)

# fix names for specific universities
corrected_df <- read.csv("data/org_names_corrected.csv") 

corrected_df %>%
  right_join(term_points) %>%
  mutate(org_name = ifelse(is.na(corrected_name) | corrected_name == "", org_name, corrected_name)) %>%
  select(-corrected_name)-> term_points

write.csv(term_points, "output/terminated_points.csv", row.names = FALSE)


sum(term_points$terminated_loss)

```

```{r Calculate Losses}
# state-specific multipliers: read in 
state_mult <- read.csv("data/state-multiplier-urm.csv") %>% 
                rename(tot_fund = 2, 
                       econ_act = 9,
                       state.name= 1,
                       jobs_per_m = 3) %>%
                mutate(econ_act = as.numeric(gsub(",", "", econ_act)),
                       tot_fund = as.numeric(gsub(",", "", tot_fund)),
                       job_val = as.numeric(gsub(",", "", jobs_per_m))/1000000) %>%
                mutate(econ_mult = econ_act/tot_fund) %>%
                select(state.name, econ_mult, job_val) %>%
                left_join(data.frame(state.name=state.name,
                                     state=state.abb)) %>%
                select(-state.name) %>%
                filter(!is.na(state))

# use state-specific multipliers
NIH_clean_fips %<>%
  left_join(state_mult)
                        
#if a 15% indirect cost rate max is imposed, this would be the loss
NIH_clean_fips$IDC_loss <- NIH_clean_fips$indirect - .15 * NIH_clean_fips$direct

#set negative values to zero; these are cases where there were no indirects before, or the IDC rate was <15% already
NIH_clean_fips$IDC_loss <- ifelse(NIH_clean_fips$IDC_loss < 0, 0, NIH_clean_fips$IDC_loss) 

#calculate estimated economic impact, using a state-specific multiplier (from a 2025 NIH report)
NIH_clean_fips$grant_funds_econ <- NIH_clean_fips$grant_funds * NIH_clean_fips$econ_mult #status quo; economic activity resulting from grants active in 2024

NIH_clean_fips$IDC_econ_loss <- NIH_clean_fips$IDC_loss * NIH_clean_fips$econ_mult #projected losses with 15% IDC

NIH_clean_fips$IDC_job_loss <- NIH_clean_fips$IDC_loss*NIH_clean_fips$job_val

# identify overlap between terminated grants and potential IDC loss (will use this to calculate total losses without repetition)
term_tab_noself <- term_tab %>% filter(cancellation_source != "Self reported")

term_tab$full_award_number <-  substring(term_tab$full_award_number, 2)

NIH_clean_fips %<>%
  mutate(
    overlap_loss = ifelse(proj.num %in% term_tab$full_award_number, IDC_loss, 0),
    overlap_econ_loss = ifelse(proj.num %in% term_tab$full_award_number, IDC_econ_loss, 0),
    overlap_job_loss = ifelse(proj.num %in% term_tab$full_award_number, IDC_job_loss, 0),
    overlap_loss_noself = ifelse(proj.num %in% term_tab_noself$full_award_number, IDC_loss, 0),
    overlap_econ_loss_noself = ifelse(proj.num %in% term_tab_noself$full_award_number, IDC_econ_loss, 0),
    overlap_job_loss_noself = ifelse(proj.num %in% term_tab_noself$full_award_number, IDC_job_loss, 0)
  ) %>%
  select(-econ_mult)

#save out long form data with FIPS codes
write.csv(NIH_clean_fips, "./output/NIH_clean_fips.csv", row.names=F)

```

```{r Commuter Weighting}
# commute is a dataframe derived from 2022 Census data with the following columns
# ORIGIN: is the FIPS code for the county of residence
# DESTINATION: is the FIPS code for the county of work
# COMMUTES: is the number of workers commuting from the specific ORIGIN to DESTINATION

# IMPORTANT: some FIPS codes have leading zeros

commute_raw <- fread("https://osf.io/8tq67/download") %>% 
  mutate(DESTINATION = str_pad(W_geocode, 11, pad = "0"),
         ORIGIN = str_pad(H_geocode, 5, pad = "0"))

# calculate proportion weights
commute <- commute_raw %>% 
  # calculate the total number of workers who work in a given county
  group_by(DESTINATION) %>%
  summarize(DESTINATION_total_workers = sum(Commuters)) %>%
  # join this back to the full commuter dataframe
  left_join(commute_raw) %>%
  # for each ORIGIN - DESTINATION pair, calculate proportion as the number
  # of workers from ORIGIN who work in DESTINATION divided by the total 
  # number of people who work in DESTINATION
  mutate(proportion = Commuters/DESTINATION_total_workers)

# Quality check: all destination weights should sum to 1
commute %>% 
  group_by(DESTINATION) %>% 
  summarize(sum = round(sum(proportion), digits=5)) %>% 
  filter(sum < 1) %>%
  nrow() 

```



```{r County Calculations}
term_tab %>%
  select("award_remaining", "award_remaining_noself", "FIPS", "state", "state_code", "tractid") %>%
  left_join(state_mult %>% rename(state_code = state)) %>%
  mutate(award_remaining = ifelse(award_remaining <0 , 0 , award_remaining),
         award_remaining_noself = ifelse(award_remaining_noself <0 , 0 , award_remaining_noself)) %>%
  rename(terminated_loss = award_remaining,
         terminated_loss_noself = award_remaining_noself) %>%
  mutate(terminated_econ_loss = terminated_loss * econ_mult,
         terminated_job_loss = terminated_loss*job_val,
         terminated_econ_loss_noself = terminated_loss_noself * econ_mult,
         terminated_job_loss_noself = terminated_loss_noself*job_val) %>%
  group_by(tractid, FIPS) %>%
  summarise(across(starts_with(c("terminated")), ~ sum(.x, na.rm = TRUE))) -> term_sum_tract

term_sum_tract %>%
  group_by(FIPS) %>%
  summarise(across(starts_with(c("terminated")), ~ sum(.x, na.rm = TRUE))) -> term_sum_FIPS

#summarize NIH data by tract id and FIPS code
NIH_sum_tract <- NIH_clean_fips %>% 
  group_by(tractid) %>% 
  summarise(across(starts_with(c("IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) 

NIH_sum_FIPS <- NIH_clean_fips %>% 
  group_by(FIPS) %>% 
  summarise(across(starts_with(c("IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) 


NIH_sum_FIPS %>% merge(term_sum_FIPS %>% select(FIPS, terminated_econ_loss), all.x=TRUE) %>%
write.csv("output/static_loss_county.csv", row.names=FALSE)

#combine commuter data with NIH data
county_commute_NIH <- merge(commute, term_sum_tract %>% select(-FIPS), by.x = "DESTINATION", by.y = "tractid", all.x = TRUE) %>%
                  merge(NIH_sum_tract, by.x = "DESTINATION", by.y = "tractid", all.x = TRUE)

#replace all NA values with zeros
county_commute_NIH <- mutate_all(county_commute_NIH , funs(replace(., is.na(.), 0)))

write.csv(county_commute_NIH, "output/county_commute_NIH.csv", row.names=FALSE)

county_commute_NIH_sum <- county_commute_NIH %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ . * proportion)) %>%
  rename(FIPS = ORIGIN) %>%
  group_by(FIPS) %>% 
  summarise(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) %>%
  mutate(state_FIPS =  substr(FIPS, 1, 2)) %>%
  # calculate total as terminated + IDC - overlap
  mutate(combined_loss = terminated_loss + IDC_loss - overlap_loss,
         combined_econ_loss = terminated_econ_loss + IDC_econ_loss - overlap_econ_loss,
         combined_job_loss = terminated_job_loss + IDC_job_loss - overlap_job_loss,
         combined_loss_noself = terminated_loss_noself + IDC_loss - overlap_loss_noself,
         combined_econ_loss = terminated_econ_loss_noself + IDC_econ_loss - overlap_econ_loss_noself,
         combined_job_loss = terminated_job_loss_noself + IDC_job_loss - overlap_job_loss_noself) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ log(.x + 1), .names =  "{.col}_log"))

# read in state FIPS to name dictionary
state_dict <- read.csv("data/state_fips_master.csv") %>%
  mutate(state_FIPS = str_pad(fips, 5, pad = "0")) %>%
  filter(endsWith(state_FIPS, "000")) %>%
  mutate(state_FIPS = substr(state_FIPS, 1, 2)) %>%
  mutate(state = str_to_title(tolower(name))) %>%
  select(-c(name, fips))

# manually fix issue with capitalization of "of" 
state_dict$state <- ifelse(state_dict$state == "District Of Columbia", "District of Columbia" , state_dict$state)

#add state name
county_commute_NIH_sum <- read.csv("data/county_pop.csv") %>%
                      mutate(FIPS = str_pad(FIPS, 5, pad = "0")) %>%
                      right_join(county_commute_NIH_sum) 

county_commute_NIH_sum <- read.csv("data/state_abbrev.csv") %>%
  rename(state_code = Abbreviation,
         state = State) %>%
  left_join(state_dict) %>%
  right_join(county_commute_NIH_sum) %>%
  select(-state_FIPS) %>%
  filter(FIPS != "0")

print(nrow(county_commute_NIH_sum))
sum(is.na(county_commute_NIH_sum$state_code))

write.csv(county_commute_NIH_sum, "output/NIH_impact_county.csv", row.names=FALSE)
```

```{r institution level calculations}
NIH_clean_fips %>%
  group_by(org_name, tractid) %>%
  summarize(across(starts_with("IDC"), ~ sum(., na.rm=TRUE)))-> IDC_inst

#combine commuter data with NIH data
county_commute_NIH_inst <- merge(commute, IDC_inst, by.x = "DESTINATION", by.y = "tractid", all.x = TRUE)

#replace all NA values with zeros
county_commute_NIH_inst <- mutate_all(county_commute_NIH_inst , funs(replace(., is.na(.), 0))) %>%
  mutate(across(starts_with(c("IDC")), ~ . * proportion)) %>%
  rename(FIPS = ORIGIN) %>%
  group_by(FIPS) %>%
  group_by(FIPS, org_name) %>% 
  summarize(across(starts_with(c("IDC")), ~ sum(., na.rm=TRUE))) %>%
  slice_max(IDC_econ_loss, n=5)

write.csv(county_commute_NIH_inst, "output/county_commute_NIH_inst.csv", row.names=FALSE)
  
```

```{r county-stayvflow}
#proportion of funds we assume stays in origin county - no commuting
county_varystatic <- map_df(seq(0, 1, by = 0.05), function(rho) {
  county_commute_NIH %>%
    mutate(prop_origin = ifelse(ORIGIN == substr(DESTINATION, 1, 5), rho + (1 - rho) * proportion, (1 - rho) * proportion)) %>%
    mutate(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ . * prop_origin)) %>%
    rename(FIPS = ORIGIN) %>%
    group_by(FIPS) %>% 
    summarise(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) %>%
    mutate(rho = rho)
})


print(nrow(county_varystatic))

write.csv(county_varystatic, "output/NIH_impact_county-sens.csv", row.names=FALSE)

```


```{r county_county-sens}

commute_raw_ctc <- fread("https://osf.io/8tq67/download") %>% 
  mutate(DESTINATION = str_pad(W_geocode, 11, pad = "0"),
         ORIGIN = str_pad(H_geocode, 5, pad = "0")) %>%
  mutate(DESTINATION = substr(1, 5, DESTINATION))

# calculate proportion weights
commute_ctc <- commute_raw_ctc %>% 
  # calculate the total number of workers who work in a given county
  group_by(DESTINATION) %>%
  summarize(DESTINATION_total_workers = sum(Commuters)) %>%
  # join this back to the full commuter dataframe
  left_join(commute_raw_ctc) %>%
  # for each ORIGIN - DESTINATION pair, calculate proportion as the number
  # of workers from ORIGIN who work in DESTINATION divided by the total 
  # number of people who work in DESTINATION
  mutate(proportion = Commuters/DESTINATION_total_workers)

#combine commuter data with NIH data
county_commute_NIH_ctc <- merge(commute_ctc, term_sum_FIPS, by.x = "DESTINATION", by.y = "FIPS", all.x = TRUE) %>%
                  merge(NIH_sum_FIPS, by.x = "DESTINATION", by.y = "FIPS", all.x = TRUE)

#replace all NA values with zeros

county_commute_NIH_ctc_sum <- county_commute_NIH_ctc %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ . * proportion)) %>%
  rename(FIPS = ORIGIN) %>%
  group_by(FIPS) %>% 
  summarise(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) %>%
  mutate(state_FIPS =  substr(FIPS, 1, 2)) %>%
  # calculate total as terminated + IDC - overlap
  mutate(combined_loss = terminated_loss + IDC_loss - overlap_loss,
         combined_econ_loss = terminated_econ_loss + IDC_econ_loss - overlap_econ_loss,
         combined_job_loss = terminated_job_loss + IDC_job_loss - overlap_job_loss,
         combined_loss_noself = terminated_loss_noself + IDC_loss - overlap_loss_noself,
         combined_econ_loss = terminated_econ_loss_noself + IDC_econ_loss - overlap_econ_loss_noself,
         combined_job_loss = terminated_job_loss_noself + IDC_job_loss - overlap_job_loss_noself) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ log(.x + 1), .names =  "{.col}_log"))

write.csv(county_commute_NIH_ctc_sum, "output/NIH_impact_county-ctc.csv", row.names=FALSE)


```


```{r natl-mult-sensitivity}
# use national multipliers
county_commute_NIH_diffmult <- county_commute_NIH %>% 
  mutate(terminated_econ_loss = terminated_loss * 2.56,
         IDC_econ_loss = IDC_loss * 2.56,
         overlap_econ_loss = overlap_loss * 2.56,
         overlap_econ_loss_noself = overlap_loss_noself * 2.56) 
                        
diff_mult_county <-  county_commute_NIH_diffmult %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ . * proportion)) %>%
  rename(FIPS = ORIGIN) %>%
  group_by(FIPS) %>% 
  summarise(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) %>%
  mutate(state_FIPS =  substr(FIPS, 1, 2)) %>%
  # calculate total as terminated + IDC - overlap
  mutate(combined_loss = terminated_loss + IDC_loss - overlap_loss,
         combined_econ_loss = terminated_econ_loss + IDC_econ_loss - overlap_econ_loss,
         combined_job_loss = terminated_job_loss + IDC_job_loss - overlap_job_loss,
         combined_loss_noself = terminated_loss_noself + IDC_loss - overlap_loss_noself,
         combined_econ_loss = terminated_econ_loss_noself + IDC_econ_loss - overlap_econ_loss_noself,
         combined_job_loss = terminated_job_loss_noself + IDC_job_loss - overlap_job_loss_noself) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ log(.x + 1), .names =  "{.col}_log"))

print(nrow(diff_mult_county))
sum(is.na(diff_mult_county$state_code))

write.csv(diff_mult_county, "output/NIH_impact_county-natlmult.csv", row.names=FALSE)

```


```{r 2022-flows}
commute_raw_22 <- read_csv("https://osf.io/tqkj7/download") %>% 
  mutate(DESTINATION = str_pad(work_tract, 11, pad = "0"),
         ORIGIN = str_pad(home_county, 5, pad = "0"))

# calculate proportion weights
commute_22 <- commute_raw_22 %>% 
  # calculate the total number of workers who work in a given county
  group_by(DESTINATION) %>%
  summarize(DESTINATION_total_workers = sum(commutes)) %>%
  # join this back to the full commuter dataframe
  left_join(commute_raw_22) %>%
  # for each ORIGIN - DESTINATION pair, calculate proportion as the number
  # of workers from ORIGIN who work in DESTINATION divided by the total 
  # number of people who work in DESTINATION
  mutate(proportion = commutes/DESTINATION_total_workers)

#combine commuter data with NIH data
county_commute_NIH_22 <- merge(commute_22, term_sum_tract %>% select(-FIPS), by.x = "DESTINATION", by.y = "tractid", all.x = TRUE) %>%
                  merge(NIH_sum_tract, by.x = "DESTINATION", by.y = "tractid", all.x = TRUE)

#replace all NA values with zeros
county_commute_NIH_22 <- mutate_all(county_commute_NIH_22 , funs(replace(., is.na(.), 0)))

county_commute_NIH_sum_22 <- county_commute_NIH_22 %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ . * proportion)) %>%
  rename(FIPS = ORIGIN) %>%
  group_by(FIPS) %>% 
  summarise(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) %>%
  mutate(state_FIPS =  substr(FIPS, 1, 2)) %>%
  # calculate total as terminated + IDC - overlap
  mutate(combined_loss = terminated_loss + IDC_loss - overlap_loss,
         combined_econ_loss = terminated_econ_loss + IDC_econ_loss - overlap_econ_loss,
         combined_job_loss = terminated_job_loss + IDC_job_loss - overlap_job_loss,
         combined_loss_noself = terminated_loss_noself + IDC_loss - overlap_loss_noself,
         combined_econ_loss = terminated_econ_loss_noself + IDC_econ_loss - overlap_econ_loss_noself,
         combined_job_loss = terminated_job_loss_noself + IDC_job_loss - overlap_job_loss_noself) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ log(.x + 1), .names =  "{.col}_log"))


#add state name
county_commute_NIH_sum_22 <- read.csv("data/county_pop.csv") %>%
                      mutate(FIPS = str_pad(FIPS, 5, pad = "0")) %>%
                      right_join(county_commute_NIH_sum_22) 

county_commute_NIH_sum_22 <- read.csv("data/state_abbrev.csv") %>%
  rename(state_code = Abbreviation,
         state = State) %>%
  left_join(state_dict) %>%
  right_join(county_commute_NIH_sum_22) %>%
  select(-state_FIPS) %>%
  filter(FIPS != "0") 

write.csv(county_commute_NIH_sum_22, "output/NIH_impact_county-22.csv", row.names=FALSE)
```

```{r rac-wac-output}
# calculate the proportion of residents in relevant sectors
rac_df <- read.csv("data/RAC_countyAll.csv") %>% 
          mutate(ORIGIN = str_pad(HomeCounty, 5, pad = "0")) %>%
          mutate(broad_sum_rac = sum.CNS03.+sum.CNS04.+sum.CNS05.+sum.CNS07.+
                             sum.CNS08.+sum.CNS12.+sum.CNS13.+sum.CNS14.+
                             sum.CNS16.,
                 narrow_sum_rac = sum.CNS12.+sum.CNS13.+sum.CNS14.+sum.CNS16.) %>%
          mutate(broad_prop_rac = broad_sum_rac/AllCommutes,
                 narrow_prop_rac = narrow_sum_rac/AllCommutes) %>%
          mutate(broad_prop_rac = ifelse(is.infinite(broad_prop_rac), 0, broad_prop_rac),
               narrow_prop_rac = ifelse(is.infinite(narrow_prop_rac), 0, narrow_prop_rac)) %>%
          select(ORIGIN, broad_prop_rac, narrow_prop_rac)

 
# calculate proportion of workers in relevant sectors
wac_df <- read.csv("data/WAC_tractAll.csv") %>% 
          mutate(DESTINATION = str_pad(WorkTract, 11, pad = "0")) %>%
          mutate(broad_sum_wac = sum.CNS03.+sum.CNS04.+sum.CNS05.+sum.CNS07.+
                             sum.CNS08.+sum.CNS12.+sum.CNS13.+sum.CNS14.+
                             sum.CNS16.,
                 narrow_sum_wac = sum.CNS12.+sum.CNS13.+sum.CNS14.+sum.CNS16.) %>%
          mutate(broad_prop_wac = broad_sum_wac/AllCommutes,
                 narrow_prop_wac = narrow_sum_wac/AllCommutes) %>%
          mutate(broad_prop_wac = ifelse(is.infinite(broad_prop_wac), 0, broad_prop_wac),
                 narrow_prop_wac = ifelse(is.infinite(narrow_prop_wac), 0, narrow_prop_wac)) %>%
          select(DESTINATION, broad_prop_wac, narrow_prop_wac)

# join rac and wac files to commuter flow
commute_raw_sect <- fread("https://osf.io/8tq67/download") %>% 
  mutate(DESTINATION = str_pad(W_geocode, 11, pad = "0"),
         ORIGIN = str_pad(H_geocode, 5, pad = "0")) %>%
  select(DESTINATION, ORIGIN, Commuters) %>%
  left_join(rac_df) %>%
  left_join(wac_df) %>%
  # calculate geometric mean of rac and wac 
  mutate(broad_prop = sqrt(broad_prop_wac*broad_prop_rac),
         narrow_prop = sqrt(narrow_prop_wac*narrow_prop_rac)) %>%
  # rescale commuters 
  select(-c(broad_prop_wac, narrow_prop_wac, broad_prop_rac, narrow_prop_rac)) %>%
  mutate(broad_commuters = broad_prop * Commuters,
         narrow_commuters = narrow_prop * Commuters) 
  
# calculate proportion weights
commute_sect <- commute_raw_sect %>% 
  # calculate the total number of workers who work in a given county
  group_by(DESTINATION) %>%
  summarize(DESTINATION_total_workers_broad = sum(broad_commuters, na.rm=TRUE),
            DESTINATION_total_workers_narrow = sum(narrow_commuters, na.rm=TRUE)) %>%
  # join this back to the full commuter dataframe
  left_join(commute_raw_sect) %>%
  # for each ORIGIN - DESTINATION pair, calculate proportion as the number
  # of workers from ORIGIN who work in DESTINATION divided by the total 
  # number of people who work in DESTINATION
  mutate(proportion_broad = broad_commuters/DESTINATION_total_workers_broad,
         proportion_narrow = narrow_commuters/DESTINATION_total_workers_narrow)

#combine commuter data with NIH data
county_commute_NIH_sect <- merge(commute_sect, term_sum_tract, by.x = "DESTINATION", by.y = "tractid", all.x = TRUE) %>%
                  merge(NIH_sum_tract, by.x = "DESTINATION", by.y = "tractid", all.x = TRUE)

#replace all NA values with zeros
county_commute_NIH_sect <- mutate_all(county_commute_NIH_sect , funs(replace(., is.na(.), 0)))

county_commute_NIH_sum_broad <- county_commute_NIH_sect %>%
  select(-FIPS) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ . * proportion_broad)) %>%
  rename(FIPS = ORIGIN) %>%
  group_by(FIPS) %>% 
  summarise(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) %>%
  mutate(state_FIPS =  substr(FIPS, 1, 2)) %>%
  # calculate total as terminated + IDC - overlap
  mutate(combined_loss = terminated_loss + IDC_loss - overlap_loss,
         combined_econ_loss = terminated_econ_loss + IDC_econ_loss - overlap_econ_loss,
         combined_job_loss = terminated_job_loss + IDC_job_loss - overlap_job_loss,
         combined_loss_noself = terminated_loss_noself + IDC_loss - overlap_loss_noself,
         combined_econ_loss = terminated_econ_loss_noself + IDC_econ_loss - overlap_econ_loss_noself,
         combined_job_loss = terminated_job_loss_noself + IDC_job_loss - overlap_job_loss_noself) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ log(.x + 1), .names =  "{.col}_log"))

#add state name
county_commute_NIH_sum_broad <- read.csv("data/county_pop.csv") %>%
                      mutate(FIPS = str_pad(FIPS, 5, pad = "0")) %>%
                      right_join(county_commute_NIH_sum_broad) 

county_commute_NIH_sum_broad <- read.csv("data/state_abbrev.csv") %>%
  rename(state_code = Abbreviation,
         state = State) %>%
  left_join(state_dict) %>%
  right_join(county_commute_NIH_sum_broad) %>%
  select(-state_FIPS) %>%
  filter(FIPS != "0")

write.csv(county_commute_NIH_sum_broad, "output/NIH_impact_county_broad.csv", row.names=FALSE)

county_commute_NIH_sum_narrow <- county_commute_NIH_sect %>%
  select(-FIPS) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ . * proportion_narrow)) %>%
  rename(FIPS = ORIGIN) %>%
  group_by(FIPS) %>% 
  summarise(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) %>%
  mutate(state_FIPS =  substr(FIPS, 1, 2)) %>%
  # calculate total as terminated + IDC - overlap
  mutate(combined_loss = terminated_loss + IDC_loss - overlap_loss,
         combined_econ_loss = terminated_econ_loss + IDC_econ_loss - overlap_econ_loss,
         combined_job_loss = terminated_job_loss + IDC_job_loss - overlap_job_loss,
         combined_loss_noself = terminated_loss_noself + IDC_loss - overlap_loss_noself,
         combined_econ_loss = terminated_econ_loss_noself + IDC_econ_loss - overlap_econ_loss_noself,
         combined_job_loss = terminated_job_loss_noself + IDC_job_loss - overlap_job_loss_noself) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ log(.x + 1), .names =  "{.col}_log"))

#add state name
county_commute_NIH_sum_narrow <- read.csv("data/county_pop.csv") %>%
                      mutate(FIPS = str_pad(FIPS, 5, pad = "0")) %>%
                      right_join(county_commute_NIH_sum_narrow) 

county_commute_NIH_sum_narrow <- read.csv("data/state_abbrev.csv") %>%
  rename(state_code = Abbreviation,
         state = State) %>%
  left_join(state_dict) %>%
  right_join(county_commute_NIH_sum_narrow) %>%
  select(-state_FIPS) %>%
  filter(FIPS != "0")

write.csv(county_commute_NIH_sum_narrow, "output/NIH_impact_county_narrow.csv", row.names=FALSE)

```

# Congressional Districts

```{r Prep District Commuter Flows}
cong_commute <- fread("https://osf.io/m5kb6/download") %>%
  mutate(Cong_DESTINATION = str_pad(w_tract, 11, pad = "0"),
         Cong_ORIGIN = str_pad(h_cd119, 4, pad = "0")) %>%
  rename(commuters = commutes) %>%
  filter(!is.na(Cong_ORIGIN))

cong_info <- read.csv("./data/us_representatives_119th_congress.csv") %>%
  mutate(GEOID = str_pad(GEOID, 4, pad = "0")) %>%
  select(GEOID, STATE_NAME, NAME, PARTY) %>%
  rename(state = STATE_NAME,
         pol_party = PARTY,
         rep_name = NAME)

#get total workers per cong district
cong_commute_sum <- cong_commute %>%
  group_by(Cong_DESTINATION) %>%
  summarise(total_workers = sum(commuters, na.rm=TRUE))

#merge in total workers
cong_commute <- left_join(cong_commute, cong_commute_sum)

#calculate proportion of workers who work in district X who live in district Y
cong_commute$proportion <- cong_commute$commuters / cong_commute$total_workers
```

```{r Merge Cong Data}
#summarize NIH data by cong district code
NIH_sum_GEOID <- NIH_clean_fips %>% 
  group_by(GEOID) %>% 
  summarise(across(starts_with(c("IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) 

NIH_sum_tract <- NIH_clean_fips %>% 
  group_by(tractid) %>% 
  summarise(across(starts_with(c("IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) 

term_tab %>%
  select("award_remaining", "award_remaining_noself", "GEOID", "state", "state_code", "tractid") %>%
  left_join(state_mult %>% rename(state_code = state)) %>%
  mutate(award_remaining = ifelse(award_remaining <0 , 0 , award_remaining),
         award_remaining_noself = ifelse(award_remaining_noself <0 , 0 , award_remaining_noself)) %>%
  rename(terminated_loss = award_remaining,
         terminated_loss_noself = award_remaining_noself) %>%
  mutate(terminated_econ_loss = terminated_loss * econ_mult,
         terminated_job_loss = terminated_loss*job_val,
         terminated_econ_loss_noself = terminated_loss_noself * econ_mult,
         terminated_job_loss_noself = terminated_loss_noself*job_val) %>%
  group_by(tractid, GEOID) %>%
  summarise(across(starts_with(c("terminated")), ~ sum(.x, na.rm = TRUE))) -> term_sum_tract

term_sum_tract %>%
  group_by(GEOID) %>%
  summarise(across(starts_with(c("terminated")), ~ sum(.x, na.rm = TRUE)))  -> term_sum_GEOID


NIH_sum_GEOID %>% merge(term_sum_GEOID %>% select(GEOID, terminated_econ_loss), all.x=TRUE) %>%
write.csv("output/static_loss_cong.csv", row.names=FALSE)

#combine commuter data with NIH data
cong_commute_NIH <- merge(cong_commute, term_sum_tract, by.x = "Cong_DESTINATION", by.y = "tractid", all = TRUE) %>%
                      merge(., NIH_sum_tract, by.x = "Cong_DESTINATION", by.y = "tractid", all = TRUE)
```


```{r District Calculations}
#replace all NA values with zeros
cong_commute_NIH <- mutate_all(cong_commute_NIH, funs(replace(., is.na(.), 0)))

#use commuter flows (e.g., proportion of people who work in district X but live in district Y) to weight the impact of NIH funding
cong_commute_NIH_sum <- cong_commute_NIH %>%
  select(-GEOID) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ . * proportion)) %>%
  rename(GEOID = Cong_ORIGIN) %>%
  group_by(GEOID) %>% 
  summarise(across(starts_with(c("terminated", "IDC", "grant", "overlap")), ~ sum(.x, na.rm = TRUE))) %>%
  mutate(state_FIPS =  substr(GEOID, 1, 2)) %>%
  mutate(combined_loss = terminated_loss + IDC_loss - overlap_loss,
       combined_econ_loss = terminated_econ_loss + IDC_econ_loss - overlap_econ_loss,
       combined_job_loss = terminated_job_loss + IDC_job_loss - overlap_job_loss,
       combined_loss_noself = terminated_loss_noself + IDC_loss - overlap_loss_noself,
       combined_econ_loss = terminated_econ_loss_noself + IDC_econ_loss - overlap_econ_loss_noself,
       combined_job_loss = terminated_job_loss_noself + IDC_job_loss - overlap_job_loss_noself) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ log(.x + 1), .names =  "{.col}_log"))

#add in rep names
cong_commute_NIH_sum <- left_join(cong_commute_NIH_sum, cong_info)

#add state name
cong_commute_NIH_sum <- read.csv("data/state_abbrev.csv") %>%
  rename(state_code = Abbreviation,
         state = State) %>%
  right_join(cong_commute_NIH_sum) %>%
  select(-state_FIPS)

print(nrow(cong_commute_NIH_sum))

write.csv(cong_commute_NIH_sum, "output/NIH_impact_cong.csv", row.names=F)

```

## State-Level Summaries

```{r State Estimates}
cong_commute_NIH_sum %>%
  select(-ends_with("log")) %>%
  group_by(state, state_code) %>%
  summarise(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ sum(.x, na.rm = TRUE))) %>%
  mutate(across(starts_with(c("terminated", "IDC", "grant", "combined")), ~ log(.x + 1), .names =  "{.col}_log")) %>%
  filter(!is.na(state_code))  -> state_cong

print(nrow(state_cong))

write.csv(state_cong, file="output/NIH_impact_state.csv", row.names=FALSE)

```