Tidycensus Markdown
Jingmiao Fei
2024-09-05
Setup
Load packages
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(tidycensus)
library(sf)
library(knitr)
library(kableExtra)
library(rmarkdown)
library(pander)Access the Census API
census_api_key("e6ed37e145b901c8c418a193d5bc174948045488", overwrite = TRUE,install=TRUE)Load census data dictionaries
acs_variable_list.2020 <- load_variables(2020, #year
"acs5", #five year American Community Service estimates
cache = TRUE)
acs_variable_list.2016 <- load_variables(2016, #year
"acs5", #five year ACS estimates
cache = TRUE)Downloading Data from Tidycensus
Create a vector of census variables
acs_vars <- c("B01001_001E", # ACS total Pop estimate
"B25002_001E", # Estimate of total housing units
"B25002_003E", # Number of vacant housing units
"B19013_001E", # Median HH Income ($)
"B02001_002E", # People describing themselves as "white alone"
"B06009_006E") # Total graduate or professional degreeCall the Census API to get tract level data for 2020 for all of Philadelphia
acsTractsPHL.2020 <- get_acs(geography = "tract",
year = 2020,
variables = acs_vars,
geometry = FALSE,
state = "PA",
county = "Philadelphia",
output = "wide") Wrangling Data with dplyr
Mutating, selecting and renaming variables
acsTractsPHL.2020 <- acsTractsPHL.2020 %>%
dplyr::select (GEOID, NAME, all_of(acs_vars))
acsTractsPHL.2020 <- acsTractsPHL.2020 %>%
rename (total_pop.2020 = B01001_001E,
total_HU.2020 = B25002_001E,
total_vacant.2020 = B25002_003E,
med_HH_Income.2020 = B19013_001E,
total_White.2020 = B02001_002E,
total_GradDeg.2020 = B06009_006E)
acsTractsPHL.2020 <- acsTractsPHL.2020 %>%
mutate(vacancyPct.2020 = total_vacant.2020/total_HU.2020,
pctWhite.2020 = total_White.2020/total_pop.2020)acsTractsPHL.2016 <- get_acs(geography = "tract",
year = 2016,
variables = acs_vars,
geometry = FALSE,
state = "PA",
county = "Philadelphia",
output = "wide") %>%
dplyr::select (GEOID, NAME, all_of(acs_vars)) %>%
rename (total_pop.2016 = B01001_001E,
total_HU.2016 = B25002_001E,
total_vacant.2016 = B25002_003E,
med_HH_Income.2016 = B19013_001E,
total_White.2016 = B02001_002E,
total_GradDeg.2016 = B06009_006E) %>%
mutate(vacancyPct.2016 = total_vacant.2016/total_HU.2016,
pctWhite.2016 = total_White.2016/total_pop.2016)Joining data
allACS <- left_join(acsTractsPHL.2016, acsTractsPHL.2020,
by= c("GEOID"))Doing column math using mutate
allACS <- allACS %>%
mutate(change_med_HH_Income = med_HH_Income.2020 - (med_HH_Income.2016 * 1.08),
change_Grad_Degree_Pct = (total_GradDeg.2020/total_pop.2020)-(total_GradDeg.2016/total_pop.2016))Summarizing Census Data
Exploring central tendancies
mean(allACS$change_med_HH_Income, na.rm = TRUE)## [1] 6012.472median(allACS$change_med_HH_Income, na.rm = TRUE)## [1] 4873.76Exploring distributions
hist(allACS$change_med_HH_Income)
ggplot(allACS)+
geom_histogram(aes(change_med_HH_Income), binwidth = 5000)+
labs(
title = "Change in Philadelphia HH median income by tract, 2016-2020",
caption = "Data: US Census Bureau, ACS 5-year estimates",
x="Change in Med HH Income (2020 dollars)",
y="Number of tracts")
Making a summary table
summaryTable <- allACS %>%
summarize(mean_change_HH_Income = mean(change_med_HH_Income, na.rm = TRUE),
med_change_HH_Income = median(change_med_HH_Income, na.rm = TRUE))Comparing geographies
myTracts <- c("42101023500",
"42101023600",
"42101023700",
"42101025300",
"42101025400",
"42101025500",
"42101025600",
"42101038800")
allACS <- allACS %>%
mutate(mtAiry = ifelse(GEOID %in% myTracts, "MT AIRY", "REST OF PHILADELPHIA"))summaryTable2 <- allACS %>%
group_by(mtAiry) %>%
summarize(mean_change_HH_Income = mean(change_med_HH_Income, na.rm = TRUE),
med_change_HH_Income = median(change_med_HH_Income, na.rm = TRUE))Graphic comparisons Using ggplot2
ggplot(allACS)+
geom_histogram(aes(change_med_HH_Income),
binwidth = 5000)+
labs(
title = "Change in Philadelphia HH median income by tract, 2016-2020",
caption = "Data: US Census Bureau, ACS 5-year estimates",
x="Change in Med HH Income (2020 dollars)",
y="Number of tracts")+
facet_wrap(~mtAiry, scales = "free")+
theme_minimal() 
ggplot(allACS)+
geom_point(aes(x =med_HH_Income.2016 * 1.08,
y = med_HH_Income.2020,
color = mtAiry))+
geom_abline(intercept = 0, slope = 1)+
labs(
title = "2020 Median HH Income as a Function of 2016 Median HH Income",
subtitle = "All figures in 2020 dollars",
caption = "Data: US Census Bureau, ACS 5-year estimates",
x="Med HH Income 2016 ($)",
y="Med HH Income 2020 ($)")
ggplot(allACS)+
geom_point(aes(x = 100* pctWhite.2020,
y = med_HH_Income.2020,
color = mtAiry))+
geom_smooth(aes(x = 100* pctWhite.2020,
y = med_HH_Income.2020),
method = "lm", se = FALSE)+
labs(
title = "2020 Median HH Income as a Function of Pct White",
subtitle = "All figures in 2020 dollars",
caption = "Data: US Census Bureau, ACS 5-year estimates",
x="Pct. Residents Identifying as 'White Only'",
y="Med HH Income 2020 ($)")
Spatial Data and Tidycensus
ggplot()+
geom_sf(data = acsTractsPHL.2020.sf, aes(fill = pctWhite.2020),
color = "transparent")+
geom_sf(data = acsTractsPHL.2020.sf %>%
filter(mtAiry == "MT AIRY") %>%
st_union(),
color = "white",
fill = "transparent")+
labs(
title = "Percentage of those identifying as 'white only' by tract",
subtitle = "",
caption = "Data: US Census Bureau, ACS 5-year estimates")
Assignment
Question 2a - Change in Raw Vacant Housing Units in Philadelphia (2016-2020)
ggplot(allACS)+
geom_point(aes(x = total_vacant.2016,
y = total_vacant.2020,
color = mtAiry))+
geom_smooth(aes(x = total_vacant.2016,
y = total_vacant.2020),
method = "lm", se = FALSE)+
geom_abline(intercept = 0, slope = 1)+
labs(
title = "Change in Raw Vacant Housing Units in Philadelphia (2016-2020)",
color = "Area",
caption = "Data: US Census Bureau, ACS 5-year estimates",
x="2016 Toal Vacant Housing Units",
y="2020 Toal Vacant Housing Units")
The plot above shows the number of vacant housing units in Philadelphia in 2016 and 2020, where the orange and blue points indicate Mt. Airy tracts and the rest of Philadelphia respectively. The blue line represents the linear fit and the black line means no change over the two years.
Question 2b - percentage of vacant housing units in Philadelphia (2016)
ggplot(acsTractsPHL.2016.sf)+
geom_sf(data = acsTractsPHL.2016.sf, aes(fill = vacancyPct.2016),
color = "transparent")+
geom_sf(data = acsTractsPHL.2016.sf %>%
filter(mtAiry == "MT AIRY") %>%
st_union(),
color = "white",
fill = "transparent")+
labs(
title = "Percentage of Vacant Housing Units in Philadelphia (2016)",
fill = "Housing Vacancy Percentage",
caption = "Data: US Census Bureau, ACS 5-year estimates"
)
The map above shows the vacant housing units percentage in Philadelphia in 2016 by tract, with Mt. Airy area highlighted in white. It shows that Mt. Airy area has a relatively low vacancy percentage compared with the rest of Philadelphia, which indicates a strong demand for housing and robust housing market.
Question 2c - Mean Number of Vacant Housing Units per Tract in Philadelphia (2020)
summaryTableVhu <- acsTractsPHL.2020.sf %>%
st_drop_geometry() %>%
group_by(mtAiry) %>%
summarize(mean_vacant_housing_units = mean(total_vacant.2020, na.rm = TRUE))%>%
rename(Area = mtAiry)%>%
rename("mean vacant housing units" = mean_vacant_housing_units)
summaryTableVhu %>%
kbl(caption = "Mean Number of Vacant Housing Units per tract in Philadelphia (2020)") %>%
kable_styling(full_width = FALSE)| Area | mean vacant housing units |
|---|---|
| MT AIRY | 156.6250 |
| REST OF PHILADELPHIA | 186.8575 |