covid-immigrant-foodworkers.Rmd

---
title: "covid-immigrant-foodworkers"
author: "Joe Yerardi"
date: "9/28/2020"
output:
  html_notebook:
    toc: true
    toc_float: true
---

```{r setup, include = F}
# We don't want to see any warnings from our code
knitr::opts_chunk$set(warning = F)

# We don't want to see any messages
knitr::opts_chunk$set(message = F)

# Set key file for the Census Bureau API key
key_file <- "key_file.txt"
```

```{r, include = F}
library(ipumsr)
library(labelled)
library(purrr)
library(rgdal)
library(sf)
library(tidycensus)
library(tidyverse)
library(tigris)
```

# Import and format data

## Import IPUMS

```{r, echo=FALSE, message=FALSE}
# IPUMS (2014-18 ACS; population = 16 and older; employed - at work,
# employed - not at work)
ipums <- read.table("data/usa_00015.dat")
ddi <- read_ipums_ddi("data/usa_00015.xml")
ipums <- read_ipums_micro(ddi) %>%
  mutate(
    state_code = as_factor(US2018C_ST, levels = "values"),
    state = as_factor(US2018C_ST, levels = "labels"),
    puma_code = paste0(state_code, as_factor(US2018C_PUMA)),
    industry_code = as_factor(US2018C_INDP, levels = "values"),
    industry = as_factor(US2018C_INDP, levels = "labels"),
    occupation_code =  as_factor(US2018C_OCCP, levels = "values"),
    occupation = as_factor(US2018C_OCCP, levels = "labels"),
    citizenship_code = as_factor(CITIZEN, levels = "values"),
    citizenship = as_factor(CITIZEN, levels = "labels"),
    ethnicity_code = as_factor(HISPAN, levels = "values"),
    ethnicity = as_factor(HISPAN, levels = "labels")
  ) %>%
  select(16:26, 5)

ipums
```

## PUMA shapefiles

```{r, echo=FALSE, message=FALSE}
pumas <- readOGR("data/ipums_puma_2010/ipums_puma_2010.shp") %>%
  # Convert shapefile to sf object
  st_as_sf(crs = 5070, remove = F) %>%
  select(puma_code = GEOID, puma = Name, geometry)
```

## County shapefiles

```{r, echo=FALSE, message=FALSE}
counties <-
  readOGR("data/tl_2019_us_county/tl_2019_us_county.shp") %>%
  # Convert shapefile to sf object
  st_as_sf(crs = 5070, remove = F) %>%
  select(county_code = GEOID, county = NAMELSAD, geometry)
```

## COVID county cases and deaths

```{r}
covid <- read_csv("data/us-counties.csv",
                  col_types = cols(.default = "?", date = "D")) %>%
  rename(county_code = fips) %>%
  # Filter to the latest data for each county
  group_by(county_code) %>%
  slice(which.max(date))

covid
```

## County populations

### Read-in Census API key

```{r}
census_key <- (readLines(key_file)[1])
```

### Create list of all US states to iterate through in the API calls

```{r}
us <- unique(fips_codes$state)[1:51]

us
```

### Download population data from the 2014-2018 ACS

```{r, echo=FALSE, message=FALSE}
population <- map_dfr(us, function(x) {
  get_acs(
    geography = "county",
    variables = "B01001_001",
    state = x,
    key = census_key
  )
}) %>%
  mutate(income_quartile = ntile(estimate, 4)) %>%
  select(county_code = GEOID,
         county = NAME,
         population = estimate)

population
```

## PUMA-to-county crosswalk

```{r}
county_to_puma_crosswalk <- read_csv("data/geocorr2018.csv", skip = 1) %>%
  mutate(puma_code = paste0(`State code`, `PUMA (2012)`)) %>%
  select(
    puma_code,
    county_code = `County code (2014)`,
    state = `State abbreviation`,
    county = `2014 county name`,
    puma = `PUMA12 name`,
    allocation_factor = `puma12 to county14 allocation factor`
  )

county_to_puma_crosswalk
```

