---
title: "Workers vulnerable to health risks"
author: "Robby"
date: "April 16, 2020"
output:
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(ipumsr)
#library(survey) # Some easy weighting functions.
library(srvyr) # Like survey but with syntax like dplyr.
library(labelled) # for working with labelled data
library(tidycensus)
library(ggrepel)
library(readxl)


file.path <- "C:/Users/Robert Habans/OneDrive - Nonprofit Knowledge Works/covid/Rproject/covid-19"
source(paste0(file.path,"/03_graphics themes.R"))

```

## Objective 

Use IPUMS to recreate similar analyses to the following about the workers vulnerable to COVID-19 exposure. 

We're basing this on a few examples:

- [CEPR](https://cepr.net/a-basic-demographic-profile-of-workers-in-frontline-industries/) -- This also has a [Github repo](https://github.com/ceprdata/frontline-workers). 
- [Boston Indicators](https://www.bostonindicators.org/article-pages/2020/april/frontline_workers)
- [NYC Comptroller](https://comptroller.nyc.gov/reports/new-york-citys-frontline-workers/)
- [ACLU](https://www.aclum.org/en/publications/data-show-covid-19-hitting-essential-workers-and-people-color-hardest)

## IPUMS setup

Load the raw IPUMS data. 

```{r, include=FALSE}

# Note that you can pass in the loaded DDI into the `read_ipums_micro()`
# NOTE: To load data, you must download both the extract's data and the DDI
# and also set the working directory to the folder with these files (or change the path below).

if (!require("ipumsr")) stop("Reading IPUMS data into R requires the ipumsr package. It can be installed using the following command: install.packages('ipumsr')")

ddi <- read_ipums_ddi("inputs/ipums/usa_00007.xml")
ipums.la <- read_ipums_micro(ddi)
```

Clean some of the variables. 

```{r, include=FALSE}
# Note this list of variables converted to factors may be incomplete.

# to_factor pulls the labels out and makes them factors.
# This is an alternative to the hybrid int+lbl default. 
ipums.la.fac <- ipums.la %>% 
  mutate_at(vars(NUMPREC:SUBSAMP,
                 HHTYPE,
                 STATEFIP:PUMA,
                 GQ:MORTGAGE,
                 TAXINCL:INSINCL,
                 FOODSTMP,
                 SSMC:MULTGEND,
                 FAMUNIT:SCHOOL,
                 CLASSWKR,
                 EDUC:LABFORCE,
                 WKSWORK2:WORKEDYR,
                 MIGRATE1:MIGPUMA1,
                 DIFFREM:GCRESPON), 
            to_factor) %>% 
  mutate(OCC2010.fac = to_factor(OCC2010), # to have factor and numeric versions of occ and ind codes
         OCCSOC.fac = to_factor(OCCSOC),
         IND.fac = to_factor(IND)) #%>% 
  #mutate_at(vars(INCTOT:INCEARN), to_character) # this takes the labels for the income variables but makes them character

# clear the original object to save memory.
rm(ipums.la)

```

Fix income variables using CPI adjustment (available [here](https://cps.ipums.org/cps/cpi99.shtml)). Convert dollar values to numeric and exclude certain codes (Note that it's strange that INCEARN doesn't have special codes). 

```{r, include=FALSE}
# Use the included CPI99, then bring to 2018
cpi99to2018 <- 1.507

ipums.la.fac <- ipums.la.fac %>% 
  mutate(INCWAGE.cpi = ifelse(INCWAGE %in% c(999999 ,999998), NA, CPI99*as.numeric(INCWAGE)*cpi99to2018)) %>% 
  mutate(INCEARN.cpi = CPI99*as.numeric(INCEARN)*cpi99to2018) %>% 
  mutate(INCTOT.cpi = ifelse(INCTOT %in% c(9999999), NA, CPI99*as.numeric(INCTOT)*cpi99to2018)) %>% 
  mutate(FTOTINC.cpi = ifelse(FTOTINC %in% c(9999999), NA, CPI99*as.numeric(FTOTINC)*cpi99to2018)) %>% 
  mutate(HHINCOME.cpi = ifelse(HHINCOME %in% c(9999999), NA, CPI99*as.numeric(HHINCOME)*cpi99to2018)) 
