WFH_geographicComparison.qmd

---
title: "Working From Home in Illinois: Who Can and Who Does?"
format:
  html:
    df-print: paged
    toc: true
---


```{r setup, warning=FALSE, message=FALSE, include=FALSE}

library(scales)
library(reldist)
library(pollster)
library(labelled)
library(weights)
library(tigris)
library(ipumsr)
library(srvyr)
library(survey)
library(tidyverse)
library(naniar)
library(gmodels)
library(gtsummary)
library(quarto)

library(huxtable) # for summ() and regression output formatting
# Create the DB connection with the default name expected by PTAXSIM functions
library(jtools)
library(modelsummary)
library(car)

knitr::opts_chunk$set(warning=FALSE, message=FALSE, tidy = TRUE)

```

```{r include=FALSE}
Alea_theme <- function(){ 
    theme_classic() %+replace%    #replace elements we want to change
    
    theme(
      
      #grid elements
      panel.grid.major = element_blank(),    #strip major gridlines
      panel.grid.minor = element_blank(),    #strip minor gridlines
      axis.ticks = element_blank(),          #strip axis ticks
      axis.text.x = element_blank(),
     axis.text.y = element_blank(),

      #since theme_minimal() already strips axis lines, 
      #we don't need to do that again
      
      #text elements
      plot.title = element_text(             #title
                   size = 14,                #set font size
                   face = 'bold',            #bold typeface
                   hjust = 0,                #left align
                   vjust = 2),               #raise slightly
      
      plot.subtitle = element_text(          #subtitle
                   size = 14),               #font size
      
      plot.caption = element_text(           #caption
                   size = 9,                 #font size
                   hjust = 1),               #right align
      
      axis.title = element_text(             #axis titles
                   size = 10),               #font size
      
      axis.text = element_text(              #axis text
                   size = 9)                #font size
    
      
      #since the legend often requires manual tweaking 
      #based on plot content, don't define it here
    )
}
```

## ACS Survey Data Notes