## Join the COVID, population and crosswalk data frames

```{r}
crosswalk <- list(covid, population, county_to_puma_crosswalk) %>%
  reduce(full_join, by = "county_code") %>%
  mutate(
    puma = paste(puma, state.x, sep = ", "),
    cases_per_100k = round(cases / population * 100000, digits = 0),
    deaths_per_100k = round(deaths / population * 100000, digits = 0)
  ) %>%
  select(
    1,
    state = state.x,
    county_code,
    county = county.y,
    puma_code,
    puma,
    population,
    cases,
    deaths,
    cases_per_100k,
    deaths_per_100k,
    allocation_factor
  )

crosswalk
```

## Filter to the food production industries we're interested in

```{r}
food_production_industries <-
  list("0170", # Crop production
       "0180", # Animal production and aquaculture
       "0290", # Support activities for agriculture and forestry
       "1070", # Animal food, grain and oilseed milling
       "1080", # Sugar and confectionery products
       "1090", # Fruit and vegetable preserving and specialty food manufacturing
       "1170", # Dairy product manufacturing
       "1180", # Animal slaughtering and processing
       "1270", # Bakeries and tortilla manufacturing, except retail bakeries
       "1280", # Seafood and other miscellaneous foods; n.e.c.
       "1290") # Not specified food industries

ipums_food_production <- ipums %>%
  filter(industry_code %in% food_production_industries)
```

## Calculate employment by industry-occupation grouping

```{r}
ipums_food_production_employment_by_ind_occ <-
  ipums_food_production %>%
  group_by(industry, occupation) %>%
  summarize(num_employees = sum(PERWT)) %>%
  pivot_wider(
    id_cols = c(industry, occupation),
    values_from = num_employees,
    names_prefix = "num_employees"
  ) %>%
  replace(is.na(.), 0) %>%
  # Exclude from our analysis ind-occ groupings
  # with fewer than 2,500 employees nationally
  filter(num_employees >= 2500)
```

## Export the data

```{r}
write_csv(
  ipums_food_production_employment_by_ind_occ,
  "data/exported/ipums_food_production_employment_by_ind_occ.csv"
)
```

## Filter the industry-occupation pairs

### Import data frame of food production industry-occupation pairs

```{r}
ipums_food_production_employment_by_ind_occ <-
  read_csv("data/ipums_food_production_total_employment_by_ind_occ.csv")
```

### Join the industry-occupation pairs to the ipums data and filter to just the industry-occupation pairs we want

```{r}
ipums_food_production_selected_occupations <-
  ipums_food_production %>%
  inner_join(ipums_food_production_employment_by_ind_occ,
             by = c("industry", "occupation")) %>%
  filter(exclude == FALSE) %>%
  select(-num_employees, -exclude)

ipums_food_production_selected_occupations
```

## And join with the pumas data frame to get the PUMA names

```{r}
ipums_food_production_selected_occupations <-
  ipums_food_production_selected_occupations %>%
  left_join(pumas, by = "puma_code") %>%
  mutate(puma = paste(puma, state, sep = ", ")) %>%
  select(1:3,
         13,
         4:12)

ipums_food_production_selected_occupations
```

```{r}
citizenship_lookup = tibble(
  citizenship_code = factor(c(
    0, 1, 2, 3, 4, 5
  )),
  citizenship = c(
    "native_citizen", # n/a
    "native_citizen", # born_abroad_citizen
    "naturalized_citizen", # naturalized_citizen
    "non_citizen", # non_citizen
    "non_citizen_with_papers", # non_citizen_with_papers
    "foreign_born_status_unknown"
  )
)

ethnicity_lookup = tibble(
  ethnicity_code = factor(c(
    0, 1, 2, 3, 4, 9
  )),
  ethnicity = c(
    "non_hispanic",
    "hispanic", # mexican
    "hispanic", # puerto_rican
    "hispanic", # cuban
    "hispanic", # other
    "not_reported"
  )
)

ipums_food_production_selected_occupations_mapped <-
  ipums_food_production_selected_occupations %>%
  select(-ethnicity,-citizenship) %>%
  left_join(ethnicity_lookup, by="ethnicity_code") %>%
  left_join(citizenship_lookup, by="citizenship_code")
```

