Using Environmentally Friendly Renormalization, we present an analytic calculation of the series for the renormalization constants that describe the equation of state for the $O(N)$ model in the whole critical region. The solution of the beta-function equation, for the running coupling to order two loops, exhibits crossover between the strong coupling fixed point, associated with the Goldstone modes, and the Wilson-Fisher fixed point. The Wilson functions $gamma_lambda$, $gamma_phi$ and $gamma_{phi^2}$, and thus the effective critical exponents associated with renormalization of the transverse vertex functions, also exhibit non-trivial crossover between these fixed points.