The scatter (${rmsigma_{text{sSFR}}}$) of the specific star formation rates (sSFRs) of galaxies is a measure of the diversity in their star formation histories (SFHs) at a given mass. In this paper we employ the EAGLE simulations to study the dependence of the ${rm sigma_{text{sSFR}}}$ of galaxies on stellar mass (${rm M_{star}}$) through the ${rm sigma_{text{sSFR}}}$-${rm M_{star}}$ relation in $ {rm z sim 0-4}$. We find that the relation evolves with time, with the dispersion depending on both stellar mass and redshift. The models point to an evolving U-shape form for the ${rm sigma_{text{sSFR}}}$-${rm M_{star}}$ relation with the scatter being minimal at a characteristic mass $M^{star}$ of ${rm 10^{9.5}}$ ${rm M_{odot}}$ and increasing both at lower and higher masses. This implication is that the diversity of SFHs increases towards both at the low- and high-mass ends. We find that active galactic nuclei feedback is important for increasing the ${rm sigma_{text{sSFR}}}$ for high mass objects. On the other hand, we suggest that SNe feedback increases the ${rm sigma_{text{sSFR}}}$ of galaxies at the low-mass end. We also find that excluding galaxies that have experienced recent mergers does not significantly affect the ${rm sigma_{text{sSFR}}}$-${rm M_{star}}$ relation. Furthermore, we employ the combination of the EAGLE simulations with the radiative transfer code SKIRT to evaluate the effect of SFR/stellar mass diagnostics in the ${rm sigma_{text{sSFR}}}$-${rm M_{star}}$ relation and find that the ${rm SFR/M_{star}}$ methodologies (e.g. SED fitting, UV+IR, UV+IRX-$beta$) widely used in the literature to obtain intrinsic properties of galaxies have a large effect on the derived shape and normalization of the ${rm sigma_{text{sSFR}}}$-${rm M_{star}}$ relation.