```{r}
pct <- function (part,whole) {
  round(part / whole * 100, digits = 0)
}
```


```{r}
calc_percents <- function (rows) {
  rows %>%
  mutate(
    total_native_citizen = (
      total_native_citizen_non_hispanic +
      total_native_citizen_hispanic
    ),
    total_naturalized_citizen = (
      total_naturalized_citizen_non_hispanic +
      total_naturalized_citizen_hispanic
    ),
    total_non_citizen = (
      total_non_citizen_non_hispanic +
      total_non_citizen_hispanic
    ),
    total_foreign_non_hispanic = (
      total_naturalized_citizen_non_hispanic +
      total_non_citizen_non_hispanic
    ),
    total_foreign_hispanic = (
      total_naturalized_citizen_hispanic +
      total_non_citizen_hispanic
    ),
    total_foreign = (
      total_foreign_non_hispanic +
      total_foreign_hispanic
    ),
    total_employees = (
        total_native_citizen +
        total_naturalized_citizen +
        total_non_citizen
    ),
    pct_native_citizen_non_hispanic = 
      pct(total_native_citizen_non_hispanic, total_employees),
    pct_native_citizen_hispanic =
      pct(total_native_citizen_hispanic, total_employees),
    pct_total_foreign_non_hispanic =
      pct(total_foreign_non_hispanic, total_employees),
    pct_total_foreign_hispanic =
      pct(total_foreign_hispanic, total_employees),
    pct_naturalized_citizen_non_hispanic =
      pct(total_naturalized_citizen_non_hispanic, total_employees),
    pct_naturalized_citizen_hispanic =
      pct(total_naturalized_citizen_hispanic, total_employees),
    pct_non_citizen_non_hispanic =
      pct(total_non_citizen_non_hispanic, total_employees),
    pct_non_citizen_hispanic =
      pct(total_non_citizen_hispanic, total_employees),
    pct_native_citizen =
      pct(total_native_citizen, total_employees),
    pct_total_foreign =
      pct(total_foreign, total_employees),
    pct_naturalized_citizen =
      pct(total_naturalized_citizen, total_employees),
    pct_non_citizen =
      pct(total_non_citizen, total_employees)
  )
}
```


```{r}
ipums_food_production_selected_occupations_mapped
```

## Calculate the national-level citizenship and ethnicity data for these industry-occupation pairs

```{r}
ipums_food_production_by_industry_occupation <-
  ipums_food_production_selected_occupations_mapped %>%
  group_by(industry_code,
           industry,
           occupation_code,
           occupation,
           citizenship,
           ethnicity) %>%
  summarize(num_employees = sum(PERWT)) %>%
  pivot_wider(
    id_cols = c(industry_code, industry, occupation_code, occupation),
    names_from = c(citizenship, ethnicity),
    values_from = num_employees,
    names_prefix = "total_"
  ) %>%
  replace(is.na(.), 0) %>%
  calc_percents()
```

```{r}
ipums_food_production_by_industry_occupation
```

### Export the data

```{r}
write_csv(
  ipums_food_production_by_industry_occupation,
  "data/exported/ipums_food_production_by_industry_occupation.csv"
)
```

### Group the data by PUMA

```{r}
ipums_food_production_by_puma <-
  ipums_food_production_selected_occupations_mapped %>%
  group_by(state_code,
           state,
           puma_code,
           puma,
           citizenship,
           ethnicity) %>%
  summarize(num_employees = sum(PERWT, na.rm = TRUE)) %>%
  pivot_wider(
    id_cols = c(state_code, state, puma_code, puma),
    names_from = c(citizenship, ethnicity),
    values_from = num_employees,
    names_prefix = "total_"
  ) %>%
  replace(is.na(.), 0) %>%
  calc_percents()
```

