The key is to realize what choroplethr actually does and returns.

As the answer above shows, creating choropleth maps with ggplot2 involves many steps. Choroplethr does all those steps for you (as well as others, such as alerting you about duplicate or missing values), and returns a ggplot2 object.

library(choroplethr)

df1 <- structure(list(region = c(6001L, 6003L, 6005L, 6007L, 6009L, 
                                 6011L, 6013L, 6015L, 6017L, 6019L, 6021L, 6023L, 6025L, 6027L, 
                                 6029L, 6031L, 6033L, 6035L, 6037L, 6039L, 6041L, 6043L, 6045L, 
                                 6047L, 6049L, 6051L, 6053L, 6055L, 6057L, 6059L, 6061L, 6063L, 
                                 6065L, 6067L, 6069L, 6071L, 6073L, 6075L, 6077L, 6079L, 6081L, 
                                 6083L, 6085L, 6087L, 6089L, 6091L, 6093L, 6095L, 6097L, 6099L, 
                                 6101L, 6103L, 6105L, 6107L, 6109L, 6111L, 6113L, 6115L), value = c(1.6, 
                                                                                                    0, 0, 0.2, 0, 0, 0.7, 0.1, 0.1, 1.2, 0.1, 0.3, 0.1, 0.1, 1.4, 
                                                                                                    0.2, 0.3, 0, 11.8, 0.2, -0.2, 0, -0.1, 0.5, 0, 0, -0.1, 0.2, 
                                                                                                    -0.2, -5, -0.3, 0.1, -1.9, 0.5, -0.1, -2.1, -2.5, -1, -0.5, -0.4, 
                                                                                                    0, 0.1, -1.8, 0, 0.3, 0, 0.2, 0.6, -0.5, 0.2, 0.2, 0.1, 0.2, 
                                                                                                    0.2, 0.2, -0.8, 0.1, 0.3)), class = "data.frame", row.names = c(NA, 
                                                                                                                                                                    -58L))
choro = county_choropleth(df1,
                      legend="PercentnOvernExpectednBased OnnPopulation",
                      num_colors = 0,
                      state_zoom=c("california"))
class(choro)
[1] "gg"     "ggplot"

A nice feature of ggplot2 is that you can override existing features of a plot by using the + operator. This means that you can override the default scale choroplethr uses with the one from the previous answer like this:

library(ggplot2)

county_choropleth(df1,
                  legend="PercentnOvernExpectednBased OnnPopulation",
                  num_colors = 0,
                  state_zoom=c("california")) + 
scale_fill_gradient2(name = "PercentnOvernExpectednBased OnnPopulation",
                       low = "blue",
                       mid = "white",
                       high = "firebrick",
                       midpoint = 0)  

Modern versions of ggplot2 even give a helpful warning saying that you are overriding the scale, and that it is using the last scale provided. When I run that code, for example, I get this message on the console:

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

IMHO I think you actually had two questions. The first was “How do I create the scale / legend I want for this map?” The second was “How do I get choroplethr to use that scale / legend?”

While choroplethr et al. are extremely useful tools for automating the mapping process, it is often difficult (or not possible) to achieve the custom results you want. I am not implying that there isn’t a way to achieve what you want using choroplethr alone, but I generally find it easier to stick to using ggplot2.

With that in mind, here is an alternative that uses data downloaded using the tigris package. This will load an sf object into R that you can then dplyr::left_join() your df1 object to. As the FIPS values are strings in the tigris data (GEOID), they will first be converted to integers so they match your region values.

Then you can easily apply any ggplot functions to your data. For instance, I have chosen theme_void() as it is the option that most closely resembles the output from your code.

library(tigris)
library(sf)
library(dplyr)
library(ggplot2)

df1 <- structure(list(region = c(6001L, 6003L, 6005L, 6007L, 6009L, 
6011L, 6013L, 6015L, 6017L, 6019L, 6021L, 6023L, 6025L, 6027L, 
6029L, 6031L, 6033L, 6035L, 6037L, 6039L, 6041L, 6043L, 6045L, 
6047L, 6049L, 6051L, 6053L, 6055L, 6057L, 6059L, 6061L, 6063L, 
6065L, 6067L, 6069L, 6071L, 6073L, 6075L, 6077L, 6079L, 6081L, 
6083L, 6085L, 6087L, 6089L, 6091L, 6093L, 6095L, 6097L, 6099L, 
6101L, 6103L, 6105L, 6107L, 6109L, 6111L, 6113L, 6115L), value = c(1.6, 
0, 0, 0.2, 0, 0, 0.7, 0.1, 0.1, 1.2, 0.1, 0.3, 0.1, 0.1, 1.4, 
0.2, 0.3, 0, 11.8, 0.2, -0.2, 0, -0.1, 0.5, 0, 0, -0.1, 0.2, 
-0.2, -5, -0.3, 0.1, -1.9, 0.5, -0.1, -2.1, -2.5, -1, -0.5, -0.4, 
0, 0.1, -1.8, 0, 0.3, 0, 0.2, 0.6, -0.5, 0.2, 0.2, 0.1, 0.2, 
0.2, 0.2, -0.8, 0.1, 0.3)), class = "data.frame", row.names = c(NA, 
-58L))

# Load counties, join df1 data to counties sf
california_counties <- counties(state = "CA", cb = TRUE, year = 2023) |>
  mutate(GEOID = as.integer(GEOID)) |>
  left_join(df1, by = c("GEOID" = "region"))
# Using FIPS code '06' for state 'CA'

# Create state outline
california <- california_counties |>
  summarise(geometry = st_union(geometry))

# Plot
ggplot() +
  geom_sf(data = california_counties,
          aes(fill = value),
          colour = "grey50") +
  geom_sf(data = california, colour = "black", fill = NA) +
  scale_fill_gradient2(name = "PercentnOvernExpectednBased OnnPopulation",
                       low = "blue",
                       mid = "white",
                       high = "firebrick",
                       midpoint = 0) +
  theme_void()