  We explore two widely used empirical models for the galaxy-halo connection,
subhalo abundance matching (SHAM) and the halo occupation distribution (HOD)
and compare their predictions with the hydrodynamical simulation IllustrisTNG
(TNG) for a range of statistics that quantify the galaxy distribution at
$n_{\rm gal}\approx1.3\times10^{-3}\,[{\rm Mpc}/h]^{-3}$. We observe that in
their most straightforward implementations, both models fail to reproduce the
two-point clustering measured in TNG. We find that SHAM models constructed
using the relaxation velocity, $V_{\rm relax}$, and the peak velocity, $V_{\rm
peak}$, perform best, and match the clustering reasonably well, although
neither model captures adequately the one-halo clustering. Splitting the total
sample into sub-populations, we discover that SHAM overpredicts the clustering
of high-mass, blue, star-forming, and late-forming galaxies and uderpredicts
that of low-mass, red, quiescent, and early-forming galaxies. We also study
various baryonic effects, finding that subhalos in the dark-mater-only
simulation have consistently higher values of their SHAM-proxy properties than
their full-physics counterparts. We then consider a two-dimensional
implementation of the HOD model augmented with a secondary parameter
(environment, velocity anisotropy, $\sigma^2R_{\rm halfmass}$, and total
potential) and tuned so as to match the two-point clustering of the
IllustrisTNG galaxies on large scales. We analyze these galaxy populations
adopting alternative statistical tools such as galaxy-galaxy lensing,
void-galaxy cross-correlations and cumulants of the smoothed density field,
finding that the hydrodynamical galaxy distribution disfavors $\sigma^2 R_{\rm
halfmass}$ and the total potential as secondary parameters, while the
environment and velocity anisotropy samples are consistent with full-physics
across all statistical probes examined.
