Now we can use splom with panel.hexbinplot and panel.loess. Besides, I have included some changes to diag.panel in order to show the univariate density of each variable (adapted from here). UPDATED: I include the use of colramp and BTC to change the color of the hexagon bins.

It is possible but (as far as I know) it’s not easy. The hexbin package includes the function hexlegendGrob which produces a grob with the legend. With it, you can draw the plot and after use update to include the legend. Take a look at the end of the code of hexbinplot.formula to see how to get it. The problem is that you have to provide the value of maxcnt to hexlegendGrob as an argument.