IPUMS [link for Survey package](https://usa.ipums.org/usa/repwt.shtml)

Survey questions for EMPSTAT & LABFORCE:

1.  Last week, did this person work for pay at a job or business? (Yes or no) -- Yes becomes coded as EMPSTAT = 1-Employed.
2.  Last week, did this person do ANY work for pay, even as little as one hour?(Yes or no) -- Yes becomes coded as LABFORCE = 2-Yes in the labor force.

Survey questions for INCEARN:

1.  INCEARN = INCWAGE + INCBUS00
    -   Total amount earned in last 12 months: Wages, salary, commissions, bonuses, tips. \[Yes --\> \_\_\_\_\_\_ \] is coded as INCWAGE value.

    -   INCEARN includes self-employment income, INCWAGE does not.

        -   INCWAGE does not include Farming income and self-employment income, but INCEARN does.

`usa_00011.xml` and `usa_00011.dat.gz` are the same as Box files named `IL_2021_1yr_ACS.dat.gz` and `IL_2021_1yr_ACS_datDDI.xml`

> original xml file references the file name that it is called in the download. Either change the XML file to reference the correct .dat.gz files OR just keep track of which extracts are the same as the box file names.

```{r message=FALSE, warning = FALSE}
# old version with less variables:
#ddi <- read_ipums_ddi("usa_00009.xml") # 45 variables
#data <- read_ipums_micro(ddi) # 126623 observations before any filtering

# larger version with 147 variables. uses same file as Box file named "IL_2021_1yr_ACS.dat.gz and IL_2021_1yr_ACS_datDDI.xml
ddi <- read_ipums_ddi("./data/usa_00011.xml") # downloaded April 10 2023

data2021 <- read_ipums_micro(ddi) # 126623 observations before any filtering

data2021 <- data2021 %>% select(YEAR, INCEARN, INCWAGE, INCTOT, TRANWORK, OCCSOC, CLASSWKR, EMPSTAT, LABFORCE, PERWT, COUNTYFIP, PUMA, PWSTATE2, AGE, STRATA, CLUSTER, RACE, HISPAN, SEX, CIHISPEED, CINETHH, MULTGEN, NCHILD, NCHLT5, MARST, FERTYR, EDUC, DEGFIELD, OCC, IND, OCC2010, METRO, CITY, HHINCOME, SERIAL,HHWT, NUMPREC, SUBSAMP,HHTYPE )

# same sample but with 150+ variables. 
# NEED TO CHANGE XML file that references the data file. currently says usa_00011.dat.gz so these two lines of code do not work. 

ddi <- read_ipums_ddi("./data/IL_2019_1yearACS_datDDI.xml") # downloaded April 10 2023
data2019 <- read_ipums_micro(ddi) # 126623 observations before any filtering
data2019 <- data2019 %>% select(YEAR, INCEARN, INCWAGE, INCTOT, TRANWORK, OCCSOC, CLASSWKR, EMPSTAT, LABFORCE, PERWT, COUNTYFIP, PUMA, PWSTATE2, AGE, STRATA, CLUSTER, RACE, HISPAN, SEX, CIHISPEED, CINETHH, MULTGEN, NCHILD, NCHLT5, MARST, FERTYR, EDUC, DEGFIELD, OCC, IND, OCC2010, METRO, CITY, HHINCOME, SERIAL,HHWT, NUMPREC, SUBSAMP,HHTYPE   )

data <- rbind(data2019, data2021) #125,007 observations before any filtering. 


# replaces 0 with NA for variables listed. Allows topline to calculate "Valid Percent" when it recognizes missing values

data <- data %>% replace_with_na(replace = list(
  EMPSTAT= c(0), 
  LABFORCE=c(0), 
  CLASSWKR = c(0),
  OCCSOC = c(0),
  CIHISPEED = c(0),
  CINETHH = c(0),
  TRANWORK = c("N/A","0"))) %>% 
  filter(LABFORCE == 2 & INCEARN > 0) # in labor force and 18 years old and up abd positive earned incomes. 


data <- data %>% mutate(age_cat = 
                              case_when(AGE < 24 ~ "16to24",
                                        AGE > 24 & AGE < 35 ~ "25to34",
                                        AGE > 34 & AGE < 45 ~ "35to44",
                                        AGE > 44 & AGE < 55 ~ "45to54",
                                        AGE > 54 & AGE < 65 ~ "55to64",
                                        AGE > 64 ~ "65+"))

data <-  data %>% mutate(white = if_else(RACE ==1, 1, 0),
         black = if_else(RACE ==2, 1, 0), 
         asian = if_else(RACE %in% c(4,5,6), 1, 0),
         otherrace = if_else(RACE %in% c(3,7,8,9),1,0)) %>%
  group_by(COUNTYFIP,PUMA) %>%
  mutate(pct_white = sum(white)/n(),
         pct_black = sum(black)/n()) %>% 
  ungroup() %>%
  mutate(race_cat = case_when(
    RACE ==1~"White",
    RACE ==2 ~ "Black", 
    RACE %in% c(4,5,6)~"Asian",
    RACE %in% c(3,7,8,9)~"Other"))


## numbers used for income breaks are calculated in Income Deciles section. 
# created now so that the variable exists in the joined dataset before creating the survey design object

data <- data %>% 
  mutate(incdecile_w = case_when(
    INCEARN < 8000 ~ 1, 
    INCEARN >= 8000 & INCEARN < 18000 ~ 2,
    INCEARN >= 18000 & INCEARN < 26000 ~ 3,
    INCEARN >= 26000 & INCEARN < 35000 ~ 4,
    INCEARN >= 35000 & INCEARN < 43000 ~ 5,
    INCEARN >= 43000 & INCEARN < 54000 ~ 6,
    INCEARN >= 54000 & INCEARN < 68000 ~ 7,
    INCEARN >= 68000 & INCEARN < 85000 ~ 8,
    INCEARN >= 85000 & INCEARN < 120000 ~ 9,
    INCEARN >= 120000 ~ 10)
  ) %>%
  ## Padding FIPS code for merging with spatial geometry later
  mutate(county_pop_type = if_else(COUNTYFIP==0, 
                                   "Rural Counties", "Urban Counties")) %>%
  mutate(PUMA = str_pad(PUMA, 5, pad="0"),
         countyFIP = str_pad(COUNTYFIP, 3, pad = "0"))


data <- data %>%
  mutate(occ_2digits = substr(OCCSOC,1,2)) %>% 
  mutate(occ_2dig_labels = case_when(
    occ_2digits %in% c(11,13,19,15,17,19,21,23,25,27,29) ~ "Management, Business, Science, Arts",
    occ_2digits %in% c(31,32,33,34,35,36,37,38,39) ~ "Service Occupations",
    occ_2digits %in% c(41,42,43) ~ "Sales & Office Jobs",
    occ_2digits %in% c(45,46,47,48,49 ) ~"Natural Resources, Construction",
    occ_2digits %in% c(51, 52, 53) ~ "Production, Transportation",
    occ_2digits == 55 ~ "Military")) 

data <- data %>%
  mutate(occ_2digits = substr(OCCSOC,1,2)) %>% 
  mutate(occ_2dig_labels_d = case_when(
    occ_2digits %in% c(11) ~ "Management",
    occ_2digits %in% c(13) ~ "Business & Finance",
    occ_2digits %in% c(15) ~ "Computer, Engineering & Science",
    occ_2digits %in% c(17) ~ "Architecture & Engineering",
    occ_2digits %in% c(19) ~ "Life/Social Sciences",
    occ_2digits == 21 ~ "Community & Social Services",
    occ_2digits == 23 ~ "Legal",
    occ_2digits == 25 ~ "Educational Instruction",
    occ_2digits == 27 ~ "Arts, Design, Entertainainment",
    occ_2digits == 29 ~ "Health Practictioners",
    occ_2digits == 31 ~ "Healthcare Support",
    occ_2digits == 33 ~ "Protective services",
    occ_2digits == 35 ~ "Food Services",
    occ_2digits == 37 ~ "Building Cleaning & Maintenance",
    occ_2digits == 41 ~ "Sales",
    occ_2digits == 43 ~"Office & Administration",
    occ_2digits == 45 ~ "Farm, Fish, Forest",
    occ_2digits == 47 ~ "Construction & Extraction",
    occ_2digits == 49 ~"Installation, Maintenance",
    occ_2digits == 51 ~"Production",
    occ_2digits == 53 ~ "Transportation & Material Moving",
    occ_2digits == 55 ~ "Military",
    TRUE~"Other") )

data <- data %>% 
  mutate(did_wfh = if_else(TRANWORK==80, 1, 0)) # 1 = wfh, 0 = did not wfh

data <- data %>% 
  mutate(
    PWSTATE2 = ifelse(PWSTATE2 == 0, NA, PWSTATE2),
    work_in_IL = ifelse(PWSTATE2 == "17", "In Illinois", "Out of IL"),
    did_wfh_labels = ifelse(did_wfh == 1, "Did WFH", "Did not WFH"),
    has_incearn = ifelse(INCEARN > 0, 1, 0), ## has earned income = 1
    has_occsoc = ifelse(OCCSOC > 0, 1, 0),# has occupation = 1
    has_incearn_labels = ifelse(INCEARN > 0, "Has EarnInc", "No IncData"), ## has earned income = 1
    has_occsoc_labels = ifelse(OCCSOC > 0, "Has Occ", "No Occ") ## OCCSOC code greater than zero coded as 1
    )

rm(ddi)
rm(data2019)
rm(data2021)
```

## Descriptive Statistics

251630 observations originally in ACS 2021 and 2019 1-year samples.

123,753 observations after removing those not in the labor force and with earned incomes less than or equal to zero.

::: {.callout-warning appearance="simple" icon="false"}
Note: Dropping groups of people using filter() from the sample will change the standard errors of estimates since it changes the sample size. Use the survey() or svy() command to drop subsets of people (like if we wanted to filter age groups). Google what commands to use to drop observations without impacting standard errors.
:::

```{r examine-didwfh}

# data %>% 
#   group_by(YEAR, did_wfh) %>%
#   dplyr::summarize(weightedcount=sum(PERWT),
#                    unweightedcount = n()) %>%  #weighted
#   mutate(pct_weighted = round(weightedcount/sum(weightedcount), digits = 3),
#          pct_noweight = round(unweightedcount/sum(unweightedcount), digits = 3))

#valid percent BEFORE final filtering of observations (income greater than zero)
data %>% filter(YEAR == 2019) %>%
topline(did_wfh_labels, weight = PERWT)

data %>% filter(YEAR == 2021) %>%
topline(did_wfh_labels, weight = PERWT)

```

18.2% of Illinois workers worked at home. 76.3% went to work using some form of transportation, 5.5% of observations were missing values.

-   **Valid percent:** 19.2% of observations with responses did WFH and 80.8% of observations with responses did not WFH. **This is the statistic we would use**
-   Increased from 5.3% of all workers that worked from home in 2019

```{r could-wfh}
# bring in the teleworkable scores based on D&N's work.
telework <- read_csv("./data/teleworkable_AWM.csv")
joined <- left_join(data, telework, by = c("OCCSOC" = "occ_codes"))

#May 22 2023, Changed 399011 occupation to 0. Was coding Nannies and Child care as teleworkable.
joined <- joined %>% mutate(teleworkable = ifelse(OCCSOC == "399011" | OCCSOC == "399010", 0, teleworkable))

#table(joined$teleworkable)
# mostly 0's and 1's.
#hist(joined$teleworkable)

joined <- joined %>% 
  mutate(CanWorkFromHome = case_when(
  teleworkable == 0 ~ "No WFH",
  teleworkable < 1 ~ "Some WFH",
  teleworkable == 1 ~ "Can WFH",
  TRUE ~ "Check Me")) %>% 
  # keeps observations that have earned income values and are in the labor force.
  filter(has_incearn == 1 & LABFORCE == 2) 


table(joined$CanWorkFromHome)


table(joined$did_wfh, joined$YEAR)


```

```{r eval=FALSE, include=FALSE}

# checks to make sure percentages calculated using joined dataframe are the same as the data dataframe. 

#valid percent
joined %>% filter(YEAR == 2019) %>%
topline(did_wfh_labels, weight = PERWT)

joined %>% filter(YEAR == 2021) %>%
topline(did_wfh_labels, weight = PERWT)




# data %>% 
#   group_by(YEAR, did_wfh) %>%
#   dplyr::summarize(weightedcount=sum(PERWT),
#                    unweightedcount = n()) %>%  #weighted
#   mutate(pct_weighted = round(weightedcount/sum(weightedcount), digits = 3),
#          pct_noweight = round(unweightedcount/sum(unweightedcount), digits = 3))
```

**Pre 5.23.2023**: 40,278 observations have occupations that indicate that they could Work from home, 67,339 responses have occupations indicating that they can not work from home, 16,136 might be able to work from home.

**5.23.2023**: 41,643 can wfh; 66111 no wfh possible, 15,999 some wfh

### Survey Design

There are three different versions of the data: dstrata has both 2019 and 2021 data together (combined using rbind() above.) Using this combined design object makes some graphs easier, but I think it changes any standard errors used in estimate. I also made the 2019 and 2021 design strata from the the separate survey data. They all should have the same variables.

> After comparing the 2021 5-year sample, it is possible to just use that 5-year sample download and use only the 2019 and 2021 years of observations. This automatically adjusts the dollar amounts to 2021 dollars. It is also possible to use Pre-COVID and post-COVID years to have even more observations. (pre-covid would be 2017,2018, and 2019, post-covid would be 2020, 2021)

Deciles were created from the strata in different ways below. Depending on age ranges kept, deciles shift slightly.

2019 Strata by itself:\
4 9000 17000 25000 32000 40000 50000 63000 80000 114000 933000

2019 Strata income decile breaks (using the 2021 5-year extract):\
4 9538 18016 26494 33912 42390 52988 66765 84781 121873 988757.\
**Already adjusted for inflation??????**

-   9000, 9400, 25000, 32000, 40000, 50000, 63000, 80000, 1140000, 949000

-   8000, 18000, 26000, 35000, 43000, 54000, 68000, 85000, 120000

    \

::: {.callout-warning appearance="simple" icon="false"}
survey::svyquantile() uses the survey design to calculate deciles. These deciles that are created are slightly different than the deciles assigned using ntile(). I trust survey::svyquantile more because I know it applies the weights to observations when creating the deciles for earned income.
:::

::: {.callout-important appearance="minimal" icon="false"}
the ntile() decile variable is what is used to graph all decile images and for all income decile tables/calculations. ~~Find a way to apply the summarized version of the svyquantiles as a new variable in the dstrata datasets.~~ - Done.

`incdecile_w` is made using the breaks returned from the svyquantile command.

-   ~~now use this variable in graphs? or create the summary tables using weights and then pass that to the graph commands.~~ - done
:::

```{r survey-design, message = FALSE, warning=FALSE}
#as_survey() from srvyr package

## both years together: calculations using this will have incorrect standard errors
# might be easier sometimes to graph together. Maybe. 
joined <- joined %>% filter(HHINCOME > 0 & HHINCOME!= 9999999 & HHINCOME != 9999998)  # 105 observations 
joined <- joined %>%
  mutate(HHincdecile_w = case_when(
    INCEARN < 34000 ~ 1,
    INCEARN >= 34000 & INCEARN < 51900 ~ 2,
    INCEARN >= 51900 & INCEARN < 68000 ~ 3,
    INCEARN >= 68000 & INCEARN < 83600 ~ 4,
    INCEARN >= 83699 & INCEARN < 100000 ~ 5,
    INCEARN >= 100000 & INCEARN < 120000 ~ 6,
    INCEARN >= 120000 & INCEARN < 142400 ~ 7,
    INCEARN >= 142400 & INCEARN < 175000 ~ 8,
    INCEARN >= 175000 & INCEARN < 235000 ~ 9,
    INCEARN >= 235000 ~ 10)
  )

dstrata <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~PERWT, data = joined) %>% 
  as_survey() %>%
  mutate(decile = ntile(INCEARN, 10))



# 2019 data turned into survey item
dstrata2019 <- joined %>% filter(YEAR==2019) 
dstrata2019 <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, 
                                 weights = ~PERWT, data = dstrata2019) %>% 
  as_survey() %>% 
  mutate(decile = ntile(INCEARN, 10))


dstrata2021 <- joined %>% filter(YEAR==2021) 

dstrata2021 <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~PERWT, data = dstrata2021) %>% as_survey() %>%
  mutate(decile = ntile(INCEARN, 10))

# deciles using ntile(). Not weighted!!  Close to income deciles from weighted suvey design though.




```

```{r export-occ-matched-csv, eval=FALSE, include=FALSE}
#Tables for Appendix Items with counts of observations by OCC types:
all_occs <- joined %>% 
  filter(YEAR == 2021)%>%
  group_by(OCCSOC) %>% 
  summarize(count_matched = n(),
            weighted_estimate = sum(PERWT)) %>%
  mutate(perc_responses = scales::percent(count_matched/sum(count_matched)),
            perc_workforce =scales::percent(weighted_estimate/sum(weighted_estimate)))

occs_23cat <-joined %>% 
  filter(YEAR == 2021)%>%
  dplyr::group_by(occ_2digits,occ_2dig_labels_d) %>% 
  dplyr::summarize(
    count_matched = n(),
            weighted_estimate = sum(PERWT),
    avg_teleworkable = mean(teleworkable)) %>%
  ungroup() %>%
  mutate(perc_responses = count_matched/sum(count_matched),
        perc_workforce = weighted_estimate/sum(weighted_estimate)) %>% arrange(-count_matched)

occ_5cat <- joined %>% 
  filter(YEAR == 2021)%>%
  dplyr::group_by(occ_2dig_labels) %>% 
  dplyr::summarize(
    count_matched = n(),
            weighted_estimate = sum(PERWT),
    avg_teleworkable = mean(teleworkable)) %>%
  ungroup() %>%
  mutate(perc_responses = count_matched/sum(count_matched),
        perc_workforce = weighted_estimate/sum(weighted_estimate)) %>% arrange(-count_matched)

#install.packages("openxlsx")
library(openxlsx)

dataset_names <- list('All Matched OCC Codes' = all_occs,
                      '5 Broadest OCC Groups' = occ_5cat, 
                      'OCC 2-digit Groups'= occs_23cat
                      )

write.xlsx(dataset_names, file = 'WFH_feasibility_comparison_AWM.xlsx')
```

## Comparison of Women w/ & w/o Kids under 5

Women under 40 who could work from home in the occupation category with large amount of people who could work from home (Management, Business, Science & Arts, all SOC 2-digit codes \< 30-0000)

Under5 is binary variable: 0 means no children under 5. 1 means at least 1 or more children under 5.

318,174 women in 2019 worked from home and 1,140,835 women in 2021 worked from home.

```{r}
childttest2021<-joined %>% 
  filter(YEAR == 2021)%>%
  filter(SEX == 2 & AGE < 40) %>%
  filter(CanWorkFromHome == "Can WFH" & occ_2digits <30)  %>%  
  mutate(under5 = ifelse(NCHLT5 > 0, 1, 0))
write.csv(childttest2021, "with_kids_comparison.csv")



childttest2019<-joined %>% 
  filter(YEAR == 2019)%>%
  filter(SEX == 2 & AGE < 40) %>%
  filter(CanWorkFromHome == "Can WFH" & occ_2digits <30)  %>%  
  mutate(under5 = ifelse(NCHLT5 > 0, 1, 0))

kidsdesign <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~PERWT, data = childttest2021)

kidsdesign2 <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~PERWT, data = childttest2019)

#withkids <- childttest %>% filter(under5 == 1)
#table(withkids$did_wfh_labels)
#nokids <- childttest %>% filter(under5==0)
#table(nokids$did_wfh_labels)


svytable(~under5+did_wfh_labels, kidsdesign) 

svytable(~under5+did_wfh_labels, kidsdesign) %>% as_tibble() %>%
 #  group_by(did_wfh_labels) %>% 
     mutate(Prop = round(n/sum(n), digits =3))

svytable(~under5+did_wfh_labels, kidsdesign) %>% summary()
# also chi square test?
svychisq(~under5+did_wfh_labels, kidsdesign, statistic = "Chisq")

# t-test between with kids and without kids group. Weighted.
# t-test probably isn't right because it isn't a continuous variable.
# svyttest(under5~did_wfh, kidsdesign)


#svytable(~did_wfh_labels+YEAR+under5, design = dstrata)
svytable(~under5+did_wfh_labels, design = kidsdesign) %>% summary()
svytable(~under5+did_wfh_labels, design = kidsdesign2) %>% summary()



 #svychisq(under5~did_wfh_labels, design = kidsdesign)
#crosstab(did_wfh_labels, under5, weight = PERWT, unwt_n = TRUE, df = childttest2021)

#crosstab(did_wfh_labels,under5, weight = PERWT, unwt_n = TRUE, df = childttest2019)


### Switch variable order.
crosstab(under5, did_wfh_labels, weight = PERWT, unwt_n = TRUE, df = childttest2021)

crosstab(under5,did_wfh_labels, weight = PERWT, unwt_n = TRUE, df = childttest2019)
#crosstab_3way(YEAR, under5, did_wfh_labels, weight = PERWT, unwt_n = TRUE, df = childttest#, pct_type = "cell")

#glm(did_wfh~under5, data = childttest2021) %>% summary()
#svyglm(did_wfh~under5, design = kidsdesign) %>% summary()
```

For men in 2019, 7.1% of men without kids worked from home. 7.6% of men with kids worked from home.

2019:\
21.8% of men who did not WFH had kid under 5.\
22.8% of men who did WFH had a kid under 5.\
7.2% of men who did not have kids under 5 worked from home.\
7.6% of men with kids under 5 worked from home.

In 2019:\
22% of women who did not WFH had a kid under 5.\
24% of women who did WFH had a kid under 5.\
5% of women without children under 5 worked from home.\
9% of women with children under 5 worked from home.

In 2021:\|\
21% of women who did not WFH had a kid under 5.\
34% of women who did WFH had a kid under 5.\
34% of women without children under 5 worked from home.\
37% of women with children under 5 worked from home.

2021:\
20.4% of men who did not work from home had kids under 5.\
20.45% of men who did work from home had kids under 5.\
40.45% of men without kids under 5 worked from home.\
40.55% of men with kids under 5 worked from home.

Not a statistically significant difference but still interesting.

**For women under 45:**

```{r}
childttest2021<-joined %>% 
  filter(YEAR == 2021)%>%
  filter(SEX == 2 & AGE < 45) %>%
  filter(CanWorkFromHome == "Can WFH" & occ_2digits <30)  %>%  
  mutate(under5 = ifelse(NCHLT5 > 0, 1, 0))

childttest2019<-joined %>% 
  filter(YEAR == 2019)%>%
  filter(SEX == 2 & AGE < 45) %>%
  filter(CanWorkFromHome == "Can WFH" & occ_2digits <30)  %>%  
  mutate(under5 = ifelse(NCHLT5 > 0, 1, 0))

kidsdesign <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~PERWT, data = childttest2021)

#withkids <- childttest %>% filter(under5 == 1)
#table(withkids$did_wfh_labels)
#nokids <- childttest %>% filter(under5==0)
#table(nokids$did_wfh_labels)


svytable(~under5+did_wfh_labels, kidsdesign)
svychisq(~under5+did_wfh_labels, kidsdesign)
svyttest(under5~did_wfh, kidsdesign)


#svytable(~did_wfh_labels+YEAR+under5, design = dstrata)
#svytable(~under5+YEAR+did_wfh_labels, design = dstrata)


#svychisq(did_wfh_labels~under5+YEAR, design = dstrata)
#crosstab(did_wfh_labels, under5, weight = PERWT, unwt_n = TRUE, df = childttest2021)

#crosstab(did_wfh_labels,under5, weight = PERWT, unwt_n = TRUE, df = childttest2019)
#crosstab_3way(YEAR, under5, did_wfh_labels, weight = PERWT, unwt_n = TRUE, df = childttest#, pct_type = "cell")

#glm(did_wfh~under5, data = childttest2021) %>% summary()

```

**For MEN under 40:**

In 2019:\
20.8% who did not WFH had a kid under 5.\
21.5% who did WFH had a kid under 5.

In 2021: 22% of men that did and did not WFH had a kid under 5.

```{r}
childttest2021<-joined %>% 
  filter(YEAR == 2021)%>%
  filter(SEX == 1 & AGE < 40) %>%
  filter(CanWorkFromHome == "Can WFH" & occ_2digits <30)  %>%  
  mutate(under5 = ifelse(NCHLT5 > 0, 1, 0))
write.csv(childttest2021, "with_kids_comparison.csv")



childttest2019<-joined %>% 
  filter(YEAR == 2019)%>%
  filter(SEX == 1 & AGE < 40) %>%
  filter(CanWorkFromHome == "Can WFH" & occ_2digits <30)  %>%  
  mutate(under5 = ifelse(NCHLT5 > 0, 1, 0))

kidsdesign <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~PERWT, data = childttest2021)

kidsdesign2 <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~PERWT, data = childttest2019)

#withkids <- childttest %>% filter(under5 == 1)
#table(withkids$did_wfh_labels)
#nokids <- childttest %>% filter(under5==0)
#table(nokids$did_wfh_labels)


svytable(~under5+did_wfh_labels, kidsdesign) 

svytable(~under5+did_wfh_labels, kidsdesign) %>% as_tibble() %>%
 #  group_by(did_wfh_labels) %>% 
     mutate(Prop = round(n/sum(n), digits =3))

svytable(~under5+did_wfh_labels, kidsdesign) %>% summary()
# also chi square test?
svychisq(~under5+did_wfh_labels, kidsdesign, statistic = "Chisq")

# t-test between with kids and without kids group. Weighted.
# t-test probably isn't right because it isn't a continuous variable.
# svyttest(under5~did_wfh, kidsdesign)


#svytable(~did_wfh_labels+YEAR+under5, design = dstrata)
svytable(~under5+did_wfh_labels, design = kidsdesign) %>% summary()
svytable(~under5+did_wfh_labels, design = kidsdesign2) %>% summary()



 #svychisq(under5~did_wfh_labels, design = kidsdesign)
#crosstab(did_wfh_labels, under5, weight = PERWT, unwt_n = TRUE, df = childttest2021)

#crosstab(did_wfh_labels,under5, weight = PERWT, unwt_n = TRUE, df = childttest2019)


### Switch variable order.
crosstab(under5, did_wfh_labels, weight = PERWT, unwt_n = TRUE, df = childttest2021)

crosstab(under5,did_wfh_labels, weight = PERWT, unwt_n = TRUE, df = childttest2019)
#crosstab_3way(YEAR, under5, did_wfh_labels, weight = PERWT, unwt_n = TRUE, df = childttest#, pct_type = "cell")

#glm(did_wfh~under5, data = childttest2021) %>% summary()
#svyglm(did_wfh~under5, design = kidsdesign) %>% summary()
```

**For MEN under 45:**

```{r}
childttest2021<-joined %>% 
  filter(YEAR == 2021)%>%
  filter(SEX == 1 & AGE < 45) %>%
  filter(CanWorkFromHome == "Can WFH" & occ_2digits <30)  %>%  
  mutate(under5 = ifelse(NCHLT5 > 0, 1, 0))

childttest2019<-joined %>% 
  filter(YEAR == 2019)%>%
  filter(SEX == 1 & AGE < 45) %>%
  filter(CanWorkFromHome == "Can WFH" & occ_2digits <30)  %>%  
  mutate(under5 = ifelse(NCHLT5 > 0, 1, 0))

kidsdesign <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~PERWT, data = childttest2021)

#withkids <- childttest %>% filter(under5 == 1)
#table(withkids$did_wfh_labels)
#nokids <- childttest %>% filter(under5==0)
#table(nokids$did_wfh_labels)


svytable(~under5+did_wfh_labels, kidsdesign)
svychisq(~under5+did_wfh_labels, kidsdesign)
svyttest(under5~did_wfh, kidsdesign)


#svytable(~did_wfh_labels+YEAR+under5, design = dstrata)
#svytable(~under5+YEAR+did_wfh_labels, design = dstrata)


#svychisq(did_wfh_labels~under5+YEAR, design = dstrata)
crosstab(did_wfh_labels, under5, weight = PERWT, unwt_n = TRUE, df = childttest2021)

crosstab(did_wfh_labels,under5, weight = PERWT, unwt_n = TRUE, df = childttest2019)
#crosstab_3way(YEAR, under5, did_wfh_labels, weight = PERWT, unwt_n = TRUE, df = childttest#, pct_type = "cell")

glm(did_wfh~under5, data = childttest2021) %>% summary()

```

## Household data

hhbreaks2019 = (32100, 50000, 65600, 80500, 97300, 115700, 137500, 168000, 225010) and max is 1671000

hhbreaks2021 = (34000, 51900, 68000, 83600, 1e+05 120,000 142,400 175,000 235,000) and max is 1,797,000.

```{r}
#| code-fold: true
#| output: hold

joined <- joined %>% filter(HHINCOME > 0 & HHINCOME!= 9999999 & HHINCOME != 9999998)  # 105 observations 

joined %>%
  ggplot() + geom_histogram(aes(x=HHINCOME, weight = HHWT))

joined %>% # 122000 observations 
ggplot() + geom_histogram(aes(x=HHINCOME))

HHdesign <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~HHWT, data = joined)

inc_quantiles <-survey::svyquantile(~HHINCOME, design=HHdesign, 
                    quantiles = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1) ,
                    na.rm=TRUE, ci = FALSE  )
inc_quantiles

# With  HH Weights:  4 32900 50000 66600 82000 99000 117000 140000 171000 230000 1797000
# With WRONG weights:  (34000, 51900, 68000, 83600, 1e+05 120,000 142,400 175,000 235,000
#Code done above when creating variables in beginning chunks.
joined <- joined %>%
  mutate(HHincdecile_w = case_when(
    INCEARN < 32900 ~ 1,
    INCEARN >= 32900 & INCEARN < 50000 ~ 2,
    INCEARN >= 50000 & INCEARN < 66600 ~ 3,
    INCEARN >= 66600 & INCEARN < 82000 ~ 4,
    INCEARN >= 82000 & INCEARN < 99000 ~ 5,
    INCEARN >= 99000 & INCEARN < 117000 ~ 6,
    INCEARN >= 117000 & INCEARN < 140000 ~ 7,
    INCEARN >= 140000 & INCEARN < 171000 ~ 8,
    INCEARN >= 171000 & INCEARN < 230000 ~ 9,
    INCEARN >= 230000 ~ 10)
  )

HHdesign <- survey::svydesign(id = ~CLUSTER, strata = ~STRATA, weights = ~HHWT, data = joined)
table <- svytable(~YEAR+HHincdecile_w+did_wfh_labels,  design = HHdesign) 

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR,HHincdecile_w)%>%
  mutate(Prop=round(n/sum(n), digits=3)) %>%
  filter(did_wfh_labels == "Did WFH")

table # has proportions calculated out of TOTAl for both years

table %>%
  ggplot(aes(factor(HHincdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, fill = YEAR, group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(stat="identity", position = "dodge")+
  #geom_col(stat = "identity", position = "stack") +   # scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+

 # facet_wrap(~YEAR)+
   coord_flip()+
    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_dodge(width = 0.8), hjust = 1.1,
              size = 4) + 
  labs(title ="Percent of each HOUSEHOLD income decile that did WFH",
  subtitle = "2019 vs 2021",
       caption = "ACS 1 year samples for 2019 and 2021. Working from home based on TRANWORK question on commuting and HHINCOME variable.
       All workers in the labor force, all ages included.
       Income based on HHINCOME for household income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of workers working from home") + 
 theme(legend.position = "bottom", legend.title = element_blank())+
  theme_classic()+
    scale_fill_manual(values = c("#a6bddb", "#2b8cbe")) +
  scale_y_continuous(labels = scales::percent)


```

## Individual Level Data

### Occupations

Combined into 6 major occupation groups. Broadest categories are made up of multiple 2-digit OCCSOC codes.

```{r}
#| label: fig-A2-Figure1
#| fig-cap: "Proportion of Occupation Types in  Illinois by Aggregated Occupation Groups used by BLS & ACS"
#| fig-subcap: "ACS 1 year samples for 2019 and 2021 used for weighted population estimates. Military occupations make up less than 0.5% of the labor force and were removed from the  graph. Occupation categories are based on broadest aggregated BLS categories used by the BLS."
#| code-fold: true

#table includes observations from BOTH years.
#table(joined$occ_2dig_labels, joined$did_wfh_labels)


#table(joined$occ_2digits)

crosstab_3way(joined, YEAR, occ_2dig_labels, did_wfh_labels, weight = PERWT)

table <- svytable(~YEAR+occ_2dig_labels, design = dstrata) 
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>%
  arrange(-n)
table

table %>% filter(occ_2dig_labels != "Military") %>% ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels)), y=n, group = YEAR)) + 
  geom_col(stat = "identity", fill="lightblue") +
  facet_wrap(~YEAR)+
        geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.02,Prop, "")), accuracy = .1),accuracy  = .1L ),position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "bottom", legend.title = element_blank(),
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.background = element_rect(fill='transparent'),
         legend.box.background = element_rect(fill='transparent')
       )+
  labs(
    title ="Proportion of Occupation Types in  Illinois",
    subtitle = "By Most Aggregated Occupation Groups used by BLS & ACS",
       x = "", y = "Estimated Number of Workers")+
 scale_y_continuous(labels = scales::comma)+
scale_x_discrete(labels = function(x) str_wrap(x, width=25))+ # makes labels better on axsis

  coord_flip()
```


```{r export-figA2-Fig1, include = FALSE}
table %>% filter(occ_2dig_labels != "Military") %>% ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels)), y=n, group = YEAR)) + 
  geom_col(stat = "identity", fill="lightblue") +
  facet_wrap(~YEAR)+
        geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.02,Prop, "")), accuracy = .1),accuracy  = .1L ),position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "bottom", legend.title = element_blank(),
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.background = element_rect(fill='transparent'),
         legend.box.background = element_rect(fill='transparent')
       )+
  labs(
  #  title ="Proportion of Occupation Types in  Illinois",
  #  subtitle = "By Most Aggregated Occupation Groups used by BLS & ACS",
       x = "", y = "Estimated Number of Workers")+
 scale_y_continuous(labels = scales::comma)+
scale_x_discrete(labels = function(x) str_wrap(x, width=25))+ # makes labels better on axsis

  coord_flip()

ggsave("./paper_figures/A2-Figure1.png", limitsize = FALSE, width = 8, height = 4, units = "in")

ggsave("./paper_figures/A2-Figure1.eps", limitsize = FALSE, width = 8, height = 4, units = "in")

```

ACS 1 year samples for 2019 and 2021 used for weighted population estimates. Military occupations make up less than 0.5% of the labor force and were removed from the graph. Occupation categories are based on broadest aggregated BLS categories used by the BLS.

```{r}
#| label: fig-A2-Figure2
#| fig-cap: "ACS 1 year samples for 2019 and 2021 used for weighted population estimates. Occupation categories based on first 2 digits of OCCSOC occupation codes. Labels for occupations that make up less than 2% of the workers were not labeled for legibility reasons."


### Proportion of All Workers in each Occupation Type ###
table <- svytable(~YEAR+occ_2dig_labels_d, design = dstrata) 
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>%
  arrange(-n)
table

table %>% ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels_d)), y=n, group = YEAR)) + 
  geom_col(stat = "identity", fill="lightblue") +
  facet_wrap(~YEAR)+
        geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.02,Prop, "")), accuracy = .1),accuracy  = .1L ),position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "bottom", legend.title = element_blank(),
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.background = element_rect(fill='transparent'),
         legend.box.background = element_rect(fill='transparent')
       )+
  labs(title ="Occupations in Illinois",
       #subtitle = "Little change between 2019 and 2021 Occurred",
       x = "", y = "Estimated Number of Workers") + scale_y_continuous(labels = scales::comma)+
  
  coord_flip()
```

```{r export-A1figure2, include = FALSE}
### Proportion of All Workers in each Occupation Type ###
table <- svytable(~YEAR+occ_2dig_labels_d, design = dstrata) 
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>%
  arrange(-n)
table

table %>% ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels_d)), y=n, group = YEAR)) + 
  geom_col(stat = "identity", fill="lightblue") +
  facet_wrap(~YEAR)+
        geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.02,Prop, "")), accuracy = .1),accuracy  = .1L ),position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "bottom", legend.title = element_blank(),
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.background = element_rect(fill='transparent'),
         legend.box.background = element_rect(fill='transparent')
       )  + 
  labs(#title ="WFH Feasibility by Occupation Type",
       x = "", y = "Estimated Number of Workers") + 
  scale_y_continuous(labels = scales::comma)+
  coord_flip()

ggsave("./paper_figures/A2-Figure2.png", limitsize = FALSE,width = 8, height = 4, units = "in")

ggsave("./paper_figures/A2-Figure2.eps", limitsize = FALSE,width = 8, height = 4, units = "in")
```

```{r}
### Percent of Workers working from home within each Broad Occupation type ###
table <- svytable(~YEAR+did_wfh_labels+occ_2dig_labels, design = dstrata) 

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, occ_2dig_labels)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>%
  arrange(did_wfh_labels, -n)
table

table %>% filter(occ_2dig_labels != "Military" ) %>%
  ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels)), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +facet_wrap(~YEAR)+
        geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.05,Prop, "")), accuracy = .1),accuracy  = .1L ),position = position_stack(vjust=.5), size=3) +
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Proportion of Workers in each Occupation Who Did WFH",
       subtitle = "Percentages add to 100% within each occupation",
       x = "", y = "Estimated # of People",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates. Military occupations were excluded from graph due to low occurance of observations.") + scale_y_continuous(labels = scales::comma)+
  scale_fill_manual(values = c("#a6bddb", "#2b8cbe")) +
  coord_flip()
```

```{r}
#|label: fig-Figure3
#|
## Proportion of all workers in each occupation cateogory.## 
table <- svytable(~YEAR+did_wfh_labels+occ_2dig_labels, design = dstrata) 

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>%
  arrange(did_wfh_labels, -n)
table
table %>% 
  filter(occ_2dig_labels != "Military") %>%
  ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels)), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
  geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.02,Prop, "")), accuracy = .1), accuracy  = .1L ),
            position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "none", 
        legend.title = element_blank(), 
        plot.title.position = "plot",
        panel.background = element_rect(fill='transparent'), #transparent panel bg
        plot.background = element_rect(fill='transparent', color=NA) #transparent plot bg
  )+
  labs(title ="Proportion of Illinois Workforce Who Worked From Home",
        subtitle = "All workers in labor force with occsoc codes in a year add to 100%",
       x = "", y = "Estimated Number of Workers") +#,
  # caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") 
scale_y_continuous(labels = scales::comma) +
  
    scale_x_discrete(labels = function(x) str_wrap(x, width=25))+ # makes labels better on axsis

  scale_fill_manual(values = c("#a6bddb",  "#2b8cbe")) + coord_flip()
```


```{r export-Figure3, include=FALSE}
table %>% 
  filter(occ_2dig_labels != "Military") %>%
  ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels)), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
  geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.02,Prop, "")), accuracy = .1), accuracy  = .1L ),
            position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "none", 
        legend.title = element_blank(), 
        plot.title.position = "plot",
        panel.background = element_rect(fill='transparent'), #transparent panel bg
        plot.background = element_rect(fill='transparent', color=NA) #transparent plot bg
  )+
  labs(# title ="Proportion of Illinois Workforce Who Worked From Home",
       # subtitle = "All workers in labor force with occsoc codes in a year add to 100%",
       x = "", y = "Estimated Number of Workers") +#,
  # caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") 
scale_y_continuous(labels = scales::comma) +
  
    scale_x_discrete(labels = function(x) str_wrap(x, width=25))+ # makes labels better on axsis

  scale_fill_manual(values = c("#a6bddb",  "#2b8cbe")) + coord_flip()

ggsave("./paper_figures/Figure3.eps", limitsize = FALSE,width = 8, height = 4, units = "in")
ggsave("./paper_figures/Figure3.png", limitsize = FALSE, width = 8, height = 4, units = "in")
```

```{r}
#| label: fig-A2-figure4


### Detailed Occuation Types ## 
table <- svytable(~YEAR+did_wfh_labels+occ_2dig_labels_d, design = dstrata) 

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>%
  arrange(did_wfh_labels, -n)
table

table %>%ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels_d)), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.01,Prop, "")), accuracy = .1),accuracy  = .1L ),position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "bottom", legend.title = element_blank(),
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.background = element_rect(fill='transparent'),
       )   + 
  labs(title ="Percent working from home by 2-digit Occupation type",
       x = "", y = "Estimated Number of Workers") + 
  scale_y_continuous(labels = scales::comma)+
  scale_fill_manual(values = c("#a6bddb", "#2b8cbe")) +
  coord_flip()


```

```{r export-a2-fig4}
table %>%ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels_d)), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.01,Prop, "")), accuracy = .1),accuracy  = .1L ),position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "bottom", legend.title = element_blank(),
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.background = element_rect(fill='transparent'),
       )   + 
  labs(#title ="Percent working from home by 2-digit Occupation type",
       x = "", y = "Estimated Number of Workers") + 
  scale_y_continuous(labels = scales::comma)+
  scale_fill_manual(values = c("#a6bddb", "#2b8cbe")) +
  coord_flip()

ggsave("./paper_figures/A2-Figure4.png", limitsize = FALSE, width = 8, height = 4, units = "in")
ggsave("./paper_figures/A2-Figure4.eps", limitsize = FALSE, width = 8, height = 4, units = "in")
```

ACS 1 year samples for 2019 and 2021 used for weighted population estimates. Graph interpretation: 3.6% of all worker in the labor force in 2021 were in Management occupations and worked from home. 8.3% of all workers were in management and did not work from home. Workers in Management occupations make up 11.9% of the entire workforce.

```{r}
#| label: fig-A2-Figure3
#|
# Both years, detailed observation types
table <- svytable(~YEAR+CanWorkFromHome+occ_2dig_labels_d, design = dstrata) 

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  arrange(CanWorkFromHome,-n) %>%
  mutate(Prop =round(n/sum(n), digits=3))
table

table %>% 
  ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels_d)), y=n, fill = CanWorkFromHome, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
  geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.01,Prop, "")), accuracy = .1), accuracy= 0.1L ),position = position_stack(vjust=.5), size=3) +
theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent of Illinois Workers that Could WFH by Occupation Type",
       x = "", y = "Estimated Number of Workers") + 
  scale_y_continuous(labels = scales::comma)+
  scale_fill_manual(values = c( "#117733","#44AA99","#D8D5C5")) +
  coord_flip()
```


```{r export-A2-Figure3, include=FALSE}
table %>% 
  ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels_d)), y=n, fill = CanWorkFromHome, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
  geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.01,Prop, "")), accuracy = .1), accuracy= 0.1L ),position = position_stack(vjust=.5), size=3) +
theme_classic()  + 
  theme(legend.position = "bottom", legend.title = element_blank(),
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.background = element_rect(fill='transparent')       ) +
  labs(#title ="Percent of Illinois Workers that Could WFH by Occupation Type",
       x = "", y = "Estimated Number of Workers") + 
  scale_y_continuous(labels = scales::comma)+
  scale_fill_manual(values = c( "#117733","#44AA99","#D8D5C5")) +
  coord_flip()

ggsave("./paper_figures/A2-Figure3.png", limitsize = FALSE, width = 8, height = 4, units = "in")
ggsave("./paper_figures/A2-Figure3.eps", limitsize = FALSE, width = 8, height = 4, units = "in")
```

OCCSOC codes and Teleworkable scores from occupation characteristics. 11.6% of all workers in Illinois had management occupations (6.6 Can WFH + 1.8 No WFH + 3.2 Some WFH in 2021). 6.6% of all workers in Illinois had management occupations and could feasibly WFH. ACS 1 year samples for 2019 and 2021 used for weighted population estimates.

```{r}
#| label: fig-wfh-feasibility-byoccupation
#| fig-cap: "Percent of Workers that Could Feasibily Work From Home in 2021"
table <- svytable(~YEAR+CanWorkFromHome+occ_2dig_labels_d, design = dstrata) 

table <- table %>% 
  as_tibble() %>% 
  filter(YEAR==2021)%>% 
  arrange(CanWorkFromHome, -n) %>%
  mutate(Prop =round(n/sum(n), digits=3))


table %>% 
  ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels_d)), y=n, fill = CanWorkFromHome, group=YEAR)) + 
  geom_col(position="stack", stat = "identity")+
  geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.01,Prop, "")), accuracy = .1), accuracy= 0.1L ), position = position_stack(vjust=.5), size=3) +
theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(x = "", y = "# of People",
      caption = "Occupation codes (OCCSOC)from the 2021 1-year ACS merged with work from home feasibility scores."
      ) + scale_y_continuous(labels = scales::comma)+
  scale_fill_manual(values = c( "#2b8cbe","#a6bddb","gray80")) +
  coord_flip()
```

```{r}
#| label: fig-Figure1
#| fig-cap: "WFH Feasibility for Workers in Illinois by Broadest Occupation Categories"

## Proportion of all workers in each occupation cateogory.## 
table <- svytable(~YEAR+CanWorkFromHome+occ_2dig_labels, design = dstrata) 

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>%
  arrange(CanWorkFromHome, -n)
table

Figure1 <- table %>% filter(occ_2dig_labels != "Military") %>%
  ggplot(aes(x=fct_rev(fct_inorder(occ_2dig_labels)), y=n, fill = CanWorkFromHome, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
  geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.02,Prop, "")), accuracy = .1),accuracy  = .1L ), 
            position = position_stack(vjust=.5), size=2.5) +
    theme_classic() +
  theme(plot.title = element_text(hjust=0), legend.position = "bottom", 
     #   legend.title = element_text(text = "WFH Feasibility")
        # legend.title = element_blank()
        )+
  labs(#title ="",
     #  subtitle = "All workers in labor force with occsoc codes in a year add to 100%",
       x = "", y = "Estimated Number of Workers",
   #   caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates"
     ) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width=25))+ # makes labels better on axis
  scale_y_continuous(labels = scales::comma) +
  scale_fill_manual(name = "WFH Feasibility", values = c( "#2b8cbe","#a6bddb", "gray89")) +
  coord_flip()  


Figure1

ggsave("./paper_figures/Figure1.eps", limitsize = FALSE,width = 8, height = 4, units = "in")#
ggsave("./paper_figures/Figure1.png", limitsize = FALSE, width = 8, height = 4, units = "in")


```

### Could WFH vs. did WFH

```{r}
#| label: fig-Figure7
#| fig-cap: "Percentages add up to 100 when adding all workers together within a year. Did work from home based on TRANWORK==80 variable from ACS surveys. Can Work from home based on teleworkable classification in Dingel & Niemen (2020). ACS 1 year samples for 2019 and 2021 used for weighted population estimates."
#| code-fold: true

## Totals add up to total number of workers in a year
table <- svytable(~CanWorkFromHome+YEAR+did_wfh_labels, design = dstrata) 
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% # will divide by all workers per year
  mutate(Prop =round(n/sum(n), digits=3)) %>%
  arrange(did_wfh_labels, -n)

table

## percentages add up to 100 when adding all workers together for a year
table %>% ggplot(aes(CanWorkFromHome, y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) +
  #scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
  theme_classic() + theme(legend.position = "none", legend.title = element_blank())+
  labs(title ="Did those who could work from home actually work from home?",subtitle = "2019 vs 2021",
       x = "", y = "Estimated Number of Workers") + 
  scale_y_continuous(labels = scales::comma) + 
  scale_fill_manual(values = c("#a6bddb",  "#2b8cbe"))
```


```{r export-Figure7, include=FALSE}
## percentages add up to 100 when adding all workers together for a year
table %>% ggplot(aes(CanWorkFromHome, y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) +
  #scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
  theme_classic() + 
  #theme(legend.position = "none", legend.title = element_blank()) +
  
    theme(legend.position = "none", legend.title = element_blank(),
         panel.background = element_rect(fill='transparent'),
         plot.background = element_rect(fill='transparent', color=NA),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.background = element_rect(fill='transparent')       ) +
  
  labs(#title ="Did those who could work from home actually work from home?",subtitle = "2019 vs 2021",
       x = "", y = "Estimated Number of Workers") + 
  scale_y_continuous(labels = scales::comma) + 
  scale_fill_manual(values = c("#a6bddb",  "#2b8cbe"))

ggsave("./paper_figures/Figure7.eps", limitsize = FALSE,width = 7, height = 5, units = "in")
ggsave("./paper_figures/Figure7.png", limitsize = FALSE, width = 7, height = 5, units = "in")

ggsave("./paper_figures/A1-Figure3.eps", limitsize = FALSE,width = 7, height = 5, units = "in")
ggsave("./paper_figures/A1-Figure3.png", limitsize = FALSE, width = 7, height = 5, units = "in")
```

```{r}
table <- svytable(~CanWorkFromHome+YEAR+did_wfh_labels, design = dstrata) 
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, CanWorkFromHome)%>% # divides by all workers per year within each categor
  mutate(Prop =round(n/sum(n), digits=3)) %>% 
  arrange(did_wfh_labels,-n)
table

xtabs(~did_wfh_labels+CanWorkFromHome+YEAR, data = dstrata)

## percentages add up to 100 when adding all workers together for a year
table %>% ggplot(aes(fct_inorder(CanWorkFromHome), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack")  +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_stack(vjust=.5), size=3) +
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Did those who COULD work from home actually work from home:", subtitle = "2019 vs 2021",
       x = "", y = "# of People",
      caption = "Comparison graph that might feel more correct.
      Percentages add up to 100 when adding all workers within each CanWorkFromHome category for each a year.
      Did work from home based on TRANWORK==80 variable from ACS surveys. Can Work from home based on teleworkable classification in Dingel & Niemen (2020).
            ACS 1 year samples for 2019 and 2021 used for weighted population estimates.") + scale_y_continuous(labels = scales::comma)+
  scale_fill_manual(values = c("#a6bddb",  "#2b8cbe")) 


table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, CanWorkFromHome)%>% # divides by all workers per year within each categor
  filter(CanWorkFromHome != "Some WFH") %>%
  mutate(Prop =round(n/sum(n), digits=3)) %>% 
  arrange(did_wfh_labels,-n)
table

xtabs(~did_wfh_labels+CanWorkFromHome+YEAR, data = dstrata)

## percentages add up to 100 when adding all workers together for a year
table %>% ggplot(aes(fct_inorder(CanWorkFromHome), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack")  +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_stack(vjust=.5), size=3) +
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Did those who COULD work from home actually work from home:", subtitle = "2019 vs 2021",
       x = "", y = "# of People",
      caption = "Comparison graph that might feel more correct.
      Percentages add up to 100 when adding all workers within each CanWorkFromHome category for each a year.
      Did work from home based on TRANWORK==80 variable from ACS surveys. Can Work from home based on teleworkable classification in Dingel & Niemen (2020).
            ACS 1 year samples for 2019 and 2021 used for weighted population estimates.") + scale_y_continuous(labels = scales::comma)+
  scale_fill_manual(values = c("#a6bddb",  "#2b8cbe")) 
# 
# table <- svytable(~CanWorkFromHome+YEAR+did_wfh_labels, design = dstrata) 
# # proportion of each respondant's sex and if they worked from home for each year in sample
# table <- table %>% 
#   as_tibble() %>% 
#   group_by(YEAR, did_wfh_labels)%>% # divides by all workers per year within each categor
#   mutate(Prop =round(n/sum(n), digits=3))


## percentages add up to 100 when adding all workers together for a year
## Don't like this version of the graph ### 
# table %>% ggplot(aes(did_wfh_labels, y=n, fill = CanWorkFromHome, group = YEAR)) + 
#   geom_col(stat = "identity", position = "stack") +
#   facet_wrap(~YEAR)+
#     geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) +
#   theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
#   labs(title ="Did those who COULD work from home actually work from home:", subtitle = "2019 vs 2021",
#        x = "", y = "# of People",
#       caption = "Comparison graph that might feel more correct.
#       Percentages add up to 100 when adding all workers within each did_wfh category for each a year.
#       Did work from home based on TRANWORK==80 variable from ACS surveys. Can Work from home based on teleworkable classification in Dingel & Niemen (2020).
#             ACS 1 year samples for 2019 and 2021 used for weighted population estimates.") + scale_y_continuous(labels = scales::comma)
```

> Ideally switch order of categories so that Did WFH is on the bottom of stack. Allows easier visual comparison.

### Creating Income Deciles

```{r figure6}
# svyquantile shows the breaks for the quantiles. hypothetically uses weights of observations for calculation of deciles.
# equal number of people should be in each decile after weights are applied 
inc_quantiles <- survey::svyquantile(~INCEARN, design=dstrata2019, 
    quantiles = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1) , na.rm=TRUE,  ci = FALSE  )
# $INCEARN for 2019
#[1,] 4 8000 16000 24000 30900 40000 50000 62000 80000 113000 933000
# values not adjusted to 2021 values.
inc_quantiles


inc_quantiles <-survey::svyquantile(~INCEARN, design=dstrata2021, 
                    quantiles = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1) ,
                    na.rm=TRUE, ci = FALSE  )
inc_quantiles
# $INCEARN
# [1,] 4 7200 16000 25000 34000 42000 52000 67000 85000 120000 949000



breaks2019 <- c(7200, 16000, 25000, 34000, 42000, 52000, 67000, 85000, 120000)

breaks2019adjusted <- c(8478, 16956, 25434, 32853, 42390, 52988, 65705, 84781, 120813) 
# from 2021 5 year sample and filtered for just 2019. 
# already adjusted for inflation.
# included for comparison and to decide to use 5 year ACS or 2019 and 2021 1 year ACS

breaks2021 <- c(8000, 18000, 26000, 35000, 43000, 54000, 68000, 85000, 120000)

# Code done above when creating variables in beginning chunks. 
# joined <- joined %>% 
#   mutate(incdecile_w = case_when(
#     INCEARN < 8000 ~ 1, 
#     INCEARN >= 8000 & INCEARN < 18000 ~ 2,
#     INCEARN >= 18000 & INCEARN < 26000 ~ 3,
#     INCEARN >= 26000 & INCEARN < 35000 ~ 4,
#     INCEARN >= 35000 & INCEARN < 43000 ~ 5,
#     INCEARN >= 43000 & INCEARN < 54000 ~ 6,
#     INCEARN >= 54000 & INCEARN < 68000 ~ 7,
#     INCEARN >= 68000 & INCEARN < 85000 ~ 8,
#     INCEARN >= 85000 & INCEARN < 120000 ~ 9,
#     INCEARN >= 120000 ~ 10)
#   )
# number of observations in each decile after weights used for creating the income deciles
#table(joined$incdecile_w)

# no major differnce between years in who COULD work from home based on teleworkable codes. Makes sense. 
# ggplot(joined, aes(teleworkable, weight = PERWT)) +
# geom_histogram()+facet_wrap(~YEAR)



table <- svytable(~YEAR+incdecile_w+did_wfh_labels,  design = dstrata) # proportion of each respondants sex in sample

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR,incdecile_w)%>%
  mutate(Prop=round(n/sum(n), digits=3)) %>%
  filter(did_wfh_labels == "Did WFH")

table # has proportions calculated out of TOTAl for both years

table %>%
  ggplot(aes(factor(incdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, fill = YEAR, group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(stat="identity", position = "dodge")+
  #geom_col(stat = "identity", position = "stack") +   # scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+

 # facet_wrap(~YEAR)+
   coord_flip()+
    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_dodge(width = 0.8), hjust = 1.1,
              size = 4) + 
  labs(title ="Working From Home by Earned Income Deciles",
  subtitle = "2019 vs 2021",
 #      caption = "ACS 1 year samples for 2019 and 2021. Working from home based on TRANWORK question on commuting.
#       All workers in the labor force, all ages included.
#       Income based on INCEARN for total earned income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of earners working from home") + 
 theme(legend.position = "none", legend.title = element_blank())+
  theme_classic()+
    scale_fill_manual(values = c("#a6bddb", "#2b8cbe")) +
  scale_y_continuous(labels = scales::percent)

decile_labels <- c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")
```


```{r export-figure6, include=FALSE}

table %>%
  ggplot(aes(factor(incdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, fill = YEAR, group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(stat="identity", position = "dodge")+
  #geom_col(stat = "identity", position = "stack") +   # scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+

 # facet_wrap(~YEAR)+
   coord_flip()+
    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_dodge(width = 0.8), hjust = 1.1,
              size = 4) + 
  labs(#title ="Working From Home by Earned Income Deciles",
  #subtitle = "2019 vs 2021",
 #      caption = "ACS 1 year samples for 2019 and 2021. Working from home based on TRANWORK question on commuting.
#       All workers in the labor force, all ages included.
#       Income based on INCEARN for total earned income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of earners working from home") + 

  theme_classic()+
   theme(legend.position = "none", legend.title = element_blank())+
    scale_fill_manual(values = c("#a6bddb", "#2b8cbe")) +
  scale_y_continuous(labels = scales::percent)

ggsave("./paper_figures/Figure6.eps", limitsize = FALSE, width = 8, height = 4, units = "in")
#ggsave("Figure6.pdf", limitsize = FALSE,width = 8, height = 4, units = "in")
ggsave("./paper_figures/Figure6.png", limitsize = FALSE, width = 8, height = 4, units = "in")



```

```{r eval=FALSE, include=FALSE}

table18up<-subset(dstrata,AGE>17)

table18up <- svytable(~YEAR+incdecile_w+did_wfh_labels,  design = table18up) # proportion of each respondants sex in sample


table18up <- table18up %>% 
  as_tibble() %>% 
  group_by(YEAR,incdecile_w)%>%
  mutate(Prop=round(n/sum(n), digits=3)) %>%
  filter(did_wfh_labels == "Did WFH")

table18up # has proportions calculated out of TOTAl for both years

table18up %>%
  ggplot(aes(factor(incdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, fill = YEAR, group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(stat="identity", position = "dodge")+
  #geom_col(stat = "identity", position = "stack") +   # scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+

 # facet_wrap(~YEAR)+
   coord_flip()+
    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_dodge(width = 0.8), hjust = 1.1,
              size = 4) + 
  labs(title ="Ages 18+: Percent of each income decile that worked from home",
  subtitle = "2019 vs 2021",
       caption = "ACS 1 year samples for 2019 and 2021. Working from home based on TRANWORK question on commuting.
       All workers in the labor force, all ages included.
       Income based on INCEARN for total earned income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of workers working from home") + 
 theme(legend.position = "bottom", legend.title = element_blank())+
  theme_classic()+
  scale_y_continuous(labels = scales::percent)



table25up<-subset(dstrata, AGE>24)
table25up <- svytable(~YEAR+incdecile_w+did_wfh_labels,  design = table25up) # proportion of each respondants sex in sample


table25up <- table25up %>% 
  as_tibble() %>% 
  group_by(YEAR,incdecile_w)%>%
  mutate(Prop=round(n/sum(n), digits=3)) %>%
  filter(did_wfh_labels == "Did WFH")

table25up # has proportions calculated out of TOTAl for both years

table25up %>%
  ggplot(aes(factor(incdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, fill = YEAR, group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(stat="identity", position = "dodge")+
  #geom_col(stat = "identity", position = "stack") +   # scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+

 # facet_wrap(~YEAR)+
   coord_flip()+
    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_dodge(width = 0.8), hjust = 1.1,
              size = 4) + 
  labs(title ="Ages 25+: Percent of each income decile that worked from home",
  subtitle = "2019 vs 2021",
       caption = "ACS 1 year samples for 2019 and 2021. Working from home based on TRANWORK question on commuting.
       All workers in the labor force, all ages included.
       Income based on INCEARN for total earned income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of workers working from home") + 
 theme(legend.position = "bottom", legend.title = element_blank())+
  theme_classic()+
  scale_y_continuous(labels = scales::percent)
```

```{r svy-tables}
#| include: FALSE

# joined %>% filter(YEAR == 2019)  %>% 
#   tbl_summary(include = c(-OCCSOC, -PUMA),
#               statistic = list(all_continuous() ~"{mean} ({sd})"))


# weights_2019table <- full_tidy_2019 %>% 
#   tbl_svysummary(include = c(INCEARN, INCWAGE, TRANWORK, EMPSTAT, LABFORCE, SEX, AGE, RACE, CLASSWKR))


# design2019 %>% 
#   tbl_svysummary(include = c(-OCCSOC, -PUMA))
# 
# design2021 %>% 
#   tbl_svysummary(include = c(-OCCSOC, -PUMA))
# 

# shows the breaks for the quantiles. hypothetically uses weights of observations for calculation of deciles.
# equal number of people should be in each decile after weights. 






# incomebywfh <- dstrata2021 %>% 
#   mutate(INCEARN = as.numeric(INCEARN) )%>%
#   svyby(formula = ~INCEARN, by = ~did_wfh_labels,# design = dstrata2021,
#       FUN = svyquantile, 
#       na.rm=TRUE,
#      keep.names = FALSE,
#         quantiles  = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9) )
# 
# incomebywfh
## ^ Income deciles by wfh status.
# The poorest 10% of WFH workers are earning almost more than the the bottom 20% of earners who did not WFH

```

> who did work from home in lower income brackets???

-   399011 in Service Occupations had 148 observations, 311122 had 86.

-   Vast Majority of those that did work form home in the bottom 10% of earners had management, business, sales, or office jobs. (specifically 399011, 253041,436014)

-   **outliers**: occupation 537062, 533030, 537065 in production and transportation, 399011 in service occupations which usually need to be in person.

-   399011 (2nd and 3rd decile), 311122 (2nd and 3rd decile)

```{r}
data %>% filter(did_wfh==1 & incdecile_w < 4) %>% group_by(incdecile_w, occ_2digits, occ_2dig_labels, OCCSOC) %>% summarize(obs = n()) %>% arrange(incdecile_w,-obs)
```

### Gender and Working from home

19.2% of those in the labor force worked from home in 2021.

-   10.7% were women, 9.8% were men.
-   In 2021, 17.6% of male workers did WFH and 21% of female workers did WFH.

52.5% of those in the labor force that have occupation codes were Men, 47.5% were Women.

```{r}

#round(prop.table(svytable(~did_wfh_labels, design=dstrata2019))*100,digits=2) 
#94.7 did not wfh, 5.29 did work from home in 2019

round(prop.table(svytable(~did_wfh_labels, design=dstrata2021))*100,digits=2) 

# # unweighted attempt at summary table using "data" dataframe
# table <- data %>% 
#  # filter(did_wfh==1) %>%  
#   mutate(total = n())  %>%
#   group_by(YEAR, did_wfh_labels, SEX) %>% 
#   dplyr::summarize(n_unweighted=n()) %>% 
#   mutate(Prop = n_unweighted/sum(n_unweighted))
# # unweighted sex proportions each year
# # for comparison
# 
# # of those that did work from home,51% were female in 2019 and 51.9% were female in 2021. 
# table 

table <- svytable(~SEX+YEAR+did_wfh_labels, design = dstrata) 
# proportion of each sex that did or did not work from home 
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, SEX)%>% 
  mutate(Prop =round(n/sum(n), digits=4))
# 5.4% of women worked at home in 2019 and 21% of women worked at home in 2021.
table


# # attempt 1, not what I wanted
# table %>% ggplot(aes(did_wfh_labels, y=Prop, group = YEAR, fill = YEAR)) + 
#   geom_col(stat = "identity", position = "dodge") +    
#   geom_text(aes(label = scales::percent(Prop)), position = position_dodge(width = .9), size=3,vjust=1.1) + 
#   theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
#   labs(title ="First Attempt at WFH graph for each Year",
#        x = "", y = "",
#       caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::percent)


# # attempt 2
# table %>% ggplot(aes(factor(SEX, labels = c("Male", "Female")), y=Prop, fill = did_wfh_labels, group = YEAR)) + 
#   geom_col(stat = "identity", position = "stack") +
#   facet_wrap(~YEAR)+
#     geom_text(aes(label = scales::percent(Prop)), position = position_fill(vjust=.5), size=3) + 
#   theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
#   scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
#   labs(title ="Percent working from home by Sex: 2019 vs 2021",
#        x = "", y = "",
#       caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::percent)


### Percentages add up to 100% for each gender for each year. 
table %>% ggplot(aes(factor(SEX, labels = c("Male", "Female")), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) +
  scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home by Sex: 2019 vs 2021",
       x = "", y = "# of People",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::comma)

#compare combined strata survey design to single year survey objects:
# MATCHES
# table <- svytable(~SEX+did_wfh, design = dstrata2019) # proportion of each respondants sex in sample
# table <- table %>% 
#   as_tibble() %>% 
#   group_by(SEX) %>%
#   mutate(Prop = round(n/sum(n), digits=4))
# table
# 
# table <- svytable(~SEX+did_wfh, design = dstrata2021, 
#                   round= TRUE) 
# table <- table %>% 
#   as_tibble() %>% group_by(SEX) %>%
#   mutate(Prop = round(n/sum(n), digits=4))
# table



## Totals add up to total number of workers in a year
table <- svytable(~SEX+YEAR+did_wfh_labels, design = dstrata) 
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% # grouped by Year and Sex!!
  mutate(Prop =round(n/sum(n), digits=4))
table

## percentages add up to 100 when adding all workers together for a year
table %>% ggplot(aes(factor(SEX, labels = c("Male", "Female")), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) +
  scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home by Sex: 2019 vs 2021",
       x = "", y = "# of People",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates.
      Percentages add up to 100 when adding all workers together each year.
      Working from home based on TRANWORK==80 variable from ACS surveys.") + scale_y_continuous(labels = scales::comma)

svyttest(formula = SEX~did_wfh_labels, design = dstrata2021)
svychisq(~did_wfh_labels+SEX, dstrata2021, statistic = "Chisq")

round(prop.table(svytable(~SEX+did_wfh_labels, design=dstrata2021))*100,digits=2)



```

```{r}

round(prop.table(svytable(~has_occsoc+SEX, design=dstrata2019))*100,digits=2)
# slight drop in women in the workforce (technically women with occupations). 
round(prop.table(svytable(~has_occsoc+SEX, design=dstrata2021))*100,digits=2)


CrossTable(data$SEX, data$YEAR)
xtabs(~SEX+YEAR, data = data) %>% summary()

xtabs(~did_wfh+SEX+YEAR, data = data) %>% summary()




## Most interesting change: 
round(prop.table(svytable(~SEX+did_wfh_labels, design=dstrata2019))*100,digits=2) 
# 1 is Male, 2 is female
round(prop.table(svytable(~SEX+did_wfh_labels, design=dstrata2021))*100,digits=2) 


# Chi square goodness of fit

# practice test to see if sex is 52% male and 48% female. For both years. 
chisq.test(table(joined$SEX), p = c(52, 48)/100)
# can't reject the null, so the sample represents the population if the known population is 52% male and 48% female.

svytotal(x = ~interaction(SEX, did_wfh_labels), design = dstrata2021, na.rm=TRUE)
svytotal(x = ~interaction(SEX, did_wfh_labels), design = dstrata2019, na.rm=TRUE)


#svymean(x = ~interaction(INCEARN, EMPSTAT), design = dstrata2021, na.rm=TRUE)
#svymean(x = ~interaction(INCEARN, LABFORCE), design = dstrata2021, na.rm=TRUE)


svyttest(formula = INCEARN~SEX, design = dstrata2021)
svyttest(formula = INCEARN~did_wfh_labels, design = dstrata2021)


table <- svytable(~SEX+YEAR+did_wfh_labels, design = dstrata) 
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, SEX)%>% # grouped by Year and Sex!!
  mutate(Prop =round(n/sum(n), digits=4))
table

```

Null hypothesis: There is no difference in working from home associated with sex.

Alt. Hyp: Working from home is associated with sex.

```{r}
svychisq(~SEX+did_wfh_labels, design = dstrata2021, statistic = "Chisq")
```

```{r}
incomebysex<- dstrata2021 %>% 
  mutate(INCEARN = as.numeric(INCEARN) )%>%
  svyby(formula = ~ INCEARN, by = ~SEX, 
      FUN = svyquantile, 
      na.rm=TRUE,
        quantiles  = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9) )

incomebysex %>% 
  pivot_longer(cols = c(-SEX),names_to = "columns", values_to = "values") %>% 
  mutate(values = round(values, digits =2))

```

```{r graphing-errorbars}
output <- svyby(formula = ~INCEARN, by = ~SEX, design = dstrata2021,
      FUN = svymean, na.rm=TRUE)
  
out_col <- mutate(output, 
                  lower = INCEARN - 2*se, 
                  upper = INCEARN + 2*se)

ggplot(out_col, aes(SEX, INCEARN, ymin=lower, ymax=upper)) + geom_col() +geom_errorbar(width = 0.7) + labs( y = "Average Earned Income", x = "Sex")+ scale_x_discrete(limits = c("Male", "Female") )


```

```{r}
table <- svytable(~SEX+NCHLT5+YEAR+did_wfh_labels, design = dstrata) 
# proportion of each sex that did or did not work from home 
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, did_wfh_labels)%>% 
  mutate(Prop =round(n/sum(n), digits=4))
# 5.4% of women worked at home in 2019 and 21% of women worked at home in 2021.
table
```

### Graphing Income Deciles

> percent includes both years for combined strata.

```{r graphing-inc-deciles-notweighted, include=FALSE, eval=FALSE}

#min= -8000, max=949000 for INCEARN 

#code from above

dstrata %>% group_by(YEAR) %>% summarize(count = n())


dstrata %>%
  as_data_frame() %>%
  dplyr::group_by(YEAR)%>% 
 # mutate(yrcount=)

  filter(did_wfh==1) %>%
  ggplot(aes(x=decile,
             y = (..count..)/(sum(..count..))*10,

           fill=factor(YEAR)
             )) + 
  geom_bar(position=position_dodge())+
  scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))  +
#    scale_y_discrete(limits = c("2019", "2021"))  +
  #scale_y_continuous(labels=percent) +
  labs(x="",y="",
       title = "Using dstrata dataframe: % of workers who worked at home by income decile in 2019 and 2021",
       subtitle = "PERCENT IS CALCULATED FROM TOTAL OBS FOR BOTH YEARS: FIX")+
  coord_flip()


# dstrata2019 %>%
#   as_data_frame() %>%
#   filter(did_wfh_labels != "NA") %>%
#   ggplot(aes(x=decile, 
#              y = (..count..)/(sum(..count..))*10,
#              fill=did_wfh_labels
#              )) + 
#   geom_bar(position="fill")+
#   theme(legend.position = "bottom", legend.title = element_blank()) +
#   scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))  +
#   labs(x="",y="",
#        title = "2019 Data: Who DID work from home?",
#        subtitle = "Percent of WFH workers by income decile")+
#   coord_flip()

dstrata2021 %>%
  as_data_frame() %>%
  filter(did_wfh_labels != "NA") %>%
  ggplot(aes(x=decile, 
             y = (..count..)/(sum(..count..))*10,
             fill=did_wfh_labels
             )) + 
  geom_bar(position="fill")+
  theme(legend.position = "bottom", legend.title = element_blank()) +
  scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))  +
  labs(x="",y="",
       title = "2021: Who DID work from home?",
       subtitle = "Percent of WFH workers by income decile")+
  coord_flip()


dstrata %>%
  as_data_frame() %>%
  filter(YEAR == 2019) %>%
  count(did_wfh, decile) %>%
  mutate(pct = round(n/sum(n)*1000, 1))%>%
  filter(did_wfh == 1)%>%
  #filter(did_wfh == "Can WFH") %>%
  #group_by(decile)%>%
  #summarize(count2 = n()) %>%
  ungroup()%>%
  ggplot(aes(x=decile,
            # group = YEAR,
             y = pct,
             )) + 
  geom_col() +
 # geom_col(aes(fill=YEAR)) +
  coord_flip()+
    theme(legend.position = "bottom") +

  scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
  geom_text(aes(label=paste0(pct,"%")), hjust=1.3)+
  labs(x="",y="% Worked at Home", title = "Individuals with higher income were more able to work from home",
       subtitle = "2019 ACS 1-year Sample with incorrectly weighted deciles",
       caption = "Did work from home is based on ACS variable TRANWORK.
       Delete this graph eventually. Keeping it for code comparison for now.")+
  theme_minimal()+
  theme(legend.position = "none")

# CanWFH2019 <- dstrata2019 %>%
# as_data_frame() %>%
#  # mutate(total = n()) %>%
#   filter(CanWorkFromHome != "Check Me") %>%
#   ggplot(aes(x=decile, 
#              y = (..count..)/sum(..count..)*10,
#              fill = CanWorkFromHome)) + 
#   geom_bar( position = "dodge") +
#   #coord_flip()+
#     theme(legend.position = "bottom", legend.title = element_blank()) +
# 
#   scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
#   scale_y_continuous(labels=scales::percent) +
#   labs(x="",y="", title = "2019 Data: Who Could work from home?", subtitle = "Based on Job Characteristics")
# 
# CanWFH2019
# 
# 
# 
# CanWFH2021 <- dstrata2021 %>%
# as_data_frame() %>%
#  # mutate(total = n()) %>%
#   filter(CanWorkFromHome != "Check Me") %>%
#   ggplot(aes(x=decile, 
#              y = (..count..)/sum(..count..)*10,
#              fill = CanWorkFromHome)) + 
#   geom_bar( position = "dodge") +
#   #coord_flip()+
#     theme(legend.position = "bottom", legend.title = element_blank()) +
# 
#   scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
#   scale_y_continuous(labels=scales::percent) +
#   labs(x="",y="", title = "2021 Data: Who Could work from home?", subtitle = "Based on Job Characteristics")
# 
# CanWFH2021

# dstrata2019 %>%
# as_data_frame() %>%
#  # mutate(total = n()) %>%
#   filter(CanWorkFromHome != "Check Me") %>%
#   ggplot(aes(x=decile, 
#              y = (..count..)/sum(..count..)*10,
#              fill = YEAR)) + 
#   geom_bar(aes(fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH"))),position = "stack") +
#   #coord_flip()+
#     theme(legend.position = "bottom", legend.title = element_blank()) +
# 
#   scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
#   scale_y_continuous(labels=scales::percent) +
#   labs(x="",y="", title = "2019 Data")+coord_flip()


# dstrata2021 %>%
# as_data_frame() %>%
#  # mutate(total = n()) %>%
#   filter(CanWorkFromHome != "Check Me") %>%
#   ggplot(aes(x=decile, 
#              y = (..count..)/sum(..count..)*10,
#              fill = YEAR)) + 
#   geom_bar(aes(fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH"))),position = "stack") +
#   #coord_flip()+
#     theme(legend.position = "bottom", legend.title = element_blank()) +
# 
#   scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
#   scale_y_continuous(labels=scales::percent) +
#   labs(x="",y="", title = "2021 Data")+
#   coord_flip()



```

```{r include=FALSE, eval=FALSE}


dstrata2019 %>%
as_data_frame() %>%
  filter(CanWorkFromHome != "Check Me") %>%
  count(decile, CanWorkFromHome) %>%
  mutate(pct = round(n/sum(n/10), 3)*100)%>%
  ggplot(aes(x=decile,  y = n,
         fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH"))
         )) + 
  geom_bar( stat = "identity") +
  coord_flip()+
  geom_text(aes(label=paste0(pct, "%")), 
            position = position_stack(vjust=0.5), size=2.5)+
    theme(legend.position = "bottom", legend.title = element_blank()) +

  scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
    guides(fill = guide_legend(reverse = TRUE))+

  #scale_y_continuous(labels=scales::percent) +
  labs(x="Earned Income Deciles", y="# Survey Responses", 
       title = "Ability to work from home based on occupational characteristics",
       subtitle = "HAS COUNTS on x axis: FIX THAT",
       caption = "Occupation data from ACS 1-year 2019 sample.")



dstrata2021 %>%
as_data_frame() %>%
  filter(CanWorkFromHome != "Check Me") %>%
  ggplot(aes(x=decile, 
        fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH"))
         )) + 
  geom_bar( position = "fill") +
  coord_flip()+
  geom_text( aes(label=paste0(signif(..count.. / tapply(..count.., ..x.., sum)[as.character(..x..)], digits=3)*100,"%")),
    stat="count", position = position_fill(vjust=0.5), size=3) +
    theme(legend.position = "bottom", legend.title = element_blank()) +

  scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
  guides(fill = guide_legend(reverse = TRUE))+
  #scale_y_continuous(labels=scales::percent) +
  labs(x="Earned Income Deciles", y="% Survey Responses", 
       title = "Ability to work from home based on occupational characteristics",
       caption = "Occupation data from ACS 1-year 2021 sample.") 


  
dstrata2021 %>%
as_data_frame() %>%
  filter(CanWorkFromHome != "Check Me") %>%
  ggplot(aes(x=decile, 
        fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH"))
         )) + 
  geom_bar(aes(fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH"))), position = "fill") +
  coord_flip()+
  geom_text( aes(label=paste0(signif(..count.. / tapply(..count.., ..x.., sum)[as.character(..x..)], digits=3)*100,"%")),
    stat="count", position = position_fill(vjust=0.5), size=3) +
    theme(legend.position = "bottom", legend.title = element_blank()) +

  scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
  guides(fill = guide_legend(reverse = TRUE))+
  #scale_y_continuous(labels=scales::percent) +
  labs(x="Earned Income Deciles", y="% Survey Responses", 
       title = "Ability to work from home based on occupational characteristics",
       caption = "Occupation data from ACS 1-year 2021 sample.")
 
```

```{r can-WFH-graphs}
#| code-fold: true


table <- svytable(~YEAR+incdecile_w+CanWorkFromHome,  design = dstrata) # proportion of each respondants sex in sample

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, incdecile_w)%>%
  mutate(Prop=round(n/sum(n), digits=3)) %>%
  mutate(CanWorkFromHome = factor(CanWorkFromHome, levels = c('No WFH', 'Some WFH',  'Can WFH')))

table 

table[rev(order(table$CanWorkFromHome)),]%>%
  ggplot(aes(factor(incdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), 
                    labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, 
            # fill = CanWorkFromHome,
            fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH")), 
             group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(aes(fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH")),stat="identity", position = "stack"))+
   coord_flip()+

    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_fill(vjust =.5), size = 2) + 
    guides(fill = guide_legend(reverse = TRUE))+

  labs(title ="Percent of each income decile that could potentially work from home",
       subtitle = "2019 vs 2021",
       caption = " Based on occupation codes from ACS 1 year samples for 2019 and 2021. Teleworkable coding based on Dingel & Neimen 2020.
       All workers in the labor force, all ages included.
       Income based on INCEARN for total earned income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of workers that can work from home based on occupation characteristics") +   theme_classic()+

 theme(legend.position = "bottom", legend.title = element_blank())+
  scale_y_continuous(labels = scales::percent) +   facet_wrap(~YEAR)

  
 



```

```{r agegroup-comparison, eval=FALSE, include=FALSE}
# Bottom 10% of earners change the most depending on the age groups incuded in sample. 

table <- svytable(~YEAR+incdecile_w+CanWorkFromHome,  design = dstrata) # proportion of each respondants sex in sample

table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, incdecile_w)%>%
  mutate(Prop=round(n/sum(n), digits=3)) %>%
  mutate(CanWorkFromHome = factor(CanWorkFromHome, levels = c('No WFH', 'Some WFH',  'Can WFH')))

table[rev(order(table$CanWorkFromHome)),]%>%
  ggplot(aes(factor(incdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), 
                    labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, 
            # fill = CanWorkFromHome,
            fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH")), 
             group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(aes(fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH")),stat="identity", position = "stack"))+
   coord_flip()+

    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_fill(vjust =.5), size = 2) + 
    guides(fill = guide_legend(reverse = TRUE))+

  labs(title ="All Ages (16+): Percent of each income decile that could potentially work from home",
       subtitle = "2019 vs 2021",
       caption = " Based on occupation codes from ACS 1 year samples for 2019 and 2021. Teleworkable coding based on Dingel & Neimen 2020.
       All workers in the labor force, all ages included.
       Income based on INCEARN for total earned income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of workers that can work from home based on occupation characteristics") +   theme_classic()+

 theme(legend.position = "bottom", legend.title = element_blank())+
  scale_y_continuous(labels = scales::percent) +   facet_wrap(~YEAR)




table18up<-subset(dstrata,AGE>17)

table18up <- svytable(~YEAR+incdecile_w+CanWorkFromHome,  design = table18up) # proportion of each respondants sex in sample


table18up <- table18up %>% 
  as_tibble() %>% 
  group_by(YEAR, incdecile_w)%>%
  mutate(Prop=round(n/sum(n), digits=3)) %>%
  mutate(CanWorkFromHome = factor(CanWorkFromHome, levels = c('No WFH', 'Some WFH',  'Can WFH')))


table18up[rev(order(table$CanWorkFromHome)),]%>%
  ggplot(aes(factor(incdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), 
                    labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, 
            # fill = CanWorkFromHome,
            fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH")), 
             group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(aes(fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH")),stat="identity", position = "stack"))+
   coord_flip()+

    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_fill(vjust =.5), size = 2) + 
    guides(fill = guide_legend(reverse = TRUE))+

  labs(title ="18&up: Percent of each income decile that could potentially work from home",
       subtitle = "2019 vs 2021",
       caption = " Based on occupation codes from ACS 1 year samples for 2019 and 2021. Teleworkable coding based on Dingel & Neimen 2020.
       All workers in the labor force, all ages included.
       Income based on INCEARN for total earned income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of workers that can work from home based on occupation characteristics") +   theme_classic()+

 theme(legend.position = "bottom", legend.title = element_blank())+
  scale_y_continuous(labels = scales::percent) +   facet_wrap(~YEAR)




table25up<-subset(dstrata,AGE>24)

table25up <- svytable(~YEAR+incdecile_w+CanWorkFromHome,  design = table25up) # proportion of each respondants sex in sample


table25up <- table25up %>% 
  as_tibble() %>% 
  group_by(YEAR, incdecile_w)%>%
  mutate(Prop=round(n/sum(n), digits=3)) %>%
  mutate(CanWorkFromHome = factor(CanWorkFromHome, levels = c('No WFH', 'Some WFH',  'Can WFH')))

table25up[rev(order(table$CanWorkFromHome)),]%>%
  ggplot(aes(factor(incdecile_w, levels = c(1,2,3,4,5,6,7,8,9,10), 
                    labels = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%")), 
                     y=Prop, 
            # fill = CanWorkFromHome,
            fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH")), 
             group = factor(YEAR, levels = "2021","2019"))) + 
  geom_col(aes(fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH")),stat="identity", position = "stack"))+
   coord_flip()+

    geom_text(aes(label = scales::percent(Prop, accuracy = 0.1L)), position = position_fill(vjust =.5), size = 2) + 
    guides(fill = guide_legend(reverse = TRUE))+

  labs(title ="25&up: Percent of each income decile that could potentially work from home",
       subtitle = "2019 vs 2021",
       caption = " Based on occupation codes from ACS 1 year samples for 2019 and 2021. Teleworkable coding based on Dingel & Neimen 2020.
       All workers in the labor force, all ages included.
       Income based on INCEARN for total earned income of survey respondents.", 
       x= "Income Deciles", 
       y = "Percent of workers that can work from home based on occupation characteristics") +   
  theme_classic()+
  theme(legend.position = "bottom", legend.title = element_blank())+
  scale_y_continuous(labels = scales::percent) +   
  facet_wrap(~YEAR)

  
```

```{r avginc-perdecile-graphs}
#| code-fold: TRUE

# summary(as.numeric(data$INCEARN))
# 
# # using decile made # of observations
# summary <- dstrata %>% 
#   group_by(YEAR,decile) %>% 
#   summarize(min = min(INCEARN),
#             max=max(INCEARN),
#             avg_income = mean(INCEARN),
#             average_income = survey_mean(INCEARN),
#             pop_represented = sum(PERWT),
#             obs_count = n())
# summary %>% 
#   ggplot(aes(x=decile, y=average_income, label=scales::dollar(average_income))) + 
#   geom_col()+
#   scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
#   scale_y_continuous(labels = scales::dollar)+labs(x="",y="", title = "Average earned income for each income decile")+
#   geom_text(vjust = -0.5, size = 3)+
#   facet_wrap(~YEAR)
# 
# 
# 
# # with weighted deciles made from # of pop represented
# summary <- dstrata %>% 
#   group_by(YEAR, incdecile_w) %>% 
#   summarize(min = min(INCEARN),
#             max=max(INCEARN),
#             avg_income = mean(INCEARN),
#             average_income = survey_mean(INCEARN),
#             pop_represented = sum(PERWT),
#             obs_count = n())
# 
# summary %>% 
#   ggplot(aes(x=incdecile_w, y=average_income, label=scales::dollar(average_income))) + 
#   geom_col()+
#   scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
#   scale_y_continuous(labels = scales::dollar)+labs(x="",y="", title = "Average earned income for each income decile")+
#   geom_text(vjust = -0.5, size = 3)+
#   facet_wrap(~YEAR)




output <- svyby(formula = ~INCEARN, by = ~YEAR+incdecile_w, design = dstrata, 
      FUN = svymean, na.rm=TRUE)
  
out_col <- mutate(output, 
                  lower = INCEARN - 2*se, 
                  upper = INCEARN + 2*se)

ggplot(out_col, aes(incdecile_w, INCEARN, ymin=lower, ymax=upper)) + 
  geom_col(stat = "identity") +
  geom_errorbar(width = 0.7) + 
  facet_wrap(~YEAR) + 
  labs(title = "Average income earned for each decile of income earners", y = "Average Earned Income for each Income Decile", x = "Deciles of Income Earners", caption = "NOT adjusted for inflation")+ coord_flip() +
  scale_x_discrete(limits = decile_labels ) + 
 # ylim(0,250000)+
  scale_y_continuous(label = scales::dollar, limits = c(0,275000))+
   geom_text(aes(label = scales::dollar(INCEARN)), size=3, hjust=-.3)


```

```{r other-teleworkable-graphs, include=FALSE, eval=FALSE}

dstrata %>% filter(YEAR==2019) %>%
as_data_frame() %>%
  filter(CanWorkFromHome != "Check Me" & did_wfh_labels != "NA") %>%

  ggplot(aes(x=decile, 
        fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH"))
         )) + 
  geom_bar( position = "fill") +
  coord_flip()+
  geom_text( aes(label=paste0(signif(..count.. / tapply(..count.., ..x.., sum)[as.character(..x..)], digits=3)*100,"%")),
    stat="count", position = position_fill(vjust=0.5), size=3) +
    theme(legend.position = "bottom", legend.title = element_blank()) +

  scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
  guides(fill = guide_legend(reverse = TRUE))+
  #scale_y_continuous(labels=scales::percent) +
  labs(x="Earned Income Deciles", y="% Survey Responses", 
       title = "Ability to work from home based on occupational characteristics",
       caption = "Uses original decile variable instead of incdecile_w. 
       `decile` is based on # of observations. 
       `incdecile_w` should be based on population estimates.
       Occupation data from ACS 1-year 2019 sample.") 

dstrata %>% filter(YEAR==2021) %>%
as_data_frame() %>%
  filter(CanWorkFromHome != "Check Me" & did_wfh_labels != "NA") %>%

  ggplot(aes(x=decile, 
        fill = factor(CanWorkFromHome, levels = c("No WFH", "Some WFH",  "Can WFH"))
         )) + 
  geom_bar( position = "fill") +
  coord_flip()+
  geom_text( aes(label=paste0(signif(..count.. / tapply(..count.., ..x.., sum)[as.character(..x..)], digits=3)*100,"%")),
    stat="count", position = position_fill(vjust=0.5), size=3) +
    theme(legend.position = "bottom", legend.title = element_blank()) +

  scale_x_discrete(limits = c("Bottom 10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "Top 10%"))+
  guides(fill = guide_legend(reverse = TRUE))+
  #scale_y_continuous(labels=scales::percent) +
  labs(x="Earned Income Deciles", y="% Survey Responses", 
       title = "Ability to work from home based on occupational characteristics",
       caption = "Uses original decile variable instead of incdecile_w. 
       `decile` is based on # of observations. `incdecile_w` should be based on population estimates. 
       Occupation data from ACS 1-year 2021 sample.") 
```

```{r}
# data has 2019 and 2021 observations TOGETHER
topline(data, did_wfh,  weight=PERWT)


# All members of the labor force could have said they either work from home (TRANWORK=80), go to work using some form of transportation, or didn't answer the question. 8.2% of the labor force did not answer the TRANWORK question and should not be included in calculations.

data %>% filter(YEAR == 2021) %>% crosstab(x=RACE, y=did_wfh_labels, weight = PERWT, pct_type = "row", unwt_n=TRUE, n=FALSE)  # matches Francis Total Row in Race output


data  %>% crosstab_3way(x=YEAR, y=did_wfh, z=race_cat, weight = PERWT)


```

```{r}
crosstab(data, x=LABFORCE, y=SEX, weight = PERWT, pct_type = "row", unwt_n=TRUE, n=FALSE) # matches Francis, 52% male, 48% female in laborforce

data %>% filter(YEAR ==2021) %>% crosstab(x=did_wfh_labels, y=SEX, weight = PERWT, pct_type = "row", unwt_n=TRUE, n=FALSE)

crosstab_3way(data, x=YEAR, y=did_wfh_labels, z=SEX, weight = PERWT, pct_type = "row", unwt_n=TRUE, n=FALSE)

crosstab(joined, x=LABFORCE, y=SEX, weight = PERWT, pct_type = "row", unwt_n=TRUE, n=FALSE) # matches Francis, 52% male, 48% female in laborforce

joined %>% filter(YEAR ==2021) %>% crosstab(x=did_wfh_labels, y=SEX, weight = PERWT, pct_type = "row", unwt_n=TRUE, n=FALSE)

crosstab_3way(joined, x=YEAR, y=did_wfh_labels, z=SEX, weight = PERWT, pct_type = "row", unwt_n=TRUE, n=FALSE)


```

```{r worker-descstats, eval=FALSE}
# Employment Status
table(data$EMPSTAT) # unweighted

data %>% 
  group_by(YEAR, as_factor(EMPSTAT)) %>%
  dplyr::summarize(weightedcount=sum(PERWT),
                   unweightedcount = n()) %>%  #weighted
  mutate(pct_weight = weightedcount/sum(weightedcount),
         pct_noweight = unweightedcount/sum(unweightedcount))

# Labor Force
data %>% 
  group_by(YEAR, as_factor(LABFORCE)) %>%
  dplyr::summarize(weightedcount=sum(PERWT),
                   unweightedcount = n()) %>%  #weighted
  mutate(pct_weight = weightedcount/sum(weightedcount),
         pct_noweight = unweightedcount/sum(unweightedcount))


# Class of Worker
data %>% 
  group_by(YEAR, as_factor(CLASSWKR)) %>%
  dplyr::summarize(weightedcount=sum(PERWT),
                   unweightedcount = n()) %>%  #weighted
  mutate(pct_weight = weightedcount/sum(weightedcount),
         pct_noweight = unweightedcount/sum(unweightedcount))



```

[***Unweighted***]{.underline} *-*

*EMPSTAT: 59,259 observations are employed, 4,194 unemployed observations, and 41,289 observations are not in the workforce (21,881 NAs)*

*LABFORCE: 63,453 are in labor force, 41,289 are not. (21,881 NAs)*

*CLASSWKR: Of these, 68,388 work for wages and 7183 people are self-employed. (51,052 NA)*

*`did_wfh`: 10,949 observations worked from home, 47,584 did not work from home. Based on TRANWORK variable: recoded as binary variable (either did wfh or did not wfh).*

[**Weighted**]{.underline} **-**

EMPSTAT: 6,102,522 people are employed (49%), 479,879 people are unemployed (3.8%), and 3,624,811 are not in the labor force (29%). There are 2,463,257 missing values; Same as LABFORCE.

LABFORCE: 6,582,401 people (52%) are in the labor force. 3,625,811 (28%) of people are not in the labor force. 2,463,257 (20%) of observations missing values.

-   employed and unemployed equal number of people in labor force - that's good

`did_wfh`: 19.2% did work from home and 80.8% did not work from home in Illinois (when not filtering for age or employment)

Location of primary workplace: 5.8 million people located in Illinois.

For the counties that can be identified in the data (populations \> 100,000 & \< 200,000. 1-Year ACS have minimum of 65,000 population), the census summary tables are close but not identical to the tables calculated with the ACS sample data. "In this way more densely populated areas, like Chicago and Cook County will contain many PUMAs within their boundaries, while multiple sparsely populated entire counties, e.g., Jackson, Perry, Franklin, and Williamson, will comprise one PUMA." - [IPUMS v other Geographies](https://iecam.illinois.edu/browse/about-public-use-microdata-areas-pumas#:~:text=What%20is%20a%20PUMA%3F,Community%20Survey%201%2Dyear%20estimates)

19 is Champaign, 31 is Cook, 37 is DeKalb, 43 is DuPage, 89 is Kane, 111 is McHenry, etc.

Economic Characteristics summary table: [link](https://data.census.gov/table?t=Employment&g=040XX00US17,17$0500000&tid=ACSDP1Y2021.DP03&moe=false&tp=false)

```{r message=FALSE, warning=FALSE}

joined %>% 
  filter(CanWorkFromHome == "Some WFH") %>% 
  distinct(OCCSOC)


```

### Race

Good example of [graphing survey data](https://rpubs.com/Onduma/surveydata). Make a summary table that has the WEIGHTED Freq and Prop of the variables of interest before passing it to graphing commands. Using svytable to make the weighted table.

```{r race}
#| code-fold: TRUE

table <- svytable(~race_cat+YEAR+did_wfh_labels, design = dstrata) 
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, race_cat)%>% # grouped by Year and Sex!!
  mutate(Prop =round(n/sum(n), digits=4)) %>%
  arrange(did_wfh_labels, -n)
table

table %>% ggplot(aes(race_cat, y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) +
 # scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home by race: 2019 vs 2021",
       x = "", y = "# of People",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") +    scale_fill_manual(values = c("#a6bddb", "#2b8cbe")) +
 scale_y_continuous(labels = scales::comma)



table <- svytable(~did_wfh_labels+YEAR+race_cat, design = dstrata) 
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, did_wfh_labels)%>% # grouped by Year and Sex!!
  mutate(Prop =round(n/sum(n), digits=4))
table
table %>% ggplot(aes(did_wfh_labels, y=n, fill = race_cat, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) +
 # scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home by race: 2019 vs 2021",
       x = "", y = "# of People",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::comma)

table <- svytable(~did_wfh_labels+YEAR+race_cat, design = dstrata) 
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% # grouped by Year and Sex!!
  mutate(Prop =round(n/sum(n), digits=4))

table %>% ggplot(aes(did_wfh_labels, y=n, fill = race_cat, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) +
 # scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home by race: 2019 vs 2021",
       x = "", y = "# of People",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::comma)

svyby(~did_wfh, by = ~RACE, design = dstrata2021, svymean, na.rm=TRUE)
svyby(~INCEARN, by = ~RACE, design = dstrata2021, FUN = svymean, na.rm=TRUE)


joined %>% crosstab_3way(x=YEAR, y=did_wfh, z=white, weight = PERWT)
joined %>% crosstab_3way(x=YEAR, y=did_wfh, z=RACE, weight = PERWT)

# joined %>% crosstab_3way(x=YEAR, y=did_wfh, z=white, weight = PERWT) %>% 
#   ggplot(aes(x=`1`, y=YEAR, fill = factor(white))) + 
#   geom_col(stat="identity", position = "fill") + labs(title = "Percent of White people that were able to work at home each year")
# 
# joined %>% crosstab_3way(x=YEAR, y=did_wfh, z=black, weight = PERWT) %>% 
#   ggplot(aes(x=`1`, y=YEAR, fill = factor(black))) + 
#   geom_col(stat="identity", position = "fill") + labs( "Percent of Black People that were able to work at home each year")


crosstab(joined, x=RACE, y=did_wfh, weight = PERWT, unwt_n=TRUE, n=FALSE)
crosstab(joined, x=RACE, y=did_wfh, weight = PERWT, pct_type = "column", unwt_n=TRUE, n=FALSE)
crosstab(joined, x=did_wfh, y=RACE, weight = PERWT, pct_type = "row", unwt_n=TRUE, n=FALSE)
crosstab(joined, x=did_wfh, y=RACE, weight = PERWT)


joined %>% filter(YEAR == 2021) %>%
  group_by(RACE) %>%
  summarize(Freq = n()) %>%
  mutate(Prop = Freq/sum(Freq))


svytable(~RACE, design=dstrata2021) %>% 
  as.data.frame() %>% # creates a frequency count by default?
  mutate(Prop =Freq/sum(Freq))


topline(dstrata2019, RACE, weight = PERWT)
topline(dstrata2021, RACE, weight = PERWT)


#svytotal(x = ~interaction(RACE, did_wfh_labels), design = dstrata2021, na.rm=TRUE)

# need to create race variables before the dstrata survey object.
#svytotal(x = ~interaction(white, did_wfh_labels), design = dstrata2021, na.rm=TRUE)

```

\

### Age

```{r eval = FALSE}
table(joined$AGE,joined$LABFORCE)

joined%>% filter(YEAR == 2021 & AGE>64 & INCEARN > 0) %>% summary()
# more then 75% of those under 18 make less than $5,000 a year.Not the kind of workers we are interested in anyways.  

# Even more extreme for less than 17

table(joined$age_cat, joined$did_wfh_labels)

```

### Internet Access

Other variables: CIHISPEED, CINETHH, MULTGEN, NCHILD, NCHLT5, MARST, FERTYR

-   CINETHH reports whether any member of the household accesses the Internet. Here, "access" refers to whether or not someone in the household uses or connects to the Internet, regardless of whether or not they pay for the service.

-   CIHISPEED reports whether the respondent or any member of their household subscribed to the Internet using broadband (high speed) Internet service such as cable, fiber optic, or DSL service.

```{r}
#10 is yes, 20 is no access. 00 is NA
table <- svytable(~CIHISPEED +did_wfh_labels+YEAR, design = dstrata) 
topline(dstrata2019, CIHISPEED, weight = PERWT)
topline(dstrata2021, CIHISPEED, weight = PERWT)

table
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR, CIHISPEED)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>% 
  arrange(did_wfh_labels, -n)
table


table %>% ggplot(aes(factor(CIHISPEED, labels = c("Has high speed", "No high speed")), y=Prop, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_fill(vjust=.5), size=3) + 
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home & Access to high speed Internet: 2019 vs 2021",
       x = "", y = "",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::percent) +scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))


table %>% ggplot(aes(factor(CIHISPEED, labels = c("Has high speed", "No high speed")), y=n, fill = did_wfh_labels, group = YEAR)) + 
  geom_col(stat = "identity", position = "stack") +
  facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) + 
  theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home & Access to high speed Internet: 2019 vs 2021",
       x = "", y = "",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") +scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))

### Access to any type of internet ####
table <- #dstrata %>% factor(CIHISPEED,)
  svytable(~CINETHH+did_wfh_labels+YEAR, design = dstrata) 
table
# proportion of each respondant's sex and if they worked from home for each year in sample
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  mutate(Prop =round(n/sum(n), digits=3))%>% 
  arrange(did_wfh_labels, -n)
table


table %>% #filter(YEAR == 2021) %>%
  ggplot(aes(factor(CINETHH, labels = c("Has Own Access", "Has Other Access", "No Access")), y=n, fill = did_wfh_labels,
                     group = YEAR
                     )) + 
  geom_col(stat = "identity", position = "stack") +
 facet_wrap(~YEAR)+
  geom_text(aes(label = scales::percent(as.numeric(ifelse(Prop>0.05,Prop, "")), accuracy = .1),accuracy  = .1L ),position = position_stack(vjust=.5), size=3) +
  theme_classic() + 
  theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home & Access to Internet: 2019 vs 2021",
       x = "", y = "# of people",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::comma)+coord_flip() + scale_fill_manual(values = c("#a6bddb", "#2b8cbe"))+
   annotate("text", label = "< 3% of labor force", x = 2, y =1400000, size = 3, colour = "black")  + annotate("text", label = "< 5% of labor force", x = 3, y =1400000, size = 3, colour = "black")






table %>% #filter(YEAR == 2021) %>%
  ggplot(aes(did_wfh_labels, y=n, fill = factor(CINETHH, labels = c("Has Own Access", "Has Other Access", "No Access")),
                     group = YEAR
                     )) + 
  geom_col(stat = "identity", position = "stack") +
 facet_wrap(~YEAR)+
    geom_text(aes(label = scales::percent(Prop)), position = position_stack(vjust=.5), size=3) + 
  theme_classic() + 
  theme(legend.position = "bottom", legend.title = element_blank())+
  labs(title ="Percent working from home & Access to Internet: 2019 vs 2021",
       x = "", y = "# of People",
      caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::comma)+coord_flip()



# table %>% ggplot(aes(factor(CINETHH, labels = c("Has Own Access", "Has Other Access", "No Access")), y=Prop, fill = did_wfh_labels, group = YEAR)) + 
#   geom_col(stat = "identity", position = "fill") +
#   facet_wrap(~YEAR)+
#     geom_text(aes(label = scales::percent(Prop)), position = position_fill(vjust=.5), size=3) + 
#   theme_classic() + theme(legend.position = "bottom", legend.title = element_blank())+
#   labs(title ="Percent working from home & Access to Internet: 2019 vs 2021",
#        x = "", y = "",
#       caption = "ACS 1 year samples for 2019 and 2021 used for weighted population estimates,") + scale_y_continuous(labels = scales::percent)


# Access at home AND access to high speed internet.
table <- 
  svytable(~CIHISPEED+ CINETHH+did_wfh_labels+YEAR, design = dstrata) 
table
table <- table %>% 
  as_tibble() %>% 
  group_by(YEAR)%>% 
  mutate(Prop =round(n/sum(n), digits=3)) %>% 
  arrange(YEAR, did_wfh_labels)
table


crosstab_3way(joined, YEAR, CIHISPEED, CINETHH, weight = PERWT, unwt_n = TRUE)
```

In 2019, an estimated 277,000 people worked from home and had highspeed internet and 29,500 did not have high speed internet .

In 2021, an estimated1,027,000 people worked from home and had high speed internet and 83,000 did not have high speed internet.

-   There is a chance that people realized that their internet was not as fast as they thought it was when they had multiple people using the internet at the same time during COVID-19.

```{r}
table(joined$NCHILD)
table(joined$NCHLT5)


```

## Geographic Data

```{r}
table(data$PWSTATE2)

svytable(~did_wfh_labels+YEAR+work_in_IL, design = dstrata) 
318174/ (318174+5632943) # percent of remote workers that worked outside of IL but lived in IL in 2019
1140835/(1140835+4604673) # 2021 percent of WFH workers that worked outside of IL


svytable(~has_incearn_labels+YEAR+work_in_IL, design = dstrata) 
svytable(~work_in_IL+YEAR+has_incearn_labels, design = dstrata) 

490937  / (490937  +5951117 ) # percent of remote workers that worked outside of IL but lived in IL in 2019
525507/(525507+5745508) # 2021 percent of WFH workers that worked outside of IL

```

7.6% of workers that earned income worked outside of IL in 2019.

8.4% of workers did so in 2021.

```{r}
joined %>% 
  group_by(YEAR, county_pop_type )%>% 
  summarize(avg_teleworkable = mean(teleworkable, weight = PERWT),
                                                         
            avg_did_WFH = mean(did_wfh, na.rm=TRUE, weight = PERWT))
```

## PUMA Level

PUMAs contain around 100,000 people.

For the counties that can be identified in the data (populations \> 100,000 & \< 200,000. 1-Year ACS have minimum of 65,000 population), the census summary tables are close but not identical to the tables calculated with the ACS sample data. "In this way more densely populated areas, like Chicago and Cook County will contain many PUMAs within their boundaries, while multiple sparsely populated entire counties, e.g., Jackson, Perry, Franklin, and Williamson, will comprise one PUMA." - [IPUMS v other Geographies](https://iecam.illinois.edu/browse/about-public-use-microdata-areas-pumas#:~:text=What%20is%20a%20PUMA%3F,Community%20Survey%201%2Dyear%20estimates)\

```{r}
obs_perPUMA<- joined %>% 
   group_by(PUMA, YEAR) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>%
  arrange(unweightedcount)
### Minimimum observations are 271 in PUMA area 03528
obs_perPUMA
```


```{r}
joined %>% 
   group_by(county_pop_type, YEAR, did_wfh_labels) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount)) %>% arrange(unweightedcount)


```

Economic Characteristics summary table: [link](https://data.census.gov/table?t=Employment&g=040XX00US17,17$0500000&tid=ACSDP1Y2021.DP03&moe=false&tp=false)

```{r warning = FALSE, message=FALSE, results='hide'}

# PUMA shapefiles
pumasIL <- pumas("IL", cb=T, year=2019)

#pumasIL %>% select(PUMACE10, GEOID10, NAME10) %>% write_csv("PUMAsIL.csv")
#county shapefiles
countyIL <- counties("IL", cb=T, year=2019)

#pumasdf <- fortify(pumasIL, region = 'PUMACE10')
```

#### Did WFH
```{r}
#| column: page

mapPUMAboth <- svytable(~YEAR+PUMA+did_wfh_labels, design = dstrata)
mapPUMAboth <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(did_wfh_labels == "Did WFH") %>%
     full_join(pumasIL, by = c("PUMA" = "PUMACE10"))


figure5a <- mapPUMAboth %>% 
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Where Did People Work from Home?", 
      # subtitle = "Percent of PUMA population in 2019 & 2021 that did WFH"
      ) +
  theme_classic() + 
  theme(axis.ticks = element_blank(), axis.text = element_blank())+
   scale_fill_steps2(
                         high = "darkblue", low = "black", 

                       limits = c(0,.6),
                       n.breaks =4,
                       show.limits=TRUE,
                       nice.breaks=FALSE,

                        name = "Percent of Workers", 
                       label = scales::percent)+
  facet_wrap(~YEAR) 
 
mapPUMAboth %>% select(YEAR, PUMA, Prop, n, NAME10)



```


```{r include = FALSE}

#ggsave("Figure5a.eps", limitsize = FALSE,width = 8, height = 5, units = "in")
#ggsave("Figure5a.pdf", limitsize = FALSE,width = 8, height = 5, units = "in")
#ggsave("Figure5aV2.png", limitsize = FALSE, width = 8, height = 5, units = "in")

```

```{r}
#| column: page
#| layout-ncol: 2
mapPUMAboth <- svytable(~YEAR+PUMA+did_wfh_labels, design = dstrata)

mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(did_wfh_labels == "Did WFH") %>%
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
    mutate(pct_change= `2021`-`2019`) %>%
       full_join(pumasIL, by = c("PUMA" = "PUMACE10"))%>%
  arrange(-pct_change) %>% select(-c(STATEFP10:GEOID10)) %>% head()

mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(did_wfh_labels == "Did WFH") %>%
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
    mutate(pct_change= `2021`-`2019`) %>%
       full_join(pumasIL, by = c("PUMA" = "PUMACE10"))%>%
  arrange(pct_change) %>% select(-c(STATEFP10:GEOID10)) %>% head()
```


```{r}
#| code-fold: true




figure5b <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(did_wfh_labels == "Did WFH") %>%
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
    mutate(pct_change= `2021`-`2019`) %>%
       full_join(pumasIL, by = c("PUMA" = "PUMACE10")) %>%
  ggplot(aes(fill = pct_change)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Change in Population that Did Work from Home", 
       subtitle = "2019 vs 2021", 
       caption = 
         "The did_WFH variable based on TRANWORK==80 from
         ACS 2019 & 2021 1-Year Survey individual level data. 
       Each geographic region has between 100K and 200K people. ") +
  theme_classic() + 
  theme(axis.ticks = element_blank(), axis.text = element_blank())+
  
     scale_fill_steps2(#colors = c("white", "darkblue"), 
                         high = "darkblue", low = "black", 

                       limits = c(0,.6),
                       n.breaks =4,
                       show.limits=TRUE,
                       nice.breaks=FALSE,

                        name = "% change in WFH", 
                       label = scales::percent)


closeupFIPs<- c("097", "111", "089","043", "197", "031"#,"093"
             #   ,"007"
                )



mapPUMAboth <- svytable(~YEAR+countyFIP+PUMA+
                          did_wfh_labels, design = dstrata)

mapPUMAboth <- mapPUMAboth %>% 
  as_tibble() %>%
  filter(countyFIP %in% closeupFIPs) %>%

  group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(did_wfh_labels == "Did WFH" & n != "0") %>%
     inner_join(pumasIL, by = c("PUMA" = "PUMACE10"))
 
#mapPUMAboth %>% select(YEAR, countyFIP, PUMA, Prop, n, NAME10)




countyborders<- countyIL %>% 
  filter(COUNTYFP %in% closeupFIPs)
#countyborders


figure5c <- ggplot(mapPUMAboth, aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
 # labs(title = "Where Did People Work from Home?: Close up View", 
  #     subtitle = "Percent of PUMA population in 2019 & 2021 that did WFH") +
  theme_classic() + 
  theme(axis.ticks = element_blank(), axis.text = element_blank())+

     scale_fill_binned(#colors = c("white", "darkblue"), 
                         high = "darkblue", low = "white", 

                       limits = c(0,.6),
                       n.breaks =4,
                       show.limits=TRUE,
                       nice.breaks=FALSE,

                        name = "% that did WFH",
                       label = scales::percent)+
  facet_wrap(~YEAR) 

figure5c <- figure5c +
  geom_sf(data = countyborders, fill=NA, color="black", lwd = 1) + theme_classic()+
  theme(axis.ticks = element_blank(), axis.text = element_blank())+ facet_wrap(~YEAR)



# as a dot graph ## 
mapPUMAboth <- svytable(~YEAR+PUMA+did_wfh_labels, design = dstrata)

mapPUMAboth <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(did_wfh_labels == "Did WFH") %>%
     full_join(pumasIL, by = c("PUMA" = "PUMACE10"))

order <- mapPUMAboth %>% as_tibble() %>%
 filter(YEAR == 2021  & did_wfh_labels == "Did WFH") %>% 
    filter(Prop > 0.25 | Prop < 0.1) %>%
  select(NAME10, Prop_2021 = Prop)

mapPUMAboth <- inner_join(mapPUMAboth, order) %>% filter()

mapPUMAboth %>% 
  ungroup() %>%
  filter(#did_wfh_labels == "Did WFH" #& NAME10 != "NA" & YEAR != "NA"
         )%>% 
 # filter(YEAR == 2021 & (Prop > 0.25 | Prop < 0.1) ) %>%
  select(NAME10, PUMA, Prop) %>% 
  distinct() %>% arrange(desc(Prop)) %>% head()

figure5d <- mapPUMAboth %>%
  filter(did_wfh_labels == "Did WFH" #& NAME10 != "NA" & YEAR != "NA"
         )%>% 
  filter(Prop > 0.25 | Prop < 0.1) %>%
  ggplot(aes(x = Prop*100, y= reorder(PUMA, Prop_2021*100)))+
  geom_line(aes(group = PUMA))+
  theme_classic()+
   geom_point(aes(color=YEAR), size = 3) +
      scale_color_brewer(palette="Paired", labels = c("% WFH in 2019", "% WFH in 2021"), direction = 1)+
  labs(title = "Change in PUMA Labor Force that did WFH",
       subtitle = "Percent of Labor Force within PUMA that Did Work from Home", 
#  caption = "Graph includes PUMAs with Proportions of WFH Workers greater 
#  than 25% or less than 10%. Filtered PUMAs to increase legibility. 
#  DuPage had around 8% of the labor force working from home in 2019 and 
#  around 27% of its labor force working from home in 2021. 
#       There was ~20 percentage point increase between 2019 and 2021 for DuPage County. 
#  Working from home is based on TRANWORK variable from ACS 1-year surveys.", 
x = "Percentage Points", y = "" 
)


mapPUMAboth <- svytable(~YEAR+PUMA+countyFIP+did_wfh_labels, design = dstrata)

figure5e <- mapPUMAboth %>% as_tibble() %>%
    filter(countyFIP %in% closeupFIPs) %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(did_wfh_labels == "Did WFH" & n != "0") %>%
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
    mutate(pct_change= `2021`-`2019`) %>%
       inner_join(pumasIL, by = c("PUMA" = "PUMACE10")) %>%

  ggplot(aes(fill = pct_change)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Change in Population that Did Work from Home", 
       subtitle = "2019 vs 2021") +
  theme_classic() + 
  theme(axis.ticks = element_blank(), axis.text = element_blank()) +
    scale_fill_steps2( high = "darkblue", low = "black", 
                       limits = c(0,.6),
                       n.breaks =4,
                       show.limits=TRUE,
                       nice.breaks=FALSE,

                        name =" Pct Pt change  that did WFH",
                       label = scales::percent)
```

```{r}
#| fig-cap: "Percent of Workers who Did Work From Home: 2019 & 2021"
#| fig-cap-location: "top"
#| fig-subcap: 
#| - "The did_WFH variable is based on TRANWORK==80 from ACS 2019 & 2021 1-Year Survey individual level data. Each geographic region has between 100K and 200K people."
#| - "Zoom in on Cook and Surrounding Counties"
#| column: page
#| layout-ncol: 2
#| echo: false 

plot(figure5a)
plot(figure5c)



```


```{r}
#| fig-cap-location: "top"
#| column: page
#| echo: false 
#| layout-ncol: 2
#| fig-cap: "Change from 2019 to 2021 in who did WFH"
#| fig-subcap: 
#| - "Graph includes PUMAs with Proportions of WFH Workers greater than 25% or less than 10%. Filtered PUMAs to increase legibility. DuPage had around 8% of the labor force working from home in 2019 and around 27% of its labor force working from home in 2021. There was ~20 percentage point increase between 2019 and 2021 for DuPage County. Working from home is based on TRANWORK variable from ACS 1-year surveys."
#| - "A dot plot of change"


plot(figure5b)
plot(figure5d)
```





```{r include=FALSE}
#ggsave("Figure5b.eps", limitsize = FALSE,width = 8, height = 5, units = "in")
#ggsave("Figure5b.pdf", limitsize = FALSE,width = 8, height = 5, units = "in")
#ggsave("Figure5bV2.png", limitsize = FALSE, width = 8, height = 5, units = "in")
```


```{r puma-dotgraph, include=FALSE}
# as a dot graph ## 
mapPUMAboth <- svytable(~YEAR+PUMA+countyFIP+did_wfh_labels, design = dstrata)

mapPUMAboth <- mapPUMAboth %>% as_tibble() %>% 
  filter(countyFIP %in% closeupFIPs)%>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(did_wfh_labels == "Did WFH" & n != "0") %>%
     inner_join(pumasIL, by = c("PUMA" = "PUMACE10")) 

order <- mapPUMAboth %>% as_tibble() %>%
 filter(YEAR == 2021  & did_wfh_labels == "Did WFH") %>% 
    filter(Prop > 0.25 | Prop < 0.1) %>%
  select(NAME10, Prop_2021 = Prop)

mapPUMAboth <- inner_join(mapPUMAboth, order) %>% filter()


mapPUMAboth %>%
  filter(did_wfh_labels == "Did WFH" #& NAME10 != "NA" & YEAR != "NA"
         )%>% 
  filter(Prop > 0.25 | Prop < 0.1) %>%
  ggplot(aes(x = Prop*100, y= reorder(PUMA, Prop_2021*100)))+
  geom_line(aes(group = PUMA))+
  theme_classic()+
   geom_point(aes(color=YEAR), size = 3) +
      scale_color_brewer(palette="Paired", labels = c("% WFH in 2019", "% WFH in 2021"), direction = 1)+
  labs(title = "Change in PUMA Labor Force that did WFH",
       subtitle = "Percent of County Population that Did Work from Home", 
  caption = "", x = "Percentage Points", y = "" )
```

### WFH Feasibility

We compute these shares using our O\*NET-derived classification of occupations that can be done at home and the occupational composition using ACS responses in each PUMA by 6-digit SOC in the BLS's 2018 Occupational Employment Statistics.

```{r pct-canWFH}
#| column: page
#| layout-ncol: 2

mapPUMAboth <- svytable(~YEAR+PUMA+CanWorkFromHome, design = dstrata)
mapPUMAboth <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     filter(CanWorkFromHome == "Can WFH") %>%
     full_join(pumasIL, by = c("PUMA" = "PUMACE10"))
 
mapPUMAboth %>%
ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry=geometry), color = "black")+ 
  labs(title = "Percent of PUMA area population that COULD work from home", 
       subtitle = " based on occupation characteristics: 2019 vs 2021", 
       caption = "Teleworkability based on D&N 2020 methodology. 
       OCCSOC codes from ACS 1 year Survey data at individual level.") +
  theme_classic() + 
  theme(axis.ticks = element_blank(), axis.text = element_blank())+
  scale_fill_steps2(#colors = c("white", "darkblue"), 
                         high = "darkblue", low = "black", 

                       limits = c(0,.6),
                       n.breaks =4,
                       show.limits=TRUE,
                       nice.breaks=FALSE,
                        name = "% Potentially Could WFH", 
                       label = scales::percent) +
  facet_wrap(~YEAR)



Figure2b <- mapPUMAboth %>%
filter(YEAR == 2021)%>%  
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  # labs(title = "Work from Home Feasibility Rate", 
  #      subtitle = " Based on Occupation Characteristics", 
  #      caption = "Teleworkability based on D&N 2020 methodology. 
  #      OCCSOC codes from ACS 1 year Survey data at individual level.") +
  theme_void() + 
  theme(axis.ticks = element_blank(), axis.text = element_blank(),   
        plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust=.5),  plot.subtitle = element_text(hjust = 0.5)) +
                           scale_fill_steps2(#colors = c("white", "darkblue"), 
                         high = "darkblue", low = "black", 

                       limits = c(0,.6),
                       n.breaks =4,
                       show.limits=TRUE,
                       nice.breaks=FALSE,
                        name = "Labor force with\nWFH Feasibility", label = scales::percent)
Figure2b

#ggsave(Figure2b, "Figure2bV2.png", width=3, height=3, units = "in")
```

```{r separate-maps-didWFH, include=FALSE, eval=FALSE}

## Graph of both  together is better because it is automatically sharing a legend scale
#### filtering from joined dataframe for comparison: ###
  ## Prepping data for the graphs:
PUMAmapdata2019 <- joined %>% 
  as_data_frame() %>%
  filter(YEAR == 2019 & did_wfh_labels != "NA" )%>% # create percentages withing missing values included. aka valid percent.
   group_by(PUMA, did_wfh_labels) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(pumasIL, by = c("PUMA" = "PUMACE10"))  


PUMAmapdata2021 <- joined %>% 
  filter(YEAR == 2021 & did_wfh_labels != "NA" )%>% # create percentages withing missing values included. aka valid percent.
   group_by(PUMA, did_wfh_labels) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(pumasIL,by = c("PUMA" = "PUMACE10")) 



## Graphs
PUMAmapdata2021 %>%
  filter(did_wfh_labels != "Did not WFH") %>%
  ggplot(aes(fill = pct_weight)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Percent of PUMA population that did WFH in 2021",
       caption = "Note: All counties with populations < 50,000 people were 
       aggregated together in ACS data extracts (COUNTYFIP == 000)") 

PUMAmapdata2019 %>%
  filter(did_wfh_labels != "Did not WFH") %>%
  ggplot(aes(fill = pct_weight)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Percent of PUMA population that did WFH in 2019",
       subtitle = "Compare to facet_wrap graph",
       caption = "Note: All areas with populations < 65,000 people are
       aggregated together to create PUMA boundaries")
```

### Internet Access

Make a subset of dstrata that only includes people who answered the TRANWORK/ did_WFH questions.

```{r include = FALSE}
#transub<- subset(dstrata,TRANWORK(x=80))
```

```{r internet-access}
#| code-fold: TRUE
#| column: page
#| layout-ncol: 2
#| layout-nrow: 2

## CINETHH is for access to internet 
# (either at home, somewhere else, or no access)
# from joined dataframe. Not using survey objects or survey package
mapPUMAboth <- joined %>% 
   group_by(YEAR, PUMA, CINETHH) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(pumasIL, by = c("PUMA" = "PUMACE10")) 
mapPUMAboth

## identical to graph below made using svytable()
# mapPUMAboth %>%
#   filter(CINETHH ==3) %>% 
#   ggplot(aes(fill = pct_weight)) +
#   geom_sf(aes(geometry = geometry), color = "black")+ 
#     scale_fill_continuous(low = "gray", high = "red", name = "% of Population", label = scales::percent)+
#     theme_minimal()+
#     theme(legend.title = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="right")+
#   labs(title = "Percent of population that did NOT have access to Internet", 
#        subtitle = "Geography = PUMAs", 
#        ) +   theme(legend.title = element_blank())+
#  scale_fill_gradientn(colors = c("white", "red4"), 
#                     #   limits = c(-.1,.1),
#                        n = 7,
#                         name = "Change in Access", label = scales::percent)+
#   facet_wrap(~YEAR)



## Access internet using survey objects
mapPUMAboth <- svytable(~YEAR+PUMA+CINETHH, design = HHdesign)

mapPUMAboth <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%   
     full_join(pumasIL, by = c("PUMA" = "PUMACE10"))

mapPUMAboth %>%
  filter(CINETHH == 3 ) %>% 
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+   
  scale_fill_binned(low = "white", high = "red4", name = "% of Population", label = scales::percent,
                    n.breaks=4,nice.breaks=FALSE,show.limits=TRUE)+
  labs(title = "Percent of population that lacked access to Internet", 
       caption = "Geography = PUMAs. Includes all dstrata responses for CINETHH.
       Does not subset for people who answered the commuting or occsoc questions.") + 
  facet_wrap(~YEAR) +  theme_minimal()+
    theme(legend.title = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="right")


# mapPUMAboth %>%   
#   select(YEAR,Prop,CINETHH,geometry) %>%
#   filter(CINETHH ==3) %>% 
#   pivot_wider( names_from ="YEAR", values_from = "Prop") %>%
#   mutate(pct_change = `2021`-`2019`) %>%
#   ggplot(aes(fill = pct_change)) +
#   geom_sf(aes(geometry = geometry), color = "black")+ 
#   theme_minimal()+
#     theme(legend.title = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="right")+
#                            scale_fill_gradientn(colors = c("green4", "white", "red4"), 
#                        limits = c(-.15,.15),
#                        n = 7,
#                         name = "Change in Access", label = scales::percent)+
#   labs(title = "Decrease in Lack of Internet Access", 
#        subtitle = "Needs to be flipped so not double negative",
#        caption = "Rural areas were more likely to increase their access to internet at home. 
#        Change in Percentage points = 2021 % - 2019 %.
#        Access to internet at home or another location is considered having Access to internet")

mapPUMAboth %>%   
  filter(CINETHH ==1) %>% 
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
  mutate(pct_change = `2021`-`2019`) %>%
  ggplot(aes(fill = pct_change)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  theme_minimal()+
  theme(legend.title = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="right")+
  scale_fill_steps2(low="red4", mid= "white", high= "forestgreen", 
                 #   limits = c(-.15,.15),
                    n.breaks=4,
                    show.limits=TRUE,
                    nice.breaks=FALSE,
                    name = "Chance in Access", label = scales::percent)+
  labs(title = "Change in Labor Force with Internet at home",
       subtitle = "2019 to 2021") 


```

```{r include=FALSE}
joined %>% filter(CINETHH %in% c(1,2,3)) %>% 

   group_by(county_pop_type, YEAR, PUMA) %>% 
   dplyr::summarize(weightedcount=sum(HHWT), #weighted 
                    unweightedcount = n())%>%
       mutate( pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount)) %>%
              full_join(pumasIL, by = c("PUMA" = "PUMACE10")) %>%
  arrange(unweightedcount)

joined %>% 
   group_by(county_pop_type, YEAR, PUMA, CINETHH) %>% 
   dplyr::summarize(weightedcount=sum(HHWT), #weighted 
                    unweightedcount = n())%>%
       mutate( pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount)) %>%
              full_join(pumasIL, by = c("PUMA" = "PUMACE10")) %>% 
 #  filter(unweightedcount > 30) %>% 
  arrange(desc(CINETHH), unweightedcount)
 # ungroup() %>%
 # arrange(unweightedcount) %>% 
 # group_by(county_pop_type, YEAR, PUMA, CINETHH, weightedcount, unweightedcount) %>% 
 #  dplyr::summarize(#weightedcount=sum(PERWT), #weighted 
                    #unweightedcount = n(),
   #                 pct_weight = weightedcount/sum(weightedcount), 
   #       pct_noweight = unweightedcount/sum(unweightedcount))%>%



## Keeps only PUMAs with more than 30 observations. 
mapPUMAboth <- joined %>% 
  group_by(county_pop_type, YEAR, PUMA, CINETHH) %>% 
   dplyr::summarize(weightedcount=sum(HHWT), #weighted 
                    unweightedcount = n())%>%
  mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))  %>%
  full_join(pumasIL, by = c("PUMA" = "PUMACE10"))




```

##### High Speed Internet

High speed internet is not the best measurement of due to them including DSL in their definition of "high speed."

FCC standards say minimum 25Mbps for download and 3Mbps for upload. "traditionally, the way to determine if a connection is high speed is to test its ability to connect multiple devices that allow streaming and access to modern applications at the same time.

High speed internet allows teleworkers the ability to live and work in locations of their choosing without commuting to work. Also important for education, especially during COVID and online classes. Many families did not have the bandwidth necessary to support someone working from home and having children in online classes at the same time (let alone multiple children in online classes).

The Infrastructure Investment and Jobs Act (IIJA) defines "underserved" broadband as an Internet speed \< 100 Mbps downstream and 20Mbps upstream.

Zoom needs between 2 and 3 Mbps for group video calls that look good.

1.0 Mbps to 1.5 Mbps for lower quality video calls that should still work.

[Source](https://broadbandusa.ntia.doc.gov/sites/default/files/2022-12/Introduction_to_Broadband_and_High_Speed_Internet_FINAL_0.pdf)

> ::: {.callout-note appearance="minimal"}
> For areas that had a decrease in access to high speed internet: Is it possible that they realized that their internet sucked when multiple people tried to use it? Most people didn't know internet speeds before the pandemic. They thought they had high speed internet precovid and then when multiple people were at home, they realized they didn't have enough bandwith for their internet demand?
> :::

```{r fast-internet-byPUMA}
#| code-fold: TRUE
#| 
### Access to Hi speed internet ###

# from joined dataframe. Not using survey objects or survey package
mapPUMAboth <- joined %>% 
   group_by(county_pop_type, YEAR, PUMA, CIHISPEED) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(pumasIL, by = c("PUMA" = "PUMACE10")) 
mapPUMAboth
```


```{r fig-hispeedinternet}
#| column: page
#| layout-ncol: 2


## Access to Hi speed internet using survey objects
mapPUMAboth <- svytable(~YEAR+PUMA+CIHISPEED, design = HHdesign)

mapPUMAboth <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, PUMA) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
     full_join(pumasIL, by = c("PUMA" = "PUMACE10"))

mapPUMAboth %>%
  filter(CIHISPEED == 10) %>% 
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Percent of population WITH High Speed Internet", 
       subtitle = "", caption = "Each PUMA has >100K & < 200K people.") + 
  facet_wrap(~YEAR)+
 scale_fill_binned(low ="white", high ="forestgreen", 
                       n.breaks = 4,
                   show.limits = TRUE,
                  # nice.breaks =FALSE,
                        name = "Population %", label = scales::percent)+Alea_theme()



mapPUMAboth %>% pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
  mutate(pct_change = `2021`-`2019`) %>%
  ggplot(aes(fill = pct_change)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Change in Access to High Speed Internet",
       subcaption = "2019 to 2021",
       caption = "Percentage point change:
       2021 PUMA% with CIHISPEED=10 - 2019 PUMA% 
       for ACS variable CIHISPEED=10.") +
  scale_fill_steps2(high = "forestgreen", mid =  "white",  low = "red4", 
                       n.breaks=4,
                  nice.breaks = FALSE,
                  show.limits = TRUE,
                        name = "Change in Access", label = scales::percent) + Alea_theme()


  

```

## County Level

19 is Champaign, 31 is Cook, 37 is DeKalb, 43 is DuPage, 89 is Kane, 111 is McHenry, etc.

#### Internet Access

```{r}
#| column: page
#| layout-ncol: 2

## Access internet using survey objects
mapPUMAboth <- svytable(~CINETHH+YEAR+countyFIP, design = HHdesign)

mapPUMAboth <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, countyFIP) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
   # filter(CINETHH == 3 ) %>% 
     full_join(countyIL, by = c("countyFIP" = "COUNTYFP"))

mapPUMAboth %>%
      filter(CINETHH == 3 ) %>% 

  filter(YEAR != "NA") %>%
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Percent of population that lacked access to Internet", 
       caption = "Geography = Counties with populations > ~60K") + 
  facet_wrap(~YEAR)  +
   theme_void() +
    theme(legend.title = element_blank())+
  scale_fill_steps2( low =  "white",  high = "red4", 
                       n.breaks=4,
                  nice.breaks = FALSE,
                  show.limits = TRUE,
                        "Population %", label = scales::percent) 

### Change in LACK OF ACCESS
mapPUMAboth %>% 
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
  mutate(pct_change = `2021`-`2019`) %>%
  ggplot(aes(fill = pct_change)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Change in access to Internet",
       caption = "Change in Percentage points = 2021 % lacking access - 2019 % lacking access.
       Uses CINETHH == 3 for not having internet access. A decrease is a good thing!" ) +
  scale_fill_steps2(low = "forestgreen", mid =  "white",  high = "red4", 
                       n.breaks=5,
                 # nice.breaks = FALSE,
                  show.limits = TRUE,
                        name = "Change in Access", label = scales::percent) + theme_void()

```

```{r}
#| label: fig-Figure8
#| column: page
#| layout-ncol: 2
#| code-fold: true


### Access to Hi speed internet ###

# from joined dataframe. Not using survey objects or survey package
# mapPUMAboth <- joined %>% 
#    group_by(YEAR, countyFIP, CINETHH) %>% 
#    dplyr::summarize(weightedcount=sum(PERWT), #weighted 
#                     unweightedcount = n()) %>% 
#    mutate(pct_weight = weightedcount/sum(weightedcount), 
#           pct_noweight = unweightedcount/sum(unweightedcount))%>%   
#   filter(CINETHH ==3) %>% 
#   full_join(countyIL, by = c("countyFIP" = "COUNTYFP")) 
# mapPUMAboth
# 
# mapPUMAboth %>%
#   ggplot(aes(fill = pct_weight)) +
#   geom_sf(aes(geometry = geometry), color = "black")+ 
#   labs(title = "Percent of population that did NOT have access to Internet", 
#        subtitle = "Geography = County", 
#        ) +  
#   Alea_theme() +theme(legend.title = element_blank())+
#   scale_fill_gradientn(colors = c( "white", "red4"), 
#                        n = 5,
#                         name = "Population %", label = scales::percent)+facet_wrap(~YEAR)





## Access to High speed internet using survey objects
mapPUMAboth <- svytable(~YEAR+countyFIP+CIHISPEED, design = HHdesign)

mapPUMAboth <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, countyFIP) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
    filter(CIHISPEED != 20 )  %>%   full_join(countyIL, by = c("countyFIP" = "COUNTYFP")) 


mapPUMAboth %>%
  filter(CIHISPEED == 10) %>%
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Percent of population with High-speed Internet"  ) + 
  facet_wrap(~YEAR)+
  scale_fill_binned(low  = "white", high = "green4", 
                    na.value = NA,
                    n.breaks = 4,
                    show.limits=TRUE,
                    nice.breaks = FALSE,
                    name = "Percent of Population", label = scales::percent)+ Alea_theme()


mapPUMAboth %>% 
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
  mutate(pct_change = `2021`-`2019`) %>%
  ggplot(aes(fill = pct_change*100)) +
  geom_sf(aes(geometry = geometry), color = "black") + 
  theme_void()+
   scale_fill_steps2(low ="red4", mid= "white", high = "forestgreen", 
                       n.breaks = 7,
                     na.value="gray",
                   #  nice.breaks=FALSE,
                     show.limits = TRUE,
             #          limits = c(-.2,.2),
                        name = "Pct Pt Change", 
            # label = scales::percent
             )
#+labs(title = "Change in access to High Speed Internet")
```


```{r export-figure8, include = FALSE}
mapPUMAboth %>%
  filter(CIHISPEED == 10) %>%
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Percent of population with High-speed Internet"  ) + 
  facet_wrap(~YEAR)+
  scale_fill_binned(low  = "white", high = "green4", 
                    na.value = NA,
                    n.breaks = 4,
                    show.limits=TRUE,
                    nice.breaks = FALSE,
                    name = "Percent of Population", label = scales::percent)+ Alea_theme()


mapPUMAboth %>% pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
  mutate(pct_change = `2021`-`2019`) %>%
  ggplot(aes(fill = pct_change)) +
  geom_sf(aes(geometry = geometry), color = "black") + 
  theme_void()+
  theme(legend.position = "none")+ 
   scale_fill_steps2(low ="red4", mid= "white", high = "forestgreen", 
                       n.breaks = 7,
                     na.value="gray",
                   #  nice.breaks=FALSE,
                     show.limits = TRUE,
             #          limits = c(-.2,.2),
                        #name = "Percent of Population", 
             #label = scales::percent
             )#+labs(title = "Change in access to High Speed Internet")
ggsave("./paper_figures/Figure8-map.png", limitsize = FALSE, width = 5, height = 4, units = "in")

ggsave("./paper_figures/Figure8-map.eps", limitsize = FALSE, width = 5, height = 4, units = "in")
```


```{r}
#| column: page
#| layout-ncol: 2
#| code-fold: true


# as a bar graph ## 
mapPUMAboth %>%
  filter(YEAR==2021) %>%
ggplot(aes(group = countyFIP)) +
  geom_bar(aes(y = Prop, x= reorder(NAME, Prop)), stat = "identity") +
  #Alea_theme() + 
  labs(title = "Proportion with High Speed Internet", subtitle = "In 2021", x = "", y = "" )+
  coord_flip()


# as a dot graph ## 

order <- mapPUMAboth %>% as_tibble() %>%
 filter(YEAR == 2019 & NAME != "NA" & CIHISPEED == "10") %>% 
  select(NAME, Prop_2019= Prop)

mapPUMAboth <- left_join(mapPUMAboth, order)

mapPUMAboth %>%
  filter(NAME != "NA" & YEAR != "NA")%>%
  ggplot(aes(x = Prop, y= reorder(NAME, Prop_2019)))+
  geom_line(aes(group = NAME))+
   geom_point(aes(color=YEAR), size = 3) +
  theme_minimal() + 
  theme(legend.position = "none", legend.title = element_blank(),
               plot.title.position = "plot",
     #   panel.background = element_rect(fill='transparent'), #transparent panel bg
    plot.background = element_rect(fill='transparent', color=NA) #transparent plot bg
   )+
        scale_color_brewer(palette="Paired", labels = c("2019", "2021"), direction = 1)+

  labs(title = "Change in Access to High-speed Internet", x = "", y = "" ) +
  scale_x_continuous(label = scales::percent)
```


```{r export-figure8dotplot}
mapPUMAboth %>%
  filter(NAME != "NA" & YEAR != "NA")%>%
  ggplot(aes(x = Prop, y= reorder(NAME, Prop_2019)))+
  geom_line(aes(group = NAME))+
   geom_point(aes(color=YEAR), size = 3) +
  theme_minimal() + 
  theme(legend.position = "none", legend.title = element_blank(),
               plot.title.position = "plot",
     #   panel.background = element_rect(fill='transparent'), #transparent panel bg
    plot.background = element_rect(fill='transparent', color=NA) #transparent plot bg
   )+
        scale_color_brewer(palette="Paired", labels = c("2019", "2021"), direction = 1) +
# 
#   labs(#title = "Change in Access to High-speed Internet", 
#     x = "", y = "" 
#     ) +
  scale_x_continuous(label = scales::percent)

ggsave("./paper_figures/Figure8-dotplot.png", limitsize = FALSE, width = 5, height = 4, units = "in")

ggsave("./paper_figures/Figure8-dotplot.eps", limitsize = FALSE,width = 5, height = 4, units = "in")

```

#### Did work from home

```{r}
#| column: page
#| layout-ncol: 2
#| code-fold: true


mapPUMAboth <- svytable(~did_wfh_labels+YEAR+countyFIP, design = dstrata)

mapPUMAboth <- mapPUMAboth %>% as_tibble() %>%
   group_by(YEAR, countyFIP) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
  full_join(countyIL, by = c("countyFIP" = "COUNTYFP"))

mapPUMAboth %>%
  filter(did_wfh_labels == "Did WFH") %>%
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Percent of Labor  Force that Worked from Home", 
       caption = "Geography = Counties with populations > ~60K") + 
  facet_wrap(~YEAR)  +
  scale_fill_binned(low= "white", high ="darkblue", 
                    na.value = "gray",
                    show.limits=TRUE,
                    n.breaks = 4,
                    name = "Population %", 
                    label = scales::percent) + 
  Alea_theme()

mapPUMAboth %>% 
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
  mutate(pct_change = `2021`-`2019`) %>%
  ggplot(aes(fill = pct_change*100)) +
  geom_sf(aes(geometry = geometry), color = "black")+ Alea_theme()+ 
  scale_fill_binned(low= "white", high ="darkblue", 
                    na.value = "gray",
                    show.limits=TRUE,
                    nice.breaks=FALSE,
                    n.breaks = 4,
                    limits = c(0,20),
                    name = "Percentage Point\nIncrease"
                    #label = scales::percent
                    ) +
  labs(title = "Percentage Point Change in Labor Force that Did WFH")

```

```{r export-figure4-map, include=FALSE, eval=FALSE}
mapPUMAboth %>% 
  pivot_wider(id_cols = -c(n), names_from ="YEAR", values_from = "Prop") %>%
  mutate(pct_change = `2021`-`2019`) %>%
  ggplot(aes(fill = pct_change*100)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  theme_void()+
  theme(legend.position = "none") +
  
  scale_fill_binned(low= "white", high ="darkblue", 
                    na.value = "gray",
                    show.limits=TRUE,
                    nice.breaks=FALSE,
                    n.breaks = 4,
                    limits = c(0,20),
                    name = "Percentage Point\nIncrease"
                   # label = scales::percent
                    ) #+
 # labs(title = "Percentage Point Change in Labor Force that Did WFH")

ggsave("./paper_figures/Figure4-map.eps", limitsize = FALSE, width = 4, height = 5, units = "in")
ggsave("./paper_figures/Figure4-map.png", limitsize = FALSE, width = 4, height = 5, units = "in")
```

```{r}
#| column: page
#| layout-ncol: 2
#| code-fold: true


# as a dot graph ## 
order <- mapPUMAboth %>% as_tibble() %>%
 filter(YEAR == 2021 & did_wfh_labels == "Did WFH") %>% 
  select(countyFIP, Prop_2021 = Prop)

mapPUMAboth <- left_join(mapPUMAboth, order)


mapPUMAboth %>%
  filter(did_wfh_labels == "Did WFH" & NAME != "NA" & YEAR != "NA")%>%
  ggplot(aes(x = Prop, y= reorder(NAME, Prop_2021)))+
  geom_line(aes(group = NAME))+
  theme_classic()+
   geom_point(aes(color=YEAR), size = 3) +
      scale_color_brewer(palette="Paired",
                         # labels = c("% WFH in 2019", "% WFH in 2021"), 
                         direction = 1)+
  labs(title = "Change in County Labor Force that Did WFH",
#  caption = "DuPage had around 8% of the labor force working from home in 2019 and 
#  around 27% of its labor force working from home in 2021. 
#       There was ~20 percentage point increase between 2019 and 2021 for DuPage County. 
#  Working from home is based on TRANWORK variable from ACS 1-year surveys.", 
x = "Labor Force Working from Home (%)", 
y = "" 
)+  scale_x_continuous(label = scales::percent)

```

```{r export-figure4-dotplot, include=FALSE, eval=FALSE}

mapPUMAboth %>%
  filter(did_wfh_labels == "Did WFH" & NAME != "NA" & YEAR != "NA")%>%
  ggplot(aes(x = Prop*100, y= reorder(NAME, Prop_2021*100)))+
  geom_line(aes(group = NAME))+
  theme_classic()+
  theme(legend.position = "none") + 
   geom_point(aes(color=YEAR), size = 3) +
      scale_color_brewer(palette="Paired",
                         # labels = c("% WFH in 2019", "% WFH in 2021"), 
                         direction = 1)+
  labs(#title = "Change in County Labor Force that did WFH",
#  caption = "DuPage had around 8% of the labor force working from home in 2019 and 
#  around 27% of its labor force working from home in 2021. 
#       There was ~20 percentage point increase between 2019 and 2021 for DuPage County. 
#  Working from home is based on TRANWORK variable from ACS 1-year surveys.", 
x = "Labor Force Working From Home (%)", 
y = "" 
)#+  scale_x_continuous(label = scales::percent)

ggsave("./paper_figures/Figure4-dotplot.eps", limitsize = FALSE,width = 8, height = 5, units = "in")
ggsave("./paper_figures/Figure4-dotplot.png", limitsize = FALSE, width = 8, height = 5, units = "in")
```

```{r wfh-bycounty, include=FALSE}
#| code-fold: TRUE

# calculates a value for rural counties to add contrast to counties that do have data. 
# would need to recalculate the "rural" average that 
# is hard coded into the code below.
mapdataboth <- joined %>% 
  filter(did_wfh_labels != "NA" )%>% # create percentages withing missing values included. aka valid percent.
   group_by(county_pop_type, YEAR, countyFIP, did_wfh_labels) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(countyIL, by = c("countyFIP" = "COUNTYFP"))  %>%
    mutate(pct_weight = ifelse(is.na(county_pop_type), 0.04486351, pct_weight))%>%
 #  mutate(pct_weight = ifelse(is.na(county_pop_type), 0.09237573, pct_weight))# %>%
  mutate(did_wfh_labels = ifelse(is.na(did_wfh_labels),"Rural", did_wfh_labels))
mapdataboth

mapdataboth %>%
  filter(did_wfh_labels != "Did not WFH") %>%
  ggplot(aes(fill = pct_weight)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
    Alea_theme()+

  labs(title = "Percent of county population that did work from home in 2021", 
  subtitle = "Countyfips != 000", 
  caption = "Note: Shows all counties with populations large enough to have their own values.
  Rural counties went from 4.5% WFH to 9.2% WFH") + 
  facet_wrap(~YEAR)

mapdata2019 <- joined %>% 
  filter(YEAR == 2019 & did_wfh_labels != "NA" )%>% # create percentages withing missing values included. aka valid percent.
   group_by(countyFIP, did_wfh_labels) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(countyIL, by = c("countyFIP" = "COUNTYFP"))  %>%
    mutate(pct_weight = ifelse( is.na(did_wfh_labels), 0.04486351, pct_weight),
           did_wfh_labels = ifelse(is.na(did_wfh_labels),"Rural", did_wfh_labels))


mapdata2021 <- joined %>% 
  filter(YEAR == 2021 & did_wfh_labels != "NA" )%>% # create percentages withing missing values included. aka valid percent.
   group_by(countyFIP, did_wfh_labels) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(countyIL, by = c("countyFIP" = "COUNTYFP"))  %>%
    mutate(pct_weight = ifelse( is.na(did_wfh_labels), 0.09237573, pct_weight),
           did_wfh_labels = ifelse(is.na(did_wfh_labels),"Rural", did_wfh_labels))

# ### excludes rural counties unintentially due to filter ###
# mapdata2021 %>%
#   filter(did_wfh_labels == "Did WFH") %>%
#   ggplot(aes(fill = pct_weight)) +
#     Alea_theme()+
#   geom_sf(aes(geometry = geometry), color = "black")+ labs(title = "Percent of county population that did work from home in 2021", subtitle = "Countyfips != 000", caption = "Note: Shows all counties with populations large enough to have their own values")

### Keeps empty counties in map ### 
mapdata2021 %>%
  filter(did_wfh_labels != "Did not WFH") %>%
  ggplot(aes(fill = pct_weight)) +
  geom_sf(aes(geometry = geometry), color="black")+ 
    Alea_theme()+

  labs(title = "Percent of population that did WFH in 2021",
       caption = "Note: All counties with populations < 50,000 people were 
       aggregated together in ACS data extracts (COUNTYFIP == 000)") 

mapdata2019 %>%
  filter(did_wfh_labels != "Did not WFH") %>%
  ggplot(aes(fill = pct_weight)) +
  geom_sf(aes(geometry = geometry), color="black")+ 
  Alea_theme()+
  labs(title = "Percent of population that did WFH in 2019",
       caption = "Note: All counties with populations < 50,000 people were 
       aggregated together in ACS data extracts (COUNTYFIP == 000)") 
```

#### Could work from home

Based on Dingel & Niemen paper. A person's occupation indicates the likelihood of if they could work from home.

```{r Figure2-couldWFH-bycounty}
#| column: page
#| layout-ncol: 2
#| code-fold: true


# teleworkable is numeric
# CanWorkFromHome is categorical variable with 3 options. 


mapdataboth <- joined %>% 
  filter(CanWorkFromHome != "NA" )%>% # create percentages withing missing values included. aka valid percent.
   group_by(county_pop_type, YEAR, countyFIP, CanWorkFromHome) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()
                    ) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(countyIL, by = c("countyFIP" = "COUNTYFP")) # %>%
    #mutate(pct_weight = ifelse(is.na(county_pop_type), 0.04486351, pct_weight))%>%
 #  mutate(pct_weight = ifelse(is.na(county_pop_type), 0.09237573, pct_weight))# %>%
 # mutate(did_wfh_labels = ifelse(is.na(did_wfh_labels),"Rural", did_wfh_labels))
mapdataboth 



mapdataboth <- svytable(~CanWorkFromHome+YEAR+countyFIP, design = dstrata)

mapdataboth <- mapdataboth %>% as_tibble() %>%
   group_by(YEAR, countyFIP) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
   # filter(CanWorkFromHome == "Can WFH") %>%
filter(YEAR == 2021 & CanWorkFromHome == "Can WFH")%>%  

     full_join(countyIL, by = c("countyFIP" = "COUNTYFP"))

mapdataboth %>%
  filter(CanWorkFromHome == "Can WFH") %>%
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs(title = "Percent of county population that COULD work from home", subtitle = "Countyfips != 000", caption = "Note: Shows all counties with populations large enough to have their own values.
                                                           Rural counties went from ...
                                                           Counties had very small changes in occupations that COULD be done at work.") + 
  facet_wrap(~YEAR)+
  Alea_theme()+ 
 scale_fill_binned(low= "white", high ="darkblue", 
                   na.value = "gray",
                   show.limits=TRUE,
                   nice.breaks=FALSE,
                       n.breaks = 4,
                        limits = c(0,.45),
      label = scales::percent)


Figure2a <- mapdataboth %>%
#filter(YEAR == 2021 & CanWorkFromHome == "Can WFH")%>%  
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  # labs(title = "Work from Home Feasibility Rate", 
  #      subtitle = " Based on Occupation Characteristics", 
  #      caption = "Teleworkability based on D&N 2020 methodology. 
  #      OCCSOC codes from ACS 1 year Survey data at individual level.") +
 theme_void() + 
  theme(axis.ticks = element_blank(), axis.text = element_blank(),   
        plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust=.5),  
        plot.subtitle = element_text(hjust = 0.5)) +
 scale_fill_binned(low= "white", high ="darkblue", 
                   na.value = "gray",
                   show.limits=TRUE,
                   nice.breaks=FALSE,
                       n.breaks = 4,
                       limits = c(0,.6),
                        name = "Labor Force with\nWFH Feasibility", 
                       label = scales::percent)

Figure2a



# ggsave("Figure2a.eps", limitsize = FALSE,width = 8, height = 4, units = "in")
# ggsave("Figure2a.pdf", limitsize = FALSE,width = 8, height = 4, units = "in")
#ggsave("Figure2a.png", limitsize = FALSE, width = 8, height = 4, units = "in")

Figure2b

# ggsave("Figure2b.eps", limitsize = FALSE,width = 8, height = 4, units = "in")
# ggsave("Figure2b.pdf", limitsize = FALSE,width = 8, height = 4, units = "in")
#ggsave("Figure2b.png", limitsize = FALSE, width = 8, height = 4, units = "in")
```

```{r export-Figure2-couldWFH-bycounty, include=FALSE}

# teleworkable is numeric
# CanWorkFromHome is categorical variable with 3 options. 


mapdataboth <- joined %>% 
  filter(CanWorkFromHome != "NA" )%>% # create percentages withing missing values included. aka valid percent.
   group_by(county_pop_type, YEAR, countyFIP, CanWorkFromHome) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n()
                    ) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount))%>%
  full_join(countyIL, by = c("countyFIP" = "COUNTYFP")) # %>%
    #mutate(pct_weight = ifelse(is.na(county_pop_type), 0.04486351, pct_weight))%>%
 #  mutate(pct_weight = ifelse(is.na(county_pop_type), 0.09237573, pct_weight))# %>%
 # mutate(did_wfh_labels = ifelse(is.na(did_wfh_labels),"Rural", did_wfh_labels))
mapdataboth 



mapdataboth <- svytable(~CanWorkFromHome+YEAR+countyFIP, design = dstrata)

mapdataboth <- mapdataboth %>% as_tibble() %>%
   group_by(YEAR, countyFIP) %>% 
     mutate(Prop = round(n/sum(n), digits =3)) %>%
   # filter(CanWorkFromHome == "Can WFH") %>%
filter(YEAR == 2021 & CanWorkFromHome == "Can WFH")%>%  

     full_join(countyIL, by = c("countyFIP" = "COUNTYFP"))

mapdataboth %>%
  filter(CanWorkFromHome == "Can WFH") %>%
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  labs() + 
  facet_wrap(~YEAR)+
  Alea_theme()+ 
  theme(legend.position="none") + 
 scale_fill_binned(low= "white", high ="darkblue", 
                   na.value = "gray",
                   show.limits=TRUE,
                   nice.breaks=FALSE,
                       n.breaks = 4,
                        limits = c(0,.45),
      label = scales::percent)


Figure2a <- mapdataboth %>%
#filter(YEAR == 2021 & CanWorkFromHome == "Can WFH")%>%  
  ggplot(aes(fill = Prop)) +
  geom_sf(aes(geometry = geometry), color = "black")+ 
  # labs(title = "Work from Home Feasibility Rate", 
  #      subtitle = " Based on Occupation Characteristics", 
  #      caption = "Teleworkability based on D&N 2020 methodology. 
  #      OCCSOC codes from ACS 1 year Survey data at individual level.") +
 theme_void() + 
  theme(axis.ticks = element_blank(), axis.text = element_blank(),   
        plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust=.5),  
        plot.subtitle = element_text(hjust = 0.5)) +
 scale_fill_binned(low= "white", high ="darkblue", 
                   na.value = "gray",
                   show.limits=TRUE,
                   nice.breaks=FALSE,
                       n.breaks = 4,
                       limits = c(0,.6),
                        name = "Labor Force with\nWFH Feasibility", 
                       label = scales::percent)

Figure2a



ggsave("./paper_figures/Figure2a.eps", limitsize = FALSE,width = 8, height = 4, units = "in")
ggsave("./paper_figures/Figure2a.png", limitsize = FALSE, width = 8, height = 4, units = "in")

Figure2b

ggsave("./paper_figures/Figure2b.eps", limitsize = FALSE,width = 8, height = 4, units = "in")
ggsave("./paper_figures/Figure2b.png", limitsize = FALSE, width = 8, height = 4, units = "in")
```


```{r teleworkable-bycounty}
mapdataboth <- joined %>% 
 # filter(CanWorkFromHome != "NA" )%>% # create percentages withing missing values included. aka valid percent.
   group_by(county_pop_type, YEAR, countyFIP) %>% 
   dplyr::summarize(weightedcount=sum(PERWT), #weighted 
                    unweightedcount = n(),
                    mean_telework = mean(teleworkable, na.rm=TRUE, weight = PERWT)
                    ) %>% 
   mutate(pct_weight = weightedcount/sum(weightedcount), 
          pct_noweight = unweightedcount/sum(unweightedcount)) %>%
  full_join(countyIL, by = c("countyFIP" = "COUNTYFP"))  %>%
    mutate(mean_telework = ifelse(is.na(county_pop_type), 0.2941470, mean_telework))
 #  mutate(pct_weight = ifelse(is.na(county_pop_type), 0.09237573, pct_weight))# %>%
 # mutate(did_wfh_labels = ifelse(is.na(did_wfh_labels),"Rural", did_wfh_labels))
mapdataboth

mapdataboth %>%
  filter(YEAR != "NA") %>%
  ggplot(aes(fill = mean_telework)) +
  geom_sf(aes(geometry = geometry), color = "black")+ labs(title = "Average teleworkable score for each County", subtitle = "Countyfips != 000", caption = "Note: Shows all counties with populations large enough to have their own values.
                                                           Uses teleworkable, a continuous variable, instead of CanWorkFromHome.
                                                           Counties had very small changes in occupations that COULD be done at work.") + 
  facet_wrap(~YEAR)+ 
  Alea_theme()+
 scale_fill_binned(low= "white", high ="darkblue", 
                   na.value = "gray",
                   show.limits=TRUE,
                   nice.breaks=FALSE,
                       n.breaks = 4, 
                   label = scales::percent)
```


```{r mixed-geography, warning = FALSE, message=FALSE, eval=FALSE, include=FALSE}
# 
pums_weighted <- joined %>%
  filter(YEAR == 2021) %>%
  mutate(county_pop_type = if_else(COUNTYFIP==0, "Rural Counties", "Urban Counties")
        ) %>%
  dplyr::group_by(PUMA, COUNTYFIP, county_pop_type) %>%
  dplyr::summarize(weighted_obs = sum(PERWT),
            mean_telework = mean(teleworkable, weight=PERWT, na.rm=TRUE),
         #   mean_canwfh = mean(CanWorkFromHome), ## not numeric variable
            mean_didwfh = mean(did_wfh, weight = PERWT, na.rm=TRUE),
            observs = n(),
            avg_inc = mean(INCEARN,na.rm=TRUE)
           # avg_inc_w = survey_mean(INCEARN, weight=PERWT
         ) %>%  #number of people the sample represents
  mutate(PUMA = str_pad(PUMA, 5, pad="0"),
         countyFIP = str_pad(COUNTYFIP, 3, pad = "0"))

pums_unweight <- joined %>% 
  filter(YEAR == 2021) %>%
    mutate(county_pop_type = if_else(COUNTYFIP==0, "Rural Counties", "Urban Counties")) %>%

  group_by(PUMA, COUNTYFIP, county_pop_type) %>% 
  summarize(unweight = n(),#unweighted number of observations
            mean_telework = mean(teleworkable,na.rm=TRUE ),
            mean_didwfh = mean(did_wfh, na.rm=TRUE),
            avg_inc = mean(INCEARN,na.rm=TRUE)) %>% 
  mutate(PUMA = str_pad(PUMA, 5, pad="0"),
         countyFIP = str_pad(COUNTYFIP, 3, pad = "0"))

plotweighted <- pumasIL %>% 
  left_join(pums_weighted, by = c("PUMACE10" = "PUMA"))

plotunweight <- pumasIL %>% 
  left_join(pums_unweight, by = c("PUMACE10" = "PUMA"))

## PUMA boundaries
plot(plotweighted["weighted_obs"])
plot(plotunweight["unweight"])


plot(plotweighted["avg_inc"])
#plot(plotweighted["avg_inc_w"])
plot(plotunweight["avg_inc"])
#plot(plotweighted["avg_inc_w"])

plot(plotweighted["mean_didwfh"])
plot(plotunweight["mean_didwfh"])



FIPweighted <- countyIL %>% left_join(pums_weighted, by = c("COUNTYFP" = "countyFIP"))
FIPunweight <- countyIL %>% left_join(pums_unweight, by = c("COUNTYFP" = "countyFIP"))
plot(FIPweighted["weighted_obs"])
plot(FIPunweight["unweight"])


plot(countyIL["COUNTYFP"])


# Observations using `CITY` variable: identifies observations from Chicago.

# 88 distinct PUMA areas.


#joined %>% group_by(CITY) %>% distinct(CITY)
joined %>% group_by(PUMA) %>% summarize(observs = n())
joined %>% filter(INCEARN > 0) %>% group_by(PUMA) %>% summarize(observs = n())

joined %>% group_by(COUNTYFIP, PUMA) %>% summarize(observs = n())

joined %>% 
  group_by(COUNTYFIP, PUMA) %>% 
  summarize(observs = n(),
            avg_inc = mean(INCEARN),
            avg_inc_w= mean(INCEARN, weight = PERWT))
```



#### Geography Variable Notes

**PUMAS vs COMMZONE vs COUNTIES**

[Link from Francis on COMZONE variable](https://eig.org/the-uneven-geography-of-remote-work/) (Commuter Zones)

[Interactive ESRI Map of all PUMA outlines](https://univofillinois.maps.arcgis.com/apps/mapviewer/index.html?layers=71eb4ef236ca42c4b577a0a105e790e2)

[Article on calculating mean income for groups of geographies with ACS data](https://atcoordinates.info/2019/05/13/calculating-mean-income-for-groups-of-geographies-with-census-acs-data/)

```{r stateworkedin, eval=FALSE}
# State worked in:
#0=NA, 17=Illinois



# ipums_var_desc(data, PWSTATE2)

joined <- joined %>% 
  mutate(PWSTATE2_clean = as_factor(lbl_na_if(PWSTATE2, ~.val %in% c(0))))

joined %>% 
 dplyr::group_by(YEAR) %>%
  dplyr::summarize(workers=sum(PERWT)) %>% #number of people that match that observation
  dplyr::ungroup()%>%
 dplyr:: group_by(YEAR,PWSTATE2_clean)%>%
  mutate(pct = n()/workers)%>% 
  arrange(desc(pct))
```


## Other Sources and Papers

Ability to work from home: evidence from two surveys and implications for the labor market in the COVID-19 pandemic. June 2020. [BLS Monthly Labor Review](https://www.bls.gov/opub/mlr/2020/article/ability-to-work-from-home.htm)

-   Authors used Current Population Survey data and O\*NET job-content data to categorize jobs as able or unable to telework. Followed Dingel & Neiman's methodology of classifying telework feasibility and merging with data from American Time Use Survey (ATUS)

-   Compares ability to work from home with actual occurance of working from home based on American Time Use Survey and Occupational Information Network (O\*NET). Also uses Current Population Survey data to look at how effects differed between occupations where telework was feasible or not.