```



For convenience, define PUMA codes for querying by geography. Maps with the exact definitions for 2010 PUMAs are on the [IPUMS website](https://usa.ipums.org/usa/volii/pumas10.shtml).

```{r, include=FALSE}
## Create filters for PUMA codes -- these probably changed in 1990 and 2000.
# ipums provides the best documentation. The fed websites are becoming unusable. 

super.region.pumacodes.2010 <- c(
  "1400",	# Capital Region Planning Commission 1--Baton Rouge Metropolitan Area (West) PUMA
  "1500",	# Capital Region Planning Commission 2--East Baton Rouge Parish (South) PUMA
  "1501",	# Capital Region Planning Commission 3--Baton Rouge City PUMA
  "1502",	# Capital Region Planning Commission 4--East Baton Rouge Parish (North) PUMA
  "1600",	# Capital Region Planning Commission 5--Ascension Parish PUMA
  "1700",	# Capital Region Planning Commission 6--Livingston & St. Helena Parishes PUMA
  "1800",	# Capital Region Planning Commission 7--Tangipahoa & Washington Parishes PUMA
  "1900",	# South Central Planning & Development Commission 1--River Parishes PUMA
  "2000",	# South Central Planning & Development Commission 2--Lafourche & Assumption Parishes PUMA
  "2100",	# South Central Planning & Development Commission 3--Terrebonne Parish PUMA
  "2200",	# Regional Planning Commission 1--St. Tammany Parish (Northwest) PUMA
  "2201",	# Regional Planning Commission 2--St. Tammany Parish (Southeast) PUMA
  "2300",	# Regional Planning Commission 3--Jefferson Parish (North) PUMA
  "2301",	# Regional Planning Commission 4--Jefferson Parish (Central) PUMA
  "2302",	# Regional Planning Commission 5--Jefferson Parish (West Bank) PUMA
  "2400",	# Regional Planning Commission 6--New Orleans City (Northeast) PUMA
  "2401",	# Regional Planning Commission 7--New Orleans City (Central) PUMA
  "2402",	# Regional Planning Commission 8--New Orleans City (South) PUMA
  "2500"	# Regional Planning Commission 9--St. Bernard, Jefferson (South) & Plaquemines Parishes PUMA
)
  
orleans.pumacodes <- c("2400","2401","2402")

# This should be double checked.
NOmsa.pumacodes <- c("2400","2401","2402", # Orleans
                     "2500", # St Bernard, Plaquemines, south Jefferson
                     "2300","2301","2302", # Jefferson
                     "2200","2201", # St Tammany
                     "1900") # river parishes (Includes St. James, I think)

# Define asian categories
asian.ipums <- c("Chinese",
                 "Japanese",
                 "Other Asian or Pacific Islander")

```

Define industry sector based on ranges of variables. For industry and industry sector labels, you could use the structure of the original census document without too much trouble. You can either use a modified custom groupings like the [six categories used by CEPR](https://cepr.net/a-basic-demographic-profile-of-workers-in-frontline-industries/) or use the [ACS codes for NAICS-equivalent sectors](https://www.census.gov/topics/employment/industry-occupation/guidance/code-lists.html).


```{r}
url <- "https://www2.census.gov/programs-surveys/demo/guidance/industry-occupation/2017-industry-code-list-with-crosswalk.xlsx"
destfile <- "X2017_industry_code_list_with_crosswalk.xlsx"
curl::curl_download(url, destfile)
ACS_industry_codes <- read_excel(destfile, 
    skip = 5)

names(ACS_industry_codes) <- c("sector",
                                                  "industry",
                                                  "DELETE",
                                                  "ACScode2017",
                                                  "NAICScode2017")

ACS_industry_codes <- ACS_industry_codes %>% select(-DELETE) 

