---
title: "Vulnerable Workers"
author: "Robby"
date: "March 27, 2020"
output:
  html_document:
    df_print: paged
  word_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 

To get up and running with IPUMS.

## Get the data

You need to go to the ipums website to get the data. Login with the credentials and click on My Data. Follow the instructions [here](https://cran.r-project.org/web/packages/ipumsr/vignettes/ipums.html).

### About the data extract

As of 3/27/2020, the data extract is as follows: 

- ACS from 2012-2018. That's seven years. Usually, we do 5-year samples but I grabbed two more in case we want to increase the sample size for some breakdowns. 
- Entire state of Louisiana. We will use PUMA codes to filter down to the region or its parishes as necessary/possible. Some PUMAs cover multiple parishes.
- A wide assortment of housing, demographic, economic, and workforce variables. 


### Load and clean

Basically, you load a "ddi" xml file that has the codes and metadata, and then load in the .dat file that actually has the data. If you open the DDI file in a browser, it's a handy codebook. For some reason, it doesn't work when you open it locally, but if you open it from IPUMS website, you can just keep it open in your browser for reference. 

```{r, eval=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_00005.xml")
ipums.la <- read_ipums_micro(ddi)
```

This pulls in *labelled* data. The raw data is just numeric codes for different characteristics, but the labels make it easier because you don't have to go back. It's a bit of a trade-off though. If memory serves, the labelled fields are a bit weird and sometimes don't play nicely with dplyr or other common packages. 

Let's clean up the labels to make them easier to use. The default is a "labelled" -- it includes both the original codes (which you can find in the DDI file data dictionary) and the labelles attached the code. For example, for HHTYPE, it will include the code "1" and the label "Married-couple family household" and all of the other code/label pairs. As a convenience I use `labelled::to_factor` to convert the labels to factors. This way, it just shows the verbose labels like "Married-couple family household". More importantly, I've had issues with labelled data causing conflicts with dplyr or srvyr.

Also, be aware that some of the continuous variables have special codes. For example, the income variables have top coding, so they're not pure double/integer variables. 

```{r, eval=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,
                 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:INCBUS00), to_character) # this takes the labels for the income variables but makes them character

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

```



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, eval=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)

```



## Examples

Here are some examples using srvyr. `as_survey_design` and `survey_mean` are from the [srvyr package](https://cran.r-project.org/web/packages/srvyr/srvyr.pdf). The command `as_survey_design` creates a tbl_svy survey object. 

[This link](https://usa.ipums.org/usa/complex_survey_vars/userNotes_variance.shtml) describes the STRATA and CLUSTER variables in IPUMS. Here is the important part:

> Three Technical Variables Are Needed for Analysis of the IPUMS Data:

> 1. A sampling weight must be chosen (PERWT, HHWT, or SLWT), based on the type of question a researcher is trying to answer. Analysts should review variable descriptions of the variables of interest and the "Sampling Weights" note for more information about which weight to use.
2. STRATA is an integrated variable that represents the impact of the sample design stratification on the estimates of variance and standard errors.
3. CLUSTER is an integrated variable which uniquely identifies each household record in a given sample.

As an alternative, replicate weights can be used as described [here](https://usa.ipums.org/usa/repwt.shtml). The current extract does not include the replicate weight variables, but these can easily be added, and the command `srvyr::as_survey_rep()` could be used to prep the data frame. 



```{r, eval=FALSE}
## Median income
# Find codes for employment status
ipums.la.fac %>% count(EMPSTAT)

