In this example, I want to comapre two variables (Housing units and Hispanic proportion) in Harris and Dallas counties. I use the data from ACS, then make Jenks break map for each variable and county.In the finall part, I make a brief summary for the variables.

Look at available ACS variables

library(tidycensus);  library(tidyverse); library(sf)
## -- Attaching packages ------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 2.2.1.9000     v purrr   0.2.4     
## v tibble  1.4.2          v dplyr   0.7.5     
## 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()
## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
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)
# table(v15_Profile$DP05_0066PE)
# table(v15_Profile$DP04_0019PE)
#Search for variables by 
v15_Profile[grep(x = v15_Profile$label, "Housing units"), c("name", "label")]
## # A tibble: 50 x 2
##    name        label                                                      
##    <chr>       <chr>                                                      
##  1 DP04_0091E  Estimate!!MORTGAGE STATUS!!Owner-occupied units!!Housing u~
##  2 DP04_0091PE Percent!!MORTGAGE STATUS!!Owner-occupied units!!Housing un~
##  3 DP04_0092E  Estimate!!MORTGAGE STATUS!!Owner-occupied units!!Housing u~
##  4 DP04_0092PE Percent!!MORTGAGE STATUS!!Owner-occupied units!!Housing un~
##  5 DP04_0093E  Estimate!!SELECTED MONTHLY OWNER COSTS (SMOC)!!Housing uni~
##  6 DP04_0093PE Percent!!SELECTED MONTHLY OWNER COSTS (SMOC)!!Housing unit~
##  7 DP04_0094E  Estimate!!SELECTED MONTHLY OWNER COSTS (SMOC)!!Housing uni~
##  8 DP04_0094PE Percent!!SELECTED MONTHLY OWNER COSTS (SMOC)!!Housing unit~
##  9 DP04_0095E  Estimate!!SELECTED MONTHLY OWNER COSTS (SMOC)!!Housing uni~
## 10 DP04_0095PE Percent!!SELECTED MONTHLY OWNER COSTS (SMOC)!!Housing unit~
## # ... with 40 more rows
v15_Profile[grep(x = v15_Profile$label, "enrolled"), c("name", "label")]
## # A tibble: 24 x 2
##    name        label                                                      
##    <chr>       <chr>                                                      
##  1 DP02_0052E  Estimate!!SCHOOL ENROLLMENT!!Population 3 years and over e~
##  2 DP02_0052PE Percent!!SCHOOL ENROLLMENT!!Population 3 years and over en~
##  3 DP02_0053E  Estimate!!SCHOOL ENROLLMENT!!Population 3 years and over e~
##  4 DP02_0053PE Percent!!SCHOOL ENROLLMENT!!Population 3 years and over en~
##  5 DP02_0054E  Estimate!!SCHOOL ENROLLMENT!!Population 3 years and over e~
##  6 DP02_0054PE Percent!!SCHOOL ENROLLMENT!!Population 3 years and over en~
##  7 DP02_0055E  Estimate!!SCHOOL ENROLLMENT!!Population 3 years and over e~
##  8 DP02_0055PE Percent!!SCHOOL ENROLLMENT!!Population 3 years and over en~
##  9 DP02_0056E  Estimate!!SCHOOL ENROLLMENT!!Population 3 years and over e~
## 10 DP02_0056PE Percent!!SCHOOL ENROLLMENT!!Population 3 years and over en~
## # ... with 14 more rows

Extract from ACS summary file data profile variables from 2015

ha_hisp<-get_acs(geography = "tract", state="TX", county = "Harris County", year = 2015,
                variables= c("DP05_0066PE","DP04_0019PE"),
                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)`.
## Using the ACS Data Profile
#create a county FIPS code - 5 digit
ha_hisp$county<-substr(ha_hisp$GEOID, 1, 5)

#rename variables and filter missing cases
ha_hisp1<-ha_hisp%>%
  mutate( phisp=DP05_0066PE,housing=DP04_0019PE ) %>%
#  st_transform(crs = 102740)%>%
  na.omit()

class(ha_hisp1)

Write data out to shapefile

#change the directory
sf::st_write(ha_hisp1,dsn="C:/Users/Ali/Desktop/topic motality/code",layer="ha_tract_dp", driver="ESRI Shapefile", delete_layer=T, update=T)

Some basic mapping of variables

#install.packages("classInt")
#install.packages("patchwork")

library(classInt)
#library(patchwork)


#install.packages("sf")
library(sf)
library(dplyr)
library(ggplot2)

phisp_map<-ha_hisp1 %>%
  mutate(pHispanic=cut(phisp,breaks = quantile(phisp, na.rm=T, p=seq(0,1,length.out = 6)),include.lowest = T),
           pHousing = cut(phisp,breaks=data.frame(classIntervals(var=ha_hisp1$phisp, n=5, style="jenks")[2])[,1], include.lowest = T))
library(ggsn)
# Harris County Hispanic Proportions built 2000-2009


p1<-ggplot(phisp_map, aes(fill = pHispanic, color = pHispanic)) + 
  geom_sf() + 
  ggtitle("Proportion of Hispanics", 
          subtitle = "Harris County Texas, 2015 - Jenks Breaks")+
  coord_sf(crs =  102009) + 
  scale_fill_brewer(palette = "Blues") + 
  scale_color_brewer(palette = "Blues")+
  theme(axis.text.x = element_blank(), axis.text.y = element_blank())
p1

#Harris County Housing Proportions built 2000-2009


housing_map<-ha_hisp1 %>%
  mutate(pHousing=cut(housing,breaks = quantile(housing, na.rm=T, p=seq(0,1,length.out = 6)),include.lowest = T),
           pHousing = cut(housing,breaks=data.frame(classIntervals(var=ha_hisp1$housing, n=5, style="jenks")[2])[,1], include.lowest = T))

  p2<-ggplot(housing_map, aes(fill = pHousing, color = pHousing)) + 
  geom_sf() + 
  ggtitle("Proportion Housing", 
          subtitle = "Harris County Texas, 2015 - Jenks Breaks")+
  coord_sf(crs =  102009) + 
  scale_fill_brewer(palette = "Blues") + 
  scale_color_brewer(palette = "Blues")+
    theme(axis.text.x = element_blank(), axis.text.y = element_blank())
p2

#p1/p2

Descriptive analysis

x<-ha_hisp1%>%
  ggplot()+geom_histogram(aes(x =housing , y=0.5*..density..))+
  facet_wrap(~county)+
  ggtitle(label = "Poverty Rate in Bexar and Brazos County")
x
ha_hisp<-get_acs(geography = "tract", state="TX", county = "Harris County", year = 2015,
                variables= c("DP05_0066PE","DP04_0019PE"),
                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)`.
