There are, to my knowledge, two functions in R to hexbin: hexBinning and geom_hex in the fMultivar and ggplot2 packages respectively. And both parametrize the centers of the hexagons according to the coordinates of hte lower-leftmost point in the sample.
It means that if you split your sample (according to a factor or, in my case, inside a mapreduce job) your hexagons become eccentric.
So I implemented my own hexbin function that assumes (0,0) as the center of the grid (i.e., if there were points around (0,0), the corresponding hexagon would be centered there) and requires just r (the radious of the hexagon) as parameter.
The implementation is here (sorry, text is in Spanish!). Moreover, my implementation has no explicit loops: it is fully vectorized.