# Total income for employed workers.
ipums.la.fac %>% 
  filter(PUMA %in% orleans.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  summarise(mean.income = survey_mean(INCTOT, vartype = "se"))

# vartype can be "se", "ci", "var", or "cv".

## Race
# Estimate race breakdown using srvyr
ipums.la.fac %>% 
  filter(PUMA %in% orleans.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  survey_count(RACE, vartype = "se") %>% 
  mutate(pct = n/sum(n))

# Same thing with normal dplyr functions -- same results, no variance.
ipums.la.fac %>% 
  filter(PUMA %in% orleans.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  count(RACE, wt=PERWT) %>% 
  mutate(pct = n/sum(n))

## Income by race
ipums.la.fac %>% 
  filter(PUMA %in% orleans.pumacodes, # filter to Orleans
         YEAR %in% as.character(2014:2018)) %>%  # filter by sample year
  filter(EMPSTAT == "Employed") %>% # employed workers only
  as_survey_design(ids = CLUSTER,
                   strata = STRATA, 
                   weights = PERWT) %>% 
  group_by(RACE) %>% 
  summarise(mean.income = survey_mean(INCTOT, vartype = "se"))

```

### Additional concerns

- Use the appropriate weight: PERWT for individuals and HHWT for households. 
- The difference between household and individual data is very important. I'm not exactly sure how this works with CLUSTER, but you wouldn't want to count a household with 5 people 5 times when calculating the median household income. 
- You might need a weighted median/weighted quantile function when calculating medians (e.g., Hmisc). survey_median and survey_quantile are probably preferable. 
- When you use weights on a five year sample, you will get a five year estimate, so you will have to divide by 5. If you only want a median income or a percentage, you can ignore this issue. 



## Load data from different sources

### ONET

First, I need to import the O*Net data from the [website](https://www.onetonline.org/). You could download the data, or try the R package. I tried signing up for the developer API, so let's see if that works -- It's being reviewed. 

```{r, eval=FALSE}
library(ONETr)
```

Otherwise, you can download the files [here](https://www.onetcenter.org/database.html#all-files) and [here for individual files](https://www.onetcenter.org/database.html#individual-files), but it's a relational database, so it might be pretty complicated. 




### EMSI data

Pull in EMSI data. You grabbed all NAICS/SOC levels, and employees and employees plus self-employed for Occupations and Industries. 

**Note**: Brookings did 2019 data from EMSI. I downloaded 2018 because that was the default. 

```{r}
Industry_Table_allLevels_NOMSA_Employees <- read_csv("inputs/emsi/Industry_Table_allLevels_NOMSA_Employees.csv", 
    col_types = cols(Level = col_character(), 
        NAICS = col_character()))

Occupation_Table_allLevels_NOMSA_Employees <- read_csv("inputs/emsi/Occupation_Table_allLevels_NOMSA_Employees.csv", 
    col_types = cols(Level = col_character(), 
        SOC = col_character()))
```

### O*Net analysis from St. Louis Fed

The St Louis Fed's analysis of [Contact-intensive occupations](https://www.stlouisfed.org/on-the-economy/2020/march/social-distancing-contact-intensive-occupations) is available for download. This includes some 

```{r}
url <- "https://www.stlouisfed.org/~/media/Files/Excel/BlogTables_SocialDistancingContactIntensiveOccupations.xlsx?la=en"
destfile <- "BlogTables_SocialDistancingContactIntensiveOccupations.xlsx"
curl::curl_download(url, destfile)
BlogTables_SocialDistancingContactIntensiveOccupations <- read_excel(destfile, sheet = "Appendix")

ContactIntensiveOccs <- BlogTables_SocialDistancingContactIntensiveOccupations %>% 
  rename_at(vars(...9:...39), ~str_replace(., "...","subocc")) %>% 
  rename_at(vars(subocc9:subocc39),
            ~str_replace(., str_extract(.,"[[:digit:]]+"),
                         as.character(as.numeric(str_extract(., "[[:digit:]]+")) - 7))) %>% 
  rename(subocc1 = `Suboccupations included (at 6-digit SOC 2010 level), scroll to the right for the full list`) 



```

Now, join this to the EMSI data. Nope, this looks whack, I'm not going to do it. 

```{r}
Occupation_Table_allLevels_NOMSA_Employees %>% 
  filter(Level == 3) %>% View()


```



### Brookings/Moody's Industries, version 1

Here are the brookings industries used in [this analysis](https://www.brookings.edu/blog/the-avenue/2020/03/17/the-places-a-covid-19-recession-will-likely-hit-hardest/). I compared their numbers to the actual data in EMSI, and these NAICS codes check out.

```{r}
BrookingsMoodysAtriskIndustries <- tibble(NAICS = c("21", 
                                     "481","483","484","485","486","487","488",
                                     "5613",
                                     "5615",
                                     "71","72")) %>% 
  mutate(BrMooSector = NA,
         BrMooSector = ifelse(NAICS ==  "21", "Mining/Oil & Gas", BrMooSector),
         BrMooSector = ifelse(NAICS %in%  c("481","483","484","485","486","487","488"), 
                              "Transportation", BrMooSector),
         BrMooSector = ifelse(NAICS ==  "5613", "Employment Services", BrMooSector),
         BrMooSector = ifelse(NAICS ==  "5615", "Travel Arrangements", BrMooSector),
         BrMooSector = ifelse(NAICS %in% c("71","72"), "Leisure & Hospitality", BrMooSector))
```

Then use this to filter on the EMSI data. They call these industries at highest risk. Note that the total at-risk industry percentage is a little less than Brookings -- maybe they exclude public sector?

```{r}
# Total jobs. This checks out with the sum row of the EMSI Industry Table. 
TotalEmployment.NOMSA <- Industry_Table_allLevels_NOMSA_Employees %>% 
  filter(Level == 2) %>% 
  summarise(total = sum(`2019 Jobs`)) %>% 
  pull()

# Show employment in the affected sectors
Industry_Table_allLevels_NOMSA_Employees %>% 
  left_join(BrookingsMoodysAtriskIndustries, by="NAICS") %>% 
  filter(!is.na(BrMooSector)) %>% 
  group_by(BrMooSector) %>% 
  summarise(jobs = sum(`2019 Jobs`)) %>% 
  mutate(jobs.p = jobs/TotalEmployment.NOMSA)

# Total share in affected sectors
Industry_Table_allLevels_NOMSA_Employees %>% 
  left_join(BrookingsMoodysAtriskIndustries, by="NAICS") %>% 
  filter(!is.na(BrMooSector)) %>% 
  group_by(BrMooSector) %>% 
  summarise(jobs = sum(`2019 Jobs`)) %>% 
  mutate(jobs.p = jobs/TotalEmployment.NOMSA) %>% 
  summarise(totaljobs = sum(jobs),
            sharejobs = sum(jobs.p))

```

Last thing, show how risk is disproportionate in GNO.

```{r}
# Download the Brookings table. 
url <- "https://www.brookings.edu/wp-content/uploads/2020/03/20200318_BrookingsMetro_Coronavirus_ImpactedMetros_Appendix.xlsx"
destfile <- "X20200318_BrookingsMetro_Coronavirus_ImpactedMetros_Appendix.xlsx"
curl::curl_download(url, destfile)
X20200318_BrookingsMetro_Coronavirus_ImpactedMetros_Appendix <- read_excel(destfile, 
    skip = 1)

# Get metro populations
census_api_key("530ce361defc2c476e5b5d5626d224d8354b9b9a")
#load_variables(2018, "acs1", cache = TRUE) %>% View() # Find vairable for 1 year, 2018

MSApopulations <- get_acs(geography = "metropolitan statistical area/micropolitan statistical area",
        variables = c(TotalPopulation = "B01003_001"),
        year = 2018) %>% 
  filter(str_detect(NAME, " Metro Area")) %>% 
  mutate(`Metro name` = str_replace(NAME, " Metro Area",""))

# Point version
X20200318_BrookingsMetro_Coronavirus_ImpactedMetros_Appendix %>% 
  left_join(MSApopulations, by="Metro name") %>% 
  filter(estimate > 750000) %>% 
  mutate(label = ifelse(`Metro name` == "New Orleans-Metairie, LA",
                        "New Orleans-Metairie, LA", "")) %>% 
  ggplot(aes(factor(1), `Share of metro jobs`)) + 
  geom_point(aes(color = label, size=label), alpha = .3) + 
  geom_text(aes(label = label, angle = 45, hjust = -.1), family="Asap") + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + 
  scale_color_manual(values = c(DCcolor.p1mediumblue, DCcolor.p1darkblue)) + 
  coord_flip() + 
  themeDC_scatter() + 
  theme(legend.position = "none") + 
  labs(x="", y="")

# Bar version -- this looks good.
brookingsAtRiskMetros_rankedbars <- X20200318_BrookingsMetro_Coronavirus_ImpactedMetros_Appendix %>% 
  left_join(MSApopulations, by="Metro name") %>% 
  filter(estimate > 1000000) %>% 
  mutate(label = ifelse(`Metro name` == "New Orleans-Metairie, LA",
                        "New Orleans-Metairie, LA", "")) %>% 
  ggplot(aes(reorder(`Metro name`, -`Share of metro jobs`), `Share of metro jobs`)) + 
  geom_col(aes(fill = label)) + 
  scale_fill_manual(values = c(DCcolor.p1skyblue, DCcolor.p1darkblue)) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + 
  annotate("text", x = "New Orleans-Metairie, LA", y = .25, label = "New Orleans metro", hjust = 0, family = "Asap", color="grey20") + 
  themeDC_scatter() + 
  theme(legend.position = "none",
        axis.text.x = element_blank()) + 
  labs(x="Metros ranked by “high-risk” share", y="",
       #y="Share of metro jobs in “most at risk” industry groups",
       title = "Share of employment in “high-risk” industries by metro (over 1 million population)",
       subtitle = "Brookings analysis")


```

Sources: 
- Employment in "At-risk" industries: [Brookings analysis](https://www.datacenterresearch.org/covid-19-data-and-information/covid-19-economic-analysis/) of EMSI data, based on a research note from Moody's Analytics. 
- Metro populations are from the American Community Survey 1-year estimates for 2018.


### Brookings Industries, version 2

Brookings released an [improved version](https://www.brookings.edu/research/who-are-the-workers-already-impacted-by-the-covid-19-recession/https://www.brookings.edu/research/who-are-the-workers-already-impacted-by-the-covid-19-recession/) of the metro analysis. It also includes demographics of workers in the affected industries. It has a lot of work done. Their tabulations are from ACS IPUMS. They used 1-year PUMS and the unweighted sample size is only 1661 for vulnerable workers. That's not so great. Luckily they provided the data, but perhaps I could replicate this analysis with a five-year sample.

You can just make graphs with their data directly. It will be easy. Come back later to recreate this with a 5-year IPUMS sample. I don't think you should use their local estimates of industries because of the small sample size, so join the categories to EMSI and use that instead. 

First, grab the data -- why these terrible formatted spreadsheets? Note that the Field name "Mining" occurs twice, so they have added "...63" or "...43". Also, there is an extraneous character in `Motor vehicle and parts dealers \r\r\n`.

```{r}
# Download spreadsheet
url <- "https://www.brookings.edu/wp-content/uploads/2020/04/Workers_covid19_appendix_handoff_20200403.xlsx"
destfile <- "Workers_covid19_appendix_handoff_20200403.xlsx"
curl::curl_download(url, destfile)

# All workers stats
Brookings.AllWorkersStats <- read_excel(destfile, 
    skip = 2,
    sheet = "T1. All workers stats") %>% 
  rename(`Motor vehicle and parts dealers` = `Motor vehicle and parts dealers \r\r\n`,
         `Metro name` = `...1`)

# Vulnerable workers stats
Brookings.VulnerableWorkersStats <- read_excel(destfile, 
    skip = 2,
    sheet = "T2. Vulnerable workers stats") %>% 
  select(-...90) %>% 
  rename(`Motor vehicle and parts dealers` = `Motor vehicle and parts dealers \r\r\n`,
         `Metro name` = `...1`)

# Other workers stats
Brookings.OtherWorkersStats <- read_excel(destfile, 
    skip = 2,
    sheet = "T3. Other workers stats") %>% 
  rename(`Motor vehicle and parts dealers` = `Motor vehicle and parts dealers \r\r\n`,
         `Metro name` = `...1`)

# Definitions of vulnerable industries (not they use ACS codes, not full NAICS.)
# You will have to manually crosswalk
Brookings.VulnerableIndustryDefinitions <- read_excel(destfile, 
    skip = 2,
    sheet = "T4. Vulnerable industries") %>% 
  rename(PUMS.NAICS = `NAICS code`) %>% 
  select(-Rationale) # Multiple sheets, so just look at the spreadsheet if you want this. 

```

Show some characteristics of immediately affected sectors. First make a crosswalk to EMSI. There are two ACS codes that could be problemmatic 4MS and 4462.

```{r}
brookingsRiskSectors_bars <- Industry_Table_allLevels_NOMSA_Employees %>% 
  left_join(Brookings.VulnerableIndustryDefinitions, by=c("NAICS"="PUMS.NAICS")) %>% 
  mutate(brookingsVulnerable = ifelse(!is.na(`Industry title`), "Immediate risk", "Other")) %>% 
  
  # recode ACS NAICS codes
  mutate(brookingsVulnerable = ifelse(brookingsVulnerable == "Other" & 
                                        NAICS %in% c("44612",
                                                   "44613",
                                                   "44619",
                                                   "4542",
                                                   "4543"), 
                                      "Immediate risk", brookingsVulnerable)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk",
                                             Description, NA)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk" & 
                                               str_sub(NAICS, 1,2) %in% c("44","45"), 
                                             "Retail", brookingsVulnerable.sector)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk" & 
                                               str_sub(NAICS, 1,2) == "48", 
                                             "Transportation", brookingsVulnerable.sector)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk" & 
                                               str_sub(NAICS, 1,2) == "71", 
                                             "Arts, Entertainment & Recreation", brookingsVulnerable.sector)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk" & 
                                               NAICS %in% c("5615","6244","812","8131"), 
                                             "Other Services", brookingsVulnerable.sector)) %>% 

  # summarise and make chart
  filter(brookingsVulnerable == "Immediate risk") %>% 
  group_by(brookingsVulnerable.sector) %>% 
  summarise(emp = sum(`2019 Jobs`),
            emp.p =  sum(`2019 Jobs`)/TotalEmployment.NOMSA) %>% 
  ggplot(aes(brookingsVulnerable.sector, emp)) +
  geom_col(fill = DCcolor.p2blue) +
  geom_text(aes(x = brookingsVulnerable.sector, y = .5*emp, label = scales::percent(emp.p, accuracy = .1)), color = "white", family = "Asap") + 
    geom_text(aes(x = brookingsVulnerable.sector, y = emp, label = scales::comma(emp)), family = "Asap", hjust = 0, nudge_y = 1000) + 

  scale_y_continuous(labels = scales::comma, limits = c(0,80000)) +
  coord_flip() +
  themeDC_vertical() + 
  labs(x="",
       y="Employment",
       title = "Employment and employment share at “immediate risk” \nby industry sector, New Orleans metro, 2019",
       subtitle = "Based on industries identified by Brookings")

ggsave(brookingsRiskSectors_bars, filename = "outputs/graphicsVulnerableWorkers/brookingsRiskSectors_bars.png", type = "cairo",
       width = 8, height = 3, units = "in")



