library(tidycensus)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.2.0 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ 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.8.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(tmap)
vachouse10<-get_acs(geography = "tract",
state="TX",
county = "Harris County",
year = 2010,
variables="DP04_0003P" ,
geometry = T,
output = "wide")
## Getting data from the 2006-2010 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
vachouse10 <- vachouse10 %>%
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(vachouse10)
tm_shape(vachouse10)+
tm_polygons(c("pvac"),
title=c("% Vacant Houses"),
palette="Reds",
style="quantile",
n=5)+
tm_scale_bar()+
tm_layout(title="Harris County Vacant Housing Estimates - Quantile Breaks",
title.size = 1.5,
legend.frame = TRUE,
title.position = c('right','top'))+
tm_compass()+
tm_format("World",
legend.position = c("left", "bottom"),
main.title.position=c("center"))
tm_shape(vachouse10)+
tm_polygons(c("pvac"),
title=c("% Vacant Housing"),
palette="Reds",
style="jenks",
n=5)+
tm_scale_bar()+
tm_layout(title="Harris County Vacant Housing Estimates - Jenks Breaks",
title.size=1.5,
legend.frame = TRUE,
title.position = c('right','top'))+
tm_compass()+
tm_format("World",
legend.position = c("left", "bottom"),
main.title.position=c("center"))
p1<- tm_shape(vachouse10)+
tm_polygons(c("pvac"),
title=c("% Vacant Housing"),
palette="Reds",
style="quantile",
n=5)+
tm_scale_bar()+
tm_layout(title="Harris County Vacant Housing Estimates",
title.size = 1.5,
legend.frame = TRUE,
title.position = c('right','top'))+
tm_compass()+
tm_format("World",
legend.position=c("left", "bottom"),
main.title.position=c("center"))
p2<-tm_shape(vachouse10)+
tm_polygons(c("pvac_cv"),
title=c("CV Vacant Houses"),
palette="Reds",
style="quantile",
n=5)+
tm_layout(title="Harris County Housing Estimate CV",
title.size=1.5,
legend.frame = TRUE,
title.position = c('right','top'))+
tm_scale_bar()+
tm_compass()+
tm_format("World", legend.position = c("left", "bottom"),
main.title.position=c("center"))
tmap_arrange(p1, p2)
plot(vachouse10$pvac, vachouse10$pvac_cv, main = "Error in Estimates vs Estimate Size")
vachouse19<-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
vachouse19 <- vachouse19 %>%
mutate(pvac19 = DP04_0003PE,
pvac19_er = DP04_0003PM/1.645,
pvac19_cv = 100*(pvac19_er/pvac19)) %>%
filter(complete.cases(pvac19))%>%
select(GEOID, pvac19, pvac19_er, pvac19_cv)
head(vachouse19)
st_geometry(vachouse19)<-NULL
mdat<- left_join(vachouse10, vachouse19, by=c("GEOID"="GEOID"))
head(mdat)
## Simple feature collection with 6 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -95.37805 ymin: 29.74486 xmax: -95.29166 ymax: 29.8197
## Geodetic CRS: NAD83
## GEOID pvac pvac_er pvac_cv pvac19 pvac19_er pvac19_cv
## 1 48201100000 32.8 7.0 21.34146 21.6 3.100304 14.35326
## 2 48201210900 5.5 6.4 116.36364 21.1 5.167173 24.48897
## 3 48201211000 26.0 10.6 40.76923 18.3 4.559271 24.91405
## 4 48201211600 9.9 5.9 59.59596 12.5 2.674772 21.39818
## 5 48201211900 15.8 6.0 37.97468 7.8 2.978723 38.18876
## 6 48201220200 11.3 8.5 75.22124 17.4 3.829787 22.01027
## geometry
## 1 MULTIPOLYGON (((-95.37348 2...
## 2 MULTIPOLYGON (((-95.32851 2...
## 3 MULTIPOLYGON (((-95.324 29....
## 4 MULTIPOLYGON (((-95.3159 29...
## 5 MULTIPOLYGON (((-95.3032 29...
## 6 MULTIPOLYGON (((-95.35169 2...
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<-1-pnorm(abs(test))
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)
}
diff1019 <- acstest(names = mdat$GEOID,
geoid = mdat$GEOID,
est1 = mdat$pvac,
est2 = mdat$pvac19,
err1 = mdat$pvac_er,
err2 = mdat$pvac19_er,
alpha = .1,
yr1 = 2010, yr2 = 2019,
span = 5)
head(diff1019)
## name geoid est1 est2 se1 se2 difference test
## 1 48201100000 48201100000 32.8 21.6 5.462129 2.419180 11.2 1.8748273
## 2 48201210900 48201210900 5.5 21.1 4.993947 4.031967 -15.6 -2.4305005
## 3 48201211000 48201211000 26.0 18.3 8.271224 3.557618 7.7 0.8551872
## 4 48201211600 48201211600 9.9 12.5 4.603794 2.087136 -2.6 -0.5143620
## 5 48201211900 48201211900 15.8 7.8 4.681825 2.324310 8.0 1.5305042
## 6 48201220200 48201220200 11.3 17.4 6.632585 2.988399 -6.1 -0.8385189
## result pval
## 1 significant decrease 0.030408244
## 2 significant increase 0.007538992
## 3 insignificant change 0.196223754
## 4 insignificant change 0.303499453
## 5 significant decrease 0.062945989
## 6 insignificant change 0.200869655
table(diff1019$result)
##
## insignificant change significant decrease significant increase
## 516 172 67
acs_merge <- left_join(mdat, diff1019, by=c("GEOID" = "geoid"))
tmap_mode("plot")
## tmap mode set to plotting
p1<- tm_shape(acs_merge)+
tm_polygons(c("pvac"),
title=c("% Vacant Housing"),
palette="Reds",
style="quantile",
n=5)+
tm_scale_bar()+
tm_layout(title = "Harris County Vacant Housing Estimates 2010",
title.size = 1.5,
legend.frame = TRUE,
title.position = c('right', 'top'))+
tm_compass()+
tm_format("World",
legend.position = c("left", "bottom"),
main.title.position = c('right', 'top'))
p2<- tm_shape(acs_merge)+
tm_polygons(c("pvac19"),
title=c("% Vacant Housing"),
palette="Reds",
style="quantile",
n=5)+
tm_layout(title="Harris County Vacant Housing Estimates 2019",
title.size = 1.5,
legend.frame = TRUE,
title.position = c('right', 'top'))+
tm_scale_bar()+
tm_compass()+
tm_format("World",
legend.position = c("left", "bottom"),
main.title.position = c("center"))
p3 <- tm_shape(acs_merge)+
tm_polygons(c("result"),
title=c("Changes"),
palette = "Set2")+
tm_layout(title = "Harris County Vacant Housing Estimate Changes",
title.size = 1.5,
legend.frame = TRUE,
title.position = c('right', 'top'))+
tm_scale_bar()+
tm_compass()
tmap_arrange(p1, p2, p3, nrow = 2)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(acs_merge)+
tm_polygons("result",
alpha=.7,
title=c("changes"),
palette="Set2")+
tm_layout(title = "Harris County Vacant Housing Estimate Changes",
title.size = 1.5)+
tm_scale_bar()