# THIS DOESN'T WORK
#ACS_industry_codes.range <- 
#   ACS_industry_codes %>% 
#   separate(ACScode2017, c("begin","end"), "-") %>% 
#   mutate(end = ifelse(industry == "Management of companies and enterprises
# ", "7570", end)) %>% 
#   mutate(end = ifelse(sector == "Construction", "0770", end)) %>% 
#   mutate(industry = ifelse(begin=="4670" & end=="5790", "Retail Trade", industry),
#          industry = ifelse(begin=="4070" & end=="4590", "Wholesale Trade", industry),
#          industry = ifelse(begin=="1070" & end=="3990", "Manufacturing", industry)) %>% 
#   filter(!is.na(end)) %>% print(n=100)


```

Or just define them manually. Use this to join to the IPUMS data set. 

```{r}
ACS_industry_codes.range <- tibble(ACScode.num = 1:9590, industry = NA) %>% 
  mutate(industry = ifelse(ACScode.num %in% 170:290, "Agriculture, Forestry, Fishing, and Hunting", industry),
         industry = ifelse(ACScode.num %in% 370:490, "Mining, Quarrying, and Oil and Gas Extraction", industry),
         industry = ifelse(ACScode.num %in% 770, "Construction", industry),
         industry = ifelse(ACScode.num %in% 1070:3990, "Manufacturing", industry),
         industry = ifelse(ACScode.num %in% 4070:4590, "Wholesale Trade", industry),
         industry = ifelse(ACScode.num %in% 4670:5790, "Retail Trade", industry),
         industry = ifelse(ACScode.num %in% 6070:6390, "Transportation and Warehousing", industry),
         industry = ifelse(ACScode.num %in% 570:690, "Utilities", industry),
         industry = ifelse(ACScode.num %in% 6470:6780, "Information", industry),
         industry = ifelse(ACScode.num %in% 6870:6992, "Finance and Insurance", industry),
         industry = ifelse(ACScode.num %in% 7071:7190, "Real Estate and Rental and Leasing", industry),
         industry = ifelse(ACScode.num %in% 7270:7490, "Professional, Scientific, and Technical Services", industry),
         industry = ifelse(ACScode.num %in% 7570, "Management of companies and enterprises", industry),
         industry = ifelse(ACScode.num %in% 7580:7790, "Administrative and support and waste management services", industry),
         industry = ifelse(ACScode.num %in% 7860:7890, "Educational Services", industry),
         industry = ifelse(ACScode.num %in% 7970:8470, "Health Care and Social Assistance", industry),
         industry = ifelse(ACScode.num %in% 8561:8590, "Arts, Entertainment, and Recreation", industry),
         industry = ifelse(ACScode.num %in% 8660:8690, "Accommodation and Food Services", industry),
         industry = ifelse(ACScode.num %in% 8770:9290, "Other Services, Except Public Administration", industry),
         industry = ifelse(ACScode.num %in% 9370:9590, "Public Administration", industry),
         industry = ifelse(ACScode.num %in% 9670:9870, "Military", industry)) %>% 
  mutate(ACScode = as.character(ACScode.num))
  #mutate(ACScode = str_pad(ACScode.num, side = "left", pad = "0", width = 4)) 

```





## Identifying industries

As we assemble it, I'm storing a version of this on my personal drive, then copying into a CSV in the inputs folder. Currently, this is just a test version using the BI/CEPR occupations. 

Right now, this is based on including the A and B industries. Childcare has also been included but should be handled separately. 

```{r}
# vulnerableEssentialWorkers <- read_csv("inputs/acsIndOcc/vulnerableEssentialWorkers.csv", 
#     col_types = cols(ACScode = col_character())) %>% 
#   filter(Include %in% c("A","B")) %>% 
#   select(-ACSnaics) #%>% 
#   #mutate(ACScode = str_pad(ACScode, side = "left", pad = "0", width = 4))

vulnerableEssentialWorkers <- read_csv("inputs/acsIndOcc/vulnerableEssentialWorkers_v2.csv",
    col_types = cols(ACScode = col_character())) %>% 
  filter(Include %in% c("A","B"))

vulnerableEssentialWorkers %>% count(Include)
vulnerableEssentialWorkers %>% count(category2)
vulnerableEssentialWorkers %>% count(category5)
vulnerableEssentialWorkers %>% count(categoryCEPR)

