We derive the macroscopic Einstein–Maxwell equations up to the second-order terms in the interaction for systems with dominating electromagnetic interactions between particles (e.g., radiation-dominated cosmological plasma in the expanding Universe before the recombination moment). The ensemble averaging of the microscopic Einstein and Maxwell equations and of the Liouville equations for the random functions of each type of particle leads to a closed system of equations consisting of the macroscopic Einstein and Maxwell equations and the kinetic equations for one-particle distribution functions for each type of particle. The macroscopic Einstein equations for a system of electromagnetically and gravitationally interacting particles differ from the classical Einstein equations in having additional terms in the left-hand side due to the interaction. These terms are given by a symmetric rank-two traceless tensor with zero divergence. Explicitly, these terms are represented as momentum-space integrals of the expressions containing one-particle distribution functions for each type of particle and have much in common with similar terms in the left-hand side of the macroscopic Einstein equations previously obtained for a system of self-gravitating particles. The macroscopic Maxwell equations for a system of electromagnetically and gravitationally interacting particles also differ from the classical Maxwell equations in having additional terms in the left-hand side due to simultaneous effects described by general relativity and the interaction effects.