We model the abundance of haloes in the ∼(3 Gpc h −1)3 volume of the MICE Grand Challenge simulation by fitting the universal mass function with an improved Jackknife error covariance estimator that matches theory predictions. We present unifying relations between different fitting models and new predictions for linear (b 1) and non-linear (c 2 and c 3) halo clustering bias. Different mass function fits show strong variations in their performance when including the low mass range (M h ≲ 3 × 1012 M⊙ h −1) in the analysis. Together with fits from the literature, we find an overall variation in the amplitudes of around 10 per cent in the low mass and up to 50 per cent in the high mass (galaxy cluster) range (M h > 1014 M⊙ h −1). These variations propagate into a 10 per cent change in b 1 predictions and a 50 per cent change in c 2 or c 3. Despite these strong variations, we find universal relations between b 1 and c 2 or c 3 for which we provide simple fits. Excluding low-mass haloes, different models fitted with reasonable goodness in this analysis, show per cent level agreement in their b 1 predictions, but are systematically 5–10 per cent lower than the bias directly measured with two-point halo-mass clustering. This result confirms previous findings derived from smaller volumes (and smaller masses). Inaccuracies in the bias predictions lead to 5–10 per cent errors in growth measurements. They also affect any halo occupation distribution fitting or (cluster) mass calibration from clustering measurements.