# Export a CSV.
Industry_Table_allLevels_NOMSA_Employees %>% 
  left_join(Brookings.VulnerableIndustryDefinitions, by=c("NAICS"="PUMS.NAICS")) %>% 
  mutate(brookingsVulnerable = ifelse(!is.na(`Industry title`), "Immediate risk", "Other")) %>% 
  
  # recode ACS NAICS codes
  mutate(brookingsVulnerable = ifelse(brookingsVulnerable == "Other" & 
                                        NAICS %in% c("44612",
                                                   "44613",
                                                   "44619",
                                                   "4542",
                                                   "4543"), 
                                      "Immediate risk", brookingsVulnerable)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk",
                                             Description, NA)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk" & 
                                               str_sub(NAICS, 1,2) %in% c("44","45"), 
                                             "Retail", brookingsVulnerable.sector)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk" & 
                                               str_sub(NAICS, 1,2) == "48", 
                                             "Transportation", brookingsVulnerable.sector)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk" & 
                                               str_sub(NAICS, 1,2) == "71", 
                                             "Arts, Entertainment & Recreation", brookingsVulnerable.sector)) %>% 
  mutate(brookingsVulnerable.sector = ifelse(brookingsVulnerable == "Immediate risk" & 
                                               NAICS %in% c("5615","6244","812","8131"), 
                                             "Other Services", brookingsVulnerable.sector)) %>% 
  select(Level, NAICS, Description, `2019 Jobs`, `Avg. Earnings Per Job`) %>%
  write.csv("outputs/graphicsVulnerableWorkers/Brookings_IndustriesAtRisk_detailedNOmetro.csv")


```

Sources: 

- Immediate risk sectors: Berube, Alan and Nicole Bateman. (2020). ["Who are the workers already impacted by the COVID-19 recession?"](https://www.brookings.edu/research/who-are-the-workers-already-impacted-by-the-covid-19-recession/) Brookings.
- Employment: The Data Center analysis of data from EMSI. Brookings originally used ACS microdata to estimate employment, and these EMSI-derived estimates differ slightly.
- Note: Employment by sector is not necessarily the total employment in the sector but only the portion of each sector that is in immediate-risk industries. For example, the Retail sector excludes the Pharmacy and Grocery Store industries. 



Make the ranked bar graph with the updated data. They estimate a total of 28 percent for New Orleans metro.

```{r}
# 
brookings.rankedRiskMetros_bars <- Brookings.AllWorkersStats %>% 
  left_join(MSApopulations, by="Metro name") %>% 
  filter(estimate > 1000000) %>% 
  mutate(label = ifelse(`Metro name` == "New Orleans-Metairie, LA",
                        "New Orleans-Metairie, LA", "")) %>% 
  ggplot(aes(reorder(`Metro name`, -`Share of workers that are vulnerable`), `Share of workers that are vulnerable`)) + 
  geom_col(aes(fill = label)) + 
  scale_fill_manual(values = c(DCcolor.p1skyblue, DCcolor.p1darkblue)) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + 
  annotate("text", x = "New Orleans-Metairie, LA", y = .32, label = "New Orleans metro", hjust = 0, family = "Asap", color="grey20") + 
  themeDC_scatter() + 
  theme(legend.position = "none",
        axis.text.x = element_blank()) + 
  labs(x="Metros ranked by share of employment in “immediate-risk” industries", y="",
       #y="Share of metro jobs in “most at risk” industry groups",
       title = "Share of employment in “immediate-risk” industries by metro (over 1 million population)",
       subtitle = "Brookings analysis of ACS microdata")

ggsave(brookings.rankedRiskMetros_bars, filename = "outputs/graphicsVulnerableWorkers/brookings.rankedRiskMetros_bars.png", type = "cairo",
       width = 8, height = 3, units = "in")


# What are the top metros?
Brookings.AllWorkersStats %>% 
  left_join(MSApopulations, by="Metro name") %>% 
  filter(estimate > 1000000) %>% 
  mutate(label = ifelse(`Metro name` == "New Orleans-Metairie, LA",
                        "New Orleans-Metairie, LA", "")) %>% 
  select(`Metro name`, `Share of workers that are vulnerable`) %>% 
  arrange(desc(`Share of workers that are vulnerable`))

```

Sources: 
- Employment in "immediate-risk" industries: [Brookings analysis](https://www.datacenterresearch.org/covid-19-data-and-information/covid-19-economic-analysis/) of ACS microdata from 2018.
- Metro populations: American Community Survey 1-year estimates for 2018.

Now do demographic breakdowns for vulnerable vs other: Age 18-24, Work part-time/part-year, Female, and No more than H.S. diploma, Black or Latino. Also do median annual wages and median family income. 

```{r}
# Income 
brookingsVulnerableIncome_bars <- Brookings.VulnerableWorkersStats %>% 
  mutate(brookingsVulnerable = "Immediate risk") %>% 
  filter(`Metro name` == "New Orleans-Metairie, LA") %>% 
  bind_rows(filter(Brookings.OtherWorkersStats, `Metro name` == "New Orleans-Metairie, LA")) %>% 
  mutate(brookingsVulnerable = ifelse(is.na(brookingsVulnerable), "Other", brookingsVulnerable),
         `black or latino` = `Latino or Hispanic` + `Black, non-Hispanic`,
         `No more than H.S. diploma` = `Less than high school diploma` + `High school diploma or equivalent`,
         Female = 1-Male) %>% 
  gather(-`Metro name`, -brookingsVulnerable, key = variable, value = value) %>% 
  filter(variable %in% c("Median annual wages",
                         "Median family income")) %>% 
  mutate(variable = fct_relevel(variable,
                                c("Median annual wages",
                                  "Median family income"))) %>% 
  ggplot(aes(variable, value, group=brookingsVulnerable)) + 
  geom_col(aes(fill = brookingsVulnerable), position = "dodge", width=0.7) + 
  geom_text(aes(x = variable, y = value, label = scales::dollar(value, accuracy=1)), family = "Asap", hjust = -.1, position=position_dodge(width = .7)) + 
  scale_fill_manual(values =  c("Other" = DCcolor.p1skyblue, "Immediate risk" = DCcolor.p2yellow)) + 
  scale_y_continuous(labels = scales::dollar_format(accuracy = 1), limits=c(0,100000)) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 13)) + 
  coord_flip() +
  themeDC_vertical() + 
  theme(legend.position = "top",
        legend.title = element_blank(),
        legend.justification = "left",
        legend.direction = "vertical") + 
  labs(x="",
       y="Dollars",
       title = "Income of workers in immediate-risk industries vs \nworkers in other industries, New Orleans metro",
       subtitle = "Brookings analysis of ACS microdata, 2018")
ggsave(brookingsVulnerableIncome_bars, filename = "outputs/graphicsVulnerableWorkers/brookingsVulnerableIncome_bars.png", type = "cairo",
       width = 8, height = 3.3, units = "in")



# Demog and economic characteristics 
brookingsVulnerableDemog_bars <- Brookings.VulnerableWorkersStats %>% 
  mutate(brookingsVulnerable = "Immediate risk") %>% 
  filter(`Metro name` == "New Orleans-Metairie, LA") %>% 
  bind_rows(filter(Brookings.OtherWorkersStats, `Metro name` == "New Orleans-Metairie, LA")) %>% 
  mutate(brookingsVulnerable = ifelse(is.na(brookingsVulnerable), "Other", brookingsVulnerable),
         `black or latino` = `Latino or Hispanic` + `Black, non-Hispanic`,
         `No more than H.S. diploma` = `Less than high school diploma` + `High school diploma or equivalent`,
         Female = 1-Male) %>% 
  gather(-`Metro name`, -brookingsVulnerable, key = variable, value = value) %>% 
  filter(variable %in% c("Ages 18-24",
                         "Bachelor's or more",
                         "Female",
                         "Rents home",
                         "Worked full time/full-year",
                         "Below 200% FPL",
                         "Black, non-Hispanic")) %>% 
  mutate(variable = fct_relevel(variable,
                                c("Worked full time/full-year",
                                  "Bachelor's or more",
                                  "Below 200% FPL",
                                  "Rents home",
                                  "Ages 18-24",
                                  "Female",
                                  "Black, non-Hispanic"))) %>% 
  ggplot(aes(variable, value, group=brookingsVulnerable)) + 
  geom_col(aes(fill = brookingsVulnerable), position = "dodge", width=0.7) + 
  geom_text(aes(x = variable, y = value, label = scales::percent(value, accuracy = 1)), family = "Asap", hjust = -.1, position=position_dodge(width = 0.7)) + 
  scale_fill_manual(values =  c("Other" = DCcolor.p1skyblue, "Immediate risk" = DCcolor.p2yellow)) + 
  scale_y_continuous(labels = scales::percent, limits=c(0,.85)) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20)) + 
  coord_flip() +
  themeDC_vertical() + 
  theme(legend.position = "top",
        legend.title = element_blank(),
        legend.justification = "left",
        legend.direction = "vertical") + 
  labs(x="",
       y="Share of workers",
       title = "Selected demographic and socioeconomic characteristics of \nworkers in immediate-risk and other industries, New Orleans metro",
       subtitle = "Brookings analysis of ACS microdata, 2018")
ggsave(brookingsVulnerableDemog_bars, filename = "outputs/graphicsVulnerableWorkers/brookingsVulnerableDemog_bars.png", type = "cairo",
       width = 8, height = 4.3, units = "in")


```

Sources:
- Characteristics of "immediate-risk" and other industries: [Brookings analysis](https://www.datacenterresearch.org/covid-19-data-and-information/covid-19-economic-analysis/) of ACS microdata from 2018.
- Note: To develop their estimates, Brookings used a 1-year sample of ACS microdata of 1,661 immediate-risk workers and 4,553 other workers. 


### Job Quality Index 

JQI published a list of industries. In theory, these should be convertable to NAICS codes. 

- I think these jobs are from CES, and it only shows Production and Non-supervisory jobs.
- There are 180 industry groups. 

Which NAICS level has 180 groups?

```{r}
Industry_Table_allLevels_NOMSA_Employees %>% 
  count(Level)
