ABSTRACT All modern galaxy formation models employ stochastic elements in their sub-grid prescriptions to discretize continuous equations across the time domain. In this paper, we investigate how the stochastic nature of these models, notably star formation, black hole accretion, and their associated feedback, that act on small (< kpc) scales, can back-react on macroscopic galaxy properties (e.g. stellar mass and size) across long (> Gyr) time-scales. We find that the scatter in scaling relations predicted by the EAGLE model implemented in the swift code can be significantly impacted by random variability between re-simulations of the same object, even when galaxies are resolved by tens of thousands of particles. We then illustrate how re-simulations of the same object can be used to better understand the underlying model, by showing how correlations between galaxy stellar mass and black hole mass disappear at the highest black hole masses (MBH > 108 M⊙), indicating that the feedback cycle may be interrupted by external processes. We find that although properties that are collected cumulatively over many objects are relatively robust against random variability (e.g. the median of a scaling relation), the properties of individual galaxies (such as galaxy stellar mass) can vary by up to 25 per cent, even far into the well-resolved regime, driven by bursty physics (black hole feedback), and mergers between galaxies. We suggest that studies of individual objects within cosmological simulations be treated with caution, and that any studies aiming to closely investigate such objects must account for random variability within their results.