We discuss various analytical approximation methods for following the evolution of cosmological density perturbations into the strong (i.e. nonlinear) clustering regime. We start by giving a thorough treatment of linear gravitational instability in cosmological models and discussing the statistics of primordial density fluctuations produced in various scenarios of structure formation, and the role of non-baryonic dark matter. We critically review various methods for dealing with the non-linear evolution of density inhomogeneities, in the context of theories of the formation of galaxies and large-scale structure. These methods can be classified into five types: (i) simple extrapolations from linear theory, such as the high-peak model and the lognormal model; (ii) dynamical approximations, including the Zel'dovich approximation and its extensions; (iii) non-linear models based on purely geometric considerations, of which the main example is the Voronoi model; (iv) statistical solutions involving scaling arguments, such as the hierarchical closure ansatz for BBGKY, fractal models and the thermodynamic model of Saslaw; (v) numerical techniques based on particles and/or hydrodynamics. We compare the results of full dynamical evolution using particle codes and the various other approximation schemes. To put the models we discuss into perspective, we give a brief review of the observed properties of galaxy clustering and the statistical methods used to quantify it, such as correlation functions, power spectra, topology and spanning trees.