### Export the data

```{r}
write_csv(
  ipums_food_production_by_puma,
  "data/exported/ipums_food_production_by_puma.csv"
)
```

## Group the data by PUMA and industry and occupation

```{r}
ipums_food_production_by_puma_ind_occ <-
  ipums_food_production_selected_occupations_mapped %>%
  group_by(
    state_code,
    state,
    puma_code,
    puma,
    industry_code,
    industry,
    occupation_code,
    occupation,
    citizenship,
    ethnicity
  ) %>%
  summarize(num_employees = sum(PERWT)) %>%
  pivot_wider(
    id_cols = c(
      state_code,
      state,
      puma_code,
      puma,
      industry_code,
      industry,
      occupation_code,
      occupation
    ),
    names_from = c(citizenship, ethnicity),
    values_from = num_employees,
    names_prefix = "total_"
  ) %>%
  replace(is.na(.), 0) %>%
  calc_percents()
```

### Export the data

```{r}
write_csv(
  ipums_food_production_by_puma_ind_occ,
  "data/exported/ipums_food_production_by_puma_ind_occ.csv"
)
```

```{r}
ipums_food_production_by_puma
```

## Join the employment data to the crosswalk data frame

```{r}
ipums_food_production_by_puma
```

```{r crosswalk the employment data from PUMAs to counties}
ipums_food_production_by_puma_crosswalked_to_county <-
  ipums_food_production_by_puma %>%
  left_join(crosswalk, by = "puma_code") %>%
  mutate(
    total_employees = round(
      total_employees * allocation_factor,
      digits = 0
    ),
    total_native_citizen_non_hispanic = round(total_native_citizen_non_hispanic *
                                          allocation_factor, digits = 0),
    total_native_citizen_hispanic = round(total_native_citizen_hispanic *
                                      allocation_factor, digits = 0),
    total_foreign_non_hispanic = round(total_foreign_non_hispanic *
                                         allocation_factor, digits = 0),
    total_foreign_hispanic = round(total_foreign_hispanic *
                                     allocation_factor, digits = 0),
    total_naturalized_citizen_non_hispanic = round(total_naturalized_citizen_non_hispanic
                                             * allocation_factor, digits = 0),
    total_naturalized_citizen_hispanic = round(total_naturalized_citizen_hispanic
                                         * allocation_factor, digits = 0),
    total_non_citizen_non_hispanic = round(total_non_citizen_non_hispanic
                                     * allocation_factor, digits = 0),
    total_non_citizen_hispanic = round(total_non_citizen_hispanic
                                 * allocation_factor, digits = 0),
    total_native_citizen = total_native_citizen_non_hispanic + total_native_citizen_hispanic,
    total_foreign = total_foreign_non_hispanic + total_foreign_hispanic,
    total_naturalized_citizen = total_naturalized_citizen_non_hispanic +
      total_naturalized_citizen_hispanic,
    total_non_citizen = total_non_citizen_non_hispanic + total_non_citizen_hispanic
  ) %>%
  select(1,
         state = state.x,
         3,
         puma = puma.x,
         32:33,
         35:40,
         5:17)
```

