Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
83,875 characters · 25 sections · 56 citation commands
Gaussian and Bootstrap Approximation for Matching-based Average Treatment Effect Estimators
Nearest neighbor matching estimators are non-parametric methods in causal inference used to estimate treatment effects by comparing treated and untreated units that are most similar in observable characteristics. This approach is commonly applied in observational studies where random assignment is not possible, necessitating statistical methods to estimate counterfactual outcomes (i.e., what would have happened to a treated unit if it had not received treatment). Matching involves pairing each treated unit with one or more untreated units that closely resemble it based on these characteristics, creating a comparison group that approximates the treated group but without the treatment. This process helps estimate the treatment effect by minimizing confounding due to observed differences. In particular, the aforementioned procedure is used to calculate the Average Treatment Effect (ATE) which reflects the treatment's effect across the entire population, including those who did and did not receive the treatment. Such ATE estimators have been widely used in various fields imbens2004nonparametric,morgan2006matching,rosenbaum2010design,stuart2010matching.
In two seminal works, abadie2006large,abadie2011bias proposed ATE estimators based on nearest neighbor matching and provided their first asymptotic analysis. In particular abadie2006large showed that their proposed estimator has non-negligible bias if the dimension is greater than one. As a remedy, abadie2011bias proposed a bias-correction procedure and established asymptotic properties (including asymptotic normality) under the crucial assumption that the number of matches is fixed. More recently, lin2023estimation established the consistency and asymptotic normality for matching-based ATE estimators allowing the number of matches to diverge with the number of observations. Through their analysis, they further showed that the estimator is doubly robust and semiparametrically efficient.
Our primary objective in this work is to derive precise Gaussian approximation results for nearest-neighbor matching-based ATE estimators. Existing asymptotic normality results, commonly used for constructing confidence intervals, have significant limitations. In particular, they provide no information on when Gaussianity “kicks in”, making the resulting confidence intervals valid only asymptotically. In addition, key parameters, such as the number of matches and the balance between treatment groups, become obscured in asymptotics. For instance, the rate at which the number of matches increases with the sample size directly impacts the accuracy of the Gaussian approximation, and consequently, the validity of the constructed confidence intervals. Establishing a fine-grained Gaussian approximation bound in this context would allow us to quantify these effects non-asymptotically, improving the reliability of the confidence intervals for the practitioner.
To this end, we introduce a novel approach for quantifying the accuracy of Gaussian approximations in matching-based ATE estimators. Our method combines stabilization theory with the Malliavin-Stein method. A key observation for our analysis is that the leading term in the bias-corrected ATE estimator (specifically, the term \( E_n \) defined at (ref)) exhibits a local geometric property termed as stabilization, as illustrated in detail in Section (ref) later. To the best of our knowledge, only shi2022berry has previously applied Stein's method in the context of causal inference problems. However, their work focused on leveraging classical results from the Stein's method literature on permutational statistics to derive Berry–Esseen bounds for design-based causal inference.
By leveraging the aforementioned insight and by refining the Gaussian approximation results in lachieze2019normal (see our Theorem (ref) specifically), we present in Theorem (ref) the first Gaussian approximation result for the bias-corrected matching-based ATE estimator. This result explicitly quantifies the approximation’s accuracy based on key parameters, including the number of matches and the treatment group balance. For example, a consequence of our result in the one-dimensional setting with balanced data, gives the Gaussian approximation upper bound $M^5n^{-\frac{1}{2}}+M^{-\frac{1}{2}}$, where $M$ is the number of nearest neighbor matches and $n$ is the number of observations; see Corollary (ref) and (ref) for details. Similarly, in Theorem (ref), we establish Gaussian approximation results for the \(\phi\)-transformed rank-based ATE estimator proposed and analyzed in cattaneo2023rosenbaum.
On a more technical note, another contribution of our work is a refinement on Theorems B.3 and B.4 in lin2023estimation. In this context, we derive a mathematically rigorous, fully non-asymptotic bound for the estimation error of the nearest-neighbor-based density ratio, as presented in Lemma (ref). In contrast, the error bounds in lin2023estimation included asymptotic simplifications tailored to their purpose without providing a fully non-asymptotic expression.
As an application of our main results, we analyze a multiplier bootstrap method to estimate the limiting distribution and establish bootstrap approximation rates for both the covariate-based and rank-based ATE estimators in Theorem (ref). Notably, our bootstrap approximation results allow the number of matches to increase with the sample size. This contrasts with the findings in abadie2008failure, which demonstrate that the naive bootstrap procedure is inconsistent when the number of matches remains fixed. Our results, on the other hand, specify the exact rate at which the number of matches can diverge with the sample size for the multiplier bootstrap method to remain consistent.
During the final stages of preparing this manuscript, we became aware of a concurrent work by lin2024consistency that establishes consistency results for the naive bootstrap procedure when the number of matches is allowed to grow with the sample size. We would like to point out three significant distinctions between this and our current work: (i) they employ the naive bootstrap procedure, similar to abadie2008failure, (ii) they do not provide rates for bootstrap approximation, and (iii) their proof techniques differ fundamentally, being more canonical, whereas our approach relies on the Malliavin-Stein method and stabilization theory.
Throughout the paper, we will use the following notation and conventions.
Following the framework by abadie2006large,abadie2011bias, we are interested in estimating the average treatment effect (ATE) on outcomes in a binary treatment experiment. Consider $(X,Y,D) \in \mathbb{X} \times \mathbb{R} \times \{0,1 \}$, where $X$ corresponds to the unit, $Y$ corresponds to the response, and $D$ is the binary variable (possibly dependent on $X$) such that $D=1$ and $D=0$ corresponds to when the unit $X$ belongs to the treatment and the control groups, respectively. Under the framework of rubin1974estimating (also see rosenbaum1995observationalstudies,imbens2009recent), $X$ has two potential outcomes, $Y(0)$ and $Y(1)$, depending on whether $D=0$ or $D=1$, but we observe only one of them. In other words, $Y=DY(1)+(1-D)Y(0)$. The central goal is to estimate the population ATE defined as
given observations $\{(X_i,Y_i,D_i)\}_{i=1}^{n}$ that are assumed to be independent and identically distributed copies of $(X,Y,D)$.
We first discuss covariate-based matching estimators for ATE. Let $n_1=\sum_{i=1}^{n}D_i$ and $n_0=\sum_{i=1}^{n}(1-D_i) = n-n_1$ denote the number of treated and control units, respectively. Note here that while $n$ is a deterministic variable, $n_0$ and $n_1$ are random variables depending on the specific instance of $\{D_i\}_{i=1}^n$.
For a point collection $\tilde{\nu}_n=\{\tilde{x}_{i} \in \widetilde{{\mathbb X}} : i \in [n]\}$, $\omega\in \{0,1\}$ and an integer $M>0$, we denote by $\mathcal{J}^{\omega}_{M}(i, \tilde{\nu}_n)$ the index set of $M$-NNs of ${x}_i$ within the set $\nu_n^\omega :=\{x \in {\mathbb X} : (x,\omega) \in \tilde{\nu}_n\}$, namely, the set of all indices $j\in[n]$ such that $d_j=\omega$ and
Since only one of the potential outcomes $Y_i(0)$ and $Y_{i}(1)$ is observed, we impute the missing potential outcomes (see abadie2006large) via nearest neighbor (NN) matching as
and
Then, the matching-based estimator for ATE is defined as the empirical counterpart of (ref) based on the $\hat Y_i(0)$ and $\hat Y_i(1)$, i.e., $$\hat \tau_M := \frac{1}{n}\sum_{i=1}^n \hat Y_i(1) - \frac{1}{n}\sum_{i=1}^n \hat Y_i(0). $$
Let $K^{\omega}_{M}(i,\tilde{\nu}_n)$ be the matched times for unit $i$ with $d_i=\omega \in \{0,1\}$, i.e.,
In other words, $K^{\omega}_{M}(i,\tilde{\nu}_n)$ denotes the total number of units $j$ in $\tilde{\nu}_n$ with $d_j=1-\omega$ (i.e. with the opposite label to $x_i$) such that $x_i$ is one of its $M$-NNs in $\nu_n^\omega$. Then $\hat{\tau}_{M}$ can be further expanded as
The estimator $\hat{\tau}_{M}$ however suffers from an asymptotically non-negligible bias when the dimension $m$ is strictly larger than $1$ (abadie2006large). To circumvent this, in a follow-up work, abadie2011bias proposed a bias-corrected version $\hat{\tau}_{M}^{bc}$ (see (ref) below), defined as follows. Consider the regression model $Y_i=\mu_{D_i}(X_i)+\bm{\varepsilon}_{i}$, $i\in[n]$, define conditional expectations
and let $\hat{\mu}_{0}(x)$ and $\hat{\mu}_{1}(x)$, respectively, be corresponding regression estimators. The estimator $\hat{\tau}_{M}^{bc}$ is then obtained by replacing the $\hat Y_i(0)$ and $\hat Y_i(1)$ in $\hat \tau_M$ by the corrected quantities $\hat Y_i(0) + \hat \mu_0(X_i) - \hat \mu_1(X_i)$ and $\hat Y_i(1) + \hat \mu_1(X_i) - \hat \mu_0(X_i),$ respectively. Writing the residuals as
and denoting the regression based estimator of the population ATE as
the bias corrected estimator of ATE can be expressed as
lin2023estimation showed that the estimator (ref) is indeed doubly robust. Consequently, it should also enjoy all the desirable properties of doubly robust estimators scharfstein1999adjusting,bang2005doubly. Moreover, following lin2023estimation, the estimator (ref) can be conveniently decomposed as
where $\bm{\varepsilon}_i$ are the errors in our regression model, and
with $j_{m}^{\omega}(i,\widetilde{\mathcal{X}}_n)$ denoting the $m$-th nearest neighbor of the point $X_i$ in $\{X_j:D_j=\omega\}_{j=1}^{n}$ for $m\in[n]$. Here, $E_n$ can be viewed as the main contributing term and $(B_{M}-\hat{B}_{M})$ as the bias term. In Theorem (ref) and Corollary (ref) in Section (ref), using stabilization theory and Malliavin-Stein method, we provide a quantitative estimate for the error in the Gaussian approximation of $\hat{\tau}_{M}^{\text{bc}}$ (appropriately centered and scaled).
The above covariate-based matching uses the Euclidean distance for determining the nearest neighbor matching. It may however exhibit sensitivity to alterations in scale and to the existence of outliers or heavy-tailed distributions. Also, in practice distance metrics are often derived from a `standardized' representation of the data, and the selection of a metric is an important factor in causal inference because different metrics can lead to different conclusions (rosenbaum2010design, chapter 9). Therefore, in two influential contributions, rosenbaum2005exact,rosenbaum2010design advocated for using the distances between component-wise ranks, instead of the original data, to measure covariate similarity when constructing matching estimators of average treatment effects. This approach is called Rosenbaum’s rank-based matching estimator for ATE.
Compared to the covariate-based matching ATE estimator, Rosenbaum’s rank-based matching estimator is obtained by replacing the original values of the $X_i$'s with their component-wise ranks when performing nearest neighbor matching. The detailed construction is as follows.
Step 1. Write $X_i=(X_{i,1},\ldots,X_{i,m})^{T}$ for $i\in [n]$. Define the vector of the marginal empirical cumulative distribution functions, $\widehat{\mathbf{F}}_{n}:\mathbb{R}^{m}\rightarrow [0,1]^{m}$, as follows: for any $x=(x_1,\ldots,x_m)^{T}\in\mathbb{R}^{m}$,
For each $i\in[n]$, define $\hat{L}_{i}:=\widehat{\mathbf{F}}_{n}(X_i)$ and note that for $k \in [m]$, the $k$-th component of $n\hat{L}_{i}$ is the corresponding rank of $X_{i,k}$ among $\{X_{j,k}\}_{j=1}^{n}$ (with ties broken arbitrarily). Also, let $\mathbf{F}:\mathbb{R}^{m}\rightarrow [0,1]^{m}$ be the vector of marginal population cumulative distribution functions, i.e. for any $x=(x_1,\ldots,x_m)^{T}\in\mathbb{R}^{m}$,
Write $L=\mathbf{F}(X)$ with $L_i:=\mathbf{F}(X_i)$, for $i\in[n]$.
Step 2. Similar to the covariate-based matching, regression adjustment is employed to correct the bias. Let $\hat{\mu}_{r,0}(\cdot)$ and $\hat{\mu}_{r,1}(\cdot)$ be estimators of the conditional means
respectively.
Step 3. The rank-based ATE estimator $\hat{\tau}_{r,M}^{bc}$ now is constructed by applying bias-correction and matching to $\{(\hat{L}_i,D_i,Y_i)\}_{i=1}^{n}$: Let $\mathcal{J}_{r,M}(i)$ denote the index set of $M$-NNs of $\hat{L}_{i}$ in $\{\hat{L}_{j}:D_{j}=1-D_{i}\}_{j=1}^{n}$ with ties broken in some arbitrary way. The rank-based ATE estimator $\hat{\tau}_{r,M}^{bc}$ is then defined as
where, for $\omega\in\{0,1\}$,
This can be further generalized by considering a functional transform of the data, as considered in cattaneo2023rosenbaum. For the sake of completeness, we also explain this general case below. For $\omega\in\{0,1\}$, consider functions $\phi_{\omega}:\mathbb{X}\rightarrow\mathbb{X}_{\phi}$ with ${X}_{\phi}\subseteq \mathbb{R}^{m'}$ for some $m'\in \mathbb{N}$. Note here that $m'$ can indeed be different from $m$. Then, for possibly unknown $\phi_{\omega}$, let $\hat{\phi}_{0}$ and $\hat{\phi}_{1}$ be generic estimators based on the sample $\{(X_i,D_i,Y_i)\}_{i=1}^{n},$ and define
Note that when $\phi_0=\phi_1=\mathbf{F}$ and $\hat{\phi}_{0}=\hat{\phi}_{1}=\hat{\mathbf{F}}_n$, it recovers $L$ and $\hat{L}_i$ from Step 1 above. Let $\mathcal{J}_{\phi,M}(i)$ represent the index set of $M$-NN matches of $\hat{L}_{\phi,1-D_i,i}$ in $\{\hat{L}_{\phi,1-D_i,j}:D_j=1-D_i\}_{j=1}^{n}$ with ties broken in an arbitrary way. In other words, for determining the nearest neighbors, this approach measures the similarity based on the Euclidean distance between transformed data points with the transformation function possibly also needing to be learned from the same data. Let $K_{\phi}(i)$ stand for the number of matched times for the unit $i$, i.e.,
Moreover, let $\hat{\mu}_{\phi,\omega}(l)$ be mappings from $\mathbb{X}_{\phi}$ to $\mathbb{R}$ that estimate the conditional means
The general $\phi$-transformed rank-based bias-corrected matching estimator $\hat{\tau}_{\phi,M}^{bc}$ is then given by
where, for $\omega\in\{0,1\}$,
Note again that when $\phi_0=\phi_1=\mathbf{F}$ and $\hat{\phi}_{0}=\hat{\phi}_{1}=\hat{\mathbf{F}}_n$, (ref) indeed recovers (ref). In the general setting, we provide a Gaussian approximation bound for $\hat{\tau}_{\phi,M}^{bc}$ in Theorem (ref) and Corollary (ref) in Section (ref).
In this paper, we work under the standard assumptions put forward in prior works by abadie2011bias, lin2023estimation and cattaneo2023rosenbaum. We refer the reader to these works for further motivations for these assumptions. We do however make a few minor modifications to the assumptions, as required for our stabilization-based Gaussian approximation techniques.
For the covariate-based matching, we assume the following two sets of conditions.
(Data Distribution)
(Regression functions)
To accommodate the changes in the rank-based matching, some of the assumptions above need appropriate adjustments for the $\phi$-transformation; see cattaneo2023rosenbaum for more details.
(Data Distribution)
(Regression functions)
Note again that similar to the case of covariate bases matching, compared to the assumptions in cattaneo2023rosenbaum, we do not require their assumption 5.5 (i) that $ \mathbb{E} [(Y(\omega)-\mu_{\phi,\omega}(U_{\phi,\omega}))^2 | U_{\phi,\omega}=u]$ is uniformly bounded away from zero for almost all $(u,\omega)$, which is needed there to apply the Lindeberg-Feller CLT, as well as the uniform boundedness of the $(2+\kappa)$-th moment in their assumption 5.5 (ii).
We now present our main results on Gaussian approximation bounds for the matching-based and rank-based ATE estimators. All our bounds are stated for the Kolmogorov distance $\mathsf{d}_{K}(\cdot,\cdot)$, which, for two real-valued random variables $Z_1,Z_2$, is defined as
For $r_{0}:=r_{\text{ratio}}\vee r_{\text{ratio},\phi}$ (with $r_{\text{ratio}}$ and $r_{\text{ratio},\phi}$ from Assumptions A.3 and C.4, respectively), define the quantities
For more details about the above terms, see the discussion following Theorem (ref). All our results hold as long as $n \geq 9$. This is due to the fact that a certain tail bound on the radius of stabilization, required for the stabilization technique, holds as long as $n \geq 9$; see Section (ref) and Lemma (ref) for additional details.
The above Gaussian approximation bound consists of three parts. The term $B_1$ corresponds to the Gaussian approximation bound for $E_n$ in (ref), centered at the true ATE and scaled by the sample variance. Similar to the classical Berry-Essen Theorem (see, for example, (ref)), the polynomials involving $M$ in the numerators are from Assumption B.2 in Section (ref), and the denominator $n^{1/2}$ in $B_1$ corresponds to a variance lower bound.
The term $B_2$ arises from the bias correction for $(B_M-\hat{B}_M)$. This can be further improved by assuming existence of higher order moments in Assumption B.4 as in Section (ref) instead of just $L^2$ moments. To ensure that the bias term $B_2\rightarrow 0$, one could pick $M\lesssim n^{\iota}$ with
where $\epsilon_{\mu}$ and $\gamma_{l}$'s are defined in Assumption B.4 in Section (ref).
The term $B_3$ relies on the convergence rate of the sample variance to its limiting variance $\sigma^2$. For the sake of generality, we keep track of the data balance parameter $\eta$ in Assumption A.3 in (ref), and have made minimal assumptions on the relationship between $\eta$ and $n$ (or $M$); see also Lemma (ref). It is obtained via a modified and rigorous non-asymptotic convergence rate argument for the density ratio estimation in contrast to the crude asymptotic arguments in lin2023estimation; see Lemma (ref). This also necessitates the addition of the non-asymptotic error terms $\delta_{H_1}$-$\delta_{H_3}$. Moreover, the assumptions that there exists positive constants $C_0,C_1$ such that $M\le C_0n\eta$ and $C_1\le n\eta^2$ are mild in the sense that they are also required for $B_3$ and $\delta_{H_1}$ to tend to zero.
The following corollary to Theorem (ref) provides a user-friendly bound under some additional mild assumptions that ensure that all the error terms involving $\delta_{H_1}$-$\delta_{H_3}$ are negligible compared to the rest of the summands.
A further simplified bound is provided in (ref) later, where we also compare covariate and rank-based ATE estimators in the univariate setting. We now make some remarks regarding Theorem (ref) and Corollary (ref).
The doubly robust estimator of ATE considered in lin2023estimation actually uses a $K$-fold random partition of the data and averages the estimation on each subset to output a final estimator. We emphasize here that the stabilization technique could also be applied in a similar way as for the bound $B_1$, since both of these estimators use nearest neighbor matching. Carrying out this exercise is left as a future work. Furthermore, although both Theorem (ref) and Corollary (ref) are stated in the context of the Euclidean setting, the stabilization technique introduced later in Section (ref) used to obtain $B_1$ is valid for general metric spaces. Particularly, it can also be applied in the $m$-dimensional manifold setting, for example, as in penrose2013limit.
Similar to Theorem (ref), we aim to keep track of $\eta$, which measures the balance of the data. Consequently, in practice one should also be careful while picking $M$ and $n$ in relation to the speed of decay of $\eta$ to $0$. Different from Theorem (ref), the Gaussian approximation bound for the $\phi$-transformation based estimator $\hat{\tau}_{\phi,M}^{bc}$ depends on the convergence of the estimator for $\phi_{\omega}$, that is $\hat{\phi}_{\omega}$, as appearing in $B_5$ and $B_6$. Comparing to the asymptotic normality result in cattaneo2023rosenbaum, for $B_5$, they assumed a Donsker-type condition: for any $\omega\in\{0,1\}$ and any $\epsilon_e>0$,
which is an asymptotic property. In our case, such a term directly appears in our non-asymptotic bound. Moreover, for $B_6$, cattaneo2023rosenbaum requires that
which does not offer any rates of convergence. We again have this difference featuring in our bound instead.
From Theorem (ref), it is also seen that the Gaussian approximation bound depends on the choice of the transformation $\phi$ (including the dimension of the embedded space, $m'$) and its estimation $\hat{\phi}$. This selection is crucial as different transformations can lead to different conclusions for ATE estimation. We refer interested readers to (rosenbaum2010design, chapter 9) and rosenbaum2005exact,rosenbaum2010design for more details on the influence of choices for the transformation $\phi$. This general issue is beyond the scope of the current paper which focuses on Gaussian approximation bounds. However, focusing on the particular choice of $\phi$-transformation as the cumulative distribution function (CDF), we present the corresponding simplified result in Corollary (ref). For the sake of simplicity, we make further mild assumptions that $M^{-1}\log n=o(1), n^{-1}M\log n=o(1)$ and $\eta$ bounded away from $0$.
It is instructive to compare the rates with covariate and CDF-base rank ATE. From Corollaries (ref) and (ref), we note that $B_1'=B_4'$ and the additional term $n^{-1/4}$ in $B_5'$ is dominated by $n^{-1/6}$ in $B_6'$, thus not affecting the overall rate. Then, the terms that actually affect the overall rate are $B_3'$ and $B_6'$. When $m\geq 3$, the last term in $B_6'$ decays fast leading to the same overall rate with or without a CDF transformation. On the contrary, for example, when $m=1$, the last summand in $B_6'$ (which equals $M^{-1/6}$) dominates $B_5'$ as well as every other summand in $B_6'$. Hence, we have (in the case of CDF-based rank ATE and $m=1$)
On the other hand, for the case without a CDF transformation and $m=1$ (i.e., covariate-based rank ATE), it holds that
We can then see that when $M^{\frac{40}{8+p}}n^{-\frac{1}{2}}\lesssim M^{-1/2}$, we expect a strictly worse rate with a CDF transformation.
In the context of matching based estimators, the asymptotic normality results from abadie2006large,lin2023estimation could be used to construct confidence intervals for the ATE parameter $\tau$. Specifically, one claims that $\hat{\tau}_{M}^{\text{bc}}\pm z_{1-\frac{\alpha}{2}}\frac{\sigma}{\sqrt{n}}$, where $z_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal distribution, provides a $1-\alpha$ confidence interval in this context. However, there are two main shortcomings of such a claim: (i) The validity of the obtained confidence intervals holds only asymptotically as $n\to\infty$, and (ii) The limiting standard deviation $\sigma$ has to be consistently estimated (see abadie2006large).
By the definition of the Kolmogorov metric, a direct application of Theorem (ref) (or Corollary (ref)) yields that, for any $0<\alpha<1$,
While this avoids the shortcoming in point (i) above, the unknown $\sigma$ makes it still impractical. To overcome both shortcomings (i) and (ii) simultaneously, we now provide an application of the Gaussian approximation results developed in the previous sections, for obtaining confidence intervals that are valid in a non-asymptotic sense.
Bootstrap serves as one of the most important inferential techniques for non-parametric statistical analysis. However, abadie2008failure provided an example showing that the naive bootstrap (i.e., resampling from the empirical distribution of the observations) fails to provide an asymptotically valid standard error and quantiles for a matching-based ATE estimators. In addition, they argue that the main reason for this failure is that the naive bootstrap fails to reproduce the distribution of the number of matched times $K^{\omega}_{M}(i,\tilde{\nu}_n)$ in (ref), with fixed $M$.
Later, otsu2017bootstrap proposed an alternative bootstrap called the weighted bootstrap to overcome this difficulty and showed validity of their procedure. They only considered the setting with fixed $M$, i.e., when the number of the nearest neighbors stays fixed. Other methodological studies include works by abadie2022robust, walsh2023nearest and kosko2024fast. In the following, we adopt a multiplier bootstrap (or wild bootstrap) to not only bootstrap the distribution of the statistic $\sqrt{n}(\hat{\tau}_{M}^{bc}-\tau)$, but also provide a bound for the approximation accuracy in terms of the Kolmogorov distance, while allowing $M$ to diverge with $n$.
We now describe the details of the multiplier bootstrap procedure and present the approximation bounds. Recall the definition of the bias-corrected ATE estimator in (ref), given by
Below, for notational ease, we denote $$\Delta\hat{\mu}(X_i):=\hat{\mu}_1(X_i)-\hat{\mu}_0(X_i),\; i\in[n]\quad \text{ and }\quad \widebar{\Delta\hat{\mu}}:=\frac{1}{n}\sum_{i=1}^{n}(\hat{\mu}_1(X_i)-\hat{\mu}_0(X_i)). $$ The multiplier bootstrap is constructed via the following steps:
Turning to the rank-based ATE estimator, note that by definition (ref), we can rewrite the estimator in (ref) as
By replacing $X_i$ with the $\phi$-transformed sample $\hat{L}_{\phi,D_i,i}$, we can analogously construct the bootstrapped version of the rank-based ATE estimators, which we denote by $\hat\tau_{\phi,M}^{\text{boot}}$. The following result provides rate of convergence for the above multiplier bootstrapping procedures.
Similarly, under the assumptions of there exists a constant $C>0$ not depending on $n,M,\eta$ or $p$ such that with probability at least
where $B_4,B_5$ and $B_6$ are given in Theorem (ref). \fi
To interpret the bounds in Theorem (ref), first note that the terms appearing in the numerators are similar to those appearing in Theorems (ref) and (ref), except $\mathsf{E}_1$ and $\mathsf{E}_2$, which estimate the quality of the approximation of $\mu_\omega, \mu_{\phi,\omega}$ and $\phi_\omega$. The quantity $L$ appearing in the denominators serve as a lower bound of the conditional variance of the bootstrap estimators. For its eventual positivity, one needs that $\mathsf{E}_1$ and $\mathsf{E}_2$ tend to zero, i.e., the construction of the regression estimator $\hat{\mu}$ and the transformation estimator $\hat{\phi}$ need to be consistent with high probability. One then should expect with high probability that $L(\mu,\hat{\mu},n),L(\mu,\hat{\mu},\phi,\hat{\phi},n)\ge L_{l}>0$ for a strictly positive constant $L_{l}>0$ for $n$ large enough, which can be explicitly determined from the convergence rates of $\hat \mu$ and $\hat \phi$. The bound, in the covariate-based case for instance, then simplifies to the following: with high probability one has
for some finite constant $C>0$.
In the proof of Theorem (ref), we use the limiting Gaussian distribution as the bridge to bound the distributions between the bootstrap and the original estimators. While showing the consistency of bootstrap does not necessarily need a Gaussian limit, our aim here is to highlight a general procedure to derive the rate of convergence of bootstrap via the rate of convergence for the Gaussian approximation using stabilization theory (see Section (ref)), which could potentially be useful for many geometric statistics having some form of local dependency structure.
\end{Theorem} \fi
In this section, we discuss briefly our approach to prove the results in Section (ref) based on the notion of stabilization and Malliavin-Stein method. Our proof of Theorem (ref) relies on Theorem (ref), that we introduce here. It is a refinement of the seminal work of lachieze2019normal, providing a quantitative bound for Gaussian approximation of Poisson functionals, and serves as the key step towards obtaining Theorems (ref) and (ref). Before we can state the result, we first explicitly introduce the setting for functionals of point processes, and the notion add-one cost operators acting on such functionals.
Recall that $\widecheck{{\mathbb X}}:={\mathbb X}\times \{0,1\}\times \mathbb{R}$, with ${\mathbb X} \subseteq \mathbb{R}^m$. As before, we will often denote an $\widecheck{{\mathbb X}}$-valued random vector $(X,D,\bm{\varepsilon})$ or a the nonrandom vector $(x,d,\varepsilon) \in \widecheck{{\mathbb X}}$ by $\widecheck{X}$ and $\check{x}$, respectively. While the concepts and results in this section can be extended to more general spaces, we will stick to the space $\widecheck{{\mathbb X}}$ related to our ATE estimation (see Section (ref)). We refer to lachieze2019normal for results in more general spaces. Let $(\widecheck{{\mathbb X}},\mathcal{F})$ be a measure space with the joint probability measure $\widecheck{\mathbb{Q}}$ of $(X,D,\bm{\varepsilon})$, with $\mathbb{Q}$ denoting the marginal distribution of $X$. We also define a semi-metric $d_S(\cdot,\cdot)$ on $\widecheck{{\mathbb X}}$ as $$ d_S(\check{x},\check{y})= d(x,y), \quad \check{x},\check{y} \in \widecheck{{\mathbb X}}, $$ where $d(\cdot,\cdot)$ denotes the Euclidean metric on $\mathbb{R}^m$. Let N be the set of $\sigma$-finite counting measures on $(\Omega,\mathcal{F})$, which can be interpreted as point configurations in $\Omega$. The set N is equipped with the smallest $\sigma$-field $\mathscr{N}$ such that the maps $m_{A}:\textbf{N}\rightarrow \mathbb{N}\cup\{0,\infty\},\mathcal{M}\mapsto\mathcal{M}(A)$ are measurable for all $A\in\mathcal{F}$. A point process is a random element in N. For $\mu \in \textbf{N}$, we write $x\in\mu$ if $\mu(\{x\})\ge 1$. Denote by $\textbf{F}(\textbf{N})$ the class of all measurable functions $f:\textbf{N}\rightarrow \mathbb{R}$, and by $L^{0}(\widecheck{{\mathbb X}}):=L^{0}(\widecheck{{\mathbb X}},\mathcal{F})$ the class of all real-valued, measurable functions $F$ on $\Omega$. Note that, as $\mathcal{F}$ is the completion of $\sigma(\mu)$, each $F\in L^{0}(\widecheck{{\mathbb X}})$ can be written as $F=f(\mu)$ for some measurable function $f\in \textbf{F}(\textbf{N})$. Such a mapping $f$, called a representative of $F$, is $\widecheck{\mathbb{Q}}\circ \mu^{-1}$-a.s. uniquely defined. In order to simplify the discussion, we make the following convention: whenever a general function $F$ is introduced, we will select one of its representatives and denote such a representative mapping by the same symbol $F$.
Notice that to apply Theorem (ref), one needs to find a lower bound for the variance of $E_n$, which we present in the following result, along with quantitative bounds for the variance approximation. Below, we denote $\sigma_{\bm{\varepsilon}}^2 = \operatorname{Var}(\bm{\varepsilon})$.
We note here that the lower bound in Lemma (ref) from our Assumptions A(1) and B(2)-(3) in Section (ref) implies that $n \operatorname{Var} E_n \ge C$ for some constant $C>0$.
While a result such as Theorem (ref) is a very powerful first step in providing (optimal) rates for Gaussian convergence of functionals of Binomial processes, often the integrals appearing in the bound involving the functions $\Gamma_n$ and $\psi_{n}$ are very difficult to directly bound. An assumption on the functional $F$ that can be very effectively used to simplify such computations, is that the $F$ can be expressed as a sum of local contributions from each point of the underlying process. This phenomenon is often referred to as stabilization in the relevant literature. There is an ever-growing literature on the application of stabilization in combination with result such as Theorem (ref) arising from the Malliavin-Stein method. We refer to works by lachieze2019normal,lachieze2022quantitative,shi2024flexible,bhattacharjee2022gaussian,shi2024multivariate for additional details and its applications in various statistical problems. In the following section, we describe how stabilization helps us to obtain the Gaussian approximation results in Section (ref) from Theorem (ref).
Let $\widecheck{X}_1,\widecheck{X}_2,...,\widecheck{X}_n$ be i.i.d.\ random variables sampled from $\widecheck{\mathbb{Q}}$. The binomial point process $\widecheck{\mathcal{X}}_n$ associated with $\{\widecheck{X}_1,\widecheck{X}_2,...,\widecheck{X}_n\}$ is defined as $\widecheck{\mathcal{X}}_n:=\sum_{i=1}^{n}\delta_{\widecheck{X}_1},$ where $\delta$ is the Dirac measure. Given its association with the i.i.d.\ sample $\{\widecheck{X}_i\}_{i=1}^{n}$, with a slight abuse of notation, we will often interchangeably use the binomial process $\widecheck{\mathcal{X}_n}$ and the i.i.d.\ sample $\{\widecheck{X}_1\}_{i=1}^{n}$. In this paper, we concern ourselves with functionals $F_{n}$ of the binomial process $\widecheck{\mathcal{X}_n}$ that can be represented as a sum of the form
where $f_{n}$ is called a {\em score function}. For $i \in [n]$, when contributions $f_n(\widecheck{X}_i,\widecheck{\mathcal{X}}_n)$ are sufficiently `local', one can expect a Gaussian limit for the sequence of functionals $F_n$ as $n \to \infty$.
For $n\ge 1$, the score function $f_n$ is said to be stabilizing if there exists an almost surely finite random variable $R_{n}:\widecheck{{\mathbb X}}\times \mathbf{N}\rightarrow \mathbb{R}_{+}$ such that for all $\mu\in\mathbf{N}$, $\check{x}\in\mu$, and all finite $\widecheck{\mathcal{A}}\subset \widecheck{{\mathbb X}}$ with $|\widecheck{\mathcal{A}}|\le 7$,
where the random variable $R_{n}$ is called the radius of stabilization.
Informally speaking, the above definition states that the value of the score function $f_{n}$ at a point $\check{x}$ is completely determined by all those point in the configuration that lie in the ball centered at $\check{x}$ with radius $R_{n}$. Note that in (ref), the score functions $f_n$ are possibly dependent on all the points $\widecheck{\mathcal{X}}_n$. Stabilization is a geometrical localization of the dependence between all the score functions.
We further need the radius of stabilization $R_n$ to satisfy certain tail decay conditions. It is said to decay exponentially if there exist constants $C_1,C_2>0$ such that for $\check{x}\in\widecheck{{\mathbb X}}$, $n\ge 9$ and $r\ge 0$,
The presence of $n-8$ in the above definition is for certain technical reasons, see lachieze2019normal for more details. An exponential decay as above for the tail probability of the radius of stabilization ensures that the dependence between score functions remain relatively local. Instead of exponentially decaying, penrose2005normal also proposed the polynomially decaying condition.
We now connect the terminology above with the matching-based ATE estimator notation introduced in Section (ref). The same, of course, also applies to the rank-based ATE estimator in (ref).
Recall the definition of $\hat{\tau}_{M}^{bc}$ in (ref) given by
Our focus is the main term $E_n.$ Write
Since $\bm{\varepsilon}$ can possibly depend on the covariate pair $(X,D)$, we are naturally led to consider a binomial point process in the aforementioned space $\widecheck{{\mathbb X}}={\mathbb X}\times\{0,1\}\times \mathbb{R}$. Also, as mentioned before, we see that $\widecheck{\mathcal{X}}_n=\{(\widecheck{X}_{i},i\in[n])\}$, as the collection of the triplets $\{(X_i,D_{i},\bm{\varepsilon}_{i})\}_{i=1}^{n}$, can be viewed as a binomial point process of size $n$ in the space $\widecheck{{\mathbb X}}$ distributed as $(X,D,\bm{\varepsilon})$. Note that the collection $\{(X_i,D_{i},\bm{\varepsilon}_{i})\}_{i=1}^{n}$ and the sample $\{(X_i,Y_i,D_i)\}_{i=1}^{n}$ are linked through the true regression functions (conditional means) $\mu_0(x)$ and $\mu_1(x)$, in particular, one has $\bm{\varepsilon}_i=Y_i - \mu_{D_i}(X_i)$ for $i \in [n]$.
Therefore, $nE_n$ can be viewed as a functional of the binomial process $\widecheck{\mathcal{X}}_n$ and be represented as a sum of score functions as
where the score function $\xi_n$ for $\check{\nu}_k =\{\check{x}_i\}_{i=1}^k \in \mathbf{N}, k \in \mathbb{N}$ is given by
Note that all terms in the score $\xi_n(\widecheck{X}_j, \widecheck{\mathcal{X}}_n)$ corresponding to the $j$-th sample $\widecheck{X}_j$ is determined by $\widecheck{X}_j$ except for $$ \sum_{i=1,D_i=1-D_j}^{n}\bm{\varepsilon}_i \mathds{1}(i\in \mathcal{J}^{D_i}_{M}(j,\widetilde{\mathcal{X}}_n)) = \sum_{i=1}^{n}\bm{\varepsilon}_i \mathds{1}(i\in \mathcal{J}^{1-D_j}_{M}(j,\widetilde{\mathcal{X}}_n)) . $$ This is a function of $\widecheck{X}_j$ and all those points in $\widecheck{\mathcal{X}}_n$ that are $M$-NNs of $\widecheck{X}_j$ (in the metric $d_S$) when considering only those points in $\widecheck{\mathcal{X}}_n^{1-D_j}$, where for a point collection $\check{\nu}_n=\{\check{x}_i\}_{i=1}^n$, we write
So, it is straightforward to see that for $j \in [n]$, the score function $\xi(\widecheck{X}_j, \widecheck{\mathcal{X}}_n)$ is stabilizing with radius of stabilization (in the metric $d_S$) given by the $M$-NN distance from $\widecheck{X}_j$ among the points in $\widecheck{\mathcal{X}}_n^{1-D_j}$. In other words, for a point collection $\check{\nu}_k = \{\check{x}_i\}_{i=1}^k$, we can take for $j \in [k]$,
The radius is also non-increasing in the point collection $\check{\nu}_k$, so the above holds true even if we add additional points $\mathcal{A}$ with $|\mathcal{A}| \le 7$ to $\check{\nu}_k$ (see the discussion in the beginning of the proof of lachieze2019normal).
We say that the score function $f_n$ satisfies the $(4+p)$-moment condition if there exists $p\in (0,1]$ such that for all $n\ge 9$, $\check{x}\in\widecheck{{\mathbb X}}$, $\widecheck{\mathcal{A}}\subset \widecheck{{\mathbb X}}$ with $|\widecheck{\mathcal{A}}|\le 7$,
In the sequel, we will usually write $M_{n}(\check{x}, \widecheck{A})$ instead for notational convenience.
The above moment condition could be motivated by connections to the classical Berry-Esseen theorem, which provides Gaussian approximation bounds under the assumption of i.i.d.\ score functions, i.e, in (ref), $f_{n}(\widecheck{X}_i,\widecheck{\mathcal{X}}_n)\equiv f_{n}(\widecheck{X}_i)$ for $i\in [n]$ with bounded third moments $\mathbb{E}|f_n(\widecheck{X}_1)|^{3}<\infty$. Specifically, it states that
where $\mathcal{N}$ is the standard normal random variable. Note here that our score functions are far from independent and can have significant local dependencies. The motivation behind (ref) is to go beyond independence. Assumptions of the radius of stabilization in Section (ref) restricts the dependence between the scores within balls with radii that decay exponentially. However, in contrast to the classical Berry-Esseen bound that require finite third moments, due to the dependence in our model, we require a slightly stronger $(4+p)$-moment condition.
We now outline the high-level ideas behind the proofs of our Gaussian and bootstrap approximation results.
Recall from (ref) the bias corrected ATE estimator $\hat{\tau}_{M}^{bc}=E_n+(B_M-\hat{B}_{M})$ for the covariate based matching. The proof of Theorem (ref) follows in three steps:
The proof of Theorem (ref) follows similar steps suitably adapted to the rank-based matching setting.
In order to bound the distributional distance between $\sqrt{n}(\tau_{M}^{\text{boot}}-\hat{\tau}_{M}^{\text{bc}})$ and $\sqrt{n}(\hat{\tau}_{M}^{\text{bc}}-\tau)$ appearing in the first assertion in Theorem (ref), we use the Gaussian limit in Theorem (ref) as the bridge. Roughly, The proof of Theorem (ref) broadly follow in three steps:
A similar procedure also applies for the bootstrap approximation in the case of the rank based ATE estimator.
The first, third and the fourth listed authors were supported by the NSF Grant DMS-2053918. The second author was supported in part by the German Research Foundation (DFG) Project 531540467.