DEM 7093 - Mapping Using the American Community Survey - Change Mapping
Jyoti Nepal, MSW
15 February, 2022
Load libraries
options(Ncores = 6)
library(tidycensus)
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(sf)## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUElibrary(tmap)
library(dplyr)
library(tmap)
library(ggplot2)
library(classInt)
library(patchwork)
library(tmaptools)
library(ggsn)## Loading required package: gridCheck the variable name of my ACS variable
v11_Profile <- load_variables(2011, "acs5/profile", cache = TRUE) #demographic profile tables
v19_Profile <- load_variables(2019, "acs5/profile", cache = TRUE) #demographic profile tables
#Search for variables by using grep()
v11_Profile[grep(x = v11_Profile$label, "VACANT", ignore.case = TRUE), c("name", "label")] #DP04_0003P## # A tibble: 2 x 2
## name label
## <chr> <chr>
## 1 DP04_0003 Estimate!!HOUSING OCCUPANCY!!Vacant housing units
## 2 DP04_0003P Percent!!HOUSING OCCUPANCY!!Vacant housing unitsv19_Profile[grep(x = v19_Profile$label, "VACANT", ignore.case = TRUE), c("name", "label")] #DP04_0003P## # A tibble: 2 x 2
## name label
## <chr> <chr>
## 1 DP04_0003 Estimate!!HOUSING OCCUPANCY!!Total housing units!!Vacant housing u~
## 2 DP04_0003P Percent!!HOUSING OCCUPANCY!!Total housing units!!Vacant housing un~Extract from ACS summary file data profile variables from 2011 and 2019 for Harris County, TX Census Tracts
# Get 2011 ACS data for Harris County
vac11<-get_acs(geography = "tract",
state = "TX",
county = "Harris County",
year = 2011,
variables = "DP04_0003P",
geometry = T,
output = "wide")## Getting data from the 2007-2011 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## Using the ACS Data Profile##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%#rename variables and filter missing cases
vac11 <- vac11 %>%
mutate(pvac = DP04_0003PE,
pvac_er = DP04_0003PM/1.645,
pvac_cv = 100*(pvac_er/pvac)) %>%
filter(complete.cases(pvac), is.finite(pvac_cv)==T) %>%
select(GEOID, pvac, pvac_er, pvac_cv)
head(vac11)## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -95.40302 ymin: 29.8129 xmax: -95.3233 ymax: 29.8869
## Geodetic CRS: NAD83
## GEOID pvac pvac_er pvac_cv geometry
## 1 48201222100 12.5 3.647416 29.17933 POLYGON ((-95.34437 29.8805...
## 2 48201221400 16.3 3.768997 23.12268 POLYGON ((-95.3931 29.84572...
## 3 48201221300 11.0 4.133739 37.57944 POLYGON ((-95.3618 29.84628...
## 4 48201530500 24.0 3.586626 14.94428 POLYGON ((-95.38838 29.8306...
## 5 48201530400 20.5 4.680851 22.83342 POLYGON ((-95.38838 29.8306...
## 6 48201530300 22.4 4.984802 22.25358 POLYGON ((-95.37563 29.8140...# get the 2019 estimates
vac19<-get_acs(geography = "tract",
state = "TX",
county = "Harris County",
year = 2019,
variables = "DP04_0003P",
geometry = T,
output = "wide")## Getting data from the 2015-2019 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## Using the ACS Data Profile##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|======= | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 54%
|
|====================================== | 54%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%#rename variables and filter missing cases
vac19<- vac19%>%
mutate(pvac19 = DP04_0003PE,
pvac19_er = DP04_0003PM/1.645,
pvac19_cv = 100*(pvac19_er/pvac19)) %>%
filter(complete.cases(pvac19), is.finite(pvac19_cv)==T)%>%
select(GEOID, pvac19, pvac19_er, pvac19_cv)
head(vac19)## Simple feature collection with 6 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -95.65378 ymin: 29.70875 xmax: -95.12716 ymax: 30.00429
## Geodetic CRS: NAD83
## GEOID pvac19 pvac19_er pvac19_cv geometry
## 1 48201311900 10.5 3.890578 37.05312 MULTIPOLYGON (((-95.33319 2...
## 2 48201450200 7.2 3.708207 51.50287 MULTIPOLYGON (((-95.59029 2...
## 3 48201450400 27.9 4.924012 17.64879 MULTIPOLYGON (((-95.62377 2...
## 4 48201554502 1.5 1.094225 72.94833 MULTIPOLYGON (((-95.65325 2...
## 5 48201550603 6.4 2.857143 44.64286 MULTIPOLYGON (((-95.48031 2...
## 6 48201252200 15.4 3.100304 20.13184 MULTIPOLYGON (((-95.18517 2...#Merge 2011 and 2019 data
st_geometry(vac19) <- NULL
mdat <- left_join(vac11, vac19, by = c("GEOID"="GEOID"))
head(mdat)## Simple feature collection with 6 features and 7 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -95.40302 ymin: 29.8129 xmax: -95.3233 ymax: 29.8869
## Geodetic CRS: NAD83
## GEOID pvac pvac_er pvac_cv pvac19 pvac19_er pvac19_cv
## 1 48201222100 12.5 3.647416 29.17933 9.3 3.221884 34.64392
## 2 48201221400 16.3 3.768997 23.12268 13.4 3.465046 25.85855
## 3 48201221300 11.0 4.133739 37.57944 3.2 2.188450 68.38906
## 4 48201530500 24.0 3.586626 14.94428 8.3 2.492401 30.02893
## 5 48201530400 20.5 4.680851 22.83342 10.6 2.978723 28.10116
## 6 48201530300 22.4 4.984802 22.25358 21.9 5.167173 23.59440
## geometry
## 1 POLYGON ((-95.34437 29.8805...
## 2 POLYGON ((-95.3931 29.84572...
## 3 POLYGON ((-95.3618 29.84628...
## 4 POLYGON ((-95.38838 29.8306...
## 5 POLYGON ((-95.38838 29.8306...
## 6 POLYGON ((-95.37563 29.8140...Calculate the differences between the two years, including the test for significance of the differences
acstest<-function(names,geoid, est1, err1, est2, err2, alpha, yr1, yr2, span){
se1<-err1/qnorm(.90)
se2<-err2/qnorm(.90)
yrs1<-seq(yr1, to=yr1-span)
yrs2<-seq(yr2, to=yr2-span)
C<-mean(yrs2%in%yrs1)
diff <- (est1-est2)
test <- (est1-est2) / (sqrt(1-C)*sqrt((se1^2+se2^2))) #this is a t-test
crit <- qnorm(1-alpha/2)
pval <- 1-pnorm(abs(test)) #2*pnorm(abs(test),lower.tail=F)
result <- NULL
result[pval > alpha] <- "insignificant change"
result[pval < alpha & test < 0] <- "significant increase"
result[pval < alpha & test > 0] <- "significant decrease"
data.frame(name=names,geoid=geoid, est1=est1, est2=est2, se1=se1, se2=se2,difference=diff, test=test, result=result, pval=pval)
}Run the function above to get comparison of percent housing units vacant and test for significant changes
diff1119 <- acstest(names = mdat$GEOID,
geoid = mdat$GEOID,
est1 = mdat$pvac,
est2 = mdat$pvac19,
err1 = mdat$pvac_er,
err2 = mdat$pvac19_er,
alpha = .1,
yr1 = 2011, yr2=2019,
span = 5)
head(diff1119)## name geoid est1 est2 se1 se2 difference test
## 1 48201222100 48201222100 12.5 9.3 2.846094 2.514050 3.2 0.8426686
## 2 48201221400 48201221400 16.3 13.4 2.940964 2.703789 2.9 0.7259137
## 3 48201221300 48201221300 11.0 3.2 3.225573 1.707656 7.8 2.1371536
## 4 48201530500 48201530500 24.0 8.3 2.798659 1.944831 15.7 4.6067287
## 5 48201530400 48201530400 20.5 10.6 3.652487 2.324310 9.9 2.2867289
## 6 48201530300 48201530300 22.4 21.9 3.889662 4.031967 0.5 0.0892484
## result pval
## 1 insignificant change 1.997069e-01
## 2 insignificant change 2.339458e-01
## 3 significant decrease 1.629275e-02
## 4 significant decrease 2.045265e-06
## 5 significant decrease 1.110583e-02
## 6 insignificant change 4.644422e-01table(diff1119$result)##
## insignificant change significant decrease significant increase
## 393 257 101Create a map using quantile breaks of percent housing units vacant for each year
#create a map using quantile breaks of percent housing units vacant for each year
acs_merge<-left_join(mdat, diff1119, by=c("GEOID"="geoid"))
tmap_mode("plot")## tmap mode set to plottingp1<-tm_shape(acs_merge)+
tm_polygons(c("pvac"), title=c("% Housing Units Vacant 2011"), palette="Blues", style="quantile", n=5)+
tm_scale_bar(10)+
tm_layout(title="Houston Estimated Percent Housing Units Vacant 2011", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_compass() +
tm_format("World",
legend.position = c("left", "bottom"),
legend.frame = FALSE,
legend.title.size = 1.5,
legend.text.size = 1,
main.title.position =c("center"))
p2<-tm_shape(acs_merge)+
tm_polygons(c("pvac19"), title=c("% Housing Units Vacant 2019"), palette="Blues", style="quantile", n=5)+
tm_layout(title="Houston Estimated Percent Housing Units Vacant 2019", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_scale_bar(10)+
tm_compass() +
tm_format("World",
legend.position = c("left", "bottom"),
legend.frame = FALSE,
legend.title.size = 1.5,
legend.text.size = 1,
main.title.position =c("center"))
tmap_arrange(p1,p2)## Warning: First scale_bar breaks value should be 0.
## Warning: First scale_bar breaks value should be 0.
## Warning: First scale_bar breaks value should be 0.
## Warning: First scale_bar breaks value should be 0.
Create a map showing the differences by census tract
p3 <- tm_shape(acs_merge)+
tm_polygons(c("result"), title=c("Change from 2011 to 2019"), palette = "Set2")+
tm_layout(title="Houston Vacant Housing Units Estimate Changes", title.size =1.5, legend.frame = TRUE, title.position = c('center', 'top'))+
tm_scale_bar(10)+
tm_compass() +
tm_format("World",
legend.position = c("left", "bottom"),
legend.frame = FALSE,
legend.title.size = 1.3,
legend.text.size = 1,
main.title.position =c("center"))
tmap_arrange(p3)## Warning: First scale_bar breaks value should be 0.
## Warning: First scale_bar breaks value should be 0.