DEM GIS_Homework2
Monica Cruz
June 10, 2018
install.packages(“censusapi”) install.packages(“tidycensus”) install.packages(“sf”) devtools::install_github(“oswaldosantos/ggsn”)
library(censusapi)##
## Attaching package: 'censusapi'## The following object is masked from 'package:methods':
##
## getFunctionlibrary(tidycensus)
library(ggplot2)
library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggsn)## Loading required package: gridlibrary(mapview)
library(RColorBrewer)
census_api_key(key ="f71931186e975910aa774fdd6aff604b31e9fe71")## To install your API key for use in future sessions, run this function with `install = TRUE`.Look at available ACS variables
library(tidycensus); library(tidyverse); library(sf)## -- Attaching packages ---------------------------------------------------------- tidyverse 1.2.1 --## v tibble 1.4.2 v purrr 0.2.4
## v tidyr 0.8.1 v stringr 1.3.1
## v readr 1.1.1 v forcats 0.3.0## -- Conflicts ------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()v15_Profile <- load_variables(2015, "acs5/profile", cache = TRUE) #demographic profile tables
v15_tables <- load_variables(2015 , "acs5", cache = TRUE) #all tables
#v10_sf1_tables <- load_variables(2010 , "sf1", cache = TRUE) #all tables for 2010 SF1
#v10_sf1_tables <- load_variables(2000 , "sf3", cache = TRUE) #all tables for 2000 SF#
View(v15_Profile)
#Search for variables by
v15_Profile[grep(x = v15_Profile$label, "Poverty"), c("name", "label")]## # A tibble: 0 x 2
## # ... with 2 variables: name <chr>, label <chr>Extract from ACS summary file data profile variables from 2010 and 2015 for Webb County, TX Census Tracts
webb2010_acs<-get_acs(geography = "tract", state="TX", county = c("Webb"), year = 2010,
variables=c("B17001_002") ,
summary_var = "B17001_001",
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)`.webb2015_acs<-get_acs(geography = "tract", state="TX", county = c("Webb"), year = 2015,
variables=c("B17001_002") ,
summary_var = "B17001_001",
geometry = T, output = "wide")## Getting data from the 2011-2015 5-year ACS
## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.#create a county FIPS code - 5 digit
webb2010_acs$county<-substr(webb2010_acs$GEOID, 1, 5)
webb2015_acs$county<-substr(webb2015_acs$GEOID, 1, 5)compare poverty rates over time
Here we take the poverty rate in tracts derived from the 2006-2010 ACS and compare it to the estimate from the 2011-2015 ACS
#v10_Profile <- load_variables(2010 ,dataset = "acs5") #demographic profile
options(tigris_use_cache = TRUE)
#2010 ACS
webb2010_acs<-get_acs(geography = "tract", state="TX", county = "Webb", year = 2010, variables="B17001_002",summary_var = "B17001_001", geometry = TRUE, output = "wide")## Getting data from the 2006-2010 5-year ACS#2015 ACS
webb2015_acs<-get_acs(geography = "tract", state="TX", county = "Webb", year = 2015, variables="B17001_002",summary_var = "B17001_001",geometry = TRUE, output = "wide")## Getting data from the 2011-2015 5-year ACS#create variables with nice names
webb2010_acs<-webb2010_acs%>%
mutate(pov10=B17001_002E/summary_est, poverr10=moe_prop(B17001_002E, denom = summary_est, moe_num = B17001_002M, moe_denom = summary_moe))%>%
select(GEOID, NAME, pov10, poverr10)
#create variables with nice names
webb2015_acs<-webb2015_acs%>%
mutate(pov15=B17001_002E/summary_est, poverr15=moe_prop(B17001_002E, denom = summary_est, moe_num = B17001_002M, moe_denom = summary_moe))%>%
select(GEOID, NAME, pov15, poverr15)
#merge the two years worth of data
mdat<-left_join(webb2010_acs, as.data.frame(webb2015_acs), by="GEOID", st_join=FALSE)
head(mdat)## Simple feature collection with 6 features and 7 fields
## Active geometry column: geometry.x
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -99.50279 ymin: 27.47402 xmax: -99.44769 ymax: 27.53058
## epsg (SRID): 4269
## proj4string: +proj=longlat +datum=NAD83 +no_defs
## GEOID NAME.x pov10 poverr10
## 1 48479000106 Census Tract 1.06, Webb County, Texas 0.3328264 0.1250222
## 2 48479000107 Census Tract 1.07, Webb County, Texas 0.6368375 0.1591927
## 3 48479000300 Census Tract 3, Webb County, Texas 0.3618396 0.1179660
## 4 48479000903 Census Tract 9.03, Webb County, Texas 0.5636975 0.1155601
## 5 48479001003 Census Tract 10.03, Webb County, Texas 0.3406822 0.1925951
## 6 48479001103 Census Tract 11.03, Webb County, Texas 0.3965812 0.1662029
## NAME.y pov15 poverr15
## 1 Census Tract 1.06, Webb County, Texas 0.5017866 0.12427429
## 2 Census Tract 1.07, Webb County, Texas 0.4048738 0.08696312
## 3 Census Tract 3, Webb County, Texas 0.7001698 0.13228184
## 4 Census Tract 9.03, Webb County, Texas 0.5237918 0.14829979
## 5 Census Tract 10.03, Webb County, Texas 0.3663872 0.15673194
## 6 Census Tract 11.03, Webb County, Texas 0.3020908 0.14394579
## geometry.x geometry.y
## 1 MULTIPOLYGON (((-99.47155 2... MULTIPOLYGON (((-99.4717 27...
## 2 MULTIPOLYGON (((-99.46845 2... MULTIPOLYGON (((-99.47458 2...
## 3 MULTIPOLYGON (((-99.49455 2... MULTIPOLYGON (((-99.50269 2...
## 4 MULTIPOLYGON (((-99.48655 2... MULTIPOLYGON (((-99.49271 2...
## 5 MULTIPOLYGON (((-99.45719 2... MULTIPOLYGON (((-99.46422 2...
## 6 MULTIPOLYGON (((-99.46285 2... MULTIPOLYGON (((-99.46799 2...
## geometry
## 1 GEOMETRYCOLLECTION EMPTY
## 2 GEOMETRYCOLLECTION EMPTY
## 3 GEOMETRYCOLLECTION EMPTY
## 4 GEOMETRYCOLLECTION EMPTY
## 5 GEOMETRYCOLLECTION EMPTY
## 6 GEOMETRYCOLLECTION EMPTYHere I create a function that implements the testing procedure used by the Census for comparing estimates across year
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)
}Here I use the function I just made to do the comparisons
diff1015<-acstest(names = mdat$NAME.x, geoid = mdat$GEOID, est1 = mdat$pov10, est2 = mdat$pov15, err1 = mdat$poverr10, err2=mdat$poverr15, alpha = .1, yr1 = 2010, yr2=2015, span = 5)
head(diff1015)## name geoid est1 est2
## 1 Census Tract 1.06, Webb County, Texas 48479000106 0.3328264 0.5017866
## 2 Census Tract 1.07, Webb County, Texas 48479000107 0.6368375 0.4048738
## 3 Census Tract 3, Webb County, Texas 48479000300 0.3618396 0.7001698
## 4 Census Tract 9.03, Webb County, Texas 48479000903 0.5636975 0.5237918
## 5 Census Tract 10.03, Webb County, Texas 48479001003 0.3406822 0.3663872
## 6 Census Tract 11.03, Webb County, Texas 48479001103 0.3965812 0.3020908
## se1 se2 difference test result
## 1 0.09755536 0.09697174 -0.16896026 -1.3455757 insignificant change
## 2 0.12421875 0.06785768 0.23196378 1.7952152 significant decrease
## 3 0.09204934 0.10322007 -0.33833017 -2.6798018 significant increase
## 4 0.09017200 0.11571894 0.03990566 0.2979787 insignificant change
## 5 0.15028275 0.12229858 -0.02570497 -0.1453282 insignificant change
## 6 0.12968880 0.11232150 0.09449035 0.6033146 insignificant change
## pval
## 1 0.178439397
## 2 0.072619418
## 3 0.007366577
## 4 0.765719455
## 5 0.884451747
## 6 0.546299430options(scipen=999)
acs_merge<-left_join(webb2010_acs, diff1015, by=c("GEOID"="geoid"))## Warning: Column `GEOID`/`geoid` joining character vector and factor,
## coercing into character vectorbrks<-quantile(c(webb2010_acs$pov10, webb2015_acs$pov15), probs = seq(0,1,.2), na.rm=T)
webb2010map<-webb2010_acs%>%
mutate(pov_q = cut(pov10,breaks = brks, include.lowest = T))
p1<-ggplot(webb2010map, aes( fill=pov_q, color=pov_q))+
geom_sf() +
ggtitle("Poverty Rate 2010 ",
subtitle = "Webb County Texas")+
scale_fill_brewer(palette = "Blues") +
scale_color_brewer(palette = "Blues")+
#scale_fill_manual(values=myfil)+
theme(axis.text.x = element_blank(), axis.text.y = element_blank())+
north(webb2010map)+
scalebar(webb2010map, dist = 5, dd2km =T, model="GRS80", st.size = 2)
webb2015map<-webb2015_acs%>%
mutate(pov_q = cut(pov15,breaks = brks, include.lowest = T))
p2<-ggplot(webb2015map, aes( fill=pov_q, color=pov_q))+
geom_sf() +
ggtitle("Poverty Rate 2015 ",
subtitle = "Webb County Texas")+
scale_fill_brewer(palette = "Blues") +
scale_color_brewer(palette = "Blues")+
#scale_fill_manual(values=myfil)+
theme(axis.text.x = element_blank(), axis.text.y = element_blank())+
north(webb2015map)+
scalebar(webb2015map, dist = 5, dd2km =T, model="GRS80", st.size = 2)
myfil<-c("grey", "blue", "red")
p3<-ggplot(acs_merge, aes( fill=result))+
geom_sf() +
ggtitle("Comparison of Poverty Rate 2010 to 2015",
subtitle = "Webb County Texas")+
#scale_fill_brewer(palette = "Accent") +
#scale_color_brewer(palette = "Accent")+
scale_fill_manual(values=myfil)+
theme(axis.text.x = element_blank(), axis.text.y = element_blank())+
north(acs_merge)+
scalebar(acs_merge, dist = 5, dd2km =T, model="GRS80", st.size = 2)
p1
p2
p3