The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
79,147 characters
Isotonic propensity score matching
\title{Isotonic propensity score matching\thanks{We are grateful to Markus Frölich, Daniel Gutknecht, Phillip Heiler, Lihua Lei, Yoshi Rai, Christoph Rothe, Carsten Trenkler, and participants at the econometrics seminar at Mannheim 2022, NASMES 2023, and IAAE 2023, for helpful comments and discussions. We also would like to thank a co-editor and anonymous referees for helpful comments to revise the paper.}}
\author{Mengshan Xu\thanks{Department of Economics, University of Mannheim, L7 3-5, 68161, Mannheim, Germany. Email: [email removed]} and Taisuke Otsu\thanks{Department of Economics, London School of Economics, Houghton Street, London, WC2A 2AE, UK. Email: [email removed]}}
\maketitle
\begin{abstract}
We propose a one-to-many matching estimator of the average treatment
effect based on propensity scores estimated by isotonic regression.
This approach is predicated on the assumption of monotonicity in the
propensity score function, a condition that can be justified in many
economic applications. We show 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 parameters,
robustness to propensity score misspecification, and bootstrap validity.
As a by-product, a uniformly consistent isotonic estimator is developed
for our proposed matching method.
\end{abstract}
\section{Introduction}
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{subsec:Uniformly-consistent-isotonic}.
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{sec:Multi}.
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 \emph{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 \emph{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
\emph{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
\begin{equation}
Y=\begin{cases}
1 & \text{if }X^{\prime}\beta_{0}>\varepsilon\\
0 & \text{if }X^{\prime}\beta_{0}\leq\varepsilon
\end{cases},\label{eq:binary}
\end{equation}
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)$, \eqref{eq:binary}
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 \eqref{eq:binary}
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 \eqref{eq:binary}.\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{subsec:Uni-rate-of=000020II},
one of the key conditions for $\sqrt{N}$-consistency and efficient
estimation of ATE is \eqref{eq:uni_rate_II} 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 \emph{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 \emph{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{subsec:DML}.
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 \emph{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{sec:Multi} 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{sec:main} shows the implementation
and asymptotic properties of the proposed isotonic matching estimator
with a univariate covariate. Section \ref{sec:Comparison} 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{sec:Multi},
where the propensity score is modeled by a semiparametric single-index
model with an unknown monotone increasing link function. In Section
\ref{sec:Boot}, we establish the validity of the nonparametric bootstrap.
Monte-Carlo simulation studies are presented in Section \ref{sec:Monte-Carlo}.
All proofs are presented in Appendix, while additional theoretical
details and simulation comparisons are provided in Supplementary Material.
\section{Main results\protect\label{sec:main}}
\subsection{Setup and isotonic propensity score\protect\label{subsec:setup}}
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{sec:Multi}. 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
\begin{equation}
\hat{p}(\cdot)=\arg\min_{p\in\mathcal{M}}\sum_{i=1}^{N}\{W_{i}-p(X_{i})\}^{2},\label{eq:iso}
\end{equation}
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 \emph{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{eq:uc_isoton})
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
\begin{equation}
\hat{p}(X_{n_{1}})<\hat{p}(X_{n_{2}})<\cdots<\hat{p}(X_{n_{K}}),\label{eq:unique}
\end{equation}
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.
\begin{asm} \label{asm:1} {[}Sampling{]} $\{Y_{i},W_{i},X_{i}\}_{i=1}^{N}$
is an independent and identically distributed (iid) sample of $(Y,W,X)\in\mathbb{R}\times\{0,1\}\times\mathcal{X}$,
where $\mathcal{X}=[x_{L},x_{U}]\in\mathbb{R}$. $X$ is continuously
distributed, and the sample $\{Y_{i},W_{i},X_{i}\}_{i=1}^{N}$ is
indexed according to $X_{1}<X_{2}<\cdots<X_{N}$. \end{asm}
\begin{prop} \label{prop:disjoint-grouping} Under Assumption \ref{asm:1},
the isotonic estimator $\hat{p}(\cdot)$ satisfying \eqref{eq:unique}
partitions the sample into $K$ disjoint groups in the sense that
for each $k=1,\ldots,K$ and $i=n_{k},\ldots,n_{k}+N_{k}-1$,
\begin{equation}
\hat{p}(X_{i})=\frac{1}{N_{k}}\sum_{j=n_{k}}^{n_{k}+N_{k}-1}W_{j}.\label{eq:iso_estimate_k}
\end{equation}
\end{prop}
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{prop:disjoint-grouping}.
\begin{prop} \label{prop:matching-score}\textbf{ }Suppose Assumption
\ref{asm:1} holds.
\begin{description}
\item [{(i)}] {[}Isotonic estimator{]} For each integer $k=1,\ldots,K$,
the isotonic estimator $\hat{p}(\cdot)$ for $p(\cdot)$ is represented
as
\begin{equation}
\hat{p}(x)=\frac{N_{k,1}}{N_{k}},\label{eq:PS_score}
\end{equation}
for each $x\in\{X_{i}\}_{i=n_{k}}^{n_{k}+N_{k}-1}$.
\item [{(ii)}] {[}Existence of matching counterparts{]} For any $i\in\{1,\ldots,N\}$
with $0<\hat{p}(X_{i})<1$, the set of its matching counterparts $\{j:W_{j}=1-W_{i},\hat{p}(X_{j})=\hat{p}(X_{i})\}$
is non-empty.
\end{description}
\end{prop}
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{asm:3} 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}$.
\subsection{Uniformly consistent isotonic estimator\protect\label{subsec:Uniformly-consistent-isotonic}}
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$,
\eqref{eq:PS_score} 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:
\begin{equation}
\tilde{W}_{i}=\begin{cases}
\frac{1}{\lfloor N^{2/3}\rfloor}\sum_{i=1}^{\lfloor N^{2/3}\rfloor}W_{i} & \text{for }i\leq\lfloor N^{2/3}\rfloor\\
W_{i} & \text{for }\lfloor N^{2/3}\rfloor<i\leq N-\lfloor N^{2/3}\rfloor\\
\frac{1}{\lfloor N^{2/3}\rfloor}\sum_{i=N-\lfloor N^{2/3}\rfloor+1}^{N}W_{i} & \text{for }i>N-\lfloor N^{2/3}\rfloor.
\end{cases}\label{eq:data_transform}
\end{equation}
Our proposed UC-isotonic estimator is obtained by implementing the
standard isotonic regression of $\tilde{W}$ on $X$:
\begin{equation}
\tilde{p}(x)=\begin{cases}
\max_{s\le i}\min_{t\ge i}\sum_{j=s}^{t}\tilde{W}_{j}/(t-s+1) & \text{for }x=X_{i}\\
\tilde{p}(X_{i}) & \text{for }X_{i-1}<x\leq X_{i}\\
\tilde{p}(X_{N}) & \text{for }x>X_{N}.
\end{cases}\label{eq:uc_isoton}
\end{equation}
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.
\begin{asm} \label{asm:2} {[}Monotonicity and continuity{]} (i)
$p(x)=\mathbb{E}[W|X=x]$ is a monotone increasing function of $x\in\mathcal{X}$,
(ii) $p(x)$ is continuously differentiable with its first derivative
$p^{(1)}(x)>0$ for all $x\in\mathcal{X}$, and (iii) $X$ has a continuous
density $f(x)$ satisfying that for some positive constants $\overline{f}$
and $\underline{f}$, it holds $\underline{f}<f(x)<\overline{f}$
all $x\in\mathcal{X}$. \end{asm}
Assumption \ref{asm:2} (i) is our main assumption, the monotonicity
of $p(\cdot)$. Assumption \ref{asm:2} (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 \eqref{eq:ebar-C} in Appendix \ref{appsub:Proof-of-T_univariate},
$|\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 \eqref{eq:ebar-C}.
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{asm:2} (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\emph{ }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{prop:matching-score} for
the UC-isotonic estimator.
\begin{prop} \label{prop:matching-score-modified} Under Assumptions
\ref{asm:1} and \ref{asm:2}, it holds
\begin{description}
\item [{(i)}] $N_{1}\geq\lfloor N^{2/3}\rfloor$ and $N_{K}\geq\lfloor N^{2/3}\rfloor$.
\item [{(ii)}] $\tilde{p}(x)=\frac{N_{k,1}}{N_{k}}$ for each $k=1,\dots,K$
and $x\in\{X_{i}\}_{i=n_{k}}^{n_{k}+N_{k}-1}.$
\item [{(iii)}] $N_{1}=O_{p}(N^{2/3})$ and $N_{K}=O_{p}(N^{2/3}).$
\end{description}
\end{prop}
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.
\begin{thm} \label{thm:uc_isoton} Under Assumptions \ref{asm:1}
and \ref{asm:2}, it holds
\[
\sup_{x\in\mathcal{X}}|\tilde{p}(x)-p(x)|=O_{p}\left(\frac{\log N}{N}\right)^{1/3}.
\]
\end{thm}
Finally, to guarantee the existence of matching counterparts by $\tilde{p}(\cdot)$,
we impose the strict overlap condition.
\begin{asm} \label{asm:3} {[}Strict overlaps{]} There exist positive
constants $\underline{p}$ and $\bar{p}$ such that $0<\underline{p}\leq p(x)\leq\bar{p}<1$
for all $x\in\mathcal{X}$. \end{asm}
Assumption \ref{asm:3} is standard in the treatment effect literature.
It is necessary for the identification and $\sqrt{N}$-consistent
estimation of the ATE. Combining Proposition \ref{prop:matching-score-modified}
and Theorem \ref{thm:uc_isoton} with Assumption \ref{asm:3}, the
existence of the matching counterparts by $\tilde{p}(\cdot)$ is obtained
as follows.
\begin{cor} \label{cor:non-empty-matched-set} Suppose Assumptions
\ref{asm:1}-\ref{asm:3} hold. For each $i=1,\dots,N$, the set of
its matching counterparts $\{j=1,\dots,N:W_{j}=1-W_{i},\tilde{p}(x_{j})=\tilde{p}(x_{i})\}$
is non-empty with probability approaching one. \end{cor}
\subsection{Isotonic propensity score matching\protect\label{subsec:Isotonic-matching}}
Based on the UC-isotonic estimator $\tilde{p}(\cdot)$, the isotonic
propensity score matching estimator for the ATE $\tau$ can be implemented
as follows.
\begin{enumerate}
\item Transform the sample $\{Y_{i},W_{i},X_{i}\}_{i=1}^{N}$ indexed by
$X_{1}<\cdots<X_{N}$ into $\{Y_{i},\tilde{W}_{i},X_{i}\}_{i=1}^{N}$
using \eqref{eq:data_transform}.
\item Compute the UC-isotonic estimator $\tilde{p}(\cdot)$ using \eqref{eq:uc_isoton}.
\item For each $i=1,\ldots,N$, compute the matching counterparts
\begin{equation}
\mathcal{J}(i)=\left\{ j=1,\dots,N:W_{j}=1-W_{i}\text{ and }\tilde{p}(X_{j})=\tilde{p}(X_{i})\right\} .\label{eq:matching}
\end{equation}
\item Calculate the matching estimator for the ATE $\tau$ by
\begin{equation}
\hat{\tau}=\frac{1}{N}\sum_{i=1}^{N}(2W_{i}-1)\left(Y_{i}-\frac{1}{M_{i}}\sum_{j\in\mathcal{J}(i)}Y_{j}\right),\label{eq:ATE=000020sample}
\end{equation}
where $M_{i}=|\mathcal{J}(i)|$ is the number of matches for $i$.
\end{enumerate}
We proceed with the following assumptions.
\begin{asm} \label{asm:4} {[}Data generating process{]} (i) $\mathbb{E}[Y(0)^{2}]<\infty$
and $\mathbb{E}[Y(1)^{2}]<\infty$, (ii) $\mathbb{E}[Y(0)|X=x]$ and
$\mathbb{E}[Y(1)|X=x]$ are continuously differentiable for all $x\in\mathcal{X}$,
(iii) for $D(Z)=\frac{WY}{p(X)^{2}}+\frac{Y(1-W)}{\{1-p(X)\}^{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{asm}
Assumption \ref{asm:4} (i)-(iii) regulates the tail behaviors of
the (conditional functions of) potential outcomes, which are necessary
for $\sqrt{N}$-consistent estimation. Assumption \ref{asm:4} (iv)
is the standard unconfoundedness assumption. Under these assumptions,
we have the following key equivalence result.
\begin{thm} \label{thm:matching_IPW} Under Assumptions \ref{asm:1}-\ref{asm:4},
the matching estimator $\hat{\tau}$ for the ATE $\tau$ using the
UC-isotonic estimator $\tilde{p}(\cdot)$ is equal to the corresponding inverse probability weighting (IPW) estimator. \end{thm}
\subsubsection*{Remark on Theorem \ref{thm:matching_IPW}}
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 \emph{et al.}, 2002; Cattaneo and Farrell, 2013). See Section
\ref{subsec:link-partition} 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{subsec:block-strat-radius} for a comparison with
these methods.
Moreover, as mentioned in the introduction, the equivalence result
in Theorem \ref{thm:matching_IPW} relies crucially on the implementation
of the UC-isotonic estimator \eqref{eq:uc_isoton}, 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
\eqref{eq:data_transform} 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{supp:thresh_theorey} and \ref{supp:thresh_sim}
of the supplementary material, respectively.
Our main result, consistency and asymptotic normality of the isotonic
propensity score matching estimator, is obtained as follows.
\begin{thm} \label{thm:ATE-matching} Under Assumptions \ref{asm:1}-\ref{asm:4},
it holds $\hat{\tau}\overset{p}{\to}\tau$ and
\[
\sqrt{N}(\hat{\tau}-\tau)\overset{d}{\to}N(0,\Omega),
\]
where $\Omega=\mathbb{V}(\mathbb{E}[Y(1)-Y(0)|X])+\mathbb{E}[\mathbb{V}(Y(1)|X)/p(X)]+\mathbb{E}[\mathbb{V}(Y(0)|X)/(1-p(X))]$.
\end{thm}
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{sec:Boot}.
\section{Comparison to related propensity score methods\protect\label{sec:Comparison}}
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{subsec:link-partition}, the comparison with propensity
score methods of blocking, stratification, and radius matching is
presented in Section \ref{subsec:block-strat-radius}, the comparison
with matching methods based on propensity score estimated by regression
trees is presented in Section \ref{subsec:trees}, and the comparison
with the double machine learning (DML) estimator for the ATE can be
found in Section \ref{subsec:DML}.
\subsection{Propensity score estimated by partitioning estimator\protect\label{subsec:link-partition}}
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 \eqref{eq:matching}. 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 \emph{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 \eqref{eq:IPW_tau_ATE}
in Appendix \ref{appsub:equivalent-result}. 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{thm:ATE-matching}.
\subsection{Propensity score methods of blocking, stratification, and radius
matching\protect\label{subsec:block-strat-radius}}
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{thm:matching_IPW}, 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.
\subsection{Propensity score estimated by regression trees\protect\label{subsec:trees}}
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{sec:Multi} 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 \emph{et al.} (2018). See Section \ref{subsec:DML}
below for more details about the comparison of our approach with the
double machine learning estimator.
\subsection{Double machine learning estimator\protect\label{subsec:DML}}
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:
\begin{eqnarray}
& & \frac{1}{N}\sum_{i=1}^{N}\{\hat{\psi}_{1}(X_{i})-\hat{\psi}_{0}(X_{i})\}+\frac{1}{N}\sum_{i=1}^{N}\left[\frac{W_{i}(Y_{i}-\hat{\psi}_{1}(X_{i}))}{\hat{p}(X_{i})}-\frac{(1-W_{i})(Y_{i}-\hat{\psi}_{0}(X_{i}))}{1-\hat{p}(X_{i})}\right]\nonumber \\
& = & \frac{1}{N}\sum_{i=1}^{N}\left[\frac{Y_{i}W_{i}}{\hat{p}(X_{i})}-\frac{Y_{i}(1-W_{i})}{1-\hat{p}(X_{i})}\right]-\frac{1}{N}\sum_{i=1}^{N}\left[\frac{W_{i}-\hat{p}(X_{i})}{\hat{p}(X_{i})}\hat{\psi}_{1}(X_{i})-\frac{W_{i}-\hat{p}(X_{i})}{1-\hat{p}(X_{i})}\hat{\psi}_{0}(X_{i})\right],\label{eq:AIPW}
\end{eqnarray}
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 \eqref{eq:AIPW} 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 \emph{et al}. (2018) combined orthogonal moment functions
– of which the formula \eqref{eq:AIPW} 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 \emph{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 \eqref{eq:uni_rate_II}
in Appendix). In contrast, the DML for the ATE represents a typical
bias correction approach. The second terms in both lines of \eqref{eq:AIPW},
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.
\section{Multivariate covariates \protect\label{sec:Multi}}
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:
\begin{equation}
W=p_{0}(X^{\prime}\alpha_{0})+\varepsilon,\qquad\mathbb{E}[\varepsilon|X]=0,\label{eq:=000020single=000020index}
\end{equation}
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 \eqref{eq:binary},
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 \eqref{eq:=000020single=000020index} 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
\begin{equation}
\mathbb{E}[X\{W-p_{0}(X^{\prime}\alpha_{0})\}]=0.\label{eq:=000020moment-single=000020index}
\end{equation}
To estimate $p_{0}$ and $\text{\ensuremath{\alpha}}_{0}$, we can
apply the method of BGH. For a fixed $\alpha$, define
\begin{equation}
\hat{p}_{\alpha}=\arg\min_{p\in\mathcal{M}}\frac{1}{N}\sum_{i=1}^{N}\{W_{i}-p(X_{i}^{\prime}\alpha)\}^{2},\label{eq:link_mono_index}
\end{equation}
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
\begin{equation}
\hat{\alpha}=\text{arg}\underset{\alpha}{\text{min}}\left\Vert \frac{1}{N}\sum_{i=1}^{N}X_{i}^{\prime}\{W_{i}-\hat{p}_{\alpha}(X_{i}^{\prime}\alpha)\}\right\Vert ^{2}.\label{eq:para_mono_index}
\end{equation}
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 \eqref{eq:para_mono_index},
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{subsec:Uniformly-consistent-isotonic},
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.
\begin{enumerate}
\item Compute $\hat{\alpha}$ by \eqref{eq:link_mono_index} and \eqref{eq:para_mono_index}.
\item Define $\tilde{\alpha}=\hat{\alpha}$, and transform the sample $\{Y_{i},W_{i},X_{i}\}_{i=1}^{N}$
indexed by $X_{1}^{\prime}\tilde{\alpha}<\cdots<X_{N}^{\prime}\tilde{\alpha}$
into $\{Y_{i},\tilde{W}_{i},X_{i}\}_{i=1}^{N}$ with \eqref{eq:data_transform}.
\item Compute the UC-iso-index estimator $\tilde{p}_{\tilde{\alpha}}$ by
\[
\tilde{p}_{\tilde{\alpha}}=\arg\min_{p\in\mathcal{M}}\frac{1}{N}\sum_{i=1}^{N}\{\tilde{W}_{i}-p(X_{i}^{\prime}\tilde{\alpha})\}^{2}.
\]
\item For each $i=1,\ldots,N$, compute the matching counterparts
\begin{equation}
\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\} .\label{eq:matching-multi}
\end{equation}
\item Calculate the matching estimator for the ATE $\tau$ by
\begin{equation}
\tilde{\tau}=\frac{1}{N}\sum_{i=1}^{N}(2W_{i}-1)\left(Y_{i}-\frac{1}{M_{i}}\sum_{j\in\mathcal{J}(i)}Y_{j}\right),\label{eq:ATE=000020sample-multi}
\end{equation}
where $M_{i}=|\mathcal{J}(i)|$ is the number of matches for $i$.
\end{enumerate}
We modify Assumptions \ref{asm:1}-\ref{asm:4} in Section \ref{sec:main}
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 \eqref{eq:link_mono_index}:
\[
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}$.\emph{ }\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 \eqref{eq:link_mono_index} and \eqref{eq:para_mono_index}
have desirable properties.
\begin{asm} \label{asm:5} For all $\alpha\ne\alpha_{0}$ such that
$\alpha\in\mathcal{B}(\alpha_{0},\delta_{0})$, the random variable
$\text{Cov}[(\alpha-\alpha_{0})^{\prime}X,p_{0}(X^{\prime}\alpha_{0})|X^{\prime}\alpha]$
is not equal to 0 almost surely. \end{asm}
\begin{asm} \label{asm:6} {[}Potential outcomes{]} Let $p_{0}^{(1)}(u)$
denote the first derivative of $p_{0}(u)$. The matrix $\mathbb{E}[p_{0}^{(1)}(X^{\prime}\alpha_{0})\text{Cov}(X|X^{\prime}\alpha_{0})]$
has rank $k-1$. \end{asm}
Based on Assumptions 1', 2', \ref{asm:5}, and \ref{asm:6}, we have
a result similar to Proposition \ref{prop:matching-score-modified},
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.
\begin{thm} \label{thm:uc_isoton-m} Under Assumptions 1', 2', \ref{asm:5},
and \ref{asm:6}, it holds
\[
\sup_{x\in\mathcal{X}}|\tilde{p}_{\tilde{\alpha}}(x^{\prime}\tilde{\alpha})-p_{0}(x^{\prime}\alpha_{0})|=O_{p}\left(\frac{\log N}{N}\right)^{1/3}.
\]
\end{thm}
The existence of matching counterparts is guaranteed by an argument
similar to Corollary \ref{cor:non-empty-matched-set}. 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.
\begin{thm} \label{thm:ATE-matching-m} Under Assumptions \ref{asm:1}'-\ref{asm:4}',
\ref{asm:5}, and \ref{asm:6}, it holds $\tilde{\tau}\overset{p}{\to}\tau$
and
\[
\sqrt{N}(\tilde{\tau}-\tau)\overset{d}{\to}N(0,\Sigma),
\]
where $\Sigma=\mathbb{E}[\{m(Z)+M(Z)+A(Z)\}\{m(Z)+M(Z)+A(Z)\}]$,
and
\begin{eqnarray}
m(Z) & = & \frac{YW}{p_{0}(X^{\prime}\alpha_{0})}-\frac{Y(1-W)}{1-p_{0}(X^{\prime}\alpha_{0})}-\tau,\qquad D(Z)=\frac{YW}{p_{0}(X^{\prime}\alpha_{0})^{2}}+\frac{Y(1-W)}{(1-p_{0}(X^{\prime}\alpha_{0}))^{2}},\nonumber \\
M(Z) & = & -\mathbb{E}_{\alpha_{0}}[D(Z)|X^{\prime}\alpha_{0}]\{W-p_{0}(X^{\prime}\alpha_{0})\},\nonumber \\
A(Z) & = & -\mathbb{E}[\mathrm{Cov}_{\alpha_{0}}(D(Z),X|X^{\prime}\alpha_{0})p_{0}^{(1)}(X^{\prime}\alpha_{0})],\nonumber \\
& & \times\mathbb{E}[p_{0}^{(1)}(X^{\prime}\alpha_{0})\mathrm{Cov}_{\alpha_{0}}(X|X^{\prime}\alpha_{0})]^{-}\{X-\mathbb{E}_{\alpha_{0}}[X|X^{\prime}\alpha_{0}]\}\{W-p_{0}(X^{\prime}\alpha_{0})\}.\label{eq:variance=000020component}
\end{eqnarray}
\end{thm}
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{subsec:Monte-Carlo-m}
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.
\section{Bootstrap inference\protect\label{sec:Boot}}
The asymptotic variances in Theorems \ref{thm:ATE-matching} and \ref{thm:ATE-matching-m}
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 \emph{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{thm:matching_IPW}.
In this section, we discuss an asymptotically valid bootstrap procedure
for the estimator $\hat{\tau}$ in Theorem \ref{thm:ATE-matching}.
This result can be similarly adapted to $\tilde{\tau}$ in Theorem
\ref{thm:ATE-matching-m}.
The nonparametric bootstrap is implemented as follows.
\begin{enumerate}
\item $\{Y_{i}^{*},W_{i}^{*},X_{i}^{*}\}_{i=1}^{N}$ is a bootstrap sample
from $\{Y_{i},W_{i},X_{i}\}_{i=1}^{N}$, and the numbering is according
to $X_{1}^{*}\leq\cdots\leq X_{N}^{*}.$
\item $\tilde{p}^{*}(\cdot)$ is the UC-isotonic estimator based on $\{Y_{i}^{*},W_{i}^{*},X_{i}^{*}\}_{i=1}^{N}$.
\item The bootstrap counterpart $\hat{\tau}^{*}$ of $\hat{\tau}$ is given
by
\begin{eqnarray*}
\hat{\tau}^{*} & = & \frac{1}{N}\sum_{i=1}^{N}(2W_{i}^{*}-1)\left(Y_{i}^{*}-\frac{1}{M_{i}^{*}}\sum_{j\in\mathcal{J^{*}}(i)}Y_{j}^{*}\right),\\
\mathcal{J^{*}}(i) & = & \left\{ j=1,\dots,N:W_{j}^{*}=1-W_{i}^{*}\text{ and }\tilde{p}^{*}(X_{j}^{*})=\tilde{p}^{*}(X_{i}^{*})\right\} ,
\end{eqnarray*}
where $M_{i}^{*}=|\mathcal{J^{*}}(i)|$ is the number of matches for
the $i$-th observation in the bootstrap sample.
\item After repeating Step (1)-(3) for $B$ times and obtaining estimator
$\hat{\tau}_{1}^{*},\hat{\tau}_{2}^{*},\dots,\hat{\tau}_{B}^{*},$
we can conduct inference for $\tau$.
\end{enumerate}
The asymptotic validity of this bootstrap approximation is obtained
as follows.
\begin{thm}\label{thm:=000020boot} Let $\mathbb{P}^{*}$ be the
bootstrap distribution conditional on the data, and $c_{1-\alpha}^{*}$
be the $(1-\alpha)$-th sample quantile of $(\sqrt{N}\left(\hat{\tau}_{1}^{*}-\hat{\tau}\right),\sqrt{N}\left(\hat{\tau}_{2}^{*}-\hat{\tau}\right),\dots,\sqrt{N}\left(\hat{\tau}_{B}^{*}-\hat{\tau}\right))$.
Under Assumptions \ref{asm:1}-\ref{asm:4}, it holds
\begin{itemize} \item[(i)] $\sup_{t\in\mathbb{R}}|\mathbb{P}^{*}\{\sqrt{N}(\hat{\tau}^{*}-\hat{\tau})\le t\}-\mathbb{P}\{\sqrt{N}(\hat{\tau}-\tau)\le t\}|\overset{p}{\to}0;$
\item[(ii)] $\mathbb{P}\{\sqrt{N}(\hat{\tau}-\tau)\le c_{1-\alpha}^{*}\}\overset{p}{\to}1-\alpha$.\end{itemize}
\end{thm}
\section{Monte-Carlo simulations \protect\label{sec:Monte-Carlo}}
In this section, we use three simulation studies to assess the finite
sample properties of our isotonic propensity score matching estimator.
\subsection{Univariate case\protect\label{subsec:sim_uni}}
Let $X=0.15+0.7Z$, where $Z$ and $\nu$ are independently uniformly
distributed on $[0,1]$, and
\begin{eqnarray}
W & = & \begin{cases}
0 & \text{if }X<\nu\\
1 & \text{if }X\geq\nu
\end{cases},\nonumber \\
Y & = & 0.5W+2X+\varepsilon,\nonumber \\
\varepsilon & \sim & N(0,1).\label{sim:uni_setup}
\end{eqnarray}
The true ATE is the coefficient of $W$, which is 0.5. The simulation
results are presented in Table \ref{tab:ATE-model}, 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.
\begin{table}[H]
\caption{\protect\label{tab:ATE-model}Matching estimators of ATE: the univariate
case}
\centering{}{\small{}
\begin{tabular}{cr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
\toprule
\multicolumn{5}{c}{with UC-isotonic} & \multicolumn{2}{c}{} & \multicolumn{6}{c}{with logit and $M=1$}\tabularnewline
\midrule
$N$ & \multicolumn{2}{c}{$\hat{\mu}_{\tau}$} & \multicolumn{2}{c}{$MSE$} & \multicolumn{2}{c}{} & \multicolumn{2}{c}{$N$} & \multicolumn{2}{c}{$\hat{\mu}_{\tau}$} & \multicolumn{2}{c}{$MSE$}\tabularnewline
\cmidrule{1-5}\cmidrule{8-13}
100 & 0&4977 & 5&2723 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{100} & 0&4997 & 7&1068\tabularnewline
1000 & 0&4934 & 5&2589 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{1000} & 0&5009 & 7&0630\tabularnewline
2000 & 0&4946 & 5&2158 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{2000} & 0&4999 & 7&0816\tabularnewline
5000 & 0&4963 & 4&9418 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{5000} & 0&4995 & 6&8376\tabularnewline
10000 & 0&4974 & 4&9785 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{10000} & 0&5000 & 6&8238\tabularnewline
\midrule
$\infty$ & 0&5 & 4&96 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{$\infty$} & 0&5 & 4&96\tabularnewline
\bottomrule
\end{tabular}}{\small\par}
\end{table}
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):
\begin{align*}
\Omega_{\text{SEB}} & =\text{Var}(\mathbb{E}[Y(1)-Y(0)|X])+\mathbb{E}[\text{Var}(Y(1)|X)/p_{0}(X)]+\mathbb{E}[\text{Var}(Y(0)|X)/(1-p_{0}(X))]\\
& =\text{Var}(0.5)+\mathbb{E}[1/p_{0}(X)]+\mathbb{E}[1/(1-p_{0}(X))]\\
& =0+\int_{0.15}^{0.85}\frac{1}{x}\frac{1}{0.7}dx+\int_{0.15}^{0.85}\frac{1}{1-x}\frac{1}{0.7}dx\approx4.96.
\end{align*}
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.
\subsection{Multivariate case\protect\label{subsec:Monte-Carlo-m}}
Consider the following setting:
\begin{eqnarray*}
Y & = & X^{\prime}\gamma_{0}+W\tau_{0}+\varepsilon,\\
W & = & \begin{cases}
0 & \text{if }X^{\prime}\alpha_{0}<\nu\\
1 & \text{if }X^{\prime}\alpha_{0}\geq\nu
\end{cases},\\
\varepsilon & \sim & N(0,1),\qquad\nu\sim N(0,1),\qquad\varepsilon\perp\nu,
\end{eqnarray*}
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{tab:ATE-model-m},
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.
\begin{table}[H]
\caption{\protect\label{tab:ATE-model-m}Matching estimators of ATE: the multivariate
case}
\centering{}{\small{}
\begin{tabular}{cr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
\toprule
\multicolumn{5}{c}{with UC-iso-index} & \multicolumn{2}{c}{} & \multicolumn{6}{c}{with probit and $M=1$}\tabularnewline
\midrule
$N$ & \multicolumn{2}{c}{$\hat{\mu}_{\tau}$} & \multicolumn{2}{c}{$MSE$} & \multicolumn{2}{c}{} & \multicolumn{2}{c}{$N$} & \multicolumn{2}{c}{$\hat{\mu}_{\tau}$} & \multicolumn{2}{c}{$MSE$}\tabularnewline
\cmidrule{1-5}\cmidrule{8-13}
100 & 0&5080 & 5&0442 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{100} & 0&5114 & 7&3459\tabularnewline
1000 & 0&5016 & 5&0014 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{1000} & 0&5030 & 6&9813\tabularnewline
2000 & 0&4991 & 5&0727 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{2000} & 0&4997 & 7&2275\tabularnewline
5000 & 0&5003 & 5&2115 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{5000} & 0&5010 & 7&2640\tabularnewline
10000 & 0&5001 & 5&0161 & \multicolumn{2}{c}{} & \multicolumn{2}{c}{10000} & 0&5002 & 7&0509\tabularnewline
\bottomrule
\end{tabular}}{\small\par}
\end{table}
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.
\subsection{Bootstrap}
Table \ref{tab:coverage} 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{tab:coverage}, we see clear trends that the bootstrap
coverage rates are converging to their theoretical limits.
\begin{table}[H]
\caption{\protect\label{tab:coverage}Bootstrap coverage rates}
\centering{}{\small{}
\begin{tabular}{ccrr}
\toprule
\multirow{1}{*}{$n$} & & 90\% CI & 95\% CI\tabularnewline
\cmidrule{1-1}\cmidrule{3-4}
100 & & 0.860 & 0.918\tabularnewline
1000 & & 0.889 & 0.938\tabularnewline
2000 & & 0.881 & 0.940\tabularnewline
5000 & & 0.901 & 0.945\tabularnewline
10000 & & 0.891 & 0.948\tabularnewline
\cmidrule{1-1}\cmidrule{3-4}
$\infty$ & & 0.90 & 0.95\tabularnewline
\bottomrule
\end{tabular}}{\small\par}
\end{table}
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{supp:add_Monte_Carlo}
in the supplementary material.
\section{Conclusion}
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.