Homework 4
Asiya
02/14/2021
#Load Libraries
library(tidycensus)
library(tidyverse)## ── Attaching packages ──────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.3 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ─────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(sf)## Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1library(tmap)
library(ggplot2)
library(classInt)
library(patchwork)
library(dplyr)
library(tmaptools)
library(ggsn)## Loading required package: gridv11_Profile <- load_variables(2011, "acs5/profile", cache = TRUE) #demographic profile tables
v19_Profile <- load_variables(2019, "acs5/profile", cache = TRUE) #demographic
head(v11_Profile [2])## # A tibble: 6 x 1
## label
## <chr>
## 1 Estimate!!HOUSEHOLDS BY TYPE!!Total households
## 2 Percent!!HOUSEHOLDS BY TYPE!!Total households
## 3 Estimate!!HOUSEHOLDS BY TYPE!!Family households (families)
## 4 Percent!!HOUSEHOLDS BY TYPE!!Family households (families)
## 5 Estimate!!HOUSEHOLDS BY TYPE!!Family households (families)!!With own children…
## 6 Percent!!HOUSEHOLDS BY TYPE!!Family households (families)!!With own children …#Search for variables by using grep()
v11_Profile[grep(x = v11_Profile$label, "HOUSING", ignore.case = TRUE), c("name", "label")]## # A tibble: 54 x 2
## name label
## <chr> <chr>
## 1 DP03_0038 Estimate!!INDUSTRY!!Transportation and warehousing, and utilities
## 2 DP03_0038P Percent!!INDUSTRY!!Transportation and warehousing, and utilities
## 3 DP04_0001 Estimate!!HOUSING OCCUPANCY!!Total housing units
## 4 DP04_0001P Percent!!HOUSING OCCUPANCY!!Total housing units
## 5 DP04_0002 Estimate!!HOUSING OCCUPANCY!!Occupied housing units
## 6 DP04_0002P Percent!!HOUSING OCCUPANCY!!Occupied housing units
## 7 DP04_0003 Estimate!!HOUSING OCCUPANCY!!Vacant housing units
## 8 DP04_0003P Percent!!HOUSING OCCUPANCY!!Vacant housing units
## 9 DP04_0004 Estimate!!HOUSING OCCUPANCY!!Homeowner vacancy rate
## 10 DP04_0004P Percent!!HOUSING OCCUPANCY!!Homeowner vacancy rate
## # … with 44 more rowsv19_Profile[grep(x = v19_Profile$label, "HOUSING", ignore.case = TRUE), c("name", "label")]## # A tibble: 232 x 2
## name label
## <chr> <chr>
## 1 DP03_0038 Estimate!!INDUSTRY!!Civilian employed population 16 years and ove…
## 2 DP03_0038P Percent!!INDUSTRY!!Civilian employed population 16 years and over…
## 3 DP04_0001 Estimate!!HOUSING OCCUPANCY!!Total housing units
## 4 DP04_0001P Percent!!HOUSING OCCUPANCY!!Total housing units
## 5 DP04_0002 Estimate!!HOUSING OCCUPANCY!!Total housing units!!Occupied housin…
## 6 DP04_0002P Percent!!HOUSING OCCUPANCY!!Total housing units!!Occupied housing…
## 7 DP04_0003 Estimate!!HOUSING OCCUPANCY!!Total housing units!!Vacant housing …
## 8 DP04_0003P Percent!!HOUSING OCCUPANCY!!Total housing units!!Vacant housing u…
## 9 DP04_0004 Estimate!!HOUSING OCCUPANCY!!Total housing units!!Homeowner vacan…
## 10 DP04_0004P Percent!!HOUSING OCCUPANCY!!Total housing units!!Homeowner vacanc…
## # … with 222 more rowsExtract from ACS summary file data profile variables from 2011 and 2019 for Dallas County, TX Census Tracts
perc_vac11<-get_acs(geography = "tract",
state="TX",
county = "Dallas",
year = 2011,
variables="DP04_0003P" , # change the variable
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## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.#rename variables and filter missing cases
perc_vac11<-perc_vac11%>%
mutate(pvac11=DP04_0003PE, # Change the code
pvac11_er=DP04_0003PM/1.65, # Change the code
pvac11_cv=100* (pvac11_er/pvac11)) %>%
filter(complete.cases(pvac11), is.finite(pvac11_cv)==T)
perc_vac19<-get_acs(geography = "tract",
state="TX",
county = "Dallas",
year = 2019,
variables="DP04_0003P" , # change the variable
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 Profileperc_vac19<-perc_vac19%>%
mutate(pvac19=DP04_0003PE, # Change the code
pvac19_er=DP04_0003PM/1.65, # Change the code
pvac19_cv=100* (pvac19_er/pvac19)) %>%
filter(complete.cases(pvac19), is.finite(pvac19_cv)==T)
#merge the two years worth of data
mdat<-tigris::geo_join(perc_vac11, as.data.frame(perc_vac19), by_sp="GEOID", by_df="GEOID")## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.Create a map using quantile breaks of this variable for each year
#2011
p1<-tm_shape(perc_vac11)+
tm_polygons(c("pvac11"), title=c("Percent Vacant Housing 2011"), palette="Blues", style="quantile", n=5)+
tm_scale_bar()+
tm_layout(title="Dallas Percent Vacant Housing Rate Estimates 2011", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_compass()
p2<-tm_shape(perc_vac11)+
tm_polygons(c("pvac11_cv"), title=c("Percent Vacant Housing Rate CV 2011"), palette="Blues", style="quantile", n=5)+
tm_layout(title="Dallas Percent Vacant Housing Rate CV 2011", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_scale_bar()+
tm_compass()
tmap_arrange(p1, p2)
#2019
p1<-tm_shape(perc_vac19)+
tm_polygons(c("pvac19"), title=c("Percent Vacant Housing 2019"), palette="Blues", style="quantile", n=5)+
tm_scale_bar()+
tm_layout(title="Dallas Percent Vacant Housing Rate Estimates 2019", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_compass()
p2<-tm_shape(perc_vac19)+
tm_polygons(c("pvac19_cv"), title=c("Percent Vacant Housing Rate CV 2019"), palette="Blues", style="quantile", n=5)+
tm_layout(title="Dallas Percent Vacant Housing Rate CV 2019", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_scale_bar()+
tm_compass()
tmap_arrange(p1, p2)
Change map between two time points
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)))
crit<-qnorm(1-alpha/2)
pval<-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)
}
diff<-acstest(names = mdat$NAME.x, geoid = mdat$GEOID, est1 = mdat$pvac11, est2 = mdat$pvac19, err1 = mdat$pvac11_er, err2=mdat$pvac19_er,alpha = .1, yr1 = 2011, yr2=2019, span = 5)Create a map showing the differences by census tract
acs_merge<-left_join(mdat, diff, by=c("GEOID"="geoid"))
tmap_mode("plot")## tmap mode set to plottingp1<-tm_shape(acs_merge)+
tm_polygons(c("pvac11"), title=c("Percent Vacant Housing 2011"), palette="Blues", style="quantile", n=5)+
tm_scale_bar()+
tm_layout(title="Dallas Percent Vacant Housing Rate Estimates 2011", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_compass()
p2<-tm_shape(acs_merge)+
tm_polygons(c("pvac19"), title=c("Percent Vacant Housing 2019"), palette="Blues", style="quantile", n=5)+
tm_layout(title="Dallas Percent Vacant Housing Rate Estimates 2019", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_scale_bar()+
tm_compass()
p3 <- tm_shape(acs_merge)+
tm_polygons(c("result"), title=c("Changes"), palette = "Dark2")+
tm_layout(title="Dallas Percent Vacant Housing Rate Estimates Changes", title.size =1.5, legend.frame = TRUE, title.position = c('right', 'top'))+
tm_scale_bar()+
tm_compass()
tmap_arrange(p1, p2, p3)