Adding loess regresion line on a hexbin plot

I would suggest using ggplot2 to build the plot.

Since you didn't include any example data, I've used the palmerpenguins package dataset for the example below.

library(palmerpenguins) # For the data
library(ggplot2)        # ggplot2 for plotting



ggplot(penguins, aes(x = body_mass_g, 
                     y = bill_length_mm)) + 
  geom_hex(bins = 40) + 
  geom_smooth(method = 'loess', se = F, color = 'red') 

^{Created on 2021-01-05 by the reprex package (v0.3.0)}

I don't have a solution for base, but it's possible to do this with ggplot. It should be possible with base too, but if you look at the documentation for ?hexbin, you can see the quote:

Note that when plotting a hexbin object, the grid package is used. You must use its graphics (or those from package lattice if you know how) to add to such plots.

I'm not familiar with how to modify these. I did try ggplotify to convert the base to ggplot and edit that way, but couldn't get the loess line added to the plot window properly.

So here is a solution with ggplot with some fake data that you can try on your Datasets:

library(hexbin) 
library(ggplot2) 

# fake data with a random walk, replace with your data
set.seed(100)
N <- 1000
x <- rnorm(N)
x <- sort(x)
y <- vector("numeric", length=N)
for(i in 2:N){
  y[i] <- y[i-1] + rnorm(1, sd=0.1)
}

# current method
# In documentation for ?hexbin it says:
# "You must use its graphics (or those from package lattice if you know how) to add to such plots."
(bin <- hexbin(x, y, xbins=40))
plot(bin)

# ggplot option. Can play around with scale_fill_gradient to
# get the colour scale similar or use other ggplot options
df <- data.frame(x=x, y=y)
d <- ggplot(df, aes(x, y)) + 
  geom_hex(bins=40) + 
  scale_fill_gradient(low = "grey90", high = "black") +
  theme_bw()
d

# easy to add a loess fit to the data
# span controls the degree of smoothing, decrease to make the line
# more "wiggly"
model <- loess(y~x, span=0.2)
fit <- predict(model)
loess_data <- data.frame(x=x, y=fit)
d + geom_line(data=loess_data, aes(x=x, y=y), col="darkorange",
              size=1.5)

Here are two options; you will need to decide if you want to smooth over the raw data or the binned data.

library(hexbin)
library(grid)

# Some data
set.seed(101)
d <- data.frame(x=rnorm(1000))
d$y <- with(d, 2*x^3 + rnorm(1000))

Method A - binned data

# plot hexbin & smoother : need to grab plot viewport
# From ?hexVP.loess : "Fit a loess line using the hexagon centers of mass 
# as the x and y coordinates and the cell counts as weights."
bin <- hexbin(d$x, d$y)
p <- plot(bin)
hexVP.loess(bin, hvp = p$plot.vp, span = 0.4, col = "red", n = 200)

Method B - raw data

# calculate loess predictions outside plot on raw data
l = loess(y ~ x, data=d, span=0.4)
xp = with(d, seq(min(x), max(x), length=200))
yp = predict(l, xp)

# plot hexbin
bin <- hexbin(d$x, d$y)
p <- plot(bin)

# add loess line
pushHexport(p$plot.vp)
grid.lines(xp, yp, gp=gpar(col="red"), default.units = "native")
upViewport()