How to get geom_hex to overlay on a map?

I am trying to overlay lunge_count averages across hexagons on a map; however, I am only able to produce two separate plots (one of the geom_hex and one of the map).

#libraries
library(dplyr) # for tidying
library(lubridate) # for working with datetimes 
library(readr) # for reading in data
library(sf) # for working with spatial data 
library(ggplot2) # for plotting/mapping 
library(ggspatial) # mapping aesthetics 
library(rnaturalearth) # for getting land data
library(rnaturalearthhires) # for getting land data
library(patchwork) #patching maps together
library(leaflet) # Has neat colour palette function
library(plotrix)

Here is the data I am using:

      ptt depth_m duration_s duration_m lunge_count year               id       lon      lat
9946 5883      51        778   12.96667           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9947 5883      53        714   11.90000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9948 5883      49        696   11.60000           1 2023 2023MX-Mno-05883 -107.2284 22.81825
9949 5883      45        756   12.60000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9950 5883      45        780   13.00000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9951 5883      45        752   12.53333           0 2023 2023MX-Mno-05883 -107.2284 22.81825

Here is the code I am using:

#ocean map information

esri_ocean <- paste0('https://services.arcgisonline.com/arcgis/rest/services/','Ocean/World_Ocean_Base/MapServer/tile/${z}/${y}/${x}.jpeg')

#create ocean basemap

basemap36 <- ggplot() + annotation_map_tile(type = esri_ocean, zoomin = 1, progress = "none") + coord_sf(crs = 4326,xlim = c(-100, -125),ylim = c(17, 36))

#plotting map and hexagon lunge count averages

basemap36 & ggplot(breeding_redo_5801, aes(lon, lat, z = lunge_count)) & stat_summary_hex(fun = mean)

#result

#other attempts

I was able to get the plot to overlay the map with the code below; however, it produced sums instead of averages across the hexagons.

basemap36 & geom_hex(data = breeding_redo_5801, aes(x = lon, y = lat, z = lunge_count)) & stat_summary_hex(fun = mean)

#result of this attempt

I believe this will work, but I don't think we have enough example data to try it:

basemap36 & stat_summary_hex(
  data = breeding_redo_5801, 
  mapping = aes(lon, lat, z = lunge_count), 
  fun = mean
)