Using $$p(\theta_1|y) \approx \frac{1}{S}\sum_{s=1}^Sp(\theta_1|\theta_2^{(s)}, y)$$ as an approximation to the marginal posterior density is called a Rao-Blackwellisation in the MCMC literature, started in Gelfand & Smith (1990) foundationalfoundational¹ paper, appropriately entitled "Sampling-Based Approaches to Calculating Marginal Densities". This density estimator is converging at a parametric speed, $\text{O}(1/\sqrt S)$.
When proposing the Gibbs sampling algorithm as a way to simulating from marginal densities and (hence) posterior distributions, Gelfand and Smith (1990) explicitly relate to the Rao-Blackwell theorem, as shown by the following quote²
...we consider the problem of calculating a final form of marginal density from the final sample produced by either the substitution or Gibbs sampling algorithms. Since for any estimated marginal the corresponding full conditional has been assumed available, efficient inference about the marginal should clearly be based on using this full conditional distribution. In the simplest case of two variables, this implies that $[X \mid Y]$ and the $y^{(i)}_j$'s $(j = 1, \ldots, m)$ should be used to make inferences about $[X]$, rather than imputing $X^{(i)}_j$ $(j = 1, \ldots, nm)$ and basing inference on these $X^{(i)}_j$'s. Intuitively, this follows, because to estimate $[X]$ using the $x^{(i)}_j$'s requires a kernel density estimate. Such an estimate ignores the known form $[X \mid Y]$ that is mixed to obtain $[X]$. The formal argument is essentially based on the Rao--Blackwell theorem. We sketch a proof in the context of the density estimator itself. If $X$ is a continuous p-dimensional random variable, consider any kernel density estimator of $[X]$ based on the $X^{(i)}_j$'s (e.g., see Devroye and Györfi, 1985) evaluated at $x_0$:$$\Delta^{(i)}_{x_0} = (1/h_m^p) \sum_{j=1}^m K[(X_0 - X^{(i)}_j)/h_m]$$ say, where $K$ is a bounded density on $\mathbb R^p$ and the sequence $\{h_m\}$ is such that as $m \rightarrow \infty$, $h_m \rightarrow 0$, whereas $mh_m \rightarrow \infty$. To simplify notation, set$$Q_{m,x_0}(X) = (1/h_m^p) K[(X - X^{(i)}_j)/h_m]$$ so that$$\Delta^{(i)}_{x_0} = (1/m) \sum_{j=1}^m Q_{m,x_0}(X_j^{(i)})$$ Define$$\gamma_{x_0}^i = (1/m) \sum_{j=1}^m \mathbb E[Q_{m,x_0}(X)\mid Y^{(i)}_j]$$ By our earlier theory, both $\Delta^{(i)}_{x_0}$ and $\gamma_{x_0}^i$ have the same expectation. By the Rao--Blackwell theorem,$$\text{var}\, \mathbb E[Q_{m,x_0}(X)] \mid Y) \le \text{var}\, Q_{m,x_0}(X)$$ and hence$$\text{MSE}(\gamma_{x_0}^i)\le\text{MSE}(\Delta^{(i)}_{x_0})$$ where MSE denotes the mean squared error of the estimate of $[X]$.
¹ "foundational" as it launched the MCMC revolution.
² The text has been retyped by me and may hence contains typos. The notations are those introduced by Gelfand and Smith (1990) and used for a while in the literature with $[X \mid Y]$ denoting the conditional density of $X$ given $Y$. The double indexation of the sequence is explained in the quote.