This work presents a comprehensive analysis of the statistical mechanics of randomly cross-linked polymer gels, starting from a microscopic model of a network made of instantaneously cross-linked Gaussian chains with excluded volume, and ending with the derivation of explicit expressions for the thermodynamic functions and for the density correlation functions which can be tested by experiments. Using replica field theory we calculate the mean field density in replica space and show that this solution contains statistical information about the behavior of individual chains in the network. The average monomer positions change affinely with macroscopic deformation and fluctuations about these positions are limited to length scales of the order of the mesh size. We prove that a given gel has a unique state of microscopic equilibrium which depends on the temperature, the solvent, the average monomer density and the imposed deformation. This state is characterized by the set of the average positions of all the monomers or, equivalently, by a unique inhomogeneous monomer density profile. Gels are thus the only known example of equilibrium solids with no long-range order. We calculate the RPA density correlation functions that describe the statistical properties of small deviations from the average density, due to both static spatial heterogeneities (which characterize the inhomogeneous equilibrium state) and thermal fluctuations (about this equilibrium). We explain how the deformation-induced anisotropy of the inhomogeneous equilibrium density profile is revealed by small angle neutron scattering and light scattering experiments, through the observation of the butterfly effect. We show that all the statistical information about the structure of polymer networks is contained in two parameters whose values are determined by the conditions of synthesis: the density of cross-links and the heterogeneity parameter. We find that the structure of instantaneously cross-linked gels becomes increasingly inhomogeneous with the approach to the cross-link saturation threshold at which the heterogeneity parameter diverges. Analytical expressions for the correlators of deformed gels are derived in both the long wavelength and the short wavelength limits and an exact expression for the total static structure factor, valid for arbitrary wavelengths, is obtained for gels in the state of preparation. We adapt the RPA results to gels permeated by free labelled chains and to gels in good solvents (in the latter case, excluded volume effects are taken into account exactly) and make predictions which can be directly tested by scattering and thermodynamic experiments. Finally, we discuss the limitations and the possible extensions of our work.