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.
79,160 characters · 18 sections · 0 citation commands
Isotonic propensity score matching
In both randomized experiments and observational studies, matching estimators are widely used to estimate treatment effects. This paper proposes a novel one-to-many propensity score matching method of the average treatment effect (ATE), where the propensity score is assumed to be monotone increasing in the exogenous covariate and is estimated by the isotonic regression. Our matching scheme is exact, i.e., for the outcome $Y$, the binary treatment $W$, the covariate $X$, and a sample of size $N$, the matched set for the $i$-th unit is defined as \[ \mathcal{J}(i)=\left\{ j=1,\dots,N:W_{j}=1-W_{i}\text{ and }\tilde{p}(X_{j})=\tilde{p}(X_{i})\right\} , \] where $\tilde{p}(\cdot)$ is a uniformly consistent isotonic estimator of the propensity score developed in Section (ref). For multi-dimensional covariates $X$, we employ a monotone index model and consider the matched set: \[ \mathcal{J}(i)=\left\{ j=1,\dots,N:W_{j}=1-W_{i}\text{ and }\tilde{p}_{\tilde{\alpha}}(X_{j}^{\prime}\tilde{\alpha})=\tilde{p}_{\tilde{\alpha}}(X_{i}^{\prime}\tilde{\alpha})\right\} , \] where $\tilde{p}_{\tilde{\alpha}}(\cdot)$ is a uniformly consistent monotone single-index estimator of the propensity score developed in Section (ref).
Remarkably, the isotonic estimator proves to be especially well-suited as the initial nonparametric estimator in a two-stage semiparametric approach to estimating the ATE. It incorporates features of both matching and weighting estimators into the second-stage ATE estimator, addressing at least five issues commonly encountered by existing matching methods in the causal inference literature.
First, it is well known that the existing matching estimators of the ATE with a fixed number of matches are inefficient (Abadie and Imbens, 2006) since they do not balance bias and variance in the second-stage estimation. In comparison, our isotonic matching estimator is more efficient. In the univariate case, our method attains the semiparametric efficiency bound; in the multivariate case, where the efficiency bound becomes more complicated, we show that our proposed estimator performs better than those based on a fixed number of matches with propensity scores derived from widely used parametric models such as probit and logit, which are prevalent in applied research.
Second, although the performance of fixed-number matching estimators can be improved by increasing the number of matches with the sample size, the efficiency gain is somewhat artificial (Imbens, 2004) since the optimal number of matches and data-dependent ways of choosing it have been open questions. However, these issues are addressed by recent papers by Armstrong and Kolesár (2021) and Lin et al. (2023). By specifying a large enough Lipschitz constant, Armstrong and Kolesár (2021) showed that the matching estimator with the number of matches set to one is minimax optimal if the conditional mean is restricted to be Lipschitz; by adding an estimated correction term. Lin et al. (2023) gave the optimal number of matches for a bias-corrected matching estimator. In this paper, we argue that the isotonic estimator can provide an alternative solution: It gives a piece-wise monotone increasing estimator, which partitions observations into different groups. Within these groups, the treated and untreated observations have the same estimated propensity scores, so they can be naturally matched to each other without the need of choosing the number of matches, weights, and relevant distance measures. (For our method, the distance is zero under any measure.) In contrast, these choice problems are unavoidable in traditional methods for both covariates matching and propensity score matching, no matter whether they are based on the inverse variance matrix (e.g., Abadie and Imbens, 2006) or the (empirical) density function (e.g., Imbens, 2004) of covariates. Surprisingly, the set of the matching counterparts adaptively selected by the isotonic estimator automatically becomes the optimal choice in the second stage, in that it achieves the semiparametric efficiency bound of ATE (Hahn, 1998) for the univariate case.
Third, compared to other semiparametric matching methods, where the first stage propensity score is estimated with kernel or series-based techniques, our method is more practical in a twofold sense. It is not only free from the choice of the optimal number of matches, as mentioned in the second point above, but also does not involve smoothing parameters of conventional nonparametric methods, such as series length or bandwidth. In general, choosing the tuning parameters of a first-stage nonparametric estimator remains a difficult open question in the semiparametric estimation literature. The MSE optimal tuning parameter is usually not a good choice since the optimal first-stage estimator of the nuisance function does not imply the optimality of the second-stage semiparametric estimation (Bickel and Ritov, 2003). To ensure the $\sqrt{N}-$consistency of a semiparametric estimator, “undersmoothed” tuning parameters should be applied (Newey, 1994). But it is difficult to find a clear standard for shrinking tuning parameters below their MSE optimal values. The non-smooth nature of the isotonic estimator, on the other hand, turns out to automatically render an adequate amount of undersmoothing. At the cost of a monotonicity assumption imposed on the nuisance function, our proposed estimator avoids this choice problem while still maintaining other desirable properties of a decent semiparametric estimator, such as $\sqrt{N}-$consistency or efficiency.
Fourth, compared to popular parametric models of propensity scores, such as probit and logit, our proposed method contains a nonparametric first stage, so it is more robust to model misspecification. We acknowledge that combined with a single index structure, the probit and logit models can also approximate many different data-generating processes. But our method will always be more robust than them since both probit and logistic functions are monotone increasing themselves. In other words, the isotonic regression can well estimate all the data generating processes that can be well approximated by probit or logit model, but not vice versa. In addition, this robustness is achieved without costing the efficiency (compared to parametric methods) of the second-stage matching estimator.
Fifth, it is well known that the nonparametric bootstrap of the fixed-number matching estimator is invalid in the presence of continuous covariates (Abadie and Imbens, 2008). In the past decade, much work has been done to solve this problem by proposing cleverly structured wild bootstrap procedures. Otsu and Rai (2017) proposed a consistent wild bootstrap for covariates matching, and their approach was extended by Bodory et al. (2016) and Adusumilli (2020) to propensity score matching estimators. In our paper, we show that all these intricate bootstraps are no longer necessary in the case of monotone increasing propensity scores since the nonparametric bootstrap inference is asymptotically valid for our isotonic matching estimator.
Our method relies on the monotonicity assumption on propensity scores. Monotonicity is a natural shape restriction that can be justified in many applications in social science, economic studies, and medical research. Well-known examples in economics include the demand function, which is usually monotone decreasing in prices, and the supply or the utility functions, which are often monotone increasing in quantities. Furthermore, many functions derived from cumulative distribution functions (CDF) inherit the monotonicity from the latter. For example, in a threshold-crossing binary choice model
the conditional expectation of $Y$ on $X$ can be written as $\mathbb{E}[Y|X]=\mathbb{P}(Y=1|X)=F_{\varepsilon}(X^{\prime}\beta_{0})$, where $F_{\varepsilon}(\cdot)$ is the CDF of an independent noise $\varepsilon$. If we assume $\varepsilon\sim N(0,1)$, (ref) becomes a probit model; if we assume $\varepsilon\sim\mathrm{Logistic}(0,\frac{\pi^{2}}{3})$, it becomes a logit model. Although both parametric models are widely applied in estimating the probability of treatments, we can relax the distributional assumptions on $\varepsilon$ and express (ref) with a semiparametric model $Y=F_{\varepsilon}(X^{\prime}\beta_{0})+\nu$, where $F_{\varepsilon}(\cdot)$ is a nonparametric link function. We emphasize that the link function is monotone increasing by construction. See Cosslett (1983, 1987, 2007), Matzkin (1992), and Klein and Spady (1993) for more discussions of the model (ref).\footnote{Although this paper explores monotonicity of the propensity score function, our isotonic regression approach may be extended to the regression-based estimators with monotonicity constraints on the expected outcome functions $\mathbb{E}[Y(1)|X]$ and $\mathbb{E}[Y(0)|X]$. However, it should be noted that if monotonicity is imposed on the link functions of index models, the regression-based approach is clearly more restrictive than the propensity-score-based approach (because monotonicity on $F_{\varepsilon}(\cdot)$ is not substantive).}
One of the main challenges of developing the asymptotic properties of the proposed estimator is the inconsistency of the isotonic estimator at its boundaries, sometimes called the “spiking” problem in the literature. If the dependent variable is binary, there is a non-trivial probability for a non-shrinking group of left-end estimates to be exactly zero even under the strict overlap condition, regardless of the sample size; the right-end estimates have a similar issue. As a result, the matched sets for observations at two ends are empty, and we cannot construct a valid sample analog of ATE. Furthermore, observations near two ends are matched according to inconsistently estimated propensity scores, which are biased towards zero or one, resulting in a detrimental effect on the ATE estimator similar to the one caused by limited overlaps (Khan and Tamer, 2010; Rothe, 2017; among others). Although truncating those observations, whose propensity scores (either estimated parametrically or nonparametrically) are closer to 0 and 1, is widely implemented in applied work, this strategy has two caveats if one works with the isotonic estimator. The first problem is the size of truncation: If too little was truncated, it might be insufficient to correct the boundary problem. A safe choice of truncation in the literature for different problems involving isotonic estimators is to truncate the first and last $\alpha_{N}$-th quantile, with $\alpha_{N}\sim N^{-1/3}$ (or up to a logarithmic factor, see Wright, 1981; Durot, Kulikov and Lopuhaä, 2013; and Babii and Kumar, 2021). However, this truncation scheme is too much for our purpose. In fact, for any $\alpha_{N}$ such that $\alpha_{N}N^{1/2}\to\infty$, the truncated ATE estimator might be no longer $\sqrt{N}$-consistent.\footnote{This problem is not universal for every semiparametric estimator. For example, for a partially linear model $Y=X\beta+\psi(Z)+\varepsilon$, we can truncate more than its $N^{-1/2}$-th quantile, and the estimator of $\beta$ maintains $\sqrt{N}$-consistency. In fact, one can get $\sqrt{N}$-rates even if $\beta$ is estimated from an arbitrary sub-sample with a size proportional to $N$ since different $X$'s are linked to the same $\beta$. However, for ATE, in general, the truncated parts directly constitute estimation bias.} Second, as discussed in Appendix (ref), one of the key conditions for $\sqrt{N}$-consistency and efficient estimation of ATE is (ref) below, but whether this condition still holds after truncation is unclear. To solve these two problems, we extend the everywhere-consistent isotonic estimator of Meyer (2006) to a uniformly consistent isotonic (hereafter, UC-isotonic) estimator, which is by design to suit our two-stage semiparametric matching estimator. The proposed estimation procedure does not involve any truncation, the above-mentioned favorable properties of the isotonic estimator remain intact, and the full set of data is utilized in both the first stage estimation of the propensity score and the second stage estimation of ATE.
Our proposed method builds on the large literature of causal inference for covariate and propensity score matching estimators, e.g., Rosenbaum and Rubin (1983, 1984), Rosenbaum (1989), Heckman, Ichimura and Todd (1997, 1998), Heckman, Ichimura, Smith and Todd (1998), Dehejia and Wahba (1999), Abadie and Imbens (2006, 2008, 2011, 2016), Imbens (2004), Frölich (2004), Frölich, Huber and Wiesenfarth (2017), Otsu and Rai (2017), Bodory, Camponovo, Huber and Lechner (2016), Adusumilli (2020), among others. The propensity score matching estimators studied in the literature mainly use parametrically estimated propensity scores, such as probit and logit. Our proposed method, in contrast, uses a special type of nonparametric estimator, the isotonic estimator, to estimate the propensity score.
The isotonic estimator has a long history. The earlier work includes Ayer et al. (1955), Grenander (1956), Rao (1969, 1970), and Barlow and Brunk (1972), among others. The isotonic estimator of a regression function can be formulated as a least square estimation with a monotonicity constraint. Suppose that the conditional expectation $\mathbb{E}[Y|X]=p_{0}(X)$ is monotone increasing. Then, for an iid random sample $\{Y_{i},X_{i}\}_{i=1}^{N}$, the isotonic estimator is the minimizer of the sum of squared errors, $\min_{p\in\mathcal{M}}\sum_{i=1}^{N}\{Y_{i}-p(X_{i})\}^{2},$ where $\mathcal{M}$ is the class of monotone increasing functions. The minimizer can be calculated with the pool adjacent violators algorithm (Barlow and Brunk, 1972), or equivalently by solving the greatest convex minorant of the cumulative sum diagram $\{(0,0),(i,\sum_{j=1}^{i}Y_{j}),i=1,\ldots,N\}$, where the corresponding $\{X_{i}\}_{i=1}^{N}$ are ordered sequence. See Groeneboom and Jongbloed (2014) for a comprehensive discussion of different aspects of isotonic regression.
Our work is linked to the vast literature on semiparametric estimation, e.g., Chamberlain (1987), Robinson (1988), Newey (1990, 1994), van der Vaart (1991), Andrews (1994), Hahn (1998), Ai and Chen (2003), Bickel and Ritov (2003), Chen, Linton and Van Keilegom (2003), Chen and Santos (2018), among others. In most of the works cited above, nonparametric methods involving smoothing parameters were applied at the initial stage, while our work uses the isotonic estimation that is non-smooth and does not involve smoothing parameters. On the other hand, the double machine learning estimators (hereafter, DML; see, e.g., Robins, Rotnitzky and Zhao, 1995; Chernozhukov et al., 2017, 2018; among others) provide efficient estimators of the ATE that do not rely on subjective choices of smoothing parameters, thereby, to some extent, sharing many advantages of our approach. We provide a detailed comparison between the isotonic propensity score matching estimator and the DML for the ATE in Section (ref).
There are some authors working on concrete semiparametric models with plug-in isotonic estimators. Huang (2002) studied the properties of the monotone partially linear model, and his work was extended by Cheng (2009) and Yu (2014) to the monotone additive model. Balabdaoui, Durot and Jankowski (2019) studied the monotone single index model with the monotone least square method, and Groeneboom and Hendrickx (2018), Balabdaoui, Groeneboom and Hendrickx (2019), and Balabdaoui and Groeneboom (2021) (the last two papers are called BGH hereafter) developed a score-type approach for the monotone single index model and show the single index parameter can be estimated at $\sqrt{N}$-rate. Building on previous works, Xu (2021) studied a general framework of semiparametric Z-estimation with plug-in isotonic estimators, monotone single-index estimators, or monotone additive estimators, and applied the generic result to inverse probability weighting (IPW) estimators of ATE. For the augmented IPW (AIPW) model, Qin et al. (2019) and Yuan, Yin and Tan (2021) applied the monotone single index model to estimate the propensity score, then plugged the estimated propensity scores with other estimates of potential outcomes into a doubly-robust moment function. Their asymptotic results rely on the consistent estimations of both propensity scores and potential outcomes, and thus differ from our approach.
In terms of applying isotonic regression to estimate the ATE, the primary difference between this paper and Chapter 3 of Xu (2021) is that we address the boundary issue inherent in the isotonic estimator, while Xu (2021) relies on a stronger assumption adapted from Assumption 5.1 in Newey (1994). In the process of writing this paper, we have gradually realized that this assumption does not automatically apply to the IPW estimator, although it straightforwardly holds for some other semiparametric models, such as the monotone partially linear model and the monotone single index model, wherein the plugged-in isotonic estimator is not in the denominator. Compared to Chapter 3 of Xu (2021), the main contributions of this paper are: (i) proposing a UC-isotonic estimator that is suitable as the first-stage estimator in a propensity score matching estimator of the ATE; (ii) revealing the equivalence between the matching estimator and the IPW estimator when the first-stage propensity score is estimated via UC-isotonic regression; and (iii) based on this equivalence, enriching the literature on propensity score matching by introducing a new approach that addresses several problems of the existing matching methods, as detailed at the beginning of this introduction.\footnote{At almost the same time, an independent work by Liu and Qin (2022) derived a similar equivalence result for the average treatment effect on treated (ATT). Recently, a revised version of Liu and Qin (2022) is published as Liu and Qin (2024). There are two main differences between our paper and their papers. First, we formally address the boundary problem of the isotonic estimator and achieve the $\sqrt{N}$-normality of the ATE estimator by proposing a uniformly consistent isotonic estimator. Second, our asymptotic analysis of the model with multivariate covariates in Section (ref) focuses on a more general case, where the influence of the estimation errors from the parametric component of the first-stage monotone single index model is maintained.}
The rest of the paper is organized as follows. After introducing the setting and notations, Section (ref) shows the implementation and asymptotic properties of the proposed isotonic matching estimator with a univariate covariate. Section (ref) compares our approach with existing matching estimators as well as the double machine learning estimator for the ATE. The univariate results are extended to the case of multivariate covariates in Section (ref), where the propensity score is modeled by a semiparametric single-index model with an unknown monotone increasing link function. In Section (ref), we establish the validity of the nonparametric bootstrap. Monte-Carlo simulation studies are presented in Section (ref). All proofs are presented in Appendix, while additional theoretical details and simulation comparisons are provided in Supplementary Material.
Suppose we observe the triple $(Y,W,X)$ drawn randomly from the product space $\mathcal{Z}=\mathbb{R}\times\{0,1\}\times\mathcal{\mathcal{X}}$. Within the triple, $W\in\{0,1\}$ is a binary treatment variable, $Y=W\cdot Y(1)+(1-W)\cdot Y(0)$ is an outcome variable with potential outcomes $Y(1)$ and $Y(0)$ for $W=1$ and $0$, respectively, and $X$ is a scalar covariate with continuous domain $\mathcal{X}=[x_{L},x_{U}]\subset\mathbb{R}$. In this section, we tentatively assume $X$ is scalar, and discuss extensions for multivariate $X$ in Section (ref). Without loss of generality, $\{Y_{i},W_{i},X_{i}\}_{i=1}^{N}$ is an iid sample of $(Y,W,X)$ and is ordered by $X$. If $X$ is continuously distributed, we should have $X_{1}<X_{2}<\cdots<X_{N}$ with probability one (i.e., no ties).
In this section, we consider the estimation of the ATE, $\tau=\mathbb{E}[Y(1)-Y(0)]$, by matching the propensity score $p(x)=\mathbb{P}(W=1|X=x)=\mathbb{E}[W|X=x]$, where $p(\cdot)$ is an unknown monotone increasing function. In particular, we estimate $p(\cdot)$ by the isotonic estimator
where $\mathcal{M}$ is the class of all monotone increasing functions defined on $\mathcal{X}$. Since Brunk (1958), this isotonic regression estimator has been extensively studied in the statistics literature (see, e.g., Barlow et al., 1972, and Groeneboom and Jongbloed, 2014, for an overview). One of the well-known features of isotonic regression is that the estimator $\hat{p}(\cdot)$ is a monotone increasing piecewise constant function with jump points at $\{X_{n_{k}}\}_{k=1}^{K}$ for some integer $K$ with $1\leq K\le N$. By these jump points, the sample is divided into $K$ disjoint groups, with $\{n_{k}\}_{k=1}^{K}$ denoting the first indices of these $K$ groups. Further, we let $N_{k}$ denote the number of observations belonging to the $k$-th group. Based on these definitions, it holds that $n_{k}+N_{k}=n_{k+1}$ for each $k=1,\ldots,K-1$, and $\sum_{k=1}^{K}N_{k}=N$. Note that the integer $K$ and the corresponding disjoint groups are automatically determined by the isotonic estimation algorithm (see the formula ((ref)) below), rather than being chosen by the user.
To avoid ambiguity caused by splitting a flat piece into several sub-pieces with the same estimated value, we impose
to ensure uniqueness of this partition (i.e., if $\hat{p}(X_{n_{k}})=\hat{p}(X_{n_{k+1}})$, we simply combine the groups $k$ and $k+1$). Then, the isotonic estimator $\hat{p}(\cdot)$ is characterized as follows.
To define our propensity score matching estimator based on $\hat{p}(\cdot)$, let $N_{k,1}$ and $N_{k,0}$ denote the numbers of treated and controlled observations within group $k$, i.e., $N_{k,1}=\sum_{i=n_{k}}^{n_{k}+N_{k}-1}W_{i}$ and $N_{k,0}=N_{k}-N_{k,1}$. Our one-to-many matching method is implemented within each of these $K$ groups, and each treated (controlled) observation in group $k$ will be matched with its $N_{k,0}$ ($N_{k,1}$) counterparts, which belong to the same group and have the same value of the estimated propensity score. The following results directly follow from Proposition (ref).
Before we proceed, we need to solve the problem of the potential lack of matching counterparts for those $i$'s with $\hat{p}(X_{i})=0$ or $1$. Under the strict overlaps (in Assumption (ref) below), the problem is essentially associated with the inconsistency of the isotonic estimator at the boundary. In the next subsection, we propose a modified isotonic estimator that is uniformly consistent on $\mathcal{X}$.
Like other nonparametric estimators, the isotonic estimator is imprecise at the boundary. If we apply the isotonic estimator to the binary dependent variable $W$, there is a non-trivial probability of $\hat{p}(X_{i})=0$ or $1$ even if the true propensity score $p(x)$ is bounded away from zero and one for all $x\in\mathcal{X}$. For example, if $\hat{p}(X_{1})=0$, (ref) implies $N_{1,1}=0$, i.e., there are no matching counterparts for the treated units.
To fix this problem, we propose a modified isotonic estimator that is uniformly consistent in its domain at a $\left(\log N\right)^{1/3}N^{-1/3}$ rate and is easy to implement. For the sample $\{W_{i},X_{i}\}_{i=1}^{N}$ with $X_{1}<\cdots<X_{N}$, we transform $\{W_{i}\}_{i=1}^{N}$ into $\{\tilde{W}_{i}\}_{i=1}^{N}$ by averaging its first and last $\lfloor N^{2/3}\rfloor$ observations:
Our proposed UC-isotonic estimator is obtained by implementing the standard isotonic regression of $\tilde{W}$ on $X$:
A similarly modified estimator was proposed by Meyer (2006), where she averaged the first and last $\lceil\log(N)\rceil$ dependent variables instead of the first and last $\lfloor N^{2/3}\rfloor$ ones. The choices are different because she focuses on the consistency of the isotonic estimator itself, while we are interested in the performance of the second-stage matching estimator. To achieve an $N^{-1/2}$ rate at the second stage, we need the isotonic estimator to be uniformly consistent at a rate faster than $N^{-1/4}$, which won't be achieved under Meyer's choice. Meyer (2006) presented a theorem regarding the consistency of the modified estimator at the boundary; however, a proof of consistency was not provided, nor was the rate of convergence discussed.
In this paper, we formally establish the uniform convergence rate of the modified isotonic estimator $\tilde{p}(\cdot)$. To this end, we impose the following assumption.
Assumption (ref) (i) is our main assumption, the monotonicity of $p(\cdot)$. Assumption (ref) (ii) is required for the $\sqrt{N}-$consistency of the second-stage matching estimator. The same assumption has been adopted by Groeneboom and Hendrickx (2018; Assumption A2) and by BGH (Assumption A3 and its accompanying remark; see also Lemma 22 in the supplementary material of BGH) in the context of the monotone single index model. If we believe that the underlying propensity score function has some flat parts where $p^{(1)}(x)=0$, we could first run an isotonic estimation of $\tilde{W}_{i}+c\cdot X_{i}$ on $X_{i}$, where $c$ is a positive constant, to obtain $\tilde{p}_{c}(x)$. Then, by subtracting the linear trend $c\cdot x$ from $\tilde{p}_{c}(x)$, we obtain a consistent estimator of $p(x)$.\footnote{Technically, if $p(\cdot)$ has some flat parts where $p^{(1)}(x)=0$, then the original estimator $\tilde{p}(\cdot)$ may not satisfy the requirement in (ref) in Appendix (ref), $|\delta(x)-\bar{\delta}_{N}(x)|\leq C_{0}|p(x)-\tilde{p}(x)|$. Flat parts in $p(\cdot)$ imply that $p(x)-\tilde{p}(x)=0$ might hold within an entire interval, potentially leading to a violation of (ref). In contrast, if $p(\cdot)$ is strictly monotone increasing, then $p(\cdot)$ and $\tilde{p}(\cdot)$ will cross at most once within each partition given by the isotonic estimator, since $\tilde{p}(\cdot)$ is a piecewise flat function. We refer to Sections 10.2-10.3 and Figure 10.1 of Groeneboom and Jongbloed (2014) for more details.} Assumption (ref) (iii) imposes an upper and lower bound for the density of $X$.
To avoid unnecessarily repeatedly defined notations, we let the same set of notations, $K,$ $N_{k,1}$, $N_{k}$, and $n_{k}$, denote the number of groups, the number of treated observations in group $k$, the number of members in group $k$, and the index of the first element of group $k$, under the grouping scheme given by the UC-isotonic estimator $\tilde{p}(\cdot)$ ($N_{k,1}$ is calculated with the original treatment variable $\{W_{i}\}_{i=1}^{N}$). We obtain an analogous result to Proposition (ref) for the UC-isotonic estimator.
Part (i) of this proposition says that all the averaged $W_{i}$'s at the beginning and end of the data are absorbed in the first and the last group. Part (ii) provides an analogous representation of the UC-isotonic estimator $\tilde{p}(\cdot)$ as $\hat{p}(\cdot)$. While Part (i) gives a lower bound of the sizes of the first and the last group, Part (iii) gives (stochastic) upper bounds of them. Based on this proposition, the uniform convergence rate of the UC-isotonic estimator is obtained as follows.
Finally, to guarantee the existence of matching counterparts by $\tilde{p}(\cdot)$, we impose the strict overlap condition.
Assumption (ref) is standard in the treatment effect literature. It is necessary for the identification and $\sqrt{N}$-consistent estimation of the ATE. Combining Proposition (ref) and Theorem (ref) with Assumption (ref), the existence of the matching counterparts by $\tilde{p}(\cdot)$ is obtained as follows.
Based on the UC-isotonic estimator $\tilde{p}(\cdot)$, the isotonic propensity score matching estimator for the ATE $\tau$ can be implemented as follows.
We proceed with the following assumptions.
Assumption (ref) (i)-(iii) regulates the tail behaviors of the (conditional functions of) potential outcomes, which are necessary for $\sqrt{N}$-consistent estimation. Assumption (ref) (iv) is the standard unconfoundedness assumption. Under these assumptions, we have the following key equivalence result.
Imbens (2004) pointed out that with $M\to\infty$ and $M/N\to0$, the matching estimator is essentially like a regression estimator. In comparison, we find out that with propensity scores estimated by the UC-isotonic estimator, the (propensity score) matching estimator is numerically equal to the weighting estimator in each finite sample. This equivalence is tightly associated with the fact that the isotonic estimator can be regarded as a type of partitioning estimator (e.g., Györfi et al., 2002; Cattaneo and Farrell, 2013). See Section (ref) below for a further comparison of isotonic and partitioning estimators within the context of a two-stage matching estimator of the ATE. Additionally, our method is related to the propensity score methods of blocking, stratification, and radius matching (Rosenbaum and Rubin, 1983, 1985; Dehejia and Wahba, 1999, 2002; among others). See Section (ref) for a comparison with these methods.
Moreover, as mentioned in the introduction, the equivalence result in Theorem (ref) relies crucially on the implementation of the UC-isotonic estimator (ref), which guarantees that both the matching and IPW estimators at the second stage are well-defined.
We notice that the threshold $\lfloor N^{2/3}\rfloor$ in the algorithm (ref) can be interpreted as an implicit tuning parameter. We would like to point out that, first, it is convenient to choose since it depends only on the sample size $N$; second, it is aimed at correcting the boundary problem, which is also faced by other semiparametric and even parametric matching methods. In practice, trimming estimated propensity scores is widely adopted, and the amount of trimming is chosen subjectively in most cases. Our proposed method provides transparent guidance for correcting this common boundary issue. Furthermore, to investigate the impact of different threshold choices, we have included both theoretical analysis and simulation evidence in Sections (ref) and (ref) of the supplementary material, respectively.
Our main result, consistency and asymptotic normality of the isotonic propensity score matching estimator, is obtained as follows.
We note that the asymptotic variance $\Omega$ is the semiparametric efficiency bound for $\tau$ (see e.g., Hahn, 1998, and Hirano, Imbens and Ridder, 2003). Although we may conduct inference based on an estimator of $\Omega$, we suggest a bootstrap inference method, which will be discussed in Section (ref).
In this section, we draw comparisons of our approach with a range of related estimators for the ATE. The comparison with matching methods based on propensity score estimated by partitioning estimator is presented in Section (ref), the comparison with propensity score methods of blocking, stratification, and radius matching is presented in Section (ref), the comparison with matching methods based on propensity score estimated by regression trees is presented in Section (ref), and the comparison with the double machine learning (DML) estimator for the ATE can be found in Section (ref).
One notable feature of the proposed isotonic propensity score matching method is that it is a one-to-many matching method that provides exact matches, as illustrated by the formula (ref). This is attributed to the isotonic estimator being considered a special type of partitioning estimator, in which the volume sizes of partitions are automatically chosen by the monotonicity constraint, and a simple average is implemented within each partition.
The partitioning estimator is a nonparametric method for estimating regression functions.\footnote{We refer to Györfi et al. (2002) and Cattaneo and Farrell (2013) for comprehensive discussions of the partitioning estimator.} It divides the domain of the running variables into disjoint partitions. Within each partition, a local estimator is implemented by the user, such as the sample mean, a linear estimator, or a series estimator. Each sample point is exclusively used in the estimation within the partition to which it belongs. This feature simplifies the complex correlation structure of a matching estimator such that it achieves equivalence with an IPW estimator. For the UC-isotonic estimator, this equivalence is presented by equation (ref) in Appendix (ref). In the resulting matching estimator of ATE, the same set of partitions serves both the first- and the second-stage nonparametric estimation. Usually, these two stages are not associated with each other since they have distinct objects, the propensity score and the potential outcomes. Certainly, a matching estimator of the ATE that utilizes propensity scores estimated with a partitioning estimator should exhibit a similar equivalence to the weighting estimator. However, the selection of the number of partitions and their sizes necessitates careful consideration, as they must meet specific undersmoothing conditions to secure the desired asymptotic properties of the second-stage ATE estimator. The challenge of selecting an appropriate undersmoothed bandwidth or volume size, as mentioned in the introduction, remains a difficult open question in the semiparametric estimation. In contrast, our proposed isotonic matching estimator automatically chooses these tuning parameters, leading to the efficient estimation of the ATE, as demonstrated in Theorem (ref).
Our proposed isotonic matching estimator is also related to some of the seminal ideas introduced at the outset of the propensity score methods: blocking, stratification (Rosenbaum and Rubin, 1983; Dehejia and Wahba, 1999, 2002), and radius matching (Rosenbaum and Rubin, 1985).
The isotonic propensity score matching shares similarities with blocking and stratification matching on propensity scores, notably: (i) they initially categorize data points into distinct groups (or strata, blocks, partitions) according to estimated propensity scores, and (ii) within each group, they calculate the conditional average treatment effect as the simple difference in means of outcomes between the treatment and comparison groups. The primary distinction lies in the grouping mechanism: for isotonic propensity score matching, the groups are determined adaptively in a data-driven manner through isotonic regression, whereas for the stratification estimator of the ATE, the strata must be explicitly specified by the user. Another distinction is that for the isotonic propensity score matching method, the same set of partitions is utilized for both the first and second stages of nonparametric estimation. As presented by Theorem (ref), this characteristic leads to the equivalence between the matching and IPW estimator, resulting in the efficient estimation of the ATE. In contrast, in the case of blocking or stratification matching methods, particularly when the propensity score is estimated using parametric models, this equivalence cannot generally be established, and efficiency cannot be assured without implementing some bias correction method.
The case for the radius matching estimator is similar to the stratified matching estimator. The difference is that for stratified matching, each unit is matched solely with units from the opposite treatment group within the same stratum, while radius matching allows each unit to be matched to several local balls, the centers of which belong to the opposite treatment group. For both radius and stratified matching estimators, the sizes of strata or the radii act as tuning parameters, which must be chosen by the users when the propensity score is estimated via parametric or nonparametric methods dependent on smoothing parameters (such as kernel or series estimation). These smoothing parameters play a key role in balancing the bias and variance, thereby significantly affecting the second-stage estimator of ATE. In contrast, isotonic regression distinguishes itself by automatically generating these partitions through the application of the monotonicity constraint.
As methods of estimating the propensity scores, the isotonic estimator and regression trees share several similarities. First, both are nonparametric estimators that do not impose restrictive parametric structures on the underlying response function. Second, both approaches partition the domain of running variables (the feature space in regression tree terminology) into several regions and use the sample average within each region as estimators. As a result, both estimators take the form of piecewise-constant functions. Third, both methods form their piecewise-constant functions in data-adaptive manners. In particular, the partitions created by both methods depend on the dependent variable (the response), which differentiates them from regular nonparametric methods, such as the kernel estimator.
On the other hand, there are notable distinctions between the two methods. First, both approaches construct their piecewise-constant functions differently: the partitions in a regression tree are obtained in a stepwise manner. In each step, a partition is chosen to achieve the maximum marginal reduction of the mean square error (MSE), without imposing any shape constraints during this process. In contrast, the isotonic estimator employs a one-step approach that determines partitions to minimize the MSE over the class of monotone functions. Second, although both approaches are data-driven, the isotonic estimator is free of smoothing parameters, whereas the regression tree depends on the user to specify the tree's length. (When the tree length is determined by cross-validation, the user must select the penalty parameter.) Third, the regression tree is inherently designed for multi-dimensional problems, whereas the canonical form of isotonic regression addresses one-dimensional issues, given that the traditional definition of monotonicity characterizes the relationship between two variables. Nevertheless, the isotonic estimation can be extended to multivariate cases by being incorporated into a partially linear model or a monotone single index model. The latter is illustrated in Section (ref) below.
To summarize, the isotonic estimator necessitates the monotonicity assumption in the underlying response function, a requirement not shared by regression trees. This assumption, however, enables the isotonic estimation algorithm to automatically regulate the trade-off between bias and variance. Conversely, when using regression trees, practitioners are tasked with the challenge of selecting an appropriate tree length to effectively manage the balance between bias and the risk of overfitting. The strength of regression trees is their natural aptitude for tackling multivariate problems. When employing regression trees in the preliminary stage of propensity score estimation as part of a two-stage approach to estimating the ATE, it is commonly combined with methods for bias correction and sample splitting, as discussed by Chernozhukov et al. (2018). See Section (ref) below for more details about the comparison of our approach with the double machine learning estimator.
The isotonic propensity score matching estimator and the DML estimator for the ATE both share the benefit of not requiring subjective choices of tuning parameters. For estimating the ATE, a typical example of a DML estimator is given by applying the sample splitting to the augmented inverse probability weighting (AIPW) estimator. In the following, we abstract from sample splitting to simplify notation:
where $\hat{\psi}_{1}$$(\cdot)$ and $\hat{\psi}_{0}(\cdot)$ are estimators of $\mathbb{E}[Y(1)|X=\cdot]$ and $\mathbb{E}[Y(0)|X=\cdot]$, respectively. The first and second lines of (ref) present two formulations of the DML estimator for the ATE. The first terms in both lines correspond to the standard regression and IPW estimators, respectively, while the subsequent terms represent their bias-correction components.
The AIPW has been extensively studied since the seminal work of Robins, Rotnitzky and Zhao (1995), Robins and Rotnitzky (1995); see also Newey, Hsieh, and Robins (1998, 2004), Scharfstein, Rotnitzky and Robins (1999), Rothe and Firpo (2019), among others. In an influential work, Chernozhukov et al. (2018) combined orthogonal moment functions – of which the formula (ref) is a specific case for the ATE – with sample splitting, accommodating a broad array of the first-stage machine learners that are prone to bias due to regularization or model selection. Recent developments by Chernozhukov et al. (2022) and Chernozhukov, Newey and Singh (2022) have proposed methods for constructing the correction term without requiring an explicit function form for the bias correction.
Both estimators have their own advantages and comparative strengths. From a practical standpoint, the isotonic propensity score matching method stands out for its simplicity and ease of implementation: it does not require the correction terms, thereby sparing the effort of estimating the conditional means of potential outcomes and sidesteps the challenges associated with their correct specification. In contrast, the DML estimator’s efficiency relies on correctly specifying and effectively estimating both the propensity score and the conditional means of potential outcomes. A misstep in either leads to a consistent yet inefficient estimator. On the other hand, the DML estimator exhibits great flexibility: through the use of sample splitting, it supports a variety of first-stage estimators, accommodating high dimensional data or highly complex function classes, such as random forest, neural networks, and other advanced machine learning technologies.
From a technical standpoint, the isotonic propensity score matching and the DML for the ATE represent two distinct pathways of semiparametric estimation: undersmoothing and bias correction. Both strategies aim for $\sqrt{N}$-consistent (or efficient in certain cases) estimators (see Newey, 1994, for a relevant discussion). The undersmoothing strategy depends on a first-stage estimator with reduced bias, achievable in nonparametric estimators by selecting smoothing parameters smaller than the MSE-optimal levels. Conversely, the bias correction method addresses bias by incorporating an estimated correction term into the second-stage sample moment function, rather than concentrating on the first stage.
The proposed isotonic propensity score matching estimator utilizes the isotonic estimator, which achieves a similar effect of “undersmoothing”, and this effect is automatically rendered by enforcing monotonicity. The isotonic estimator does not really shrink its bias to a level lower than $N^{-1/2}$. However, when combined with the monotonicity, it eventually achieves a deviation from the efficient influence function that decays at a rate faster than $N^{-1/2}$ (see (ref) in Appendix). In contrast, the DML for the ATE represents a typical bias correction approach. The second terms in both lines of (ref), while achieving the “doubly robust” effect, also serve as bias-correction components. At the cost of computing additional correction terms and some efficiency loss due to sample splitting, the DML approach manages to mitigate potential bias and prevent overfitting risks, while being less restrictive on the first-stage estimation. It is not only less sensitive to the choice of the smoothing parameter for the traditional first-stage nonparametric estimator but can also accommodate many black-box machine learning methods, whose asymptotic properties remain to be fully understood. Consequently, the theoretical development of the isotonic propensity score matching and the DML for the ATE differs substantially. The DML approach significantly reduces the effort needed to address issues arising from the complexity of function classes, which is associated either with the correlation brought by plug-in estimators or with the choice of smoothing parameter. In contrast, this paper needs to address the impact of the plug-in estimator in the theoretical development of the isotonic propensity score matching estimator.
Finally, we would like to emphasize that our proposed method represents a targeted advancement within the matching estimation literature, specifically addressing several limitations present in existing matching techniques for estimating the ATE. In contrast, the DML is a versatile tool designed for broader semiparametric estimation tasks, which include a wide array of econometric problems such as average derivatives, partially linear models, and parameters of economic structural models. Our approach, therefore, complements rather than competes with the expansive toolkit that DML offers, by providing subtle yet significant improvements in the specialized area of matching estimation.
Certainly, researchers are more interested in models with multivariate covariates $X$. One way to balance the robustness and the curse of dimensionality is to estimate the propensity score with the monotone single-index model:
where $p_{0}(\cdot)$ is a monotone increasing link function of its index $X^{\prime}\alpha_{0}$ and $X\in\mathbb{R}^{k}$. For identification, $\text{\ensuremath{\alpha}}_{0}$ is a $k$-dimensional vector normalized with $||\text{\ensuremath{\alpha}}_{0}||=1$.\footnote{In the estimation, the constraint $||\text{\ensuremath{\alpha}}_{0}||=1$ can be dealt with reparametrization or the augmented Lagrange method by Balabdaoui and Groeneboom (2021). In this section, we study our model without discussing those technical details. See BGH for more details.}
For a binary dependent variable, this model can be derived from (ref), and $p_{0}(\cdot)$ is by nature monotone increasing. It was studied by Cosslett (1983, 1987, 2007), Han (1987), Matzkin (1992), Sherman (1993), Klein and Spady (1993), among others. In the case where $p_{0}(\cdot)$ is estimated with isotonic regression, Balabdaoui, Durot and Jankowski (2019) studied (ref) with the monotone least square method, and Groeneboom and Hendrickx (2018), Balabdaoui, Groeneboom and Hendrickx (2019), and Balabdaoui and Groeneboom (2021) (BGH) estimated $\alpha_{0}$ and $p_{0}(\cdot)$ by solving a score-type sample moment condition of
To estimate $p_{0}$ and $\text{\ensuremath{\alpha}}_{0}$, we can apply the method of BGH. For a fixed $\alpha$, define
where $\mathcal{M}$ is the set of monotone increasing functions defined on $\mathbb{R}$. Note that $\hat{p}_{\alpha}(u)$ can be solved with isotonic regression of $W_{i}$ on the data points $\{X_{i}^{\prime}\alpha\}_{i=1}^{N}$. Then, $\alpha_{0}$ can be estimated by minimizing the squared sum of a score function. For example, the simple score estimator in Balabdaoui and Groeneboom (2021) is given by solving
BGH showed that under certain assumptions, $\hat{\alpha}$ is a $\sqrt{N}$-consistent estimator for $\alpha_{0}$,\footnote{BGH proposed solving a “zero-crossing” root of $\frac{1}{N}\sum_{i=1}^{N}X\{W_{i}-\hat{p}_{\alpha}(X_{i}^{\prime}\alpha)\}=0$. Then they realized that there is an issue with the existence of the zero-crossing root for a finite sample (due to the discreteness of $\hat{p}_{\alpha}$). To fix this problem, Balabdaoui and Groeneboom (2021) replaced this objective function with (ref), where a minimizer always exists. If there are multiple minimizers, any of them is a $\sqrt{N}$-consistent estimator for $\alpha_{0}$. (See a discussion on p.1426 of Groeneboom and Hendrickx, 2018). BGH also proposed an efficient estimator of $\alpha_{0}$ by solving a kernel-adjusted score function. Since our aim is the second-stage ATE $\tau$ instead of the first-stage propensity score $p$, we do not apply BGH's efficient estimator. It will introduce additional tuning parameters without improving the second-stage ATE.} and $\mathbb{E}[\hat{p}_{\hat{\alpha}}(X^{\prime}\hat{\alpha})-p_{0}(X^{\prime}\text{\ensuremath{\alpha}}_{0})]=O_{P}((\log N)N^{-2/3})$. We apply their method to estimate the propensity score with multi-dimensional control variables $X$.
In this section, $\tilde{\tau}$ denotes the ATE estimator based on the multi-dimensional covariates $X$. Similarly to Section (ref), to solve the boundary problem of the isotonic estimator to ensure that each observation has a non-empty matched set, we develop a uniformly consistent monotone single-index (hereafter, UC-iso-index) estimator, which is denoted by $\tilde{p}_{\tilde{\alpha}}$. The matching procedure can be implemented as follows.
We modify Assumptions (ref)-(ref) in Section (ref) as follows.
\begin{asm1'} {[}Sampling{]} $\{Y_{i},W_{i},X_{i}\}_{i=1}^{N}$ is an iid sample of $(Y,W,X)\in\mathbb{R}\times\{0,1\}\times\mathcal{X}$, where the space $\mathcal{X}$ is a convex subset of $\mathbb{R}^{k}$ with a nonempty interior. There exists $R>0$ such that $\mathcal{X}\subset\mathcal{B}(0,R)=\{x:||x||\leq R\}$. \end{asm1'}
Given $\alpha$, we define the true link function of (ref): \[ p_{\alpha}(u)=\mathbb{E}[W|X^{\prime}\alpha=u]. \] Obviously, $p_{\alpha_{0}}=p_{0}$. Let $a_{0}$ and $b_{0}$ be the minimum and the maximum of the interval $I_{\alpha_{0}}=\{x^{\prime}\alpha_{0}:x\in\mathcal{X}\}$, respectively.
\begin{asm2'}{[}Monotonicity and continuity{]} (i) There exists $\delta_{0}>0$ such that for each $\alpha\in\mathcal{B}(\alpha_{0},\delta_{0})$, the function $u\mapsto\mathbb{E}[W|X^{\prime}\alpha=u]$ is monotone increasing in $u$ and differentiable in $\alpha$; (ii) $p_{0}(\cdot)$ is continuously differentiable with its first derivative $p^{(1)}(u)>0$ on $u\in(a_{0}-\delta_{0}R,b_{0}+\delta_{0}R)$, and (iii) $X$ has a continuous density $f(x)$ satisfying that for some positive constants $\underline{f}$ and $\overline{f}$ , it holds $\underline{f}<f(x)<\overline{f}$ all $x\in\mathcal{X}$. \end{asm2'}
\begin{asm3'} {[}Strict overlaps{]} There exist positive constants $\underline{p}$ and $\bar{p}$ such that $0<\underline{p}\leq p_{0}(x^{\prime}\alpha_{0})\leq\bar{p}<1$ for all $x\in\mathcal{X}$. \end{asm3'}
\begin{asm4'} {[}Data generating process{]} (i) $\mathbb{E}[Y(0)^{2}]<\infty$ and $\mathbb{E}[Y(1)^{2}]<\infty$, (ii) $u\mapsto\mathbb{E}[Y(1)|X=x]$ are continuously differentiable for all $x\in\mathcal{X}$ and $\alpha\in\mathcal{B}(\alpha_{0},\delta_{0})$, (iii) for $D(Z)=\frac{WY}{p_{0}(X^{\prime}\text{\ensuremath{\alpha}}_{0})^{2}}+\frac{Y(1-W)}{\{1-p_{0}(X^{\prime}\text{\ensuremath{\alpha}}_{0})\}^{2}}$, there exist positive constants $c_{0}$ and $M_{0}$ such that $\mathbb{E}[|D(Z)|^{m}|X=x]\leq m!M_{0}^{m-2}c_{0}$ holds for all integers $m\geq2$ and every $x$, and (iv) $Y(1),Y(0)\perp W|X$ almost surely. \end{asm4'}
Let $Z$ denote the triple $(Y,W,X)$, and $\mathcal{Z}$ denote the space of the random vector $Z$. For each $\alpha\in\mathcal{B}(\alpha_{0},\delta_{0})$, $u\in I_{\alpha}=\{x^{\prime}\alpha:x\in\mathcal{X}\}$, and a function $f(\cdot)$ defined on $\mathcal{Z}$, we define $\mathbb{E}_{\alpha}[f(Z)|u]=\mathbb{E}[f(Z)|X^{\prime}\alpha=u]$. Similarly, we define the conditional covariance $\text{Cov}_{\alpha_{0}}(f(Z),X|u)$. The following two assumptions are adapted from BGH, which ensure that the score estimators (ref) and (ref) have desirable properties.
Based on Assumptions 1', 2', (ref), and (ref), we have a result similar to Proposition (ref), but the numbering is according to $X_{1}^{\prime}\tilde{\alpha}<\cdots<X_{N}^{\prime}\tilde{\alpha}$. The uniform convergence rate of the UC-iso-index estimator is obtained as follows.
The existence of matching counterparts is guaranteed by an argument similar to Corollary (ref). Finally, let $\mathbf{B}^{-}$ denote the Moore-Penrose inverse of a square matrix $\mathbf{B}$. The asymptotic properties of the isotonic propensity score matching estimator are obtained as follows.
Note that the semiparametric efficiency bound for estimating $\tau$ with known $\alpha_{0}$ is given by $\mathbb{E}[\{m(Z)+M(Z)\}\{m(Z)+M(Z)\}^{\prime}]$ (see, e.g., Newey, 1994). The additional term $A(Z)$ can be interpreted as the influence of estimating the index coefficients $\alpha_{0}$. This influence is also faced by parametric matching estimators. In general, our proposed method uses the matched sets, in which the number of matches increases to infinite, so it better balances the variance and bias in the second stage and should asymptotically outperform any matching method with fixed numbers of matches. In Section (ref) below, we present simulation results to illustrate that the proposed ATE estimator $\tilde{\tau}$ outperforms the probit matching estimator in every sample size, even in the case that the true propensity score is a probit (the correct specification).
Theoretically, the additional term $A(Z)$ can be avoided by using a semiparametric weighting estimator. However, the costs are strong assumptions on the smoothness of the propensity scores (typically, $7\cdot\dim(X)$-th continuous differentiability; see Hirano, Imbens and Ridder, 2003) and a proper choice of smoothing parameters. Our proposed method only requires the propensity score to be once continuously differentiable, and it does not involve smoothing parameters, such as bandwidths or series lengths.
The asymptotic variances in Theorems (ref) and (ref) contain conditional mean and variance functions, such as $\mathbb{V}(Y(1)|X)$ and $\mathbb{E}[X|X^{\prime}\alpha_{0}]$, which need to be estimated. If we use nonparametric methods to estimate them, we still have to choose some smoothing parameters even though the point estimators are free from smoothing. To avoid the estimation of such nonparametric components, we employ a bootstrap method to approximate the asymptotic distribution of the proposed isotonic propensity score matching estimator.
After Abadie and Imbens (2008) showed that the nonparametric bootstrap of the fixed-number matching estimator is invalid in the presence of continuous covariates, much work tried to solve this problem by proposing modified wild bootstraps, including Otsu and Rai (2017) for covariates matching estimators, and Bodory et al. (2016) and Adusumilli (2020) for propensity score matching estimators. In contrast, the nonparametric bootstrap of our one-to-many matching method is valid, which is an interesting implication of Theorem (ref). In this section, we discuss an asymptotically valid bootstrap procedure for the estimator $\hat{\tau}$ in Theorem (ref). This result can be similarly adapted to $\tilde{\tau}$ in Theorem (ref).
The nonparametric bootstrap is implemented as follows.
The asymptotic validity of this bootstrap approximation is obtained as follows.
In this section, we use three simulation studies to assess the finite sample properties of our isotonic propensity score matching estimator.
Let $X=0.15+0.7Z$, where $Z$ and $\nu$ are independently uniformly distributed on $[0,1]$, and
The true ATE is the coefficient of $W$, which is 0.5. The simulation results are presented in Table (ref), where $\hat{\mu}_{\tau}$ is the Monte-Carlo mean, and the mean square errors (MSE) are rescaled by $N$. The number of Monte-Carlo simulations is 5000 for each sample size.
The left panel shows the simulation results of the proposed matching method based on propensity scores estimated by the UC-isotonic estimator, and the right panel shows those of the one-to-one matching estimator based on propensity scores estimated with the logit model $\mathbb{P}(W=1|X=x)=\frac{\text{exp}(a+bx)}{\text{exp}(a+bx)+1}.$ The last row shows the true value of ATE and the semiparametric efficiency bound of this problem calculated according to Hahn (1998):
In comparison, the logit matching estimator has a slightly smaller bias, and it seems that both estimators are asymptotically unbiased. The MSEs of the isotonic propensity score matching estimator are considerably smaller than those of the logit matching estimator in every sample size. With the sample size growing, the MSEs of isotonic propensity score matching estimator approaches to the semiparametric efficiency bound.
Consider the following setting:
where $X\sim U[-1,1]^{3}$, and the true parameters are set as $\alpha_{0}=(1,1,1)^{\prime}/\sqrt{3}$, and $\gamma_{0}=(0.1,0.2,0.3)^{\prime}$, and the ATE is $\tau_{0}=0.5$. Under this setting, we have $\mathbb{P}(W=1|X=x)=p_{0}(x)=\Phi(x^{\prime}\alpha_{0})$, where $\Phi$ is the CDF of the standard normal distribution, i.e., the propensity score is correctly specified in probit estimation.
The simulation results are presented in Table (ref), where $\hat{\mu}_{\tau}$ is the Monte-Carlo mean, and the MSEs are rescaled by $N$. The number of Monte-Carlo simulations is 5000 for each sample size. The left panel shows the simulation results of the proposed matching method based on propensity scores estimated by the UC-iso-index estimator, and the right panel shows those of the one-to-one matching estimator based on propensity scores estimated with the correctly specified probit model.
The pattern is similar to the univariate case. The biases of both estimators are small and converge to zero. The isotonic matching estimator outperforms the probit matching estimator in every sample size in terms of MSE.
Table (ref) shows the bootstrap coverage rates. We draw 2000 Monte-Carlo simulations, and for each simulation, we draw 500 bootstrap samples. The coverage rates are calculated with these 2000 sets of confidence intervals for both 90% and 95% confidence levels. From Table (ref), we see clear trends that the bootstrap coverage rates are converging to their theoretical limits.
Overall, the simulation outcomes of the univariate case, the multivariate case, and the bootstrap encourage the proposed isotonic propensity score matching method. Additionally, for further simulation comparisons of our approach with propensity score methods of one-to-many matching and radius matching, as well as the impact of thresholds for averaging treatment variables at boundaries, see Section (ref) in the supplementary material.
We develop a one-to-many matching estimator of ATE based on propensity scores estimated by modified isotonic regression. We reveal that the nature of the isotonic estimator can help us to fix many problems of existing matching methods, including efficiency, choice of the number of matches, choice of tuning parameter, robustness to the propensity score misspecification, and bootstrap validity. As by-products, a uniformly consistent isotonic estimator and a uniformly consistent monotone single-index estimator, for both univariate and multivariate cases, are designed for our proposed isotonic matching estimator, and we study their asymptotic properties. The method can be further extended to other causal estimators based on propensity scores, such as blocking on propensity scores and regression on propensity scores.