  We perform a detailed study of the cosmological bias of gravitational gave
(GW) events produced by binary black hole mergers (BBHM). We start from a BBHM
distribution modeled inside the EAGLE hydrodyamical simulation using the
population synthesis code MOBSE. We then compare our findings with predictions
from different Halo Occupation Distribution (HOD) prescriptions and find
overall agreement, provided that the modeled properties of host galaxies and
halos in the semi-analytical treatment match those in the simulations. By
highlighting the sources of these discrepancies, we provide the stepping stone
to build future more robust models that prevent the shortcoming of both
simulation-based and analytical models. Finally, we train a neural network to
build a simulation-based HOD and perform feature importance analysis to gain
intuition on which host halo/galaxy parameters are the most relevant in
determining the actual distribution and power spectrum of BBHM. We find that
the distribution of BBHM in a galaxy does not only depend on its size, star
formation rate and metallicity, but also by its kinetic state.
