We examine the efficacy of the halo mass function as a probe of $f(R)$, symmetron and DGP gravity. We develop an excursion-set method to generalise a range of popular HMF fitting functions from GR to screened MG, considering the HMF dependence on the critical density parameter $\delta_{c}$. In particular, we propose a variety of new methods to account for the environmental dependence of chameleon screening and compare their accuracy to existing ones via $N$-body data. Finally, using the nested sampling routine MultiNest, we examine two possible interpretations of the MG $N$-body results: whether they can be interpreted as $\Lambda$CDM HMFs and whether the MG HMFs display any universality. We find that the effects of MG can be interpreted as a change in best-fit parameters in the $\Lambda$CDM HMF---i.e. can be mistaken for GR---for all of the fitting functions. Alternatively, the relation can be inverted to judge the universality of the HMF---not in terms of its purported invariance to $z$, $H_{0}$, $\Omega_{m}$, $\Omega_{\Lambda}$ as in GR, but its independence on the underlying theory of gravity. Although we found no completely universal HMF, the parameter values did cluster according to the type of screening mechanism. The results suggest that a single, best-fit HMF might be used for each type of screening, independent of the parameters in the MG model. This demonstrates that the additional complexity of the gravitational collapse in screened MG theories cannot always be accounted for using the techniques developed in GR. However, the variations are small enough that it is unnecessary to develop new fitting functions and calibrate them on a case-by-case basis.

An important indicator of modified gravity is the effect of the local environment on halo properties. This paper examines the influence of the local tidal structure on the halo mass function, the halo orientation, spin and the concentration-mass relation. We generalise the excursion set formalism to produce a halo mass function conditional on large-scale structure. Our model agrees well with simulations on large scales at which the density field is linear or weakly non-linear. Beyond this, our principal result is that $f(R)$ does affect halo abundances, the halo spin parameter and the concentration-mass relationship in an environment-independent way, whereas we find no appreciable deviation from $\Lambda$CDM for the mass function with fixed environment density, nor the alignment of the orientation and spin vectors of the halo to the eigenvectors of the local cosmic web. There is a general trend for greater deviation from $\Lambda$CDM in underdense environments and for high-mass haloes, as expected from chameleon screening.

Relativistic scalar fields are ubiquitous in modified theories of gravity. An important tool in understanding their impact on structure formation, especially in the context of N-body simulations, is the quasi-static approximation in which the time evolution of perturbations in the scalar fields is discarded. We show that this approximation must be used with some care by studying linearly perturbed scalar field cosmologies and quantifying the errors that arise from taking the quasi-static limit. We focus on f(R) and chameleon models and link the accuracy of the quasi-static approximation to the fast/slow-roll behaviour of the background and its proximity to {\Lambda}CDM. Investigating a large range of scales, from super- to sub-horizon, we find that slow-rolling ({\Lambda}CDM-like) backgrounds generically result in good quasi-static behaviour, even on (super-)horizon scales. We also discuss how the approximation might affect studying the non-linear growth of structure in numerical N-body simulations.