```

I made a separate spreadsheet that included a crosswalk for the industries. You can import this and join it to the 

```{r}
# Pull in your re-creation of the JQI data
JQItoNAICS <- read_excel("inputs/JQI/JQItoNAICS.xlsx", 
    sheet = "JCQtoNAICS") %>% 
  mutate(industryNAICS = as.character(industryNAICS))
```

Then do the analysis with the JQI industries. Note that this is higher than the JQI numbers because they only include production and non-supervisory from CES.

```{r}
# By industry
Industry_Table_allLevels_NOMSA_Employees %>% 
  left_join(JQItoNAICS, by=c("NAICS"="industryNAICS")) %>% 
  mutate(JQIatRisk = ifelse(!is.na(industryJQI), "At risk", "")) %>% 
  filter(!is.na(industryJQI)) %>% 
  group_by(industryJQI) %>% 
  summarise(jobs = sum(`2019 Jobs`)) %>% 
  mutate(jobs.p = jobs/TotalEmployment.NOMSA)
  
# Total
Industry_Table_allLevels_NOMSA_Employees %>% 
  left_join(JQItoNAICS, by=c("NAICS"="industryNAICS")) %>% 
  mutate(JQIatRisk = ifelse(!is.na(industryJQI), "At risk", "")) %>% 
  filter(!is.na(industryJQI)) %>% 
  group_by(industryJQI) %>% 
  summarise(jobs = sum(`2019 Jobs`)) %>% 
  mutate(jobs.p = jobs/TotalEmployment.NOMSA) %>% 
  summarise(jobs = sum(jobs, na.rm = TRUE),
            jobs.p = jobs/TotalEmployment.NOMSA)


```

### Eyeballing occupations

First, try copying Boston's scatterplot. Try to show wages, race, and gender. You could do this with IPUMS, too. 

```{r}
# Define at risk occupations
OccupationsHighRisk.2dig <- c("Installation, Maintenance, and Repair Occupations",
                            "Transportation and Material Moving Occupations",
                            "Sales and Related Occupations",
                            "Personal Care and Service Occupations",
                            "Food Preparation and Serving Related Occupations",
                            "Arts, Design, Entertainment, Sports, and Media Occupations")

OccupationsMidRisk.2dig <- c("Construction and Extraction Occupations",
                            "Office and Administrative Support Occupations",
                            "Production Occupations",
                            "Building and Grounds Cleaning and Maintenance Occupations")


# Make a chart by occupation
occGroupWageAtrisk_scatter <- Occupation_Table_allLevels_NOMSA_Employees %>% 
  filter(Level == 2) %>% 
  mutate(risk = "moderate risk",
         risk = ifelse(Description %in% OccupationsMidRisk.2dig, "higher risk", risk),
         risk = ifelse(Description %in% OccupationsHighRisk.2dig, "immediate risk", risk),
         risk = fct_relevel(as.factor(risk), c("moderate risk","higher risk","immediate risk"))) %>% 
  mutate(Description.short = str_replace(Description," Occupations","")) %>% 
  ggplot(aes(`2019 Jobs`, `Median Hourly Earnings`)) + 
  geom_point(aes(color = risk), size = 3) + 
  geom_text_repel(aes(label = Description.short), size=3, family = "Asap", color="grey50") + 
  scale_color_manual(values = c(DCcolor.p2limegreen, DCcolor.p1mediumblue, DCcolor.p2orangered)) + 
  scale_y_continuous(labels = scales::dollar) +
  scale_x_continuous(label = scales::comma) + 
  themeDC_scatter() + 
  theme(legend.position = "top",
        legend.title = element_blank(), 
        legend.justification = "left") + 
  labs(title = "Major occupation groups and risk of job loss",
       subtitle = "New Orleans metro")
  
```

Demographics of risky occupations? Looks like Workers of color are probably disproportionate in high-risk industries. 

```{r}
Occupation_Table_allLevels_NOMSA_Employees %>% 
  filter(Level == 2) %>% 
  mutate(risk = "moderate risk",
         risk = ifelse(Description %in% OccupationsMidRisk.2dig, "higher risk", risk),
         risk = ifelse(Description %in% OccupationsHighRisk.2dig, "immediate risk", risk),
         risk = fct_relevel(as.factor(risk), c("moderate risk","higher risk","immediate risk"))) %>% 
  mutate(Description.short = str_replace(Description," Occupations","")) %>% 
  mutate(nonwhite = `2019 Jobs` - White,
         nonwhite.p = 1- `White % of Occupation`) %>% 
  group_by(risk) %>% 
  summarise(female = sum(Females, na.rm = TRUE),
            male = sum(Males, na.rm = TRUE),
            nonwhite = sum(nonwhite, na.rm = TRUE),
            white = sum(White, na.rm = TRUE)) %>% 
  gather(-risk, key = key, value = jobs) %>% 
  ungroup() %>%
  mutate(var.type = ifelse(key %in% c("female","male"),"sex","race")) %>% 
  arrange(risk) %>% 
  group_by(risk, var.type) %>% 
  mutate(jobs.p = jobs/sum(jobs))



```

### ATUS analysis

Prepare the ATUS data and ACS data. 

```{r}

tidycensus::census_api_key("530ce361defc2c476e5b5d5626d224d8354b9b9a")

# Select variables
AllACSVariables <- load_variables(2018, "acs5", cache = TRUE)
#View(AllACSVariables)

EmpByIndustryACSvars <- AllACSVariables %>% 
  filter(str_sub(name, 1,6) == "C24030")


# Download ACS employment by industry for NO and NO MSA
ACS.OrleansParish.Employment <- get_acs(geography = "county",
        variables = EmpByIndustryACSvars$name,
        state = "LA",
        county = "071",
        year = 2018)

ACS.AllMSA.Employment <- get_acs(geography = "metropolitan statistical area/micropolitan statistical area",
        variables = EmpByIndustryACSvars$name,
        year = 2018)

ACS.USA.Employment <- get_acs(geography = "us",
        variables = EmpByIndustryACSvars$name,
        year = 2018)

# Industry sectors -- to clean up the ACS labels and make consistent with ATUS
ACS.ATUScomparable.IndustrySectorLabels <- tibble(
  label = c(
  "Estimate!!Total",
  "Estimate!!Total!!Male!!Agriculture, forestry, fishing and hunting, and mining!!Agriculture, forestry, fishing and hunting",
  "Estimate!!Total!!Male!!Agriculture, forestry, fishing and hunting, and mining!!Mining, quarrying, and oil and gas extraction",
  "Estimate!!Total!!Male!!Construction",
  "Estimate!!Total!!Male!!Manufacturing",
  "Estimate!!Total!!Male!!Wholesale trade",
  "Estimate!!Total!!Male!!Retail trade",
  "Estimate!!Total!!Male!!Transportation and warehousing, and utilities",
  "Estimate!!Total!!Male!!Information",
  "Estimate!!Total!!Male!!Finance and insurance, and real estate, and rental and leasing",
  "Estimate!!Total!!Male!!Professional, scientific, and management, and administrative, and waste management services",
  "Estimate!!Total!!Male!!Educational services, and health care and social assistance",
  "Estimate!!Total!!Male!!Arts, entertainment, and recreation, and accommodation and food services",
  "Estimate!!Total!!Male!!Other services, except public administration",
  "Estimate!!Total!!Male!!Public administration",
  "Estimate!!Total!!Female!!Agriculture, forestry, fishing and hunting, and mining!!Agriculture, forestry, fishing and hunting",
  "Estimate!!Total!!Female!!Agriculture, forestry, fishing and hunting, and mining!!Mining, quarrying, and oil and gas extraction",
  "Estimate!!Total!!Female!!Construction",
  "Estimate!!Total!!Female!!Manufacturing",
  "Estimate!!Total!!Female!!Wholesale trade",
  "Estimate!!Total!!Female!!Retail trade",
  "Estimate!!Total!!Female!!Transportation and warehousing, and utilities",
  "Estimate!!Total!!Female!!Information",
  "Estimate!!Total!!Female!!Finance and insurance, and real estate, and rental and leasing",
  "Estimate!!Total!!Female!!Professional, scientific, and management, and administrative, and waste management services",
  "Estimate!!Total!!Female!!Educational services, and health care and social assistance",
  "Estimate!!Total!!Female!!Arts, entertainment, and recreation, and accommodation and food services",
  "Estimate!!Total!!Female!!Other services, except public administration",
  "Estimate!!Total!!Female!!Public administration"),
  sector = c(
    "Total",
    "Agriculture, forestry, fishing, and hunting",
    "Mining, quarrying, and oil and gas extraction",
    "Construction",
    "Manufacturing",
    "Wholesale and retail trade",
    "Wholesale and retail trade",
    "Transportation and utilities",
    "Information",
    "Financial activities",
    "Professional and business services",
    "Education and health services",
    "Leisure and hospitality",
    "Other services",
    "Public administration",
    "Agriculture, forestry, fishing, and hunting",
    "Mining, quarrying, and oil and gas extraction",
    "Construction",
    "Manufacturing",
    "Wholesale and retail trade",
    "Wholesale and retail trade",
    "Transportation and utilities",
    "Information",
    "Financial activities",
    "Professional and business services",
    "Education and health services",
    "Leisure and hospitality",
    "Other services",
    "Public administration"
  )
)

# Import ATUS tables

ATUStables_2017_2018.paidLeave <- read_excel("inputs/excelwork/ATUStables_2017-2018.xlsx", 
    sheet = "cleanPaidLeave", col_types = c("text", 
        "text", "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric"))

