When is Schwartz's method for "integrating out" a field valid? In Schwartz's QFT book, heavy fields are often "integrated out" by simply solving their equations of motion formally (i.e. allowing things like $\Box^{-1}$) and plugging them back into the Lagrangian. 
This procedure doesn't make sense to me for several reasons. First, you're not allowed to do this even in classical mechanics -- if you just plug equations of motion back into the Lagrangian, you get incorrect dynamics for the remaining degrees of freedom! The exception is if the degrees of freedom eliminated are auxiliary, but the fields Schwartz integrates out aren't. Perhaps they're "almost auxiliary" because they're so heavy, but I don't know how to make that precise.
Second, I already know two separate methods for properly integrating out the field, i.e. computing a Wilsonian effective action from integrating over it in the path integral, and doing a matching-and-running continuum EFT calculation. It's not at all clear that either one will give the same result as Schwartz's procedure. 
What is the justification for Schwartz's procedure? Does it approximately match either of the "proper" procedures above, and if so, where does it differ?
 A: I think you are probably overthinking this. "Integrating out" a field in this context of Schwartz's book is just a way to eliminate constrained variables. It is most certainly not related to the methods of effective field theory you mentioned.
When there is redundancy in the Lagrangian, the equations of motion for some variables are not physical but are just equations of constraint. To apply canonical quantization, one must remove the constrained varaibles to obtain a Lagrangian (and Hamiltonian) in terms of only the physical degrees of freedom.
In a path integral approach, that can be done by first integrating over the constrained variables (hence the term "integrating out"). As you know from doing Gaussian integrations, the result of that is given by setting the exponential (the action) to its stationary value with respect to the constrained variables, which is equivalent to plugging in the equations of motion. 
Moreover, when there is redundancy, one has the freedom to choose which variables to remove. In the example in Schwartz's page 131, you can remove $\pi$, keeping the three degrees of freedom in $A_\mu$ or choose to keep $\pi$ and removing a degree of freedom from $A_\mu$.
For a more transparent example, consider the following Lagrangian:
$$
\mathcal{L}_0=\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}F_{\mu\nu}\big(\partial^\mu A^\nu-\partial^\nu A^\mu\big)+\frac{1}{2}m^2A_\mu A^\mu
$$
with $F_{\mu\nu}$ and $A_\mu$ independent variables. There are certainly redundancy and constrained variables in this Lagrangian. The equations of motion for the $F_{\mu\nu}$ , to no one's surprise, are:
$$
F_{\mu\nu}=\partial_\mu A_\nu-\partial_\nu A_\mu
$$
We can "integrate out" the $F_{\mu\nu}$ field by plugging in this EoM to obtain the familiar Lagrangian of a free massive vector field:
$$
\mathcal{L}_1=-\frac{1}{4}\big(\partial_\mu A_\nu-\partial_\nu A_\mu\big)^2+\frac{1}{2}m^2A_\mu A^\mu
$$
Alternatively, we can "integrate out" $A_0$ and $F_{ij}$ to obtain a Lagrangian (and the action in the path integral) in terms of the canonical variables $A_i$ and $F^{0i}$.
I believe that is all there is behind Schwartz's method for “integrating out” a field. (But full disclaimer, I haven't worked out the implications of Schwartz's example in page 131, so can't say further on that particular example).
A: The key is to do a systematic Taylor expansion of the full action in fluctuations $\eta$ of the fields $\phi$ that you want to integrate out, cf. above comment by Wakabaloola. [This recipe also applies to auxiliary fields, i.e. fields with no spacetime derivatives in the action.] 
As long as the quadratic part in $\eta$ is non-degenerate, one can formally do perturbation theory with the limitations that implies. This may lead to quantum corrections and one may loose (manifest) locality. 
The justification for leaving out quantum corrections will depend on the specific situation, e.g. a heavy mode would have a propagator with a large denominator, and so forth.
The 2 examples that OP explicitly mentions in the comments [(i) integrating out the Stueckelberg field on p. 131 and (ii) integrating out the $A$-field in QED in problem 3.9] are just Gaussian integrations [if we assume that the gauge-fixing terms are quadratic as well], which can be done exactly.
A: I find that Schwartz is using the laziest way to discuss the integration of a field. What is going on can be made much more transparent rather easily.
To simplify the discussion, I will assume that the theory is described by two scalar fields $\psi$ and $\phi$, each described by a (possibly interacting) Lagrangian $\mathcal L_\psi$ and  $\mathcal L_\phi$, which are coupled linearly, $\mathcal L_c=-g \,\psi\phi$. (I use Euclidian notations.  The sign in front of $g$ is introduced for later convenience. The discussion can be adapted easily for other coupling and other kinds of fields).
The partition function is then given by
$$ Z=\int \mathcal D \psi \mathcal D \phi e^{-S_\psi-S_\phi-S_c},$$
with $S_i$ the action corresponding to $\mathcal L_i$. In the present context, integrating out a field just means formally doing the corresponding functional integral.
Say we want to integrate out the field $\phi$. In all generality, the generating function of the connected correlation function of the field $\phi$ (not coupled to the field $\psi$) is defined by
$$
W_\phi[J]=\log\int \mathcal D \phi\,  e^{-S_\phi+\int dx J \phi},
$$
with $J$ a source. Then, we have
$$
Z=\int \mathcal D \psi e^{-S_{eff}},
$$
with $S_{eff}=S_\psi-W_\phi[g\,\psi]$. This is really (and only) what integrating out a field means.
Now, several remarks are in order.
1) If $S_\phi$ is quadratic, for instance $\mathcal L =\frac12 \phi(-\nabla^2+m^2)\phi$, one has $W_\phi[J]=\frac12 \int dx dy J(x) \frac1{-\nabla^2+m^2}J(y)$, up to a constant. This could also be obtained by replacing $\phi$ by the solution of the equation of motion $\phi_c$, which is such that $(-\nabla^2+m^2)\phi_c=J$. This is because for quadratic action, the semi-classical analysis (or saddle point approximation + gaussian fluctuations) is exact. We then recover Schwartz's discussion.
2) If $S_\phi$ is not quadratic, we can formally expand $W_\phi[g\, \psi]$ in $\psi$, which gives rise to new interaction terms. If $\phi$ is massive, its correlation functions are short-range, and, as long as we are interested in the physics of $\psi$ on distances much larger than the inverse mass, these interactions can be approximated as local (it corresponds to a gradient expansion of the correlation functions, for momenta much smaller than the mass). This is pretty much the effective field theory point of view (which allows for instance to record Fermi's theory of weak interaction from the standard model).
3) In a more Wilsonian point of view, when the RG scale (effective cut-off) gets smaller than the mass of $\phi$, this field decouples (since its contribution to the flow is of order cut-off/mass to some power), and it has been in essence integrated out in much the same way than what is done above.
4) If $\phi$ is massless (either it is free and massless, or interacting and massless, as it happens close at second-order phase transitions), then the effective interactions are highly non-local (in momentum space, they are non-analytic), which in practice makes $S_{eff}$ rather useless for any real calculations. The general rule of thumb is that integrating out massless fields is a very bad idea, they should be kept.
