Soon after the mathematical properties of renmormalized relativistic quantum fields were investigated it became clear that the field equations derived from a Lagrangian couldn't make sense as operator equations. The principal reason is that the product of fields at the same point is ill-defined, due to the distributional nature of the fields. Some exactly solvable models in 2 spacetime dimensions (where normal ordering is the only kind of renormalization needed) revealed that in place of the usual field equations one had equations involving the corresponding normally ordered products, but in three dimensions it was found that even these could not be correct.
Gradually it was found that the correct form of the field equations could always be written in a limiting form involving products of fields at infinitesimally different arguments. This lead to Wilson's concept of operator product expansions (OPEs), in which such limits of products of fields are systematically explored and utilized. (Wilson's influential original paper from 1964 was faulty and therefore never published, and the first published papers on the OPE were two 1967 papers by Brandt and by Zimmermann; see the review of Wilson's second paper from 1969.)
For massive $\phi^4$ theory, the correct renormalized field equations are described in W. Zimmermann Local field equation for $A^4$-coupling in renormalized perturbation theory, Comm. Math. Phys. 6 (1967), 161-188. Note that the mass and the coupling constant appearing in the final equations have their physical, renormalized values. For massive $\phi^4$ one needs an OPE for $\phi(x)\phi(y)$ which introduces a renormalized version $\phi_2(x)$ of $\phi(x)^2$, and an OPE for $\phi(x)\phi_2(x)$, which turns out to introduce not a new renormalized field but the d'Alembertian of the original field.
The massless case is significantly more complicated since anomalous exponents must be accounted for correctly. The correct formulation is described in Brandt and Ng Wing-Chiu, Quasicanonical quantum field theory, Phys. Rev. D 9 (1974), 373-386; see also the follow-up paper in Phys. Rev. D 19 (1979), 503-510.
After the usefulness of OPEs was generally recognized (and flourishes till today), the interest in field equations in itself waned, and work on it for other field theories seems virtually inexistent. In terms of OPEs, a field equation is just a nontrivial linear relation between the local fields of sufficiently low scaling dimension and their derivatives. This point of view is discussed (for the case of wick-rotated, Euclidean QFT) in a lecture by E. Witten, Composite operators and operator product expansion (Lecture 3, pp. 445ff of Volume 1 of Quantum Fields and Strings: A Course for Mathematicians); see in particular Section 3.7.
In general, one can define the family of local operators of a theory directly as the collection of all operators satisfying causal commutation relations among themselves and with an assumed family of Wightman fields (e.g., those that figure in the Lagrangian). That this gives a valid definition was proved by Borchers, hence the resulting family is called the Borchers class. Under natural conditions, the smeared fields from the Borchers class form a Hilbert space carrying a unitary representation of the Poincare group augmented by the dilation group. Decomposing it into its irreducible constituents defines the fields figuring in the OPE. They are naturally ordered by the scaling dimension, which is determined by the action of the dilation group. The interesting part of the OPE describes how the tensor product of two fields (in practice, of low scaling dimension) decompose into irreducibles.
In this global view, the complete field content is assumed given at the start, and the family of OPEs is essentially the multiplication law, Just as in a Lie algebra one assumes the generators as given and postulates the commutation laws. On the other hand, in a recursive construction, one starts with the most fundamental fields and obtains the others from the OPEs, just as one can start in so(3) with the two ladder operators $J_\pm$, define $[J_+,J_-]=:J_3$ (corresponding in the analogy to a particular OPE), and then require closure under commutation (which corresponds to the remaining OPEs).
The advantage of the global view is that it always works, while the recursive view requires a Lagrangian as a starting point. In addition, it requires lots of hocus pocus to turn the Lagrangian formalism into a well-defined renormalized perturbation theory - and even then it works only on the perturbative level. But perturbation theory works in the textbook fashion only in the absence of bound states and other infrared obstacles. Under these restrictions, there is (perturbatively, as Witten shows) a 1-1 correspondence between normally ordered products of the free field and interacting local field operators.
The latter can be obtained by composite operator renormalization. This process uses only a small part of the information in the OPEs. For example, to renormalize the composite operator $:\phi(x)^2:$, one needs to find a corresponding local field in the OPE for $\phi(x)\phi(y)$. Now, perturbatively, this OPE contains just one new field with a singular prefactor - which is declared to be the renormalization of $:\phi(x)^2:$ (defined only up to certain renormalization constants, as explained by Witten).
Note that many CFTs have no Lagrangian, so that composite operator renormalization makes no sense! (Even the notion of compositeness becomes problematic.) The global OPE view is also preferable in curved space; see work by Hollands, e.g., http://arxiv.org/abs/1105.3375 who starts with the OPE rather with a Lagrangian.
Also note that all papers and results mentioned (except for the properties of the Borchers class, whcih is a rigorous result) work at the relaxed level of rigor conventionally used in QFT. Rigorous results on field equations are available only in 2 and 3 space-time dimensions (and even there only sparingly).