ATUStables_2017_2018.workedHome <- read_excel("inputs/excelwork/ATUStables_2017-2018.xlsx", 
    sheet = "cleanWorkedFromHome", col_types = c("text", 
        "text", "numeric", "numeric", "numeric", 
        "numeric", "numeric", "text", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric"))



# Combine sexes and industries and join with ATUS percentages.
ACS.ATUS.PaidLeave.OrleansParish.Employment <- 
  ACS.OrleansParish.Employment %>% 
  left_join(ACS.USA.Employment, by="variable") %>% 
  rename(estimate.NO = estimate.x,
         estimate.USA = estimate.y,
         moe.NO = moe.x,
         moe.USA = moe.y) %>% 
  select(variable, estimate.NO, moe.NO, estimate.USA, moe.USA) %>% 
  left_join(EmpByIndustryACSvars, by=c("variable" = "name")) %>% 
  filter(label %in% ACS.ATUScomparable.IndustrySectorLabels$label | 
           label == "") %>% 
  left_join(ACS.ATUScomparable.IndustrySectorLabels, by="label") %>% 
  group_by(sector) %>% 
  summarise(local.acs.emp = sum(estimate.NO, na.rm = TRUE),
            us.acs.emp = sum(estimate.USA, na.rm = TRUE)) %>% 
  left_join(filter(ATUStables_2017_2018.paidLeave, VariableGroup == "Industry"), 
            by=c("sector"="Variable")) %>% 
  mutate(local.emp.pct = local.acs.emp/sum(local.acs.emp),
         us.emp.pct = us.acs.emp/sum(us.acs.emp),
         local.paidLeave.est = .01*PaidLeaveYes*local.acs.emp,
         local.nopaidleave.est = .01*PaidLeaveNo*local.acs.emp)

ACS.ATUS.workedHome.OrleansParish.Employment <- 
  ACS.OrleansParish.Employment %>% 
  left_join(ACS.USA.Employment, by="variable") %>% 
  rename(estimate.NO = estimate.x,
         estimate.USA = estimate.y,
         moe.NO = moe.x,
         moe.USA = moe.y) %>% 
  select(variable, estimate.NO, moe.NO, estimate.USA, moe.USA) %>% 
  left_join(EmpByIndustryACSvars, by=c("variable" = "name")) %>% 
  filter(label %in% ACS.ATUScomparable.IndustrySectorLabels$label | 
           label == "") %>% 
  left_join(ACS.ATUScomparable.IndustrySectorLabels, by="label") %>% 
  group_by(sector) %>% 
  summarise(local.acs.emp = sum(estimate.NO, na.rm = TRUE),
            us.acs.emp = sum(estimate.USA, na.rm = TRUE)) %>% 
  left_join(filter(ATUStables_2017_2018.workedHome, VariableGroup == "Industry"), 
            by=c("sector"="Variable")) %>% 
  mutate(local.emp.pct = local.acs.emp/sum(local.acs.emp),
         us.emp.pct = us.acs.emp/sum(us.acs.emp),
         local.workedHome.est = .01*PctWorkedFromHome*local.acs.emp,
         local.noworkedHome.est = (1-.01*PctWorkedFromHome)*local.acs.emp)

  
  
```

Now make the charts.

These three tell a story together, but think about how best to show all of this. 

```{r}
# Paid leave percentage
paidLeaveByIndustry <- ACS.ATUS.PaidLeave.OrleansParish.Employment %>% 
  filter(sector != "Total" & sector != "Mining, quarrying, and oil and gas extraction") %>% 
  mutate(sector = fct_reorder(sector, local.emp.pct)) %>% 
  ggplot(aes(sector, .01*PaidLeaveYes)) + 
  geom_col(fill = DCcolor.p1grayblue) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  coord_flip() + 
  themeDC_vertical() + 
  theme(legend.position = "top",
        legend.title = element_blank()) +
  labs(x="",
       y="Percent of workers",
       title = "Share who reported access to paid leave by industry group, US",
       subtitle = "2017-2018 American Time Use Survey")

# Worked from home percentage
workedFromHomebyIndustry <- ACS.ATUS.workedHome.OrleansParish.Employment %>% 
  filter(sector != "Total" & sector != "Mining, quarrying, and oil and gas extraction") %>% 
  mutate(sector = fct_reorder(sector, local.emp.pct)) %>% 
  ggplot(aes(sector, .01*PctWorkedFromHome)) + 
  geom_col(fill = DCcolor.p2green) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  coord_flip() + 
  themeDC_vertical() + 
  theme(legend.position = "top",
        legend.title = element_blank()) +
  labs(x="",
       y="Percent of workers",
       title = "Share who worked from home by industry group, US",
       subtitle = "2017-2018 American Time Use Survey")

# Local vs national share of employment -- this is cleaner.
NOMSAandUSAemploymentShareByIndustry<- ACS.ATUS.PaidLeave.OrleansParish.Employment %>% 
  filter(sector != "Total" & sector != "Mining, quarrying, and oil and gas extraction") %>% 
  select(sector, local.emp.pct, us.emp.pct) %>% 
  gather(-sector, key=locnat, value = emp.p) %>% 
  mutate(locnat = ifelse(locnat == "local.emp.pct", "New Orleans metro", "National"),
         locnat = fct_relevel(locnat, c("New Orleans metro", "National"))) %>% 
  ggplot(aes(sector, emp.p)) + 
  geom_col(aes(fill = locnat), position = "dodge") + 
  scale_fill_manual(values =  c("New Orleans metro" = DCcolor.p1skyblue, "National" = DCcolor.p2yellow)) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0,.15)) +
  coord_flip() + 
  themeDC_vertical() + 
  theme(legend.position = "top",
        legend.title = element_blank()) +
  labs(x="",
       y="Percent of workers",
       title = "Share of total employment by industry group, US and New Orleans metro",
       subtitle = "2014-2018 American Community Survey")

  
```

### More detailed St Louis Fed occupations

This comes from this [blog post](https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-unemployment-risk-state-msa-differences). It has full classifications of all of the SOC codes. Let's apply it by doing the following:

- Show where New Orleans ranks relative to other metros.
- Join to EMSI or OES and do some breakdowns by wage level, maybe by demographics. 
- Crosswalk to IPUMS. 

```{r}
# Load MSA data (state should be the same)
url <- "https://www.stlouisfed.org/~/media/Files/Excel/ote_employment_risk_by_location.xlsx?la=en"
destfile <- "ote_employment_risk_by_location.xlsx"
curl::curl_download(url, destfile)
ote_employment_risk_by_location.metro <- read_excel(destfile, 
                                                    sheet = "metropolitan areas",
                                                    skip = 1) %>% 
  rename(MSA = ...1,
         HighRiskRank = ...2,
         TotalEmp = ...3,
         LowRiskTotal = Total...4,
         HighRiskTotal = Total...8)

# Load occupational coding
ote_employment_risk_by_location.SOCcodes <- read_excel(destfile, 
                                                       sheet = "classification by occupation") 

```

This could be graphed in essentially the same manner as the Brookings piece. Our MSA is still near the top of large metros. 

```{r}
StLFed.Gascon.rankedRiskMetros_bars <- ote_employment_risk_by_location.metro %>% 
  left_join(MSApopulations, by=c("MSA" = "Metro name")) %>% 
  filter(estimate > 1000000) %>% 
  mutate(label = ifelse(MSA == "New Orleans-Metairie, LA",
                        "New Orleans-Metairie, LA", "")) %>% 
  ggplot(aes(reorder(`MSA`, -HighRiskTotal), .01*HighRiskTotal)) + 
  geom_col(aes(fill = label)) + 
  scale_fill_manual(values = c(DCcolor.p1skyblue, DCcolor.p1darkblue)) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0,.6)) + 
  annotate("text", x = "New Orleans-Metairie, LA", y = .52, label = "New Orleans metro", hjust = 0, family = "Asap", color="grey20") + 
  themeDC_scatter() + 
  theme(legend.position = "none",
        axis.text.x = element_blank()) + 
  labs(x="Metros ranked by share of employment in “high-risk” occupations", y="",
       #y="Share of metro jobs in “most at risk” industry groups",
       title = "Share of employment in “high-risk” occupations by metro (over 1 million population)",
       subtitle = "St. Louis Fed Analysis")
ggsave(StLFed.Gascon.rankedRiskMetros_bars, filename = "outputs/graphicsVulnerableWorkers/StLFed.Gascon.rankedRiskMetros_bars.png", type = "cairo",
       width = 8, height = 3, units = "in")

```

Source:

- Gascon, Charles. (2019). ["COVID-19 and Unemployment Risk: State and MSA Differences"](https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-unemployment-risk-state-msa-differences) *On the Economy Blog*. The St. Louis Federal Reserve.


Now, let's join it to EMSI. The St Louis Fed data has 808 SOC codes. EMSI data only has 775 detailed occupations. This means that you can join, but you should verify that no important occupations are falling off.

Join fails:

- Some of the EMSI SOC codes are different. For example EMSI SOC for First-Line Supervisors of Gaming Workers is 39-1018. For the St Louis Fed data, it's 39-1010. This 8/0 discrepancy is common, though there are some others where a number of occupations seem to have been combined in the EMSI data. EMSI appears to be the weird one in this case. [This document](https://kb.economicmodeling.com/wp-content/uploads/2019/09/Emsi-SOC-Changes-2018.3-9.pdf) might allow you to write a crosswalk easily. 
- This issue can probably be ignored if all of the unmatched industries are low-risk or if the number of jobs is small. 


```{r, eval=FALSE}
# Join to EMSI
Occupation_Table_allLevels_NOMSA_Employees %>% 
  filter(Level == "5") %>% 
  left_join(ote_employment_risk_by_location.SOCcodes, by=c("SOC"="Occupation code")) %>% 
  View()

