update bivariate map, fix variable names for relative housing cost va…

…riables

01_document/example_census_race_ethnicity_calculation.qmd

@@ -62,6 +62,7 @@ library(mapview)
 library(biscale)
 library(cowplot)
 library(glue)
+library(ggtext)
 
 # conflicts ----
 library(conflicted)
@@ -871,12 +872,12 @@ And, compute the portion of households paying more than 30% / 50% of their incom
 
 ```{r}
 water_system_demographics <- water_system_demographics %>%
-    mutate(households_percent_housing_costs_above_30_percent = 
+    mutate(households_housing_costs_over30pct_percent = 
                100 * (households_mortgage_over30pct_count + 
                           households_no_mortgage_over30pct_count +
                           households_rent_over30pct_count) / 
                households_count) %>% 
-    mutate(households_percent_housing_costs_above_50_percent = 
+    mutate(households_housing_costs_over50pct_percent = 
                100 * (households_mortgage_over50pct_count + 
                           households_no_mortgage_over50pct_count +
                           households_rent_over50pct_count) / 
@@ -1107,12 +1108,12 @@ Prepare interpolation data
 ```{r}
 # census_data_interpolate  census_data_clip
 census_data_interpolate <- census_data_interpolate %>% 
-    mutate(households_percent_housing_costs_above_30_percent = 
+    mutate(households_housing_costs_over30pct_percent = 
                100 * (households_mortgage_over30pct_count + 
                           households_no_mortgage_over30pct_count +
                           households_rent_over30pct_count) / 
                households_count) %>% 
-    mutate(households_percent_housing_costs_above_50_percent = 
+    mutate(households_housing_costs_over30pct_percent = 
                100 * (households_mortgage_over50pct_count + 
                           households_no_mortgage_over50pct_count +
                           households_rent_over50pct_count) / 
@@ -1130,27 +1131,35 @@ systems_top_n <- water_system_demographics %>%
 
 ### Bivariate Map
 
-All water systems:
+This section uses the [@biscale] R package to create bivariate choropleth maps that show how two variables vary together spatially.
+
+@fig-bivariate-all shows the relationship between relative housing costs and income – using the estimated data – for the top `{r} n_systems` systems by estimated population in Sacramento.
 
 ```{r}
-# see: https://chris-prener.github.io/biscale/articles/biscale.html
+#| message: false
+#| warning: false
+#| fig-width: 8 
+#| fig-height: 6
+#| label: fig-bivariate-all
+#| code-fold: true
 
-# households_percent_housing_costs_above_30_percent
-# households_percent_housing_costs_above_50_percent
+# Table B25140 - Housing Costs as a Percentage of Household Income in the past 12 months.
+# Shows the count of households paying more than 30% or 50% of their income towards housing costs broken out by three tenure categories (owned with a mortgage, owned without a mortgage, and rented).
 
-# Housing Costs as a Percentage of Household Income in the past 12 months. Table B25140 shows the count of households paying more than 30% of their income towards housing costs broken out by three tenure categories (owned with a mortgage, owned without a mortgage, and rented). The table also shows the number of households paying more than 50% of their income toward housing costs
+# set defaults
 biscale_pal <- 'BlueOr' # 'GrPink' # 'DkViolet2'
 biscale_dim <- 3
+
 # create classes
 biscale_data <- bi_class(water_system_demographics %>% 
                              filter(WATER_SY_1 %in% systems_top_n) %>% 
                              filter(!is.na(median_household_income_hh_weighted)), 
-                         x = households_percent_housing_costs_above_30_percent, 
+                         x = households_housing_costs_over30pct_percent, 
                          y = median_household_income_hh_weighted, 
                          style = "quantile", 
                          dim = biscale_dim)
 
-# create map - all systems
+# create map
 biscale_map <- ggplot() +
     geom_sf(data = biscale_data, 
             mapping = aes(fill = bi_class), 
@@ -1160,48 +1169,70 @@ biscale_map <- ggplot() +
     bi_scale_fill(pal = biscale_pal, 
                   dim = biscale_dim) + 
     labs(
-        title = "Housing Cost and Income in Sacramento \nWater Systems", 
-        caption = "% Housing cost shows the percent of households paying more than 30% of their income towards housing costs \nIncome shows median household income (yellow = missing)",
-        subtitle = glue("Top {n_systems} systems by population")
+        title = "Estimated % of Households Paying More Than 30% of Income Towards Housing Costs \nand Estimated Median Household Income in Sacramento Water Systems",
+        subtitle = glue("Top {n_systems} systems by population"),
+        caption = glue("Data estimated from {acs_year} 5-year ACS Block Groups")
+        # title = "Estimated Housing Cost as % of Household Income and \nEstimated Median Household Income in Sacramento Water Systems", 
+        # caption = "% Housing cost shows the percent of households paying more than 30% of their income towards housing costs \nIncome shows median household income (yellow = missing)"
     ) +
-    # add the missing polygons back in
+    #   labs(
+    #   title = "Housing Cost<sup>1</sup> and Income<sup>2</sup> in Sacramento Water Systems",
+    #   caption = "<sup>1</sup>% of households paying more than 30% of their income towards housing costs<br><sup>2</sup>Median household income (yellow = missing)",
+    #   subtitle = glue("Top {n_systems} systems by population")
+    # ) +
+    # add missing polygons back in
     geom_sf(data = water_system_demographics %>% 
-                filter(WATER_SY_1 == system_plot) %>% 
+                filter(WATER_SY_1 %in% systems_top_n) %>% 
                 filter(is.na(median_household_income_hh_weighted)),
             color = "white",
             fill = 'gold'
     ) +
     geom_sf(data = counties_ca %>% filter(NAME == 'Sacramento'), 
             color = 'grey',
             fill = NA) +
-    bi_theme()
+    bi_theme() + 
+    theme(plot.title = element_text(size=12), # element_markdown(size=12)
+          plot.subtitle = element_text(size=10),
+          plot.caption = element_text(size=8, hjust = 1)) # element_markdown(size=8, hjust = 1))
 
+# create legend
 biscale_legend <- bi_legend(pal = biscale_pal,
                             dim = biscale_dim,
                             xlab = "% Housing Costs ",
                             ylab = "Income ",
                             size = 8)
 
+# construct map
 biscale_plot <- ggdraw() +
     draw_plot(biscale_map, 0, 0, 1, 1) +
     draw_plot(biscale_legend, 0.1, .65, 0.2, 0.2)
+
 biscale_plot
 ```
 
-Plot of block groups overlapping `{r} str_to_title(system_plot)`:
+@fig-bivariate-system shows the same variables (relative housing costs and income) for the portions block groups overlapping `{r} str_to_title(system_plot)` – this illustrates the data underlying the interpolation process.
 
 ```{r}
-# create map - single system
+#| message: false
+#| warning: false
+#| fig-width: 8 
+#| fig-height: 6
+#| label: fig-bivariate-system
+#| code-fold: true
+
+# set defaults
 biscale_pal_system <- 'BlueOr' # 'GrPink' # 'DkViolet2'
 biscale_dim_system <- 3
+
+# create classes
 biscale_data_system <- bi_class(census_data_interpolate %>% 
                                     filter(WATER_SY_1 == system_plot) %>% 
                                     filter(!is.na(median_household_income)), 
-                                x = households_percent_housing_costs_above_30_percent, 
+                                x = households_housing_costs_over30pct_percent, 
                                 y = median_household_income, 
                                 style = "quantile", 
                                 dim = biscale_dim_system)
-
+# create map
 biscale_map_system  <- ggplot() +
     geom_sf(data = biscale_data_system , 
             mapping = aes(fill = bi_class), 
@@ -1211,28 +1242,36 @@ biscale_map_system  <- ggplot() +
     bi_scale_fill(pal = biscale_pal_system, 
                   dim = biscale_dim_system) + 
     labs(
-        title = glue("Housing Cost and Income \nin {str_to_title(system_plot)}"), 
-                     caption = "% Housing cost shows the percent of households paying more than 30% of their income towards housing costs \nIncome shows median household income (yellow = missing)"#,
-        # subtitle = ""
+        title = glue("Estimated % of Households Paying More Than 30% of Income Towards Housing Costs \nand Estimated Median Household Income in {str_to_title(system_plot)}"),
+        # subtitle = glue(""),
+        caption = glue("Data from {acs_year} 5-year ACS Block Groups (Yellow = Missing Data)")#,
+        # title = glue("Housing Cost and Income \nin {str_to_title(system_plot)}"), 
+        # caption = "% Housing cost shows the percent of households paying more than 30% of their income towards housing costs \nIncome shows median household income (yellow = missing)"#,
     ) +
     # add the missing polygons back in
     geom_sf(data = census_data_interpolate %>% 
-                                    filter(WATER_SY_1 == system_plot) %>% 
-                                    filter(is.na(median_household_income)),
+                filter(WATER_SY_1 == system_plot) %>% 
+                filter(is.na(median_household_income)),
             color = "white",
             fill = 'gold'
-            ) +
-    bi_theme()
+    ) +
+    bi_theme() + 
+    theme(plot.title = element_text(size=12), # element_markdown(size=12)
+          plot.subtitle = element_text(size=10),
+          plot.caption = element_text(size=8, hjust = 1)) # element_markdown(size=8, hjust = 1))
 
+# create legend
 biscale_legend <- bi_legend(pal = biscale_pal_system,
                             dim = biscale_dim_system,
                             xlab = "% Housing Costs ",
                             ylab = "Income ",
                             size = 8)
 
+# construct map
 biscale_plot_system <- ggdraw() +
   draw_plot(biscale_map_system, 0, 0, 1, 1) +
   draw_plot(biscale_legend, 0.1, .55, 0.2, 0.2)
+
 biscale_plot_system
 ```
 

01_document/references.bib

@@ -70,3 +70,11 @@ @article{tigris
 	date = {2023},
 	url = {https://CRAN.R-project.org/package=tigris}
 }
+
+@article{biscale,
+	title = {biscale: Tools and Palettes for Bivariate Thematic Mapping},
+	author = {Prener, Christopher and Grossenbacher, Timo and Zehr, Angelo},
+	year = {2022},
+	date = {2022},
+	url = {https://CRAN.R-project.org/package=biscale}
+}