A modification of the effective field method that allows improving its predictions in the region of high volume fractions of inclusions is considered. In this version, the local exciting field acting on the inclusions is assumed to be a random function. The method allows deriving an infinite chain of integral equations for n-point statistical moments of the effective field. Additional assumptions allow cutting this chain at the n-th equation for the calculation of the multipoint statistical moments of the order 1, 2, ..., n. Then, the method allows the construction of the effective properties of the composites, taking into account the n-point statistical moments of the effective field. An example of the conductive medium with spherical homogeneous inclusions is considered. The two-point Percus-Yevick correlation functions of the centers of non-penetrating spheres are used in the calculations. The predictions of the method are compared with numerical simulations.