I combine the 2019 and 2021 datasets when creating new variables and then separate them by year again before creating the survey design object for each sample year.
Last week, did this person work for pay at a job or business? (Yes or no) – Yes becomes coded as EMPSTAT = 1-Employed.
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:
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.
# 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.xmlddi <-read_ipums_ddi("usa_00011.xml") # downloaded April 10 2023data2021 <-read_ipums_micro(ddi) # 126623 observations before any filtering
Use of data from IPUMS USA is subject to conditions including that users should
cite the data appropriately. Use command `ipums_conditions()` for more details.
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 referneces the data file. currently says usa_00011.dat.gz so these two lines of code do not work. ddi <-read_ipums_ddi("C:/Users/aleaw/Box/Fiscal Futures/FY22_Working/WFH/Data/IL_2019_1yearACS_datDDI.xml") # downloaded April 10 2023data2019 <-read_ipums_micro(ddi) # 126623 observations before any filtering
Use of data from IPUMS USA is subject to conditions including that users should
cite the data appropriately. Use command `ipums_conditions()` for more details.
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 valuesdata <- 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 objectdata <- 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 latermutate(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 & Extractraction",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 wfhdata <- 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 = 1has_occsoc =ifelse(OCCSOC >0, 1, 0),# has occupation = 1has_incearn_labels =ifelse(INCEARN >0, "Has EarnInc", "No IncData"), ## has earned income = 1has_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.
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.
# 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
# bring in the teleworkable scores based on D&N's work.telework <-read_csv("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)
Can WFH No WFH Some WFH
41643 66111 15999
table(joined$did_wfh, joined$YEAR)
2019 2021
0 56782 47136
1 3197 10949
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.
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.
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
#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 itemdstrata2019 <- 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.
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.
# 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()
did_wfh_labels
under5 Did not WFH Did WFH
0 174904 91658
1 49846 29442
Pearson's X^2: Rao & Scott adjustment
data: svychisq(~under5 + did_wfh_labels, design = kidsdesign, statistic = "F")
F = 1.0905, ndf = 1, ddf = 3182, p-value = 0.2964
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.
Design-based t-test
data: under5 ~ did_wfh
t = 1.6909, df = 3987, p-value = 0.09094
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
-0.004807002 0.065089151
sample estimates:
difference in mean
0.03014107
# 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()
did_wfh_labels
under5 Did not WFH Did WFH
0 149045 101265
1 38192 26053
Pearson's X^2: Rao & Scott adjustment
data: svychisq(~under5 + did_wfh_labels, design = kidsdesign, statistic = "F")
F = 0.0012511, ndf = 1, ddf = 2729, p-value = 0.9718
Design-based t-test
data: under5 ~ did_wfh
t = 0.42896, df = 3490, p-value = 0.668
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
-0.02531256 0.03949063
sample estimates:
difference in mean
0.007089036
Call:
glm(formula = did_wfh ~ under5, data = childttest2021)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.3924 -0.3924 -0.3838 0.6076 0.6162
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.392423 0.009225 42.540 <2e-16 ***
under5 -0.008596 0.018875 -0.455 0.649
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.2380974)
Null deviance: 874.82 on 3675 degrees of freedom
Residual deviance: 874.77 on 3674 degrees of freedom
(103 observations deleted due to missingness)
AIC: 5160.7
Number of Fisher Scoring iterations: 2
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.
Code
joined <- joined %>%filter(HHINCOME >0& HHINCOME!=9999999& HHINCOME !=9999998) # 105 observations joined %>%ggplot() +geom_histogram(aes(x=HHINCOME, weight = HHWT))
Code
joined %>%# 122000 observations ggplot() +geom_histogram(aes(x=HHINCOME))
Code
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
Code
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.
Code
#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)
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.
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
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())+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 axsiscoord_flip()
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.
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.
### 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
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.
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())+labs(title ="WFH Feasibility by Occupation Type",#subtitle = "Little change between 2019 and 2021 Occurred",x ="", y ="Estimated Number of Workers") +scale_y_continuous(labels = scales::comma)+coord_flip()
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.
### Percent of Workersworking 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
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.
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()
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 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 bgplot.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 axsisscale_fill_manual(values =c("#a6bddb", "#2b8cbe")) +coord_flip()
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())+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()
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.
# Both years, detailed observation typestable <-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()
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.
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(title ="Percent of Workers that Could Feasibily Work From Home in 2021", 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()
## 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_blank())+labs(title ="Proportion of All WFH Workers by Broad Occupation",# subtitle = "All workers in labor force with occsoc codes in a year add to 100%",x ="", y ="Number of Workers in Illinois",# 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 axsisscale_y_continuous(labels = scales::comma) +scale_fill_manual(values =c( "#2b8cbe","#a6bddb", "gray89")) +coord_flip() # = element_text(hjust = 0, vjust=2.12))Figure1
## Totals add up to total number of workers in a yeartable <-svytable(~CanWorkFromHome+YEAR+did_wfh_labels, design = dstrata) # proportion of each respondant's sex and if they worked from home for each year in sampletable <- table %>%as_tibble() %>%group_by(YEAR)%>%# will divide by all workers per yearmutate(Prop =round(n/sum(n), digits=3)) %>%arrange(did_wfh_labels, -n)table
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
## percentages add up to 100 when adding all workers together for a yeartable %>%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"))
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.
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 sampletable <- table %>%as_tibble() %>%group_by(YEAR, CanWorkFromHome)%>%# divides by all workers per year within each categormutate(Prop =round(n/sum(n), digits=3)) %>%arrange(did_wfh_labels,-n)table
xtabs(~did_wfh_labels+CanWorkFromHome+YEAR, data = dstrata)
, , YEAR = 2019
CanWorkFromHome
did_wfh_labels Can WFH No WFH Some WFH
Did not WFH 18433 30664 7094
Did WFH 1484 971 596
, , YEAR = 2021
CanWorkFromHome
did_wfh_labels Can WFH No WFH Some WFH
Did not WFH 13594 27599 5380
Did WFH 6217 2225 2370
## percentages add up to 100 when adding all workers together for a yeartable %>%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 categorfilter(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)
, , YEAR = 2019
CanWorkFromHome
did_wfh_labels Can WFH No WFH Some WFH
Did not WFH 18433 30664 7094
Did WFH 1484 971 596
, , YEAR = 2021
CanWorkFromHome
did_wfh_labels Can WFH No WFH Some WFH
Did not WFH 13594 27599 5380
Did WFH 6217 2225 2370
## percentages add up to 100 when adding all workers together for a yeartable %>%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
# 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
# $INCEARN# [1,] 4 7200 16000 25000 34000 42000 52000 67000 85000 120000 949000breaks2019 <-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 ACSbreaks2021 <-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 sampletable <- 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)
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)
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.
#round(prop.table(svytable(~did_wfh_labels, design=dstrata2019))*100,digits=2) #94.7 did not wfh, 5.29 did work from home in 2019round(prop.table(svytable(~did_wfh_labels, design=dstrata2021))*100,digits=2)
did_wfh_labels
Did not WFH Did WFH
80.76 19.24
# # 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 yeartable <-svytable(~SEX+YEAR+did_wfh_labels, design = dstrata) # proportion of each respondant's sex and if they worked from home for each year in sampletable <- 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 yeartable %>%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)
Design-based t-test
data: SEX ~ did_wfh_labels
t = 8.3, df = 35419, p-value < 2.2e-16
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
0.04071947 0.06589675
sample estimates:
difference in mean
0.05330811
# slight drop in women in the workforce (technically women with occupations). round(prop.table(svytable(~has_occsoc+SEX, design=dstrata2021))*100,digits=2)
Call: xtabs(formula = ~SEX + YEAR, data = data)
Number of cases in table: 123753
Number of factors: 2
Test for independence of all factors:
Chisq = 1.113, df = 1, p-value = 0.2914
xtabs(~did_wfh+SEX+YEAR, data = data) %>%summary()
Call: xtabs(formula = ~did_wfh + SEX + YEAR, data = data)
Number of cases in table: 118064
Number of factors: 3
Test for independence of all factors:
Chisq = 5246, df = 4, p-value = 0
## Most interesting change: round(prop.table(svytable(~SEX+did_wfh_labels, design=dstrata2019))*100,digits=2)
did_wfh_labels
SEX Did not WFH Did WFH
1 49.27 2.58
2 45.57 2.58
# 1 is Male, 2 is femaleround(prop.table(svytable(~SEX+did_wfh_labels, design=dstrata2021))*100,digits=2)
did_wfh_labels
SEX Did not WFH Did WFH
1 43.21 9.27
2 37.55 9.97
# 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)
Chi-squared test for given probabilities
data: table(joined$SEX)
X-squared = 0.0049174, df = 1, p-value = 0.9441
# 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)
total SE
interaction(SEX, did_wfh_labels)1.Did not WFH 2562380 20128.2
interaction(SEX, did_wfh_labels)2.Did not WFH 2226558 18405.0
interaction(SEX, did_wfh_labels)1.Did WFH 549602 9349.2
interaction(SEX, did_wfh_labels)2.Did WFH 591233 9978.3
Design-based t-test
data: INCEARN ~ SEX
t = -32.274, df = 36475, p-value < 2.2e-16
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
-23794.51 -21069.87
sample estimates:
difference in mean
-22432.19
Design-based t-test
data: INCEARN ~ did_wfh_labels
t = 27.63, df = 35419, p-value < 2.2e-16
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
32054.25 36949.32
sample estimates:
difference in mean
34501.79
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 sampletable <- 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.
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.
Code
table <-svytable(~YEAR+incdecile_w+CanWorkFromHome, design = dstrata) # proportion of each respondants sex in sampletable <- 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
Code
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)
Code
# 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)
# data has 2019 and 2021 observatins TOGETHERtopline(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)
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-
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
19 is Champaign, 31 is Cook, 37 is DeKalb, 43 is DuPage, 89 is Kane, 111 is McHenry, etc.
joined %>%filter(CanWorkFromHome =="Some WFH") %>%distinct(OCCSOC)
Race
Good example of graphing survey data. 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.
Code
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 sampletable <- 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
Code
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)
Code
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 sampletable <- table %>%as_tibble() %>%group_by(YEAR, did_wfh_labels)%>%# grouped by Year and Sex!!mutate(Prop =round(n/sum(n), digits=4))table
Code
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)
Code
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 sampletable <- 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)
Code
svyby(~did_wfh, by =~RACE, design = dstrata2021, svymean, na.rm=TRUE)
Code
svyby(~INCEARN, by =~RACE, design = dstrata2021, FUN = svymean, na.rm=TRUE)
Code
joined %>%crosstab_3way(x=YEAR, y=did_wfh, z=white, weight = PERWT)
Code
joined %>%crosstab_3way(x=YEAR, y=did_wfh, z=RACE, weight = PERWT)
Code
# 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)
joined %>%filter(YEAR ==2021) %>%group_by(RACE) %>%summarize(Freq =n()) %>%mutate(Prop = Freq/sum(Freq))
Code
svytable(~RACE, design=dstrata2021) %>%as.data.frame() %>%# creates a frequency count by default?mutate(Prop =Freq/sum(Freq))
Code
topline(dstrata2019, RACE, weight = PERWT)
Code
topline(dstrata2021, RACE, weight = PERWT)
Code
#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
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 17table(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.
#10 is yes, 20 is no access. 00 is NAtable <-svytable(~CIHISPEED +did_wfh_labels+YEAR, design = dstrata) topline(dstrata2019, CIHISPEED, weight = PERWT)
topline(dstrata2021, CIHISPEED, weight = PERWT)
table
, , YEAR = 2019
did_wfh_labels
CIHISPEED Did not WFH Did WFH
10 4613288 277027
20 878321 29514
, , YEAR = 2021
did_wfh_labels
CIHISPEED Did not WFH Did WFH
10 3871147 1027194
20 683274 82666
# proportion of each respondant's sex and if they worked from home for each year in sampletable <- 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
, , YEAR = 2019
did_wfh_labels
CINETHH Did not WFH Did WFH
1 5491609 306541
2 109151 3768
3 249879 7865
, , YEAR = 2021
did_wfh_labels
CINETHH Did not WFH Did WFH
1 4554421 1109860
2 94608 19191
3 139909 11784
# proportion of each respondant's sex and if they worked from home for each year in sampletable <- 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
, , did_wfh_labels = Did not WFH, YEAR = 2019
CINETHH
CIHISPEED 1
10 4613288
20 878321
, , did_wfh_labels = Did WFH, YEAR = 2019
CINETHH
CIHISPEED 1
10 277027
20 29514
, , did_wfh_labels = Did not WFH, YEAR = 2021
CINETHH
CIHISPEED 1
10 3871147
20 683274
, , did_wfh_labels = Did WFH, YEAR = 2021
CINETHH
CIHISPEED 1
10 1027194
20 82666
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.
, , work_in_IL = In Illinois
YEAR
did_wfh_labels 2019 2021
Did not WFH 5632943 4604673
Did WFH 318174 1140835
, , work_in_IL = Out of IL
YEAR
did_wfh_labels 2019 2021
Did not WFH 217696 184265
Did WFH 0 0
318174/ (318174+5632943) # percent of remote workers that worked outside of IL but lived in IL in 2019
[1] 0.05346458
1140835/(1140835+4604673) # 2021 percent of WFH workers that worked outside of IL
, , work_in_IL = In Illinois
YEAR
has_incearn_labels 2019 2021
Has EarnInc 5951117 5745508
, , work_in_IL = Out of IL
YEAR
has_incearn_labels 2019 2021
Has EarnInc 217696 184265
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
obs_perPUMA<- joined %>%group_by(PUMA, YEAR) %>% dplyr::summarize(weightedcount=sum(PERWT), #weighted unweightedcount =n()) %>%arrange(unweightedcount)### Minimimum observations are 271 in PUMA area 03528obs_perPUMA
# first attempt to graph data. Does appear to match the 2nd attempt below.# mapPUMAboth <- joined %>% # filter(did_wfh_labels != "NA" )%>% # create percentages withing missing values included. aka valid percent.# group_by(county_pop_type, YEAR, 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")) # mapPUMAboth# # mapPUMAboth %>%# filter(did_wfh_labels != "Did not WFH") %>%# ggplot(aes(fill = pct_weight)) +# geom_sf(aes(geometry = geometry), color = "black")+ # labs(title = "Percent of PUMA area population that did work from home in 2021", # subtitle = "Each PUMA has >100K & < 200K people", # caption = "Facet_wrapped by year") + # facet_wrap(~YEAR)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"))mapPUMAboth %>%select(YEAR, PUMA, Prop, n, NAME10)
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.
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_gradientn(colors =c("white", "darkblue"), limits =c(0,.6),n =6,name ="Percent of Workers", label = scales::percent)+# geom_sf(data = countyIL, fill=NA, color="dark gray") +facet_wrap(~YEAR)
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.
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_gradientn(colors =c("white", "darkblue"), limits =c(0,.6),n =6,name ="% change in WFH", label = scales::percent)
# 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 %>%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 County 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 ="")
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.
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.
figure5b <-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_gradientn(colors =c("white", "darkblue"), limits =c(0,.6),n =6,name ="% that did WFH", label = scales::percent)+facet_wrap(~YEAR) figure5b
mapPUMAboth <-svytable(~YEAR+PUMA+countyFIP+did_wfh_labels, design = dstrata)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_gradientn(colors =c("white", "darkblue"), limits =c(0,.6),n =6,name ="% change in WFH", label = scales::percent)
# 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.
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_gradientn(colors =c("white", "darkblue"), limits =c(0,.6),n =5,name ="% Potentially Could WFH", label = scales::percent) +geom_sf(data = countyIL, fill=NA, color ="gray")+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_gradientn(colors =c("white", "darkblue"), limits =c(0,.6),n =5,name ="Labor force with\nWFH Feasibility", label = scales::percent)Figure2b
#ggsave(Figure2b, "Figure2b.png", width=3, height=3, units = "in")
Internet Access
Make a subset of dstrata that only includes people who answered the TRANWORK/ did_WFH questions.
#transub<- subset(dstrata,TRANWORK(x=80))
Code
## CINETHH is for access to internet # (either at home, somewhere else, or no access)# from joined dataframe. Not using survey objects or survey packagemapPUMAboth <- 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
Code
## 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 objectsmapPUMAboth <-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_continuous(low ="white", high ="red4", name ="% of Population", label = scales::percent)+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")
Code
# 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_gradientn(colors =c("red4", "white", "forestgreen"), limits =c(-.15,.15),n =7,name ="Chance in Access", label = scales::percent)+labs(title ="Change in Labor Force with Internet at home",subtitle ="2019 to 2021")
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.
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?
Code
### Access to Hi speed internet #### from joined dataframe. Not using survey objects or survey packagemapPUMAboth <- 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
Code
## Access to Hi speed internet using survey objectsmapPUMAboth <-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_gradientn(colors =c("white", "forestgreen"), n =5,name ="Population %", label = scales::percent)+Alea_theme()
Code
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_gradientn(colors =c("forestgreen", "white", "red4"), n =5,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
### 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 objectsmapPUMAboth <-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", # subtitle = "Geograph = Counties with More than 60,000 People" ) +facet_wrap(~YEAR)+scale_fill_gradientn(colors =c( "white", "green4"), n =5,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()+scale_fill_gradientn(colors =c("red4", "white", "forestgreen"), n =7,limits =c(-.2,.2),name ="Percent of Population", label = scales::percent)#+labs(title = "Change in access to High Speed Internet")
ggsave("Figure7b.png", limitsize =FALSE, width =5, height =4, units ="in")#ggsave("Figure7b.eps", limitsize = FALSE,width = 8, height = 4, units = "in")#ggsave("Figure7b.pdf", limitsize = FALSE,width = 8, height = 4, units = "in")# 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 bgplot.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)
# labs(title = "Percentage Point Change in Labor Force that Did WFH")ggsave("Figure4b.eps", limitsize =FALSE, width =8, height =5, units ="in")ggsave("Figure4b.pdf", limitsize =FALSE,width =8, height =5, units ="in")ggsave("Figure4b.png", limitsize =FALSE, width =8, height =5, units ="in")# 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*100, y=reorder(NAME, Prop_2021*100)))+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 that WFH (%)", y ="")#+ scale_x_continuous(label = scales::percent)
ggsave("Figure4a.eps", limitsize =FALSE,width =8, height =5, units ="in")ggsave("Figure4a.pdf", limitsize =FALSE,width =8, height =5, units ="in")ggsave("Figure4a.png", limitsize =FALSE, width =8, height =5, units ="in")
Could work from home
Based on Dingel & Niemen paper. A person’s occupation indicates the likelihood of if they could work from home.
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_gradientn(colors =c("white", "darkblue"), na.value ="gray",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_gradientn(colors =c("white", "darkblue"), na.value ="gray",limits =c(0,.6),name ="Labor Force with\nWFH Feasibility", label = scales::percent)Figure2a
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_continuous(label = scales::percent, low ="white", high ="darkblue")
Mixed Geography
# 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 variablemean_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 representsmutate(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 observationsmean_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 boundariesplot(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.
# 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))
Regression attempt
Regression for 2019 using survey object dstrata2019.
Very simple interpretation: Compared to Business and Finance occupations who Can Work From Home in Rural Counties without any children under 5
Increases likelihood of working from home: Number of Children under 5 increases odds of working from home, living in an “Urban” County (higher pop. density),
Decreases liklihood of woring from home: Lacking access to internet, WFH not feasible & occupation labels within Production and Transportation,
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
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.