# Full join to see what's been left out. 
Occupation_Table_allLevels_NOMSA_Employees %>% 
  filter(Level == "5") %>% 
  full_join(ote_employment_risk_by_location.SOCcodes, by=c("SOC"="Occupation code")) %>% 
  select(SOC, Description, `Occupation title`, `2019 Jobs`, `Total employment (US)`) %>% 
  View()

# Full by occ title. 
Occupation_Table_allLevels_NOMSA_Employees %>% 
  #filter(Level == "5") %>% 
  full_join(ote_employment_risk_by_location.SOCcodes, by=c("SOC"="Occupation code")) %>% 
  select(SOC, Description, `Occupation title`, `2019 Jobs`, `Total employment (US)`) %>% 
  View()


ote_employment_risk_by_location.SOCcodes %>% View()


```

Let's evaluate the way the SOCs are coded:

- 150 are essential. Looks pretty legit. Material moving, food processing, government, health care. 
- 262 are work from home and therefore low-risk. Looks mostly like white collar. 
- 441 are low-risk.
- 367 are high-risk. 

High risk includes some that might not be low-risk. A lot of industrial occupations.:

- Dredge operators.
- Laborers and Freight, Stock, and Material Movers, Hand
- Semiconductor Processors
- Medical Appliance Technicians
- Nuclear Power Reactor Operators
- Assemblers and Fabricators, All Other, Including Team Assemblers -- This is 1.3 M jobs.
- Maintenance and Repair Workers, General -- this is 1.4 M.
- Telecommunications Line Installers and Repairers
- Telecommunications Equipment Installers and Repairers, Except Line Installers
- Some medical, like Orthodontists and Podiatrists.
- Foresters
- Marine Engineers and Naval Architects
- Mechanical Engineers


Office and Administrative Support Occupations that are high-risk look right. Hotel and bank-type jobs. 

I think the post-secondary are the only problem. The rest is solved by leaving the Level open-ended. Essentially, the St Louis Fed data actually uses multiple levels. All but 6 were joined on Level 5, so check out those for issues: Can you combine the 4s and 5s with no double counting?

- 13-1020 Buyers and purchasing agents -- No issues.
- 15-2090 Miscellaneous Mathematical Science Occupations -- No issues.
- 29-2010 Clinical Laboratory Technologists and Technicians -- No issues
- 39-1010 First-Line Supervisors of Gaming Workers -- No issues
- 39-7010 Tour and Travel Guides -- No issues
- 47-4090 Miscellaneous Construction and Related Workers -- No issues

```{r}
StLFed.Gascon.OccGroup_scatter <- Occupation_Table_allLevels_NOMSA_Employees %>% 
  full_join(ote_employment_risk_by_location.SOCcodes, by=c("SOC"="Occupation code")) %>% 
  select(Level, SOC, Description, `Occupation title`, `2019 Jobs`, `Total employment (US)`, `Low risk`,`High risk`) %>% 
  filter(!is.na(Description)) %>% 
  mutate(`Low risk` = ifelse(SOC == "25-1099", 1,`Low risk`),
         `High risk` = ifelse(SOC == "25-1099", 0,`High risk`)) %>% 
  filter(!is.na(Description) & !is.na(`Occupation title`)) %>% 
  mutate(SOC.2dig = str_sub(SOC, 1,2)) %>% 
  group_by(SOC.2dig, `High risk`) %>% 
  summarise(jobs = sum(`2019 Jobs`)) %>% 
  spread(key = `High risk`, value=jobs) %>% 
  rename(lowRisk = `0`, highRisk= `1`) %>% 
  mutate(lowRisk = ifelse(is.na(lowRisk), 0, lowRisk),
         highRisk = ifelse(is.na(highRisk), 0, highRisk),
         riskShare = highRisk/(lowRisk + highRisk)) %>%
  right_join(filter(mutate(Occupation_Table_allLevels_NOMSA_Employees,
                           SOC.2dig = str_sub(SOC,1,2)), 
                    Level=="2"),
             by="SOC.2dig") %>% 
  mutate(Description.short = str_replace(Description," Occupations","")) %>% 
  filter(Description.short != "Unclassified Occupation",
         Description.short != "Military-only occupations") %>% 

  ggplot(aes(`2019 Jobs`, `Median Hourly Earnings`)) + 
  geom_point(aes(color = riskShare), size = 3) + 
  geom_text_repel(aes(label = Description.short), size=2.5, family = "Asap", color="grey50",
                  hjust = 1, vjust=1,
                  box.padding = .4, point.padding = .1,
                  segment.color="gray80", segment.size=.2) + 
  scale_color_distiller(palette = "Spectral", labels = scales::percent, guide = "legend") + 
  #scale_color_viridis_c(labels = scales::percent, guide = "legend") + 
  scale_y_continuous(labels = scales::dollar) +
  scale_x_continuous(label = scales::comma) + 
  themeDC_scatter() + 
  theme(legend.position = "top",
        legend.justification = "left") + 
  labs(title = "Employment in major occupation groups by share at risk and earnings, New Orleans metro",
       subtitle = "Based on detailed occupations at risk identified by St. Louis Federal Reserve",
       color = "Share in at-risk occupations")
ggsave(StLFed.Gascon.OccGroup_scatter, filename = "outputs/graphicsVulnerableWorkers/StLFed.Gascon.OccGroup_scatter.png", type = "cairo",
       width = 8, height = 6, units = "in")


# Export a CSV.
Occupation_Table_allLevels_NOMSA_Employees %>% 
  full_join(ote_employment_risk_by_location.SOCcodes, by=c("SOC"="Occupation code")) %>% 
  select(Level, SOC, Description, `Occupation title`, 
         `Low risk`,`High risk`, Essential, `Work from home`, `Salaried/Likely paid`,
         `2019 Jobs`,`Median Hourly Earnings`) %>% 
  filter(!is.na(Description)) %>% 
  mutate(`Low risk` = ifelse(SOC == "25-1099", 1,`Low risk`),
         `High risk` = ifelse(SOC == "25-1099", 0,`High risk`)) %>% 
  filter(!is.na(Description) & !is.na(`Occupation title`)) %>% 
  mutate(SOC.2dig = str_sub(SOC, 1,2)) %>% 
  select(Level, SOC, Description, 
         `Low risk`,`High risk`, Essential, `Work from home`, `Salaried/Likely paid`,
         `2019 Jobs`,`Median Hourly Earnings`) %>% 
  write.csv("outputs/graphicsVulnerableWorkers/StLFed_occupationsAtRisk_detailedNOmetro.csv")

```

Source:
- Gascon, Charles. (2019). ["COVID-19 and Unemployment Risk: State and MSA Differences"](https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-unemployment-risk-state-msa-differences) *On the Economy Blog*. The St. Louis Federal Reserve.




And let's throw in a few basic breakdowns. 

```{r}
# Low risk
StLFed.Gascon.lowRisk_bars <- ote_employment_risk_by_location.metro %>% 
  filter(MSA == "New Orleans-Metairie, LA") %>% 
  select(Essential,
         `WFH, not essential`,
         `Paid, not essential or WFH`,
         `Food preparation & serving-related occupations`,
         `Sales & related occupations`,
         `Production occupations`,
         `Installation, maintenance, and repair occupations`,
         `All other high-risk occupations`) %>% 
  gather() %>% 
  mutate(risk=ifelse(key %in% c("Essential",
                                "WFH, not essential",
                                "Paid, not essential or WFH"),
                     "Low-risk", "High-risk")) %>% 
  filter(risk == "Low-risk") %>% 
  mutate(key = ifelse(key == "WFH, not essential", "Work from home likely possible", key),
         key = ifelse(key == "Paid, not essential or WFH", "Likely to be salaried", key)) %>% 
  ggplot(aes(key, .01*value)) + 
  geom_col(fill = DCcolor.p1grayblue) +
  geom_text(aes(x = key, y = .01*value, label = scales::percent(.01*value, accuracy = 1)), family = "Asap", hjust = -.1) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20)) + 
  scale_y_continuous(labels = scales::percent, limits = c(0,.3)) +

  #facet_wrap(~risk, nrow = 2, scales = "free_y", strip.position = "top") +
  #facet_grid(risk~., scales = "free_y", space = "free", switch = "both") + 
  coord_flip() + 
  themeDC_vertical() + 
  labs(x="",
       y="Share of total employment",
       title = "Employment share in low-risk occupations \nin the New Orleans metro",
       subtitle = "St. Louis Federal Reserve analysis of Occupational \nEmployment Statistics, 2018")
ggsave(StLFed.Gascon.lowRisk_bars, filename = "outputs/graphicsVulnerableWorkers/StLFed.Gascon.lowRisk_bars.png", type = "cairo",
       width = 5, height = 2.8, units = "in")


# High risk
StLFed.Gascon.highRisk_bars <- ote_employment_risk_by_location.metro %>% 
  filter(MSA == "New Orleans-Metairie, LA") %>% 
  select(Essential,
         `WFH, not essential`,
         `Paid, not essential or WFH`,
         `Food preparation & serving-related occupations`,
         `Sales & related occupations`,
         `Production occupations`,
         `Installation, maintenance, and repair occupations`,
         `All other high-risk occupations`) %>% 
  gather() %>% 
  mutate(risk=ifelse(key %in% c("Essential",
                                "WFH, not essential",
                                "Paid, not essential or WFH"),
                     "Low-risk", "High-risk")) %>% 
  filter(risk == "High-risk") %>% 
  mutate(key = ifelse(key == "WFH, not essential", "Work from home possible", key),
         key = ifelse(key == "Paid, not essential or WFH", "Likely to be salaried", key)) %>% 
  ggplot(aes(key, .01*value)) + 
  geom_col(fill = DCcolor.p1grayblue) +
  geom_text(aes(x = key, y = .01*value, label = scales::percent(.01*value, accuracy = 1)), family = "Asap", hjust = -.1) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20)) + 
  scale_y_continuous(labels = scales::percent, limits = c(0,.3)) +

  #facet_wrap(~risk, nrow = 2, scales = "free_y", strip.position = "top") +
  #facet_grid(risk~., scales = "free_y", space = "free", switch = "both") + 
  coord_flip() + 
  themeDC_vertical() + 
  labs(x="",
       y="Share of total employment",
       title = "Employment share in high-risk occupations \nin the New Orleans metro",
       subtitle = "St. Louis Federal Reserve analysis of Occupational \nEmployment Statistics, 2018")