```

To identify the appropriate categories, let's view earnings by the different ways of slicing industries. 

```{r}
# 2 categories
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  filter(CLASSWKR == "Works for wages") %>% 
  left_join(vulnerableEssentialWorkers, by=c("IND.fac"="ACScode")) %>% 
  mutate(essential = ifelse(Include %in% c("A","B"), "Yes", "No"),
         category2 = ifelse(is.na(category2), "Other industries", category2)) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(essential, category2, SEX, vartype = "se") %>% 
  group_by(category2) %>% 
  mutate(pct = n/sum(n))

ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  filter(CLASSWKR == "Works for wages") %>% 
  left_join(vulnerableEssentialWorkers, by=c("IND.fac"="ACScode")) %>% 
  mutate(essential = ifelse(Include %in% c("A","B"), "Yes", "No"),
         category2 = ifelse(is.na(category2), "Other industries", category2)) %>% 
  mutate(race.WOC = ifelse(RACE == "White" & HISPAN == "Not Hispanic", "White", "WOC")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(essential, category2, race.WOC, vartype = "se") %>% 
  group_by(category2) %>% 
  mutate(pct = n/sum(n))

# 5 categories
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  filter(CLASSWKR == "Works for wages") %>% 
  left_join(vulnerableEssentialWorkers, by=c("IND.fac"="ACScode")) %>% 
  mutate(essential = ifelse(Include %in% c("A","B"), "Yes", "No"),
         category5 = ifelse(is.na(category5), "Other industries", category5)) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(essential, category5, SEX, vartype = "se") %>% 
  group_by(category5) %>% 
  mutate(pct = n/sum(n))

ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  filter(CLASSWKR == "Works for wages") %>% 
  left_join(vulnerableEssentialWorkers, by=c("IND.fac"="ACScode")) %>% 
  mutate(essential = ifelse(Include %in% c("A","B"), "Yes", "No"),
         category5 = ifelse(is.na(category5), "Other industries", category5)) %>% 
  mutate(race.WOC = ifelse(RACE == "White" & HISPAN == "Not Hispanic", "White", "WOC")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(essential, category5, race.WOC, vartype = "se") %>% 
  group_by(category5) %>%
  mutate(pct = n/sum(n))


```




## Perform queries and make figures. 

First, let's see what the sample size is.

- Should I include employed and W/S workers, or also self-employed?
- Should I use worked last week WRKLSTWK? Or employed EMPSTAT? [IPUMS documentation](https://usa.ipums.org/usa-action/variables/WRKLSTWK#description_section) suggests that the two variables are closely related, though CEPR specifies that they use weeks worked last week. CEPR might have simply replicated the EMPSTAT variable from IPUMS. See [line 59](https://github.com/ceprdata/frontline-workers/blob/master/frontline-workers_acs-recode.do). I think EMPSTAT is fine. 
- Should I exclude military? This could be done with EMPSTATD or Industry.

```{r, echo=FALSE}
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  filter(IND %in% vulnerableEssentialWorkers$ACScode)

```

What is the overall percentage of the workforce when using the A and B groups?

```{r}
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  left_join(vulnerableEssentialWorkers, by=c("IND.fac"="ACScode")) %>% 
  mutate(include.group = ifelse(!is.na(Include), Include, "other")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(include.group, vartype = "se") %>% 
  mutate(pct = n/sum(n))

```



Industry sector of essential workers?

```{r}
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(CLASSWKR == "Works for wages") %>% 
  #filter(WRKLSTWK == "Worked") %>% 
  filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(IND.char = as.character(IND)) %>% 
  left_join(ACS_industry_codes.range, by=c("IND.char"="ACScode")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(industry, vartype = "se") %>% 
  mutate(pct = n/sum(n))


```

Median annual w/s earnings, family earnings, and below poverty level. Did Brookings only use w/s earnings? This would exclude self-employed. 

```{r}
# Poverty
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(CLASSWKR == "Works for wages") %>% 
  #filter(WRKLSTWK == "Worked") %>% 
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  filter(POVERTY != 0) %>% 
  mutate(pov200 = ifelse(POVERTY < 200, "Below 200%", "200% or higher")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  survey_count(pov200, vartype = "se") %>% 
  mutate(pct = n/sum(n))

# Earnings -- note that this gives a weird vcov error.
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed",
         CLASSWKR == "Works for wages") %>% # employed workers only
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  summarise(median.incwage = survey_median(INCEARN.cpi, vartype = "se"))

# Total family income -- need to restrict to family?
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed",
         CLASSWKR == "Works for wages") %>% # employed workers only
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  filter(FTOTINC != 9999999) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  summarise(median.faminc = survey_median(FTOTINC.cpi, vartype = "se"))

```


Rent burden of the essential workers. 

```{r}
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  filter(OWNERSHP == "Rented") %>% 
  filter(HHINCOME != 9999999) %>% 
  mutate(rent.burden = ifelse(RENTGRS/(HHINCOME/12) > .3, "Yes", "No")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  survey_count(rent.burden, vartype = "se") %>% 
  mutate(pct = n/sum(n))

  
```

Race and gender breakdowns of essential workers. 

```{r}
# Race
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(race.simple = "Other",
         race.simple = ifelse(HISPAN != "Not Hispanic", "Latino", race.simple),
         race.simple = ifelse(RACE == "White" & HISPAN == "Not Hispanic", "White", race.simple),
         race.simple = ifelse(RACE == "Black/African American/Negro" & HISPAN == "Not Hispanic", "Black", race.simple),
         race.simple = ifelse(RACE %in% asian.ipums & HISPAN == "Not Hispanic", "AAPI", race.simple)) %>% 

  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(race.simple, vartype = "se") %>% 
  mutate(pct = n/sum(n))

# Gender
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(SEX, vartype = "se") %>% 
  mutate(pct = n/sum(n))

```

Children and elderly in household and age of worker. 

- What is the appropriate cutoff for children? Should I use NCHILD or NCHLT5 instead of the presence of children?
- What is the appropriate cutoff for old people? Jenna used 60 on the demographics page.

```{r}
# Children
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  mutate(children = "No children", 
         children = ifelse(as.numeric(YNGCH) < 10, 
                           "Children younger than 10", 
                           children),
         children = ifelse(as.numeric(YNGCH) >= 10 & as.numeric(YNGCH) <18, 
                           "Children between 10 and 18", 
                           children)) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  survey_count(children, vartype = "se") %>% 
  mutate(pct = n/sum(n))

# Older in household
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  group_by(YEAR, SERIAL) %>% 
  mutate(max.age.in.hh = max(as.numeric(AGE)),
         old.in.hh = ifelse(max.age.in.hh >= 60, "Yes", "No")) %>% 
  ungroup() %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  survey_count(old.in.hh, vartype = "se") %>% 
  mutate(pct = n/sum(n))

  
# Older workers
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  mutate(old = "Not old",
         old = ifelse(as.numeric(AGE) >= 60, "old", old)) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  survey_count(old, vartype = "se") %>% 
  mutate(pct = n/sum(n))

# Are the workers differently aged?
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  summarise(mean.age = survey_mean(as.numeric(AGE), vartype = "se"),
            median.age = survey_median(as.numeric(AGE), vartype = "se")) 



```

Health insurance.

```{r}
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  survey_count(HCOVANY, vartype = "se") %>% 
  mutate(pct = n/sum(n))


```

Use of public transportation. 

```{r}
ipums.la.fac %>% 
  filter(PUMA %in% NOmsa.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  #filter(IND %in% vulnerableEssentialWorkers$ACScode) %>% 
  mutate(essential = ifelse(IND %in% vulnerableEssentialWorkers$ACScode, "Essential", "Other")) %>% 
  filter(TRANWORK != "N/A") %>% 
  mutate(takes.transit = ifelse(TRANWORK %in% c("Bus or trolley bus",
                                                 "Streetcar or trolley car",
                                                 "Subway or elevated",
                                                 "Railroad",
                                                 "Ferryboat"),
                                 "Yes", "No")) %>% 
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(essential) %>% 
  survey_count(takes.transit, vartype = "se") %>% 
  mutate(pct = n/sum(n))

```