## Using the ACS Data Profile
#create a county FIPS code - 5 digit
ha_hisp$county<-substr(ha_hisp$GEOID, 1, 5)

#rename variables and filter missing cases
ha_hisp1<-ha_hisp%>%
  mutate( phisp=DP05_0066PE,housing=DP04_0019PE ) %>%
#  st_transform(crs = 102740)%>%
  na.omit()

class(ha_hisp1)
#table(ha_hisp1$phisp)

Write data out to shapefile

#change the directory
sf::st_write(ha_hisp1,dsn="C:/Users/Ali/Desktop/topic motality/code",layer="ha_tract_dp", driver="ESRI Shapefile", delete_layer=T, update=T)

Some basic mapping of variables

da_hisp<-get_acs(geography = "tract", state="TX", county = "Dallas County", year = 2015,
                variables= c("DP05_0066PE","DP04_0019PE"),
                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)`.
## Using the ACS Data Profile
#create a county FIPS code - 5 digit
da_hisp$county<-substr(da_hisp$GEOID, 1, 5)

#rename variables and filter missing cases
da_hisp1<-da_hisp%>%
  mutate( phisp=DP05_0066PE,housing=DP04_0019PE ) %>%
#  st_transform(crs = 102740)%>%
  na.omit()

class(da_hisp1)
## [1] "sf"         "data.frame"

Dallas County (Hispanic and Housing Proportion )

phisp_map_da<-da_hisp1 %>%
  mutate(pHispanic=cut(phisp,breaks = quantile(phisp, na.rm=T, p=seq(0,1,length.out = 6)),include.lowest = T),
           pHousing = cut(phisp,breaks=data.frame(classIntervals(var=da_hisp1$phisp, n=5, style="jenks")[2])[,1], include.lowest = T))
library(ggsn)
# Dallas County Hispanic Proportions built 2000-2009


p3<-ggplot(phisp_map_da, aes(fill = pHispanic, color = pHispanic)) + 
  geom_sf() + 
  ggtitle("Proportion of Hispanics", 
          subtitle = "Dallas County Texas, 2015 - Jenks Breaks")+
  coord_sf(crs =  102009) + 
  scale_fill_brewer(palette = "Blues") + 
  scale_color_brewer(palette = "Blues")+
  theme(axis.text.x = element_blank(), axis.text.y = element_blank())
p3

#Dallas County Housing Proportions built 2000-2009


housing_map_da<-da_hisp1 %>%
  mutate(pHousing=cut(housing,breaks = quantile(housing, na.rm=T, p=seq(0,1,length.out = 6)),include.lowest = T),
           pHousing = cut(housing,breaks=data.frame(classIntervals(var=da_hisp1$housing, n=5, style="jenks")[2])[,1], include.lowest = T))

  p4<-ggplot(housing_map_da, aes(fill = pHousing, color = pHousing)) + 
  geom_sf() + 
  ggtitle("Proportion Housing", 
          subtitle = "Dallas County Texas, 2015 - Jenks Breaks")+
  coord_sf(crs =  102009) + 
  scale_fill_brewer(palette = "Blues") + 
  scale_color_brewer(palette = "Blues")+
    theme(axis.text.x = element_blank(), axis.text.y = element_blank())
p4

#p1/p2
## Summaries
# Histograms
hist(ha_hisp1$housing)

hist(da_hisp1$housing)

# Summaries of Harris County
summary(ha_hisp1$phisp)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   20.68   35.70   41.50   60.40   97.40
summary(ha_hisp1$housing)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    3.50   11.65   17.65   26.12   90.60
#Summaries of Dallas County
summary(da_hisp1$phisp)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   16.20   33.50   37.65   57.00   96.20
summary(da_hisp1$housing)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    2.15    6.80   12.75   18.50   83.50
# Results
# In comparison between to counties, the highest density can be seen in central and southern part in Harris (more than 70 percent) and the northern part has lower proportion of Hispanics. Dallas county, as well, has this pattern, when in the central parts of the county (both Eastern and Western), we can see high density of Hispanics, while there is a lower proportion of Hispanics in the Downtown. Again, in the central northern part, the density very low, but in the southern part this population has a higher density.
# In terms of housing, the central parts of Harris have very low density of housing, while northern part has greater proportion. Dallas county, in general has lower density than Harris county. The southern part has a higher proportion of housing than northern part. 
 #In the summary part, I calculate the summaries for each variable separately that show the mean of housing in  Harris County is close to 18, and the mean of housing in Dallas County 12.75. For the Hispanic population, the mean in Harris County in 41.50, and 37.65 in Dallas County.