ggsave(StLFed.Gascon.highRisk_bars, filename = "outputs/graphicsVulnerableWorkers/StLFed.Gascon.highRisk_bars.png", type = "cairo",
       width = 5, height = 4, units = "in")



```

Source:
- Gascon, Charles. (2019). ["COVID-19 and Unemployment Risk: State and MSA Differences"](https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-unemployment-risk-state-msa-differences) *On the Economy Blog*. The St. Louis Federal Reserve.






The crosswalks for ACS appear [here](https://usa.ipums.org/usa/volii/acs_occtooccsoc10.shtml). It will probably be tricky because the OCCSOC variable combines some occupations into a code with XX in the label. How would you combine multiple SOC codes into 1110XX when 111021 also exists? See "Postsecondary" teachers for an example, I think these are all rolled into one in EMSI. 

### Pew analysis of industries

This is what [Pew did to define high-risk industries](https://www.pewresearch.org/fact-tank/2020/03/27/young-workers-likely-to-be-hard-hit-as-covid-19-strikes-a-blow-to-restaurants-and-other-service-sector-jobs/) in an analysis of CPS and CES.

> Food services and drinking places (9.7 million); accommodation (1.5 million); arts, entertainment, and recreation (3.4 million); child day care services (1.5 million); personal and laundry services (2.7 million); retail trade (16.2 million); and selected transportation industries (3.1 million). The selected transportation industries include air, rail and water transportation, bus service and urban transit, taxi and limousine services, scenic and sightseeing transportation, and services incidental to transportation.

I could replicate this pretty easily for comparison with the Brookings analysis and JQI. When you do it, maybe you should modify retail:

- Exclude 446110 Pharmacies and Drug Stores  from 446.
- Exclude 445110 Supermarkets and Other Grocery (except Convenience) Stores from 445.

This is pretty similar to the second Brookings piece.

## Writeup of content

### First draft: Which workers might be most at-risk (preliminary)?

Detailed employment data is released slowly,[^slowData] so we won't have a better picture of how COVID-19 is affecting our economy until updated data is released. In the meantime, a number of national and local analysts have offered back-of-the-envelope estimates of the number of jobs that might be most at risk, and their characteristics. The data presented below adapts a few of these approaches to our region. The Data Center will continue to refine this analysis, especially as more data on the post-COVID-19 becomes available. 

While these approaches vary, all are rough calculations based on old data. Together, they suggest a few takeaways with confidence. First, impacts will be widespread but concentrated in certain sectors of the economy. Second, the workers most at risk are less likely to have the ability to weather an income shock or employment disruption. Third, the New Orleans metro economy may have more risk due to the composition of its industries. 

[^slowData]: McCracken, Michael. (2020). "COVID-19: Forecasting with Slow and Fast Data." St. Louis Fed: *On the Economy Blog*. https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-forecasting-slow-fast-data; Bureau of Labor Statistics Release Calendar. https://www.bls.gov/schedule/

Brookings published data on the portion of each metro area's workforce in industry groups deemed "most at risk" by a Moody's Analytics research note. These industry groups are mining/oil and gas, transportation, employment services, travel arrangements, and leisure and hospitality. In the New Orleans metro, this adds up to about 128 thousand jobs, or 22 percent of metro employment. At 91 thousand, the vast majority of these at-risk jobs are in the leisure and hospitality "super sector." 

> maybe figure: Brookings industries

The leisure and hospitality super sector encompasses large industries most associated with tourism, such as hotels and full-service (sit down) restaurants. As reported in The Data Center's 2018 brief on the tourism workforce, these industries are over-represented in our region. Given the metro's concentration of tourism-related leisure and hospitality jobs and oil and gas activities, Brookings' analysis suggests that the New Orleans metro economy is particularly at risk. Among large metros with over 1 million population, only two had a higher concentration in at risk industry groups.

```{r}
brookingsAtRiskMetros_rankedbars
```


The Job Quality Index team at Cornell University took similar approach but with more detailed industry categories. According to their national estimates, these industries employ 37 million production and non-supervisory workers who are vulnerable to layoff in the short-term.[^JQIstatement] Applying the analysis to our regional economy suggests that the most "at risk" jobs are low-wage, low-hour jobs for front-line workers in industries that combined make up 32 percent of our regional economy. Unlike their national analysis, we are unable to break out local estimates front-line workers from supervisors other less-vulnerable positions due to data limitations.

> maybe figure: Job quality index industries


[^JQIstatement]: Statement from the U.S. Private Sector Job Quality Index (“JQI”) Team on
Vulnerabilities of Jobs in Certain Sectors to the Covid-19 Economic Shutdown. https://d3n8a8pro7vhmx.cloudfront.net/prosperousamerica/pages/5561/attachments/original/1584703152/JQI_Team_Statement_on_COVID-19_Economic_Shutdown_Job_Impact_031920.pdf?1584703152

> Figure JQI

National surveys shed light on other dimensions of vulnerability. Leisure and hospitality workers are less likely to report access to paid leave than workers in other sectors. If we assume that workers who have worked from home will have an easier transition to working from home during social isolation, 93 percent of leisure and hospitality workers nationally reported that they did not work from home, unlike in other sectors like Professional and Business Services where nearly half reported working from home. Our region's disproportionate share of leisure and hospitality workers suggests...
```{r}
paidLeaveByIndustry
workedFromHomebyIndustry
NOMSAandUSAemploymentShareByIndustry