### Join the employment data to the crosswalk data frame
```{r crosswalk the employment data from PUMAs and ind-occ pairs to counties and ind-occ pairs}
ipums_food_production_by_puma_ind_occ_crosswalked_to_county <-
  ipums_food_production_by_puma_ind_occ %>%
  left_join(crosswalk, by = "puma_code") %>%
  mutate(
    total_employees = round(
      total_employees * allocation_factor,
      digits = 0
    ),
    total_native_citizen_non_hispanic = round(total_native_citizen_non_hispanic
                                        * allocation_factor, digits = 0),
    total_native_citizen_hispanic = round(total_native_citizen_hispanic
                                    * allocation_factor, digits = 0),
    total_foreign_non_hispanic = round(total_foreign_non_hispanic
                                       * allocation_factor, digits = 0),
    total_foreign_hispanic = round(total_foreign_hispanic
                                   * allocation_factor, digits = 0),
    total_naturalized_citizen_non_hispanic = round(total_naturalized_citizen_non_hispanic
                                             * allocation_factor, digits = 0),
    total_naturalized_citizen_hispanic = round(total_naturalized_citizen_hispanic
                                         * allocation_factor, digits = 0),
    total_non_citizen_non_hispanic = round(total_non_citizen_non_hispanic
                                     * allocation_factor, digits = 0),
    total_non_citizen_hispanic = round(total_non_citizen_hispanic
                                 * allocation_factor, digits = 0),
    total_native_citizen = total_native_citizen_non_hispanic + total_native_citizen_hispanic,
    total_foreign = total_foreign_non_hispanic + total_foreign_hispanic,
    total_naturalized_citizen = total_naturalized_citizen_non_hispanic +
      total_naturalized_citizen_hispanic,
    total_non_citizen = total_non_citizen_non_hispanic + total_non_citizen_hispanic
  ) %>%
  select(1,
         state = state.x,
         3,
         puma = puma.x,
         36:37,
         5:8,
         39:44,
         9:21)
```

## Group the data by county

```{r group the data by county}
# This step is required to deal with situations
# where multiple PUMAs split a single county
ipums_food_production_by_county <-
  ipums_food_production_by_puma_crosswalked_to_county %>%
  group_by(state_code, state, county_code, county) %>%
  calc_percents()
```

### Export the data

```{r}
write_csv(
  ipums_food_production_by_county,
  "data/exported/ipums_food_production_by_county.csv"
)
```

## Group the data by industry and occupation and county

```{r}
ipums_food_production_by_puma_ind_occ_crosswalked_to_county
```

```{r group the data by industry and occupation and county}
# This step is required to deal with situations
# where multiple PUMAs split a single county
ipums_food_production_by_county_ind_occ <-
  ipums_food_production_by_puma_ind_occ_crosswalked_to_county %>%
  group_by(
    state_code,
    state,
    county_code,
    county,
    industry_code,
    industry,
    occupation_code,
    occupation
  ) %>%
  calc_percents()
```


```{r}
ipums_food_production_by_county_ind_occ
```

### Export the data
```{r}
write_csv(
  ipums_food_production_by_county_ind_occ,
  "data/exported/ipums_food_production_by_county_ind_occ.csv"
)
```

# Analyze the data

## What is the total number of workers employed in these industries?

```{r analyze employment by total foreign}
sum(ipums_food_production_by_puma$total_employees)
```

## What is the total number of foreigners employed in these industries?

```{r}
sum(ipums_food_production_by_puma$total_foreign)
```

## What's that as a percent?

```{r}
sum(ipums_food_production_by_puma$total_foreign) /
  sum(ipums_food_production_by_puma$total_employees) * 100
```

## How many of the jobs have a higher proportion of foreigners employed than the national average of 17.1%?

```{r}
ipums_food_production_by_industry_occupation %>%
  filter(pct_total_foreign > 17.1) %>%
  select(
    industry,
    occupation,
    total_employees,
    pct_total_foreign,
    pct_total_foreign_hispanic
  ) %>%
  arrange(desc(total_employees))
```

## What is the total number of Hispanic foreigners employed in these industries?

```{r analyze employment by foreign and Hispanic}
sum(ipums_food_production_by_puma$total_foreign_hispanic)
```

## What's that as a percent of all foreign workers in these industries?

```{r}
sum(ipums_food_production_by_puma$total_foreign_hispanic) /
  sum(ipums_food_production_by_puma$total_foreign) * 100
```

## Which PUMAs have the highest number of immigrant foodworkers?

