Create a bivariate choropleth map using Tidycensus?

It is quite easy to create a choropleth using the tidycensus` R package as discussed in Kyle Walker's "Spatial Analysis with US Census Data".

How might one create a bivariate choropleth using Kyle's tidycensus, which acquires Census data through the Census.gov API?

Choropleth with tidycensus => Easy

miami_income <- get_acs(geography = "tract", year = 2020,
                        variables  = c(
                         Median_Income = "B19013_001"
                        ),
                        state = "FL", county = "Miami-Dade County",
                        geometry = TRUE)
plot(miami_income["estimate"])

Bivariate choropleth?

There is a recently developed biscale package by Christopher Prener, Ph.D. that might be useful. Example from https://cran.r-project.org/web/packages/biscale/vignettes/biscale.html.

## stl_race_income example data in biscale package
stlouis_race_income_bivar <- 
  biscale::bi_class(stl_race_income, x = pctWhite, y = medInc, style = "quantile", dim = 3)
# create map
ggplot() +
  geom_sf(data = stlouis_race_income_bivar, 
          mapping = aes(fill = bi_class), color = "white", 
          size = 0.1, show.legend = FALSE) +
  bi_scale_fill(pal = "GrPink", dim = 3) +
  labs(
    title = "Race and Income in St. Louis, MO",
    subtitle = "Gray Pink (GrPink) Palette"
  ) +
  bi_theme()

This was easier than I thought. Posting because I did not see any search hits combining 'tidycensus' and 'bivariate choropleth' (two-variable heat map) that depicted images. Please note that Kyle Walker discusses concerns about margins-of-error at "Spatial Analysis with US Census Data"

#### Census API key
## Only need this line because API key stored in .Renviron
## If not stored, get Census API key at http://api.census.gov/data/key_signup.html then run
## census_api_key("PLACE CENSUS KEY HERE", install = TRUE)
## per https://rdrr.io/cran/tidycensus/man/census_api_key.html
Sys.getenv("CENSUS_API_KEY")
## remotes::install_github("chris-prener/biscale")

library(tidycensus)
library(cowplot)
library(sf)
library(biscale)
library(dplyr)

miami_rent_income <- get_acs(geography = "tract", year = 2020,
                             variables  = c(
                              Median_Income = "B19013_001",
                              Median_Rent = "B25031_001"
                             ),
                             state = "FL", county = "Miami-Dade County",
                             geometry = TRUE,
                             output = "wide")

miami_rent_income <- miami_rent_income %>% 
  rename(Median_Income = Median_IncomeE,
         Median_Rent = Median_RentE)

miami_rent_income_bivar <- 
  biscale::bi_class(miami_rent_income, 
                    x = Median_Income, 
                    y = Median_Rent, style = "quantile", dim = 3)
# create map
miami_rent_income_bivar_map <- 
  ggplot() +
  geom_sf(data = miami_rent_income_bivar, 
          mapping = aes(fill = bi_class), color = "white", 
          size = 0.1, show.legend = FALSE) +
  bi_scale_fill(pal = "GrPink", dim = 3) +
  labs(
    title = "Median Rent and Income in Miami, FL",
    subtitle = "Census American Community Survey (2020)"
  ) +
  bi_theme()

miami_rent_income_bivar_legend <- bi_legend(pal = "GrPink",
                    dim = 3,
                    xlab = "Higher Median Income",
                    ylab = "Higher Median Rent ",
                    size = 8)

# combine map with legend
miami_rent_income_bivar_map2 <- ggdraw() +
  draw_plot(miami_rent_income_bivar_map, 0, 0, 1, 1) +
  draw_plot(miami_rent_income_bivar_legend, 0.2, .65, 0.2, 0.2)

miami_rent_income_bivar_map2