```



Breaking down vulnerability by occupation, or the type of position, shows how vulnerability relates to earnings. As a rough estimate, we have classified occupation groups by the degree of risk for job loss. The occupation groups with the most immediate risk tend to include lower-earning jobs. Sales, Food Preparation and Serving, and Personal Care Occupations may be particularly vulnerable. 

```{r}
occGroupWageAtrisk_scatter
```

### Second draft: Which workers might be most at-risk of job loss (preliminary)?

Detailed employment data arrives slowly,[^slowData] so the picture of how COVID-19 is affecting our economy will grow clearer as updated data is released. Until then, prior data can provide clues into which parts of the economy might be most vulnerable. Several national and local analysts have offered back-of-the-envelope estimates of the number of workers at the most immediate risk of job loss, temporary furlough, or hour reduction.[^clarifyVulnerability] We've rounded up a few of these analyses, extended them, and highlighted data relevant to our region. At present, the analysis foregrounds immediate job loss risks rather than continued health risks to workers, impacts on small businesses and the public sector, or any of the possible longer-term economic ripple effects of COVID-19. The Data Center will refine and extend the analysis as we learn more and as more data on the economy after COVID-19 becomes available. 


[^slowData]: McCracken, Michael. (2020). "COVID-19: Forecasting with Slow and Fast Data." St. Louis Fed: *On the Economy Blog*. https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-forecasting-slow-fast-data; Bureau of Labor Statistics Release Calendar. https://www.bls.gov/schedule/

[^clarifyVulnerability]: These analyses focus on estimating immediate economic and job loss risks due to the broad economic shutdown of non-essential industries and stay at home orders, not the types of risks that frontline and essential employees experience as they keep going to work. Upcoming analysis will focus on this group and their vulnerabilities.   

Below, we present an industry approach and an occupation approach, both based on assumptions of how social distancing, travel restriction, and other public- and private-sector measures to slow the spread of coronavirus may disrupt some kinds of work more than others. To be clear, *risk* of job loss does not mean that every worker in these occupations and industries *has* or *will* experience job loss. Nor does it mean that other sectors will be entirely spared. These are rough but expert- and data-driven approximations based on comparing work that is essential or can be accomplished from home with work that is not -- restaurants, hotels, retail, and personal care. 

> BLUE BOX Industries versus occupations: Industries are types of business establishments that produce similar or related products or services. Hospitals, hotels, gas stations, and grocery stores are industries. Occupations are types of positions or professions defined by specific tasks or skills. Architects, economists, musicians, and cashiers are occupations. In general terms, data on industries more closely reflects how an economy is structured, and data on occupations more closely reflects the skills, working conditions, and quality of jobs.


Estimates differ, but a few key takeaways are consistent. First, impacts will be widespread but concentrated in certain sectors of the economy. Second, many of the workers most at risk may be less equipped to weather an income shock or employment disruption. Third, the New Orleans metro may face more risk than other metros due to the composition of its industries. 



#### Risk of job loss by industry

Brookings recently published metro-level data on industries at "immediate risk" of employment disruption due to COVID-19.[^brookingsSecondReport] Using 2018 data, they found that immediate-risk industries employed 37 million people, 23 percent of the US workforce. Other approaches using different methods have arrived at similarly large estimates.[^otherIndustryEstimates] 

[^brookingsSecondReport]: Berube, Alan and Nicole Bateman. (2020). ["Who are the workers already impacted by the COVID-19 recession?"](https://www.brookings.edu/research/who-are-the-workers-already-impacted-by-the-covid-19-recession/) Brookings.

[^otherIndustryEstimates]: The Job Quality Index team at Cornell University took similar approach with a different but overlapping set of industry categories. According to their national estimates, these industries employ 37 million production and non-supervisory workers who are vulnerable to layoff in the short-term. See ["Statement from the U.S. Private Sector Job Quality Index (“JQI”) Team on
Vulnerabilities of Jobs in Certain Sectors to the Covid-19 Economic Shutdown."](https://d3n8a8pro7vhmx.cloudfront.net/prosperousamerica/pages/5561/attachments/original/1584703152/JQI_Team_Statement_on_COVID-19_Economic_Shutdown_Job_Impact_031920.pdf?1584703152) Pew Research Center arrived at a similar number of 38.1 million. See ["Young workers likely to be hard hit as COVID-19 strikes a blow to restaurants and other service sector jobs."](https://www.pewresearch.org/fact-tank/2020/03/27/young-workers-likely-to-be-hard-hit-as-covid-19-strikes-a-blow-to-restaurants-and-other-service-sector-jobs/) 

For the New Orleans metro, 172,000 workers are in immediate-risk industries as defined by Brookings, making up 30 percent of the total workforce.[^brookingsDiscrepancy] The industry sectors with the largest risk of job loss are Retail and Food Services and Drinking Places, which account for about 19 percent of jobs in the region. Food Services and Drinking Places and Accomodation in particular are associated with tourism and [over-represented in our regional workforce, especially in the city of New Orleans](https://www.datacenterresearch.org/reports_analysis/benchmarking-new-orleans-tourism-economy-hotel-and-full-service-restaurant-jobs/). Aside from Retail, many of these jobs fall within the broader Leisure and Hospitality sector. While COVID-19 has curtailed demand, previous data also suggests that few workers in these industries are positioned to work from home. According to 2017-2018 national survey, only 9 percent of leisure and hospitality workers reported that they could work at home, the lowest of any industry sector.[^ATUSreference] In contrast, one-quarter of all workers and over one-half of in some predominately white-collar sectors reported that they could work from home. 

[^brookingsDiscrepancy]: These number differ slightly different from Brookings' estimate. Our analysis applies Brookings' industry definitions to a different source, EMSI data for 2019. Brookings's estimates use ACS microdata for 2018.

[^ATUSreference]: Bureau of Labor Statistics. (2019). ["Job Flexibilities and Work Schedules -- 2017-2018 Data from the American Time Use Survey"](https://www.bls.gov/news.release/pdf/flex2.pdf). Bureau of Labor Statistics News Release.


```{r}
brookingsRiskSectors_bars
```

Sources: 

- Immediate risk sectors: Berube, Alan and Nicole Bateman. (2020). ["Who are the workers already impacted by the COVID-19 recession?"](https://www.brookings.edu/research/who-are-the-workers-already-impacted-by-the-covid-19-recession/) Brookings.
- Employment: The Data Center analysis of data from EMSI. Brookings originally used ACS microdata to estimate employment, and these EMSI-derived estimates differ slightly. Note that employment by sector is not necessarily the total employment in the sector but only the portion of each sector that is in immediate-risk industries. For example, the Retail sector excludes the Pharmacy and Grocery Store industries. 



Given the sizeable role of tourism and oil and gas in our region relative to other regional economies, Brookings' analysis suggests that the New Orleans metro economy may face steeper job losses than other large metros. Among metros with population over 1 million, only three had a higher concentration in "immediate-risk" industries. 


```{r}
brookings.rankedRiskMetros_bars
```

Sources: 

- Employment in "immediate-risk" industries: [Brookings analysis](https://www.datacenterresearch.org/covid-19-data-and-information/covid-19-economic-analysis/) of ACS microdata from 2018. To develop their estimates, Brookings used a 1-year sample of ACS microdata of 1,661 immediate-risk workers and 4,553 other workers. 
- Metro populations: American Community Survey 1-year estimates for 2018.



The Brookings analysis also sheds light on the characteristics of vulnerable workers. Workers in "immediate-risk" industries earn lower wages and salaries and have lower family incomes than workers in other industries with less immediate risk. They are also more likely to be renters, more likely to live in households below 200 percent of the federal poverty level, less likely to have completed a four-year degree, and less likely to work on a full-time and year-round basis. The vulnerable industries are more than twice as likely to employ young workers (age 18-24). Black workers and women are also over-represented in vulnerable industries. 


```{r}
brookingsVulnerableIncome_bars
brookingsVulnerableDemog_bars

```

Source:
- Berube, Alan and Nicole Bateman. (2020). ["Who are the workers already impacted by the COVID-19 recession?"](https://www.brookings.edu/research/who-are-the-workers-already-impacted-by-the-covid-19-recession/) Brookings. To develop their estimates, Brookings used a 1-year sample of 2018 ACS microdata of 1,661 immediate-risk workers and 4,553 other workers. 


#### Risk of job loss by occupation

The St. Louis Federal Reserve published an occupation-based estimate of employment at high risk of layoff.[^StLFedOccupationsReference] Occupations essential to public health and safety, capable of completing work off-site or at home, or likely to be salaried were deemed "low risk" for job loss. With all other occupations deemed "high risk," this approach led to an even greater estimate of 67 million jobs at risk nationally or 46 percent of the national workforce. 

[^StLFedOccupationsReference]: Gascon, Charles. (2019). ["COVID-19: Which Workers Face the Highest Unemployment Risk?"](https://www.stlouisfed.org/on-the-economy/2020/march/covid-19-workers-highest-unemployment-risk) *On the Economy Blog*. The St. Louis Federal Reserve.; Gascon, Charles. (2019). ["COVID-19 and Unemployment Risk: State and MSA Differences"](https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-unemployment-risk-state-msa-differences) *On the Economy Blog*. The St. Louis Federal Reserve.

For greater New Orleans, occupations that were essential accounted for about 18 percent of employment, and work from home is likely possible for 27 percent of employment. High-risk occupations encompass 49 percent of the regional workforce.[^discrepancyWithStLouisFed] Two of the largest groups are occupations related to sales and the preparation and serving of food. 


```{r}
StLFed.Gascon.lowRisk_bars
StLFed.Gascon.highRisk_bars

```

Source: Gascon, Charles. (2019). ["COVID-19 and Unemployment Risk: State and MSA Differences"](https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-unemployment-risk-state-msa-differences) *On the Economy Blog*. The St. Louis Federal Reserve.


Echoing the finding of the industry-driven analysis published by Brookings, the share of employment in high-risk occupations in the New Orleans metro ranked near the top among large metros. 

[^discrepancyWithStLouisFed]: The Data Center replicated the St. Louis Fed analysis with EMSI data instead of the original's use of Occupational Employment Statistics. Both resulted in an estimate of 48.7 percent. 




```{r}
StLFed.Gascon.rankedRiskMetros_bars

```

Source: Gascon, Charles. (2019). ["COVID-19 and Unemployment Risk: State and MSA Differences"](https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-unemployment-risk-state-msa-differences) *On the Economy Blog*. The St. Louis Federal Reserve.


The occupation-based estimate also suggests that the most at-risk occupations tend to be lower earning. The scatterplot shows different hourly earning levels (median wages or hourly salary equivalents) for broad groups of occupations. Each occupation group is colored by the share of its employment that is classified as high-risk in the St. Louis Fed analysis. The occupation groups with the most extensive risk of job loss tend to include lower-earning jobs. Sales, Food Preparation and Serving, Personal Care and Building and Grounds Maintenance Occupations are low-wage occupations that may be particularly vulnerable to the combination of limited financial resources and the most immediate threat of job loss. 



```{r}
StLFed.Gascon.OccGroup_scatter

```

Source:

- At-risk industry definitions: Gascon, Charles. (2019). ["COVID-19 and Unemployment Risk: State and MSA Differences"](https://www.stlouisfed.org/on-the-economy/2020/april/covid-19-unemployment-risk-state-msa-differences) *On the Economy Blog*. The St. Louis Federal Reserve.
- Employment and wages: The Data Center analysis of EMSI occupation data. 



Though speculative, back-of-the-envelope estimates suggest that the risk of job loss is distributed unevenly to lower-earning workers. The shock of job loss can acutely undermine the ability of families with limited liquid assets to afford basic expenses and debt payments, especially when experiencing prolonged unemployment.[^chaseIncomeVolatility] Back-of-the-envelope estimates of jobs at risk expose how precarious finances and well-being could become for workers who may lose jobs and who likely already have lower earnings. The New Orleans metro's relatively large share of at-risk industries and occupations also points to the potential ripple effects in the regional economy of a prolonged contraction of spending by unemployed workers. During the public health fight against COVID-19, a major benchmark for economic policy responses will be to maintain the well-being of displaced workers. The employer payroll, unemployment insurance, and direct relief payment provisions of the CARES Act should alleviate some of the current and impending hardships, but further steps at the state and local steps might be necessary, both as a short-term "bridge" where support is late or insufficient and in the longer-term transition of vulnerable workers to more resilient employment opportunities.[^brookingsBridgeReference] 

[^brookingsBridgeReference]: Berube, Alan and Nicole Bateman. (2020). ["Who are the workers already impacted by the COVID-19 recession?"](https://www.brookings.edu/research/who-are-the-workers-already-impacted-by-the-covid-19-recession/) Brookings.

[^chaseIncomeVolatility]: Farrell, Diana, Fiona Grieg, and Chris Wheat. (2019). ["The potential economic impacts of COVID-19 on families, small businesses, and communities: Insights from five years of big data research."](https://institute.jpmorganchase.com/institute/research/household-income-spending/potential-economic-impacts-of-covid-19-on-families-small-businesses-communities). JPMorgan Chase & Co. Institue. 