```{r analyze employment by puma}
ipums_food_production_by_puma %>%
  arrange(desc(total_foreign)) %>%
  select(
    state,
    puma,
    total_employees,
    total_foreign,
    total_foreign_hispanic
  ) %>%
  head(23) # 2,341 PUMAs / 100 = 23.41
```

## And how many immigrant workers are in these PUMAs?

```{r}
ipums_food_production_by_puma %>%
  arrange(desc(total_foreign)) %>%
  select(puma,
         total_employees,
         total_foreign,
         total_foreign_hispanic) %>%
  head(23) %>%
  ungroup() %>%
  summarize(sum(total_foreign))
```

## And how many of these foreign-born workers are Hispanic?

```{r}
ipums_food_production_by_puma %>%
  arrange(desc(total_foreign)) %>%
  select(puma,
         total_employees,
         total_foreign,
         total_foreign_hispanic) %>%
  head(23) %>%
  ungroup() %>%
  summarize(sum(total_foreign_hispanic))
```

## Which counties have the highest number of immigrant foodworkers?

```{r analyze employment by county}
ipums_food_production_by_county %>%
  arrange(desc(total_foreign)) %>%
  select(
    state,
    county,
    total_employees,
    total_foreign,
    total_foreign_hispanic
  ) %>%
  head(31) # 3,142 counties / 100 = 31.42
```

## And how many immigrant workers are in these counties?

```{r}
ipums_food_production_by_county %>%
  arrange(desc(total_foreign)) %>%
  select(
    county,
    total_employees,
    total_foreign,
    total_foreign_hispanic
  ) %>%
  head(31) %>%
  ungroup() %>%
  summarize(sum(total_foreign))
```

## And how many of these foreign-born workers are Hispanic?

```{r}
ipums_food_production_by_county %>%
  arrange(desc(total_foreign)) %>%
  select(
    county,
    total_employees,
    total_foreign,
    total_foreign_hispanic
  ) %>%
  head(31) %>%
  ungroup() %>%
  summarize(sum(total_foreign_hispanic))
```

## What is the distribution of foreign workers in these industries by state?

```{r analysis employment and COVID status by state}
foodworkers_by_state <- ipums_food_production_by_puma %>%
  group_by(state) %>%
  summarize(
    state_total_employment_selected_industries =
      sum(total_employees),
    state_total_foreign = sum(total_foreign),
    state_pct_total_foreign = round(
      sum(total_foreign) / sum(total_employees) * 100
    ),
    2,
    state_total_foreign_hispanic = sum(total_foreign_hispanic),
    state_pct_total_foreign_hispanic = round(
      sum(total_foreign_hispanic) /
        sum(total_employees) * 100
    ),
    2
  ) %>%
  mutate(
    rank_total_employment = dense_rank(desc(
      state_total_employment_selected_industries
    )),
    rank_total_foreign = dense_rank(desc(state_total_foreign)),
    rank_total_foreign_hispanic = dense_rank(desc(state_total_foreign_hispanic))
  ) %>%
  select(
    state,
    state_total_employment_selected_industries,
    state_total_foreign,
    state_pct_total_foreign,
    state_total_foreign_hispanic,
    state_pct_total_foreign_hispanic,
    rank_total_employment,
    rank_total_foreign,
    rank_total_foreign_hispanic
  ) %>%
  arrange(desc(state_total_employment_selected_industries))

foodworkers_by_state
```

## What proportion of total foodworker employment do the top 10 states account for?

```{r}
foodworkers_by_state %>%
  arrange(desc(state_total_employment_selected_industries)) %>%
  head(10) %>%
  summarize(sum(state_total_employment_selected_industries))
```

### What is the distribution of foreign workers in these industries by state?

```{r}
foodworkers_by_state %>%
  summarize(sum(state_total_employment_selected_industries))
```

## Export the states

```{r}
write_csv(foodworkers_by_state,
          "data/exported/foodworkers_by_state.csv")
```