Covariate Distribution Balance via Propensity Scores

Covariate Distribution Balance via Propensity Scores

{4pt} {4pt} {4pt} {4pt}

Introduction

Identifying and estimating the effect of a policy, treatment or intervention on an outcome of interest is one of the main goals in applied research. Although a randomized control trial (RCT) is the gold standard to identify causal effects, many times its implementation is infeasible and researchers have to rely on observational data. In such settings, the propensity score (PS), which is defined as the probability of being treated given observed covariates, plays a prominent role. Statistical methods using the PS include matching, inverse probability weighting (IPW), regression, as well as combinations thereof; for review, see, e.g., Imbens2015.

To use these methods in practice, one has to acknowledge that the PS is usually unknown and has to be estimated from the observed data. Given the moderate or high dimensionality of available covariates, researchers are usually coerced to adopt a parametric model for the PS. A popular approach is to assume a linear logistic model, estimate the unknown parameters by maximum likelihood (ML), check if the resulting PS estimates balance specific moments of covariates, and in case they do not, refit the PS model including higher-order and interaction terms and repeat the procedure until covariate balancing is achieved, see, e.g., Rosenbaum1984 and Dehejia2002. On top of involving ad hoc choices of model refinements, such model selection procedures may result in distorted inference about the parameters of interest, see, e.g., Leeb2005. An additional challenge faced by PS estimators based on ML is that the likelihood loss function does not take into account the covariate balancing property of the PS ( Rosenbaum1983 ), and, as a result, treatment effect estimators based of ML PS estimates can be very sensitive to model misspecifications, see, e.g., Kang2007.

In light of these practical issues, alternative estimation procedures that are able to resemble randomization in a closer fashion have been proposed. For instance, Graham2012, Hainmueller2012, Imai2014, Zubizarreta2015, and Zhao2018 propose alternative estimation procedures that attempt to directly balance covariates among the treated, untreated and, combined sample. Although such methods usually lead to treatment effect estimators with improved finite sample properties, they only aim to balance some specific functions of covariates. However, the covariate balancing property of the PS is considerably more powerful as it implies balance not only for some particular moments but for all measurable, integrable functions of the covariates. Indeed, the balancing property of the propensity score resembles randomization: when the data come from a randomized control trial (RCT) with perfect compliance, the entire covariate distributions among different treatment groups are balanced and, therefore, all measurable, integrable functions of the covariates are indeed balanced.

In this paper, we propose an alternative framework for estimating the PS that is arguably more suitable for causal inference, as it fully exploits the covariate balancing property of the PS. We call the resulting PS estimator the integrated propensity score (IPS). At a conceptual level, the IPS builds on the observation that the covariate balancing property of the PS can be equivalently characterized by balancing covariate distributions, namely, by an infinite, but tractable, number of unconditional moment restrictions. Upon such an observation, we consider Cram\'{e}r-von Mises-type distances between these infinite balancing conditions and zero, and show that their minima are uniquely achieved at the true PS parameters. These results, in turn, suggest that we can estimate the unknown PS parameters within the minimum distance framework, as in, for example, Dominguez2004 and Escanciano2006b, Escanciano2018. We emphasize that the IPS can be used under different \textquotedblleft research designs\textquotedblright, including not only the unconfounded treatment assignment setup, see, e.g., Rosenbaum1983, Hirano2003, Firpo2007, and Chen2008, but also the \textquotedblleft local treatment effect\textquotedblright\ setup, where selection into treatment is possibly endogenous but a binary instrumental variable is available, see, e.g., Abadie2003, and Frolich2013. In this latter case, the IPS aims to balance the covariates among the treated, non-treated, and overall complier subpopulations.

At the practical level, one can think of the IPS as an estimation procedure that attempts to estimate the unknown finite dimensional parameters of a PS model by making the underlying entire covariate distribution of different treatment groups as close to each other as possible. The IPS framework also acknowledges that, in practice, there are different ways to compare covariate distribution functions depending on how covariate distribution balance is measured and the norm chosen. We explicitly consider three natural ways to characterize covariate distribution balance: 1) using the covariates' joint cumulative distribution, 2) their joint characteristic function, or 3) exploiting the Cram\'{e}r--Wold theorem to focus on the cumulative distribution of the one-dimensional projections of the covariates. In terms of the norm, we focus on Cram\'{e}r-von Mises-type distances as they can lead to smooth criteria functions that admit a closed-form representation, allowing us to avoid using computationally heavy numerical integration procedures. In fact, our proposed method is computationally simple and easy to use as currently implemented in the new package IPS for R, available at \url{https://github.com/pedrohcgs/IPS}.

The proposed IPS enjoys several appealing properties. First, the IPS procedure guarantees that the unknown PS parameters are globally identified. This is in contrast to the traditional generalized method of moments approach based on only finitely many balancing conditions, see, e.g., Hellerstein1999 and Dominguez2004. Second, even though we aim to balance an infinite number of balancing conditions, the IPS estimator does not rely on tuning parameters. Third, the IPS does not rely on outcome data and separates the design stage (where one estimates the propensity score) from the analysis stage (where one estimates different treatment effect measures). As advocated by Rubin2007, Rubin2008, this separation is useful as it simultaneously mimics RCTs and avoids potential data snooping. Another direct consequence of this clear separation is that one can use the IPS to estimate a variety of causal effect parameters in a relatively straightforward manner. We illustrate this flexibility by deriving the asymptotic properties of inverse probability weighted (IPW) estimators for average, distributional and quantile treatment effects based on the IPS, under both the unconfoundedness and the local treatment effects setups.

Related literature: Our proposal builds on different branches of the econometrics literature. For instance, this paper is related to Shaikh2009 and SantAnna2018, who exploit the covariate balancing of the PS to propose specification tests for a given PS model. Here, instead of checking if a given PS estimator balances the covariate distribution among different treatment groups, we propose to estimate the PS unknown parameters by maximizing the covariate balancing. The IPS estimators also build on Dominguez2004 and Escanciano2006b, Escanciano2018, who propose generic estimation procedures for finite-dimensional parameters defined via an infinite number of unconditional moment restrictions. Upon characterizing the covariate balancing property of the PS as an infinite number of unconditional moment restrictions, we are able to adapt their proposals to our causal inference context.

Our proposal is also related to the growing literature on weighting-based covariate balancing methods. Among this branch of the literature, the closest papers to ours are Graham2012, Imai2014, Diaz2015 and Fan2016a. An important difference between our proposal and theirs is that all these papers focus exclusively on average treatment effects under unconfoundedness, whereas we show that one can directly use the IPS to estimate a variety of causal parameters of interest such as average, quantile and distributional treatment effects, not only under unconfoundedness but also in settings with endogenous treatment. It is also worth stressing that Graham2012 and Imai2014 propose estimating PS by balancing some specific pre-determined moments of the covariates, and that their procedure requires one to assume that the propensity score parameters are uniquely (globally) identified, see, e.g., Assumption 2.1(i) in Graham2012. In practice, it is hard to verify such important condition, and when such assumption is not satisfied, inference procedures based on their proposal will in general not be valid, see, e.g. Dominguez2004. Our proposed IPS procedure, on the other hand, does not suffer from this drawback as it aims to balance the entire covariate distribution, i.e., our proposal is based on an infinite number of balancing conditions that fully characterize the propensity score.

In a recent working paper, Fan2016a consider the case where the number of balancing moments grows with the sample size at an appropriate rate. Although this proposal bypass the identification challenge mentioned above (see, e.g., Ai2003 and Donald2003), to implement their proposal one needs to carefully choose tuning parameters and select basis functions such that the resulting balancing moments are guaranteed to be finite. Their proposal also (implicitly) relies on covariates having compact support. Our proposal avoids these practical complications.

Organization of the paper: Section (ref) introduces the framework of balancing weights and explains the estimation problem of the IPS. Section (ref) presents the large sample properties of the IPS estimator. This section also discusses how one can use the IPS to estimate and make inference about average, distributional and quantile treatment effects under the unconfoundedness assumption. In Section (ref), we discuss how one can use the IPS in the empirically relevant situation where treatment adoption is endogenous and one has access to a binary instrumental variable. Section (ref) illustrates the comparative performance of the proposed method through simulations. Section (ref) presents two empirical applications. Section (ref) concludes. Proofs, as well as additional results, are reported in the Supplemental Appendix\footnote{The Supplemental Appendix is available at \url{https://pedrohcgs.github.io/files/IPS-supplementary.pdf}}.

Covariate balancing via propensity score

Background

Let $D$ be a binary random variable that indicates participation in the program, i.e., $D=1$ if the individual participates in the treatment and $D=0$ otherwise. Define $Y\left( 1\right) $ and $Y\left( 0\right) $ as the potential outcomes under treatment and untreated, respectively. The realized outcome of interest is $Y=DY\left( 1\right) +\left( 1-D\right) Y\left( 0\right) $, and $\mathbf{X}$\ is an observable $k\times1$ vector of pre-treatment covariates.\ Denote the support of $\mathbf{X}$ by $\mathcal{X\subset}\mathbb{R}^{k}$ and the propensity score $p\left( \mathbf{x}\right) =\mathbb{P}\left( D=1|\mathbf{X}=\mathbf{x}\right) $. For $d\in\left\{ 0,1\right\} $, denote the distribution and quantile of the potential outcome $Y\left( d\right) $ by $F_{Y\left( d\right) }\left( y\right) =\mathbb{P}\left( Y\left( d\right) \leq y\right) $, and $q_{Y\left( d\right) }\left( \tau\right) =\inf\left\{ y:F_{Y\left( d\right) }\left( y\right) \geq\tau\right\} $, respectively, where $y\in\mathbb{R}$ and $\tau\in\left( 0,1\right) $. Henceforth, assume that we have a random sample $\left\{ \left( Y_{i},D_{i},\mathbf{X}_{i}^{\prime }\right) ^{\prime}\right\} _{i=1}^{n}$ from $\left( Y,D,\mathbf{X}^{\prime }\right) ^{\prime},$ where $n\geq1$ is the sample size, and all random variables are defined on a common probability space $\left( \Omega ,\mathcal{A},\mathbb{P}\right) .$ For a generic random variable $Z$, denote $\mathbb{E}_{n}\left[ Z\right] =n^{-1}\sum_{i=1}^{n}Z_{i}$.

The main goal in causal inference is to assess the effect of a treatment $D$ on the outcome of interest $Y$. Perhaps the most popular causal parameter of interest is the overall average treatment effect, $ATE=\mathbb{E}\left[ Y\left( 1\right) -Y\left( 0\right) \right] $. Despite its popularity, the ATE can mask important treatment effect heterogeneity across different subpopulations, see, e.g., Bitler2006. Thus, in order to uncover potential treatment effect heterogeneity, one usually focuses on different treatment effect parameters beyond the mean. Leading examples include the overall distributional treatment effect, $DTE\left( y\right) =F_{Y\left( 1\right) }\left( y\right) -F_{Y\left( 0\right) }\left( y\right) $, and the overall quantile treatment effect, $QTE\left( \tau\right) =q_{Y\left( 1\right) }\left( \tau\right) -q_{Y\left( 0\right) }\left( \tau\right) $. Given that these causal parameters depend on potential outcomes that are not jointly observed for the same individual, one cannot directly rely on the analogy principle to identify and estimate such functionals.

A commonly used identification strategy in policy evaluation to bypass this difficulty is to assume that selection into treatment is based on observable characteristics, and that all individuals have a positive probability of being in either the treatment or the untreated group --- the so-called unconfoundedness setup, see, e.g., Rosenbaum1983. Formally, unconfoundedness requires the following assumption.

Rosenbaum1987a shows that, under Assumption (ref), the ATE is identified by \[ ATE=\mathbb{E}\left[ \left( \frac{D}{p\left( \mathbf{X}\right) } -\frac{\left( 1-D\right) }{1-p\left( \mathbf{X}\right) }\right) Y\right] . \] Analogously, for $d\in\left\{ 0,1\right\} $, $F_{Y\left( d\right) }\left( y\right) $ is identified by \[ F_{Y\left( d\right) }\left( y\right) =\mathbb{E}\left[ \frac{1\left\{ D=d\right\} }{dp\left( \mathbf{X}\right) +(1-d)\left( 1-p\left( \mathbf{X}\right) \right) }1\left\{ Y\leq y\right\} \right] , \] with $1\left\{ \cdot\right\} $ the indicator function, implying that both DTE$\left( y\right) $ and QTE$\left( \tau\right) $ can also be written as functionals of the observed data; see, e.g., Firpo2007, and Chen2008.

These identification results suggest that, if the PS were known, one could get consistent estimators by using the sample analogue of such estimands. For instance, one can estimate the ATE using the Hajek1971-type estimator \[ \widetilde{ATE}_{n}=\mathbb{E}_{n}\left[ \left( \varpi_{n,1}^{ps}\left( D,\mathbf{X}\right) -\varpi_{n,0}^{ps}\left( D,\mathbf{X}\right) \right) Y\right] , \] where \[ \varpi_{n,1}^{ps}\left( D,\mathbf{X}\right) =\left. \frac{D}{p\left( \mathbf{X}\right) }\right/ \mathbb{E}_{n}\left[ \frac{D}{p\left( \mathbf{X}\right) }\right] ,\text{ and }\varpi_{n,0}^{ps}\left( D,\mathbf{X}\right) =\left. \frac{1-D}{1-p\left( \mathbf{X}\right) }\right/ \mathbb{E}_{n}\left[ \frac{1-D}{1-p\left( \mathbf{X}\right) }\right] . \] Estimators for $F_{Y\left( d\right) }\left( y\right) $, $d\in\left\{ 0,1\right\} $, and DTE$\left( y\right) $ are formed using an analogous strategy. For the QTE$\left( \tau\right) ,$ one can simply invert the estimator of $F_{Y\left( d\right) }\left( y\right) $ to estimate $q_{Y\left( d\right) }\left( \tau\right) $; see, e.g., Firpo2007 and Chen2008. Of course, estimators for other treatment effect measures such as the difference of Theil indexes and/or Gini coefficients can also be formed using a similar strategy, see, e.g., Firpo2016.

In observational studies, however, the propensity score $p\left( \mathbf{X}\right) $ is usually unknown, and has to be estimated. Given that $\mathbf{X}$ is usually of moderate or high dimensionality, researchers routinely adopt a parametric approach. A popular choice among practitioners is to use the logistic model, where \[ p\left( \mathbf{X}\right) =p\left( \mathbf{X};\boldsymbol{\beta} _{0}\right) =\frac{\exp{(\mathbf{X}^{\prime}}\boldsymbol{\beta}_{0}{)} }{1+\exp{(\mathbf{X}^{\prime}}\boldsymbol{\beta}_{0}{)}}, \] with $\boldsymbol{\beta}_{0}\in\Theta\subset\mathbb{R}^{k}.$ Next, one usually proceeds to estimate $\boldsymbol{\beta}_{0}$ within the maximum likelihood paradigm, i.e., \[ \boldsymbol{\widehat{\beta}}_{n}^{mle}=\arg\max_{\boldsymbol{\beta}\in\Theta }\mathbb{E}_{n}\left[ D\ln\left( p\left( \mathbf{X};\boldsymbol{\beta }\right) \right) +\left( 1-D\right) \ln\left( 1-p\left( \mathbf{X} ;\boldsymbol{\beta}\right) \right) \right] , \] and uses the resulting PS fitted values $p\left( \mathbf{X} ;\boldsymbol{\widehat{\beta}}_{n}^{mle}\right) $ to construct different treatment effect estimators. Despite the popularity of this procedure, it has been shown that it can lead to significant instabilities under mild PS misspecifications, particularly when some PS estimates are relatively close to zero or one, see e.g. Kang2007.

In light of these challenges, alternative methods to estimate the PS have emerged. A particularly fruitful direction is to exploit the covariate balancing property of the PS, that is, to exploit the fact that, for all measurable and integrable function $f\left( \mathbf{X}\right) $ of the covariates $\mathbf{X}$,

for a unique value $\boldsymbol{\beta}_{0}\in\Theta$. For example, Imai2014 propose estimating the PS parameters $\boldsymbol{\beta}_{0}$ within the generalized method of moments framework where, for a finite vector of user-chosen functions $f\left( \mathbf{X}\right) $ (e.g. $f\left( \mathbf{X}\right) =\mathbf{X}$),

Graham2012, on the other hand, propose estimating $\boldsymbol{\beta }_{0}$ as the solution to a globally concave programming problem such that \[ \mathbb{E}\left[ \left( \frac{D}{p\left( \mathbf{X};\boldsymbol{\beta} _{0}\right) }-1\right) \mathbf{X}\right] =\mathbf{0}. \] Note that both procedures rely on choosing a finite number of functions $f\left( \mathbf{X}\right) $, though there is little to no theoretical guidance on how to choose such functions.

While estimators that balance low-order moments of covariates usually enjoy more attractive finite sample properties than those based on the ML paradigm, it is important to emphasize that the aforementioned proposals do not fully exploit the covariate balancing property characterized in ((ref)). Furthermore, as emphasized by Dominguez2004, the global identification condition for $\boldsymbol{\beta}_{0}$ can fail when one adopts the generalized method of moment approach, and only attempts to balance finitely many covariate moments.

In this paper we aim to estimate the PS parameters $\boldsymbol{\beta}_{0}$ by taking advantage of all the information contained in ((ref)). Our proposed estimators do not rely on tuning parameters such as bandwidth, do not consult the outcome data, and can be implemented in a data-driven manner. Our estimation procedure also guarantees that the unknown PS parameters are globally identified.

The integrated propensity score

In this section, we discuss how we operationalize our proposal. The crucial step is to reexpress the infinite number of covariate balancing conditions ((ref)) in terms of a more tractable set of moment restrictions, and then characterize $\boldsymbol{\beta}_{0}$ as the unique minimizer of a (population) minimum distance function. We then leverage on this characterization, and make use of the analogy principle to suggest a natural estimator for $\boldsymbol{\beta}_{0}$. In what follows, we present a step-by-step description of how we achieve this.

First, note that by using the definition of conditional expectation, ((ref)) can be expressed as

where $\mathbf{h}\left( D,\mathbf{X};\boldsymbol{\beta}\right) =\left( h_{1}\left( D,\mathbf{X};\boldsymbol{\beta}\right) ,h_{0}\left( D,\mathbf{X};\boldsymbol{\beta}\right) \right) ^{\prime}$, $h_{d}\left( D,\mathbf{X};\boldsymbol{\beta}\right) =\varpi_{d}^{ps}\left( D,\mathbf{X} ;\boldsymbol{\beta}\right) -1$, $d\in\left\{ 0,1\right\} $, and \[ \varpi_{1}^{ps}\left( D,\mathbf{X};\boldsymbol{\beta}\right) =\left. \dfrac{D}{p\left( \mathbf{X};\boldsymbol{\beta}\right) }\right/ \mathbb{E}\left[ \dfrac{D}{p\left( \mathbf{X};\boldsymbol{\beta}\right) }\right] ,\text{ }\varpi_{0}^{ps}\left( D,\mathbf{X};\boldsymbol{\beta }\right) =\left. \dfrac{1-D}{1-p\left( \mathbf{X};\boldsymbol{\beta }\right) }\right/ \mathbb{E}\left[ \dfrac{1-D}{1-p\left( \mathbf{X} ;\boldsymbol{\beta}\right) }\right] . \] That is, one can express the covariate balancing conditions ((ref)) in terms of stabilized conditional moment restrictions.

Next, by exploiting the \textquotedblleft integrated conditional moment approach\textquotedblright\ commonly adopted in the specification testing literature (Gonzalez-Manteiga2013 contains a comprehensive review), one can express ((ref)) as an infinite number of unconditional covariate balancing restrictions. That is, by appropriately choosing a parametric family of functions $\mathcal{W=}\left\{ w(\mathbf{X} ;\mathbf{u}):\mathbf{u}\in\Pi\right\} $, one can equivalently characterize ((ref)) as

see, e.g., Lemma 1 of Escanciano2006a for primitive conditions on the family $\mathcal{W}$ such that the equivalence between ((ref)) and ((ref)) holds. Choices of weight $w$ satisfying this equivalence include $(a)$ $w(\mathbf{X};\mathbf{u})=1\left\{ \mathbf{X}\leq \mathbf{u}\right\} $, where $\mathbf{u}\in\left[ -\infty,\infty\right] ^{k}$, $1\left\{ A\right\} $ denotes the indicator function of the event $A$ and $\mathbf{X}\leq\mathbf{u}$ is understood coordinate-wise (see, e.g., Stute1997 and Dominguez2004, Dominguez2015), $(b)$ $w(\mathbf{X};\mathbf{u})=\exp(i\mathbf{u}^{\prime}\Phi\left( \mathbf{X} \right) \mathbf{)}$, where $\mathbf{u}$ $\in$ $\mathbb{R}^{k}$, $\Phi\left( \cdot\right) $ is a vector of bounded one-to-one maps from $\mathbb{R}^{k}$ to $\mathbb{R}^{k}$ and $i=\sqrt{-1}$ is the imaginary unit (see, e.g., Bierens1982 and Escanciano2018), and $\left( c\right) $ $w(\mathbf{X};\mathbf{u})=1\left\{ \boldsymbol{\gamma}^{\prime}\mathbf{X}\leq u\right\} $,$~$where $\mathbf{u=}\left( \boldsymbol{\gamma},u\right) \in\mathbb{S}_{k}\times\left[ -\infty,\infty\right] $, $\mathbb{S} _{k}=\left\{ \boldsymbol{\gamma}\in\mathbb{R}^{k}:\left\Vert \boldsymbol{\gamma}\right\Vert =1\right\} $, and $\left\Vert \boldsymbol{\gamma}\right\Vert $ is the Euclidean norm of real-valued vector $\boldsymbol{\gamma}$ (see, e.g., Escanciano2006b). We call ((ref)) the \textquotedblleft integrated covariate balancing condition\textquotedblright\ because it uses the integrated (cumulative) measure of covariate balancing.

Finally, let

where $\mathbf{H}_{w}(\boldsymbol{\beta},\mathbf{u})=\mathbb{E}\left[ \mathbf{h}\left( D,\mathbf{X};\boldsymbol{\beta}\right) w(\mathbf{X} ;\mathbf{u})\right] $, $\left\Vert A\right\Vert ^{2}=A^{c}A$, $A^{c}$ denotes the conjugate transpose of the column vector $A$, and $\Psi(\mathbf{u})$ is an integrating probability measure that is absolutely continuous with respect to a dominating measure on $\Pi$.

With these results in hand, in the following lemma we show that

and $\boldsymbol{\beta}_{0}$ is the unique value such that the covariate balancing condition ((ref)) is satisfied.

Lemma (ref) is a global identification result that characterizes $\boldsymbol{\beta}_{0}$ as the unique minimizer of a population minimum distance function, $Q_{w}(\boldsymbol{\beta})$. That is, from Lemma (ref) we have that $\boldsymbol{\beta}_{0}$ is the unique PS parameter that minimizes the imbalances of all measurable and integrable functions $f\left( \mathbf{X}\right) $ between the treated, untreated and the combined group. Here, it is worth mentioning that neither Graham2012 nor Imai2014 covariate balancing approach guarantee global identification of the propensity score parameters. Instead, they directly assume that the vector of user-selected balancing conditions uniquely identify the propensity score parameters; see, e.g., Assumption 2.1 (i) of Graham2012. In practice, however, it is hard if not impossible to verify if such condition indeed holds. In cases it does not hold, inference procedures that rely on their proposed propensity score estimator, in general, will not be valid; see, e.g., Dominguez2004. Lemma (ref) shows that our propose IPS procedure completely avoids this important drawback.

Another important implication of Lemma (ref) is that it suggests a natural estimator for $\boldsymbol{\beta}_{0}$ based on the sample analogue of ((ref)), namely,

where $Q_{n,w}(\boldsymbol{\beta})=\int_{\Pi}\left\Vert \mathbf{H} _{n,w}(\boldsymbol{\beta},\mathbf{u})\right\Vert ^{2}\,\Psi_{n}(d\mathbf{u})$, $\Psi_{n}$ is a uniformly consistent estimator of $\Psi$, $\mathbf{H} _{n,w}(\boldsymbol{\beta},\mathbf{u})=\mathbb{E}_{n}\left[ \mathbf{h} _{n}\left( D,\mathbf{X};\boldsymbol{\beta}\right) w(\mathbf{X} ;\mathbf{u})\right] $, with $\mathbf{h}_{n}\left( D,\mathbf{X} ;\boldsymbol{\beta}\right) =\left( h_{n,1}\left( D,\mathbf{X} ;\boldsymbol{\beta}\right) ,h_{n,0}\left( D,\mathbf{X};\right. \right. $ $\left. \left. \boldsymbol{\beta}\right) \right) ^{\prime}$, $h_{n,d}\left( D,\mathbf{X};\boldsymbol{\beta}\right) =\varpi_{n,d} ^{ps}\left( D,\mathbf{X};\boldsymbol{\beta}\right) -1$, $d\in\left\{ 0,1\right\} $, and

We call $\widehat{\boldsymbol{\beta}}_{n,w}^{ips}$ the integrated propensity score estimator of $\boldsymbol{\beta}_{0}$ because it is based on the integrated covariate balancing conditions ((ref)).

From ((ref)), one can conclude that different PS estimators that fully exploit the covariate balancing property ((ref)) can be constructed by choosing different $w$ and $\Psi_{n}$. In this article, we focus on three different combinations that are intuitive, computationally simple, and that perform well in practice:

The estimators ((ref))-((ref)) build on Dominguez2004 and Escanciano2006b, Escanciano2018, respectively. Despite the apparent differences, they all aim to minimize covariate distribution imbalances: ((ref)) aims to directly minimize imbalances of the joint distribution of covariates; ((ref)) exploits the Cram\'{e}r-Wold theorem and focuses on minimizing imbalances of the distribution of all one-dimensional projections of covariates; and ((ref)) focuses on minimizing imbalances of the (transformed) covariates' joint characteristic function. From the Cram\'{e}r-Wold theorem and the fact that the characteristic function completely defines the distribution function (and vice-versa), ((ref) )-((ref)) are indeed intrinsically related. Furthermore, we emphasize that our estimators are data-driven, and neither $w$ nor $\Psi_{n}$ plays the role of a bandwidth as they do not affect the convergence rate of the IPS estimator.

From the computational perspective, ((ref))-((ref)) are easy to estimate because they do not involve matrix inversion nor nonparametric estimation. In the supplemental Appendix (ref), we show that the objective functions in ((ref))-((ref)) can be written in closed form, which, in turn, implies a more straightforward implementation. In practice, the IPS is easy to use as it is already implemented in the new package IPS for R, available at \url{https://github.com/pedrohcgs/IPS}.

Large sample properties

In this section, we first derive the asymptotic properties of the IPS estimators, namely the consistency, asymptotic linear representation, and asymptotic normality of $\widehat{\boldsymbol{\beta}}_{n,w}^{ips}$ . We then discuss how one can build on these results to conduct asymptotically valid inference for overall average, distributional and quantile treatment effects, using inverse probability weighted estimators. Although our proposal can also be used to estimate other treatment effects of interest such as those discussed in Firpo2016, we omit such a discussion for the sake of brevity.

Asymptotic theory for IPS estimator

Here we derive the asymptotic properties of the IPS estimator. Let the score of $\mathbf{H}_{w}(\boldsymbol{\beta},\mathbf{u})$ be defined as $\mathbf{\dot{H}}_{w}(\boldsymbol{\beta},\mathbf{u})=\left( \mathbf{\dot{H} }_{1,w}^{^{\prime}}(\boldsymbol{\beta},\mathbf{u}),\mathbf{\dot{H}} _{0,w}^{^{\prime}}(\boldsymbol{\beta},\mathbf{u})\right) ^{\prime},$ a $2\times k$ matrix, where, for $d\in\left\{ 0,1\right\} ,$ $\mathbf{\dot{H} }_{d,w}(\boldsymbol{\beta},\mathbf{u})=\mathbb{E}\left[ \mathbf{\dot{h}} _{d}\left( D,\mathbf{X};\boldsymbol{\beta}\right) w(\mathbf{X} ;\mathbf{u})\right] $, with $\mathbf{\dot{h}}_{1}$ and $\mathbf{\dot{h}}_{0}$ being the $1\times k$ vectors defined as

and $\dot{p}\left( \mathbf{\cdot;}\boldsymbol{\beta}\right) =\left. \left. \partial p\left( \mathbf{\cdot;}\boldsymbol{b}\right) \right/ \partial\boldsymbol{b}\right\vert _{\boldsymbol{b=\beta}},$ the $k\times1$ vector of scores of the PS model $p\left( \cdot,\boldsymbol{\beta}\right) $. We make the following set of assumptions.

Assumption (ref) is standard in the literature, see, e.g., Theorems 2.6 and 3.4 of Newey1994c, Example 5.40 of VanderVaart1998, and Graham2012. Assumption (ref)$(i)$ states that the true PS is known up to finite dimensional parameters $\boldsymbol{\beta}_{0}$, that is, we are in a parametric setup. Assumption (ref)$\left( ii\right) $ imposes that the parametric PS is bounded from above and from below. This assumption can be relaxed by assuming that $\left( D/p\left( \mathbf{X} ;\boldsymbol{\beta}\right) ,\left( 1-D\right) /\left( 1-p\left( \mathbf{X};\boldsymbol{\beta}\right) \right) \right) ^{\prime} \leq\mathbf{b}\left( \mathbf{X}\right) $ such that $\mathbb{E}\left[ \left\Vert \mathbf{b}\left( \mathbf{X}\right) \right\Vert ^{2}\right] <\infty$. Assumptions (ref)$(iii)$-$\left( iv\right) $ impose additional smoothness conditions on the PS, whereas Assumption (ref)$(v)$ (together with Assumption (ref)) implies that, in a small neighborhood of $\boldsymbol{\beta}_{0}$ and for all $u\in\Pi$, the score $\mathbf{\dot{H}}_{w}(\boldsymbol{\beta},\mathbf{u})$ is uniformly bounded by an integrable function.

Assumption (ref) restricts our attention to the IPS estimators ((ref))-((ref)). As mentioned before, we focus on such estimators because of their computational simplicity and transparency. Nonetheless, other types of IPS estimators can also be formed, provided that the weighting function $w$ and integrating measure $\Psi_{n}$ satisfy some high-level regularity conditions.

The next theorem characterizes the asymptotic properties of the IPS estimators $\widehat{\boldsymbol{\beta}}_{n,w}^{ips}$. Define the $k\times k$ matrix \[ C_{w,\Psi}=\int_{\Pi}\left( \mathbf{\dot{H}}_{w}(\boldsymbol{\beta} _{0},\mathbf{u})^{c}~\mathbf{\dot{H}}_{w}(\boldsymbol{\beta}_{0} ,\mathbf{u})+\mathbf{\dot{H}}_{w}(\boldsymbol{\beta}_{0},\mathbf{u})^{\prime }\left( \mathbf{\dot{H}}_{w}(\boldsymbol{\beta}_{0},\mathbf{u})^{\prime }\right) ^{c}\right) \Psi(d\mathbf{u}), \] and the $k\times1$ vector

From Theorem (ref), we conclude that the proposed IPS estimator is consistent, admits an asymptotic linear representation with influence function $l_{w,\Psi}\left( D,\mathbf{X};\boldsymbol{\beta}_{0}\right) $, and converges to a normal distribution. The asymptotic linear representation ((ref)) plays a major role in establishing the asymptotic properties of causal parameters such as average, distributional, and quantile treatment effects; see Section (ref).

Estimating treatment effects under unconfoundedness

In this section, we illustrate how one can estimate and make asymptotically valid inference about average, distributional, and quantile treatment effects under the unconfoundedness assumption (ref) using IPW estimators based on the IPS estimator $\widehat{\boldsymbol{\beta}}_{n,w}^{ips}.$

Based on the discussion in Section (ref), the IPW estimators for ATE, DTE and QTE are respectively:

where, for $d\in\left\{ 0,1\right\} $, \[ \widehat{q}_{n,Y\left( d\right) }^{ips}=\arg\min_{q\in\mathbb{R}} \mathbb{E}_{n}\left[ \varpi_{n,d}^{ps}\left( D,\mathbf{X} ;\widehat{\boldsymbol{\beta}}_{n,w}^{ips}\right) \cdot\rho_{\tau}\left( Y-q\right) \right] , \] with $\rho_{\tau}\left( a\right) =a\cdot\left( \tau-1\left\{ a\leq0\right\} \right) $ the check function as in Koenker1978, and the weights $\varpi_{n,1}^{ps}$ and $\varpi_{n,0}^{ps}$ are as in ((ref))-((ref)).

To derive the asymptotic properties of ((ref))-((ref)), we need to make an additional assumption about the underlying distributions of the potential outcomes $Y\left( 1\right) $ and $Y\left( 0\right) $.

Assumption (ref)$(i)$ requires potential outcomes to be square-integrable, whereas Assumption (ref)$(ii)$ is a mild regularity condition which guarantees that, in a small neighborhood of $\boldsymbol{\beta}_{0}$, the score of the IPW estimator for the ATE is bounded by an integrable function. Assumption (ref)$(iii)$ requires potential outcomes to be continuously distributed and only plays a role in the analysis of quantile treatment effects. In principle, Assumption (ref)$(iii)$ can be relaxed at the cost of using more complex arguments, see Chernozhukov2017c for details.

Before stating the results as a theorem, let us define some important quantities. Let

where, for $j\in\left\{ ate,dte,qte\right\} $, $g^{j}\left( Y,D,\mathbf{X} \right) =g_{1}^{j}\left( Y,D,\mathbf{X}\right) -g_{0}^{j}\left( Y,D,\mathbf{X}\right) $, with

and

The functions $g^{ate}$, $g^{dte}$ and $g^{qte}$ would be the influence functions of the ATE, DTE and QTE estimators, respectively, if the PS parameters $\boldsymbol{\beta}_{0}$ were known. With some abuse of notation, denote $\Omega_{w,\Psi}^{ate}=\mathbb{E}\left[ \psi_{w,\Psi}^{ate}\left( Y,D,\mathbf{X}\right) ^{2}\right] $, $\Omega_{w,\Psi,y}^{dte}=\mathbb{E} \left[ \psi_{w,\Psi}^{dte}\left( Y,D,\mathbf{X};y\right) ^{2}\right] $, and $\Omega_{w,\Psi,\tau}^{qte}=\mathbb{E}\left[ \psi_{w,\Psi}^{qte}\left( Y,D,\mathbf{X};\tau\right) ^{2}\right] $.

Theorem (ref) indicates that one can use our proposed IPS estimator to estimate a variety of causal parameters that are able to highlight treatment effect heterogeneity\footnote{Although the results stated in Theorem (ref) for distribution and quantile treatment effects are pointwise, in Appendix (ref) we prove their uniform counterpart using empirical process techniques. We omit the details in the main text only to avoid additional cumbersome notation. We refer interested readers to the proof of Theorem (ref) in Appendix (ref) for additional details.}. Furthermore, Theorem (ref) also suggests that to conduct asymptotically valid inference for these causal parameters, one simply needs to estimate the asymptotic variance $\Omega_{w,\Psi}^{ate}$, $\Omega _{w,\Psi,y}^{dte}$, and $\Omega_{w,\Psi,\tau}^{qte}$. Under additional smoothness conditions (for instance, the PS being twice continuously differentiable with bounded second derivatives), one can show that their sample analogues are consistent using standard arguments. We omit the details for the sake of brevity.

The IPS when treatment is endogenous

In many important applications, the assumption that treatment adoption is exogenous may be too restrictive. For instance, when individuals do not comply with their treatment assignment, or more generally when they sort into treatment based on expected gains, Assumption (ref) is likely to be violated. Imbens1994 and Angrist1996 point out that when this is the case and a binary instrument ($Z)$ for the selection into treatment is available, one can only nonparametrically identify treatment effect measures for the subpopulation of compliers, that is, individuals who comply with their actual assignment of treatment, and would have complied with the alternative assignment. As shown by Abadie2003, Frolich2007, and Frolich2013, the instrument propensity score $q\left( \mathbf{X} \right) \equiv$ $\mathbb{P}(Z=1|X)$ plays a prominent role in this local treatment effect (LTE) setup. In this section, we show that one can use the IPS approach to estimate the instrument propensity score $q\left( \mathbf{X}\right) $, by maximizing covariate distribution balancing among different instrument-by-treatment subgroups.

Before providing the details about how we apply the IPS approach to estimate $q\left( \mathbf{X}\right) $ under the LTE setup, we introduce a brief description of the LTE setup. Let $Z$ be a binary instrumental variable $Z$ for the treatment assignment. Denote $D\left( 0\right) $ and $D\left( 1\right) $ the value that $D$ would have taken if $Z$ is equal to zero or one, respectively. The realized treatment is $D=ZD\left( 1\right) +\left( 1-Z\right) D\left( 0\right) $. Thus, the observed sample in the LTE setup consists of independent and identically distributed copies $\left\{ \left( Y_{i},D_{i},Z_{i},\mathbf{X}_{i}^{^{\prime}}\right) ^{\prime}\right\} _{i=1}^{n}$. To identify the average, distributional and quantile treatment effects for the compliers, we follow Abadie2003 and make the following assumption.

Assumption (ref)($i$) imposes that once we condition on $\mathbf{X}$, $Z$ is \textquotedblleft as good as randomly assigned\textquotedblright. Assumption (ref)($ii$) imposes a common support condition, and guarantees that, conditional on $\mathbf{X}$, $Z$ is a relevant instrument for $D$. Finally, Assumption (ref)($iii$) is a monotonicity condition that rules out the existence of defiers.

From Abadie2003 and Frolich2013, we have that under Assumption (ref), the average, distributional and quantile treatment effects for compliers are nonparametrically identified, i.e.,

where $\mathcal{C}$ denotes the complier subpopulation, and, for $d\in\left\{ 0,1\right\} $,

and \[ \kappa_{d}\left( q\right) \equiv\mathbb{E}\left[ \frac{1\left\{ D=d\right\} Z}{q\left( \mathbf{X}\right) }-\frac{1\left\{ D=d\right\} \left( 1-Z\right) }{1\mathbb{-}q\left( \mathbf{X}\right) }\right] , \] and $F_{\varpi_{d}^{lte}\cdot Y}^{-1}\left( \tau\right) =\inf\left\{ y:F_{\varpi_{d}^{lte}\cdot Y}\left( y\right) \geq\tau\right\} .$ From the above results, it is clear that the instrument PS plays a prominent role in the LTE setup, and that once we have an estimator for $q$ available, it is relatively straightforward to construct estimators for the LATE, LDTE, and LQTE.

To estimate the instrument PS $q$, we adopt a parametric approach, i.e., we assume that $q\left( \mathbf{X}\right) =q\left( \mathbf{X;} \boldsymbol{\beta}_{0}^{lte}\right) $, where $q$ is known up to the finite-dimensional parameters $\boldsymbol{\beta}_{0}^{lte}$. Here, as we are interested in treatment effects for the (latent) subpopulation of compliers, we will attempt to estimate $\boldsymbol{\beta}_{0}^{lte}$ by maximizing the covariate distribution balance among compliers. To do so, we build on Theorem 3.1 of Abadie2003, which establishes that, for every measurable and integrable function $f\left( \mathbf{X}\right) $ of the covariates $\mathbf{X}$,

where $\varpi_{d}^{lte}\left( D,Z,\mathbf{X;}\boldsymbol{\beta}_{0} ^{lte}\right) $ is defined as in ((ref)) but with $\boldsymbol{\beta} _{0}^{lte}$ playing the role of $q$, as we assume that $q$ is a parametric model, and \[ \varpi^{lte}\left( D,Z,\mathbf{X;}\boldsymbol{\beta}\right) =\frac{1} {\kappa\left( \boldsymbol{\beta}\right) }\left( 1-\frac{\left( 1-D\right) Z}{q\left( \mathbf{X;}\boldsymbol{\beta}\right) }-\frac{D\left( 1-Z\right) }{1-q\left( \mathbf{X;}\boldsymbol{\beta}\right) }\right) , \] with \[ \kappa\left( \boldsymbol{\beta}\right) \equiv\mathbb{E}\left[ 1-\frac{\left( 1-D\right) Z}{q\left( \mathbf{X;}\boldsymbol{\beta}\right) }-\frac{D\left( 1-Z\right) }{1-q\left( \mathbf{X;}\boldsymbol{\beta }\right) }\right] . \] As noted in Theorem 3.1 of Abadie2003, under Assumption (ref), $\mathbb{E}\left[ \varpi^{lte}\left( D,Z,\mathbf{X;}\boldsymbol{\beta} _{0}^{lte}\right) \cdot f\left( \mathbf{X}\right) \right] =\mathbb{E} \left[ f\left( \mathbf{X}\right) |\mathcal{C}\right] $, implying that ((ref)) are indeed balancing conditions for the complier subpopulation.

Next and analogously to the discussion in Section (ref), we rewrite ((ref)) as

where $\mathbf{H}_{w}^{lte}(\boldsymbol{\beta},\mathbf{u})=\mathbb{E}\left[ \mathbf{h}^{lte}\left( D,Z,\mathbf{X};\boldsymbol{\beta}\right) w(\mathbf{X};\mathbf{u})\right] $, with $\mathbf{h}^{lte}\left( D,Z,\mathbf{X};\boldsymbol{\beta}\right) =(h_{1}^{lte}\left( D,Z,\mathbf{X} ;\boldsymbol{\beta}\right) $,\allowbreak\ $h_{0}^{lte}\left( D,Z,\mathbf{X} ;\boldsymbol{\beta}\right) )^{\prime}$, and, for $d\in\left\{ 0,1\right\} ,$ $h_{d}^{lte}\left( D,Z,\mathbf{X};\boldsymbol{\beta}\right) =\varpi _{d}^{lte}\left( D,Z,\mathbf{X;}\boldsymbol{\beta}\right) -\varpi ^{lte}\left( D,Z,\mathbf{X;}\boldsymbol{\beta}\right) $.

Based on ((ref)), we then show in Lemma (ref) in the Supplemental Appendix that $\boldsymbol{\beta}_{0}^{lte}$ is be globally identified, i.e., $\boldsymbol{\beta}_{0}^{lte}$ is the unique minimizer of the population minimum distance criteria $Q_{w}^{lte}(\boldsymbol{\beta} )=\int_{\Pi}\left\Vert \mathbf{H}_{w}^{lte}(\boldsymbol{\beta},\mathbf{u} )\right\Vert \,\Psi(d\mathbf{u})$. Thus, like in the case where treatment is exogenous, we can fully exploit the balancing conditions ((ref)) and estimate $\boldsymbol{\beta}_{0}^{lte}$ by

where $\Psi_{n}$ is a uniformly consistent estimator of $\Psi$, $\mathbf{H} _{n,w}^{lte}(\boldsymbol{\beta},\mathbf{u})=\mathbb{E}\left[ \mathbf{h} _{n}^{lte}\left( D,Z,\mathbf{X};\boldsymbol{\beta}\right) w(\mathbf{X} ;\mathbf{u})\right] $, $\mathbf{h}_{n}^{lte}\left( D,Z,\mathbf{X} ;\boldsymbol{\beta}\right) =(h_{n,1}^{lte}\left( D,Z,\mathbf{X} ;\boldsymbol{\beta}\right) $, $h_{n,0}^{lte}\left( D,Z,\mathbf{X} ;\boldsymbol{\beta}\right) )^{\prime}$, $h_{n,d}^{lte}\left( D,Z,\mathbf{X} ;\boldsymbol{\beta}\right) =\varpi_{n,d}^{lte}\left( D,Z,\mathbf{X;} \boldsymbol{\beta}\right) -\varpi_{n}^{lte}\left( D,Z,\mathbf{X;} \boldsymbol{\beta}\right) $, and

As before, we focus our attention on the three weighting functions described in Assumption (ref). We call ((ref)) the local integrated propensity score (LIPS) estimator.

In what follows, we derive the asymptotic properties of the instrument IPS estimator $\widehat{\boldsymbol{\beta}}_{n,w}^{lips}$. Let the score of $\mathbf{H}_{w}^{lte}(\boldsymbol{\beta},\mathbf{u})$ be defined as $\mathbf{\dot{H}}_{w}^{lte}(\boldsymbol{\beta},\mathbf{u})=\left( \mathbf{\dot{H}}_{1,w}^{lte^{\prime}}(\boldsymbol{\beta},\mathbf{u} ),\mathbf{\dot{H}}_{0,w}^{lte^{\prime}}(\boldsymbol{\beta},\mathbf{u})\right) ^{\prime}$ where, for $d\in\left\{ 0,1\right\} ,$ $\mathbf{\dot{H}} _{d,w}^{lte}(\boldsymbol{\beta},\mathbf{u})=\mathbb{E}\left[ \mathbf{\dot{h} }_{d}^{lte}\left( D,Z,\mathbf{X};\boldsymbol{\beta}\right) w(\mathbf{X} ;\mathbf{u})\right] $, \[ \mathbf{\dot{h}}_{d}^{lte}\left( D,Z,\mathbf{X};\boldsymbol{\beta}\right) =\boldsymbol{\dot{\varpi}}_{d}^{lte}\left( D,Z,\mathbf{X;}\boldsymbol{\beta }\right) -\boldsymbol{\dot{\varpi}}^{lte}\left( D,Z,\mathbf{X;} \boldsymbol{\beta}\right) , \] with

and

and $\dot{q}\left( \mathbf{\cdot;}\boldsymbol{\beta}\right) =\left. \left. \partial q\left( \mathbf{\cdot;}\boldsymbol{b}\right) \right/ \partial\boldsymbol{b}\right\vert _{\boldsymbol{b=\beta}}$. We make the following set of assumptions, which are the analogue of Assumption (ref).

The next theorem characterizes the asymptotic properties of the instrument IPS estimators $\widehat{\boldsymbol{\beta}}_{n,w}^{lips}$. Define the $k\times k$ matrix \[ C_{w,\Psi}^{lte}=\int_{\Pi}\left( \mathbf{\dot{H}}_{w}^{lte} (\boldsymbol{\beta}_{0}^{lte},\mathbf{u})^{c}~\mathbf{\dot{H}}_{w} ^{lte}(\boldsymbol{\beta}_{0}^{lte},\mathbf{u})+\mathbf{\dot{H}}_{w} ^{lte}(\boldsymbol{\beta}_{0}^{lte},\mathbf{u})^{\prime}\left( \mathbf{\dot {H}}_{w}^{lte}(\boldsymbol{\beta}_{0}^{lte},\mathbf{u})^{\prime}\right) ^{c}\right) \Psi(d\mathbf{u}), \] and the $k\times1$ vector

With the results of Theorem (ref) at hand, we can estimate the LATE, LDTE, and LQTE by using the instrument IPS estimators:

where, for $d\in\left\{ 0,1\right\} $, $\widehat{F^{r}}_{n,\varpi_{d} ^{lte}\cdot Y}$ $\left( \cdot\right) $ denotes the rearrangement of $\widehat{F}_{n,~\varpi_{d}^{lte}\cdot Y}\left( \cdot\right) $, \[ \widehat{F}_{n,\varpi_{d}^{lte}\cdot Y}\left( \cdot\right) =\mathbb{E} _{n}\left[ \varpi_{n,d}^{lte}\left( D,Z,\mathbf{X} ;\widehat{\boldsymbol{\beta}}_{n,w}^{lips}\right) 1\left\{ Y\leq \cdot\right\} \right] \text{, } \] if $\widehat{F}_{n,\varpi_{d}^{lte}\cdot Y}$ is not monotone, see, e.g., Chernozhukov2010, and Wuthrich2019\footnote{Lack of monotonicity may appear in finite samples because the weights $w_{n,d}^{lte}$ can be negative. This poses problems for the inversion of the weighted cumulative distribution functions to obtain the quantile functions. On the other hand, under Assumption (ref), the population weights $w_{d}^{lte}$ are non-negative, implying that these potential problems disappear, asymptotically. As discussed in detail in Chernozhukov2010, we can bypass such challenges by monotonizing $\widehat{F}_{n,~w_{d}^{lte}\cdot Y}$ via rearrangements.}. Importantly, these rearrangements do not change the asymptotic properties of the estimators.

To derive the asymptotic properties of ((ref))-((ref)), we impose the following regularity conditions, which are the analogue of Assumption (ref).

Monte Carlo simulations

Unconfoundedness setup

In this section, we conduct a series of Monte Carlo experiments to study the finite sample properties of our proposed treatment effect estimators based on the IPS. We first compare the performance of different IPW estimators for the ATE and the QTE$\left( \tau\right) $, $\tau\in\left\{ 0.25,0.5,0.75\right\} $ when one estimates the PS using our proposed IPS estimators ((ref))-((ref)), the classical maximum likelihood (ML) approach, Imai2014's just-identified covariate balancing propensity score (CBPS) as in ((ref)) with $f\left( \mathbf{X}\right) =\mathbf{X}$, and Imai2014's overidentified CBPS ((ref)) with $f\left( \mathbf{X}\right) =\left( \mathbf{X}^{\prime},\dot{p}\left( \mathbf{X};\boldsymbol{\beta}\right) ^{\prime}\right) ^{\prime}$, i.e., on top of balancing the means, one also makes use of the likelihood score equation. In all cases, we consider a logistic PS model where all available covariates enter linearly. All treatment effect estimators use stabilized weights ((ref)) and ((ref)).

We consider sample size $n$ equal to $500$\footnote{Simulation results with $n=200$ and $n=1000$ lead to similar conclusions and are available on request.}. For each design, we conduct $1,000$ Monte Carlo simulations. We compare the various IPW estimators in terms of average bias, root mean square error (RMSE), relative mean square error (relMSE), empirical 95% coverage probability, the median length of a 95% confidence interval, and the asymptotic relative efficiency (ARE)\footnote{For any parameter $\eta$ of a distribution $F$, and for estimators $\widehat{\eta}_{1}$ and $\widehat{\eta }_{2}$ approximately $N\left( \eta,V_{1}/n\right) $ and $N\left( \eta ,V_{2}/n\right) $, respectively, the asymptotic relative efficiency of $\widehat{\eta}_{2}$ with respect to $\widehat{\eta}_{1}$ is given by $V_{1}/V_{2}$; see, e.g., Section 8.2 in VanderVaart1998. Thus, to compute the ARE for our estimators, we build on Theorem (ref) and replace the asymptotic variances with their sample analogues.}. For the relative measures of performance, relMSE and ARE, we treat estimators based on the overidentified CBPS as the benchmark. The confidence intervals are based on the normal approximation in Theorem (ref), with the asymptotic variances being estimated by their sample analogues. For the variance of QTE$\left( \tau\right) $ estimators, we estimate the potential outcome densities using the Gaussian kernel coupled with Silverman's rule-of-thumb bandwidth - these are the default choices of the density function in the stats package in R. We use the CBPS package in R to estimate both CBPS estimators. Finally, we emphasize that our measures of performance highlight not only the behavior of IPW point estimates but also the accuracy of their associated inference procedures.

Our simulation design is largely based on Kang2007. Let $\mathbf{X} =\left( X_{1},X_{2},\right. $ $\left. X_{3},X_{4}\right) ^{\prime}$ be distributed as $N\left( 0,I_{4}\right) $, and $I_{4}$ be the $4\times4$ identity matrix. The true PS is given by \[ p\left( \mathbf{X}\right) =\frac{\exp\left( -X_{1}+0.5X_{2}-0.25X_{3} -0.1X_{4}\right) }{1+\exp\left( -X_{1}+0.5X_{2}-0.25X_{3}-0.1X_{4}\right) }, \] and the treatment status $D$ is generated as $D=1\left\{ p\left( \mathbf{X}\right) >U\right\} $, where $U$ follows a uniform $\left( 0,1\right) $ distribution. The potential outcomes $Y\left( 1\right) $ and $Y\left( 0\right) $ are given by

where $m\left( \mathbf{X}\right) =27.4X_{1}+13.7X_{2}+13.7X_{3}+13.7X_{4}$, $\varepsilon\left( 1\right) $ and $\varepsilon\left( 0\right) $ are independent $N\left( 0,1\right) $ random variables. The ATE and the QTE$\left( \tau\right) $ are equal to 10, for all $\tau\in\left( 0,1\right) $.

We consider two different scenarios to assess the sensibility of the proposed estimators under misspecified models that are \textquotedblleft nearly correct\textquotedblright. In the first experiment, the observed data is $\left\{ \left( Y_{i},D_{i},\mathbf{X}_{i}^{^{\prime}}\right) ^{\prime }\right\} _{i=1}^{n}$, and, therefore, all IPW estimators are correctly specified. In the second experiment the observed data is $\left\{ \left( Y_{i},D_{i},\mathbf{W}_{i}^{^{\prime}}\right) ^{\prime}\right\} _{i=1}^{n}$, where $\mathbf{W}=\left( W_{1},W_{2},W_{3},W_{4}\right) ^{\prime}$ with $W_{1}=\exp\left( X_{1}/2\right) $, $W_{2}=X_{2}/\left( 1+\exp\left( X_{1}\right) \right) ,$ $W_{3}=\left( X_{1}X_{3}/25+0.6\right) ^{3}$, and $W_{4}=\left( X_{2}+X_{4}+20\right) ^{2}$. In this second scenario, the IPW estimators for ATE and QTE$\left( \tau\right) $ are misspecified.

Table (ref) displays the simulation results for both scenarios. When the PS model is correctly specified, all estimators perform well in terms of bias and coverage probability, i.e., all estimators are essentially unbiased and their associated confidence intervals have correct coverage. Comparing ML-based with CBPS-based estimators, we note that IPW estimators based on ML tend to have higher mean square error, longer confidence intervals, and lower ARE. Thus, it is clear that CBPS-based IPW estimators can improve upon those based on ML. However, our simulation results under correct specification suggest that we can improve further the performance of the CBPS estimator by fully exploiting the covariate balancing of the propensity score. For instance, the relative mean square error of estimators based on the IPS with either projection or exponential weight function tend to be at least 10% smaller than those based on the CBPS, with the exception of the QTE$\left( 0.25\right) $. The gains in terms of ARE also tend to be large. For example, the ARE of the ATE estimator based on the IPS with projection weight function with respect to the one based on the overidentified CBPS is 1.26. This implies that the ATE estimator based on the overidentified CBPS would require $1.26\times n$ observations to perform equivalently to the ATE estimator based on IPS with projection weight. IPS estimators based on the exponential weight also tend to dominate CBPS estimators in terms of mean square errors and ARE. Finally, we note that IPW estimators based on the IPS with the indicator function tend to give slightly larger confidence intervals than when using other IPS estimators, perhaps because there are multiple covariates (four in our simulation design), implying that many $1\left\{ \mathbf{X}_{i} \leq\mathbf{u}\right\} $ are equal to zero when $\mathbf{u}$ is evaluated at the sample observations.

When the PS model is misspecified, our Monte Carlo results suggest that the potential gains of using the IPS can also be pronounced. In this scenario, we note that estimators based on ML tend to be substantially biased, have relatively high RMSE, and inference tends to be misleading. These findings are in line with the results in Kang2007. Overall, estimators based on just-identified CBPS improve on ML, though under-coverage is still an unresolved issue when one focuses on the ATE and QTE$\left( 0.75\right) $. Estimators based on the overidentified CBPS tend to have better coverage than those based on the just-identified CBPS, but under-coverage of QTE$\left( 0.75\right) $ is still severe, perhaps because of the large biases. Finally, we note that our proposed IPS estimators tend to further improve upon CBPS. In particular, estimators based on the IPS with the projection weight function have the lowest bias and RMSE, and their confidence intervals are close to the nominal coverage --- the only exception is when one focuses on QTE$\left( 0.25\right) $, where estimators based on CBPS tends to perform slightly better than our proposed IPS procedure. On the other hand, we note that, in terms of mean square error, the gains of adopting the IPS estimator with either projection or exponential weighting function tend to be large in all other considered causal measures, especially for ATE and QTE$\left( 0.75\right) $.

Local Treatment Effect Setup

We now consider the setup where treatment is endogenous but one has access to a binary instrument $Z$, as described in Section (ref). Here, we compare the performance of different IPW estimators for the LATE and the LQTE$\left( \tau\right) $, $\tau\in\left\{ 0.25,0.5,0.75\right\} $ when one estimates the instrument PS $q\left( \cdot\right) $ using our proposed instrument IPS estimator ((ref)) with exponential, indicator and projection-based weights, the classical ML approach, Imai2014's just-identified and overidentified CBPS with $Z$ playing the role of $D$. In all cases, we consider a logistic instrument PS model where all available covariates enter linearly. As in the unconfoundedness case, we consider sample size $n$ equal to $500$, and conduct $1,000$ Monte Carlo simulations for each design.

The simulation design is similar to the one in Section (ref). Let $\mathbf{X}$, $\mathbf{W,}$ $Y\left( 1\right) $, and $Y\left( 0\right) $ be defined as before. The true instrument PS is given by \[ q\left( \mathbf{X}\right) =\frac{\exp\left( -X_{1}+0.5X_{2}-0.25X_{3} -0.1X_{4}\right) }{1+\exp\left( -X_{1}+0.5X_{2}-0.25X_{3}-0.1X_{4}\right) }, \] the instrument $Z$ is generated as $Z=1\left\{ q\left( \mathbf{X}\right) >U_{1}\right\} $, where $U_{1}$ follows a uniform $\left( 0,1\right) $ distribution. The potential treatments $D\left( 1\right) $ and $D\left( 0\right) $ are generated as $D\left( 1\right) =1\{p^{\ast}\left( Y\left( 1\right) -Y\left( 0\right) \right) >U_{2}\}$ and $D\left( 0\right) =0,$ where $U_{2}$ follows a uniform $\left( 0,1\right) $ distribution, and \[ p^{\ast}\left( Y\left( 1\right) -Y\left( 0\right) \right) =\frac {\exp\left( 2+0.05\cdot\left( Y\left( 1\right) -Y\left( 0\right) \right) \right) }{1+\exp\left( 2+0.05\cdot\left( Y\left( 1\right) -Y\left( 0\right) \right) \right) }. \] Finally, the realized treatment is $D=Z\cdot D\left( 1\right) +\left( 1-Z\right) \cdot D\left( 0\right) $, and the realized outcome is $Y=D\cdot Y\left( 1\right) +\left( 1-D\right) \cdot Y\left( 0\right) $. The LATE, LQTE$\left( 0.25\right) $, LQTE$\left( 0.5\right) $, and LQTE$\left( 0.75\right) $ are approximately equal to $39.25,$ 42.94, 35, and 42.94, respectively. This design is consistent with a generalized Roy model, under which individuals with higher treatment effects are more likely to be treated if they are eligible for treatment. We also emphasize that, given the one-sided non-compliance, LATE and LQTE are equal to the ATT and QTT in this scenario.

As before, we consider two scenarios. On the first one, the observed data is $\left\{ \left( Y_{i},D_{i},Z_{i},\mathbf{X}_{i}^{^{\prime}}\right) ^{\prime}\right\} _{i=1}^{n}$, and, therefore, all IPW estimators are correctly specified. In the second scenario, the observed data is $\left\{ \left( Y_{i},D_{i},Z_{i},\mathbf{W}_{i}^{^{\prime}}\right) ^{\prime }\right\} _{i=1}^{n}$, and all considered IPW estimators for LATE and LQTE$\left( \tau\right) $ are misspecified.

Table (ref) displays the simulation results for both scenarios. When the instrument PS model is correctly specified, all estimators perform well in terms of bias and coverage probability, except the estimators based on the LIPS estimator ((ref)) with the indicator weighting function --- the bias of the local treatment effect estimators based on LIPS with indicator function is non-negligible when $n=500$, and such biases distort the confidence intervals. In additional simulations, we note that the bias associated with estimators based on the LIPS with the indicator weighting function converges to zero when sample size grows, though the rate of convergence is rather slow. As such, we recommend that, in practice, one should favor the other PS estimators with respect to the LIPS with the indicator weighting function. Like in the unconfoundedness setup, we note that IPW estimators based on ML tend to have higher mean square error, longer confidence intervals, and lower ARE than the IPW estimators based on the just-identified CBPS estimator; the performance of the overidentified CBPS is, in general, worse than MLE, specially for LATE. The results in Table (ref) also show that, when the instrument propensity score is correctly specified, the LIPS estimators with the exponential or projection weighting function tend to outperform the other methods, particularly when estimating the LATE and LQTE$\left( 0.75\right) $.

When the instrument PS model is misspecified, our Monte Carlo results suggest that using the LIPS can also be attractive. In this setup, we note that estimators based on ML tend to have higher biases, RMSE and misleading confidence intervals. Local treatment effect estimators based on the (instrumented) CBPS improve upon those based on ML, with the just-identified\ CBPS estimator performing better than the overidentified CBPS. However, under-coverage is still an issue, except when one focuses on LQTE$\left( 0.25\right) $. On the other and, our simulation results suggest that our proposed LIPS estimators lead to local treatment effect estimators with even better statistical properties than those based on the (instrumented) CBPS --- such gains are specially pronounced when estimating the local treatment effect parameters based on the LIPS with the exponential or projection weighting functions.

Overall, our Monte Carlo simulations illustrate that, by fully exploiting the covariate balancing property of the (instrument) PS, we can get treatment effect estimators with improved finite sample properties. Our simulation results also point out that treatment effect estimators based on the IPS and LIPS estimators with either exponential or projection weighting functions tend to perform better than when one uses the indicator weighting function. As such, we recommend that, in practice, one should favor these weighting functions with respect to the indicator weighting function, especially when the dimension of the covariates included in the PS model is moderate or high\footnote{In unreported additional simulations, we also have found that the numerical performance of $IPS_{ind}$ and $LIPS_{ind}$ is sometimes sensitive to initial values used in the optimization procedure when the number of included covariates is moderate. We argue that this is additional reason to favor the other weighting functions with respect to the indicator one.}.

Empirical illustrations

In this section, we apply our proposed tools to two different datasets. First, we revisit Ichino2008 and use Italian data from the early 2000s to study if temporary work agency (TWA) assignment affects the probability of finding a stable job later on. Second, we study the effect of 401(k) retirement plan on asset accumulation using data from the Survey of Income and Program Participation, as in Benjamin2003a, Abadie2003, and Chernozhukov2004.

Effect of temporary work assignment on future stable employment

In temporary agency work, a company that needs employees signs a contract with a TWA, which, in turn, is in charge of hiring and subsequently leasing these workers to the company. In contrast to \textquotedblleft traditional\textquotedblright\ jobs, the TWA is in charge of paying the workers salary and fringe benefits, whereas the company's responsibility is to train and guide the workers. One of the main arguments of introducing temporary agency work is that it helps workers facing barriers to employment find a stable job later on.

To evaluate whether TWA assignment has a positive impact on employment, Ichino2008 collected data for two Italian regions, Tuscany and Sicily, in the early 2000s. The dataset contains 2030 individuals, 511 of them treated and 1519 untreated. Here, the treated group consists of individuals who were on a TWA assignment during the first 6 months of 2001, whereas the untreated group contains individuals aged 18 - 40, who belonged to the labor force but did not have a stable job on January 2001, and who did not have a TWA assignment during the first semester of 2001. Thus, both treatment groups were drawn from the same local labor market. The outcome of interest is having a permanent job at the end of 2002. A rich set of variables related to demographic characteristics, family background, educational achievements, and work experience before the treatment period were collected to adjust for potential confounding (see Table 1 in Ichino2008). Using PS matching, Ichino2008 find evidence that TWA assignment has a positive effect on permanent employment, especially in Tuscany. The results for Sicily are sensitive to small violations of the strong ignorability assumption. Therefore, in what follows, we focus on the Tuscany sub-sample\footnote{The data are publicly available at http://qed.econ.queensu.ca/jae/2008-v23.3/ichino-mealli-nannicini/.}.

We use the results in Sections (ref) to estimate the ATE. We compare different IPW estimators based on the same PS estimation methods as in the simulation studies in Section (ref), except the IPS coupled with the indicator weighting function as it tends to be numerically unstable when dimension of covariates is moderate. Table (ref) shows the point estimates and standard errors (in parentheses) for the whole Tuscany sample, and presents some heterogeneity results based on gender. The PS specification we use is the one adopted by Ichino2008, which includes all the pre-treatment variables mentioned in Table 1 of Ichino2008, squared distance, and an interaction between self-employment and one of the provinces.

The results in Table (ref) suggest that the ATE is positive, and statistically significant at the conventional levels, regardless of the estimation procedure adopted. The overall average effect of TWA assignment on the probability of having a permanent job ranges from 18 to 21, 14 to 23, and 15 to 19 percentage points when using the whole sample, the male subpopulation, and the female subpopulation, respectively. Interestingly, the IPS estimators can provide gains of efficiency when compared to both the MLE and CBPS estimators. For instance, for the subsample of females, the asymptotic relative efficiency (ARE) of the ATE estimator based on the IPS with exponential, and projection weights with respect to the one based on MLE are 1.70, and 1.58, respectively, while the ARE for the ATE based on the just and overidentified CBPS with respect to the one based on MLE are, respectively, 0.9 and 0.8. These findings suggest that the IPS can indeed lead to improved treatment effect estimators in relevant settings.

Effect of 401(k) retirement plans on asset accumulation

As discussed in Benjamin2003a, Abadie2003, Chernozhukov2004, and many others, tax-deferred retirement plans have been popular in the US since the 1980s. A main goal of these programs is to increase individual saving for retirement. Amongst the most popular tax-deferred programs is the 401(k) plan. Interestingly, 401(k) plans are provided by employers, and, therefore, only workers in firms that offer such programs are eligible. On the other hand, we emphasize that eligible employees choose whether to participate (i.e., make a contribution) or not, making the evaluation of the effectiveness of 401(k) plans on accumulated assets more challenging as a result of endogeneity concerns --- individuals who participate in 401(k) programs have stronger preferences for savings and would have saved more even in the absence of these programs.

To bypass the endogeneity challenge, Benjamin2003a uses data from the 1991 Survey of Income and Program Participation (SIPP) and compares households that are eligible with those who are non-eligible for 401(k) plans to assess the effect of eligibility on accumulated assets. He argues that since 401(k) eligibility is determined by the employers, household preference for savings plays a negligible role in determining eligibility once one controls for observed household characteristics. Using PS matching, Benjamin2003a finds evidence that 401(k) eligibility has a positive effect on asset accumulation.

Abadie2003, Chernozhukov2004 and Wuthrich2019, on the other hand, study the effect of 401(k) participation on asset accumulation, using 401(k) eligibility as an instrument for the actual participation status. Similarly to Benjamin2003a, they argue that 401(k) eligibility is exogenous after controlling for a vector of observed household characteristics. Abadie2003, using a semiparametric IPW estimator for the LATE, finds that the effect of 401(k) participation on net financial assets is significant and positive. Chernozhukov2004 and Wuthrich2019, using an IV quantile regression model, also find positive and significant effects of 401(k) participation on net financial assets.

In what follows, we apply the methodology discussed in Sections (ref) and (ref) to study the effects of eligibility and participation in 401(k) programs on saving behavior. As suggested by Benjamin2003a, Abadie2003, and Chernozhukov2004, eligibility is assumed to be exogenous after controlling for covariates. Also note that, because only eligible individuals can enroll in 401(k) plans, the monotonicity condition in Assumption (ref)$\left( iii\right) $ holds trivially, and the LATE and LQTE estimators presented in Section (ref) approximate the average and quantile treatment effect for the treated (i.e., for 401(k) participants).

We use the same dataset as Benjamin2003a, Chernozhukov2004 and Wuthrich2019. The data consists of a sample of 9,910 households from the 1991 SIPP\footnote{The original data have 9,915 households, but we follow Benjamin2003a and delete the five observations with zero or negative income. Descriptive statistics are available in Table 1 in Benjamin2003a and in Tables 1 and 2 in Chernozhukov2004.}. The outcomes of interest are net financial assets, and total wealth. For the (instrument) propensity score estimation, we adopt a logistic specification, and use all two-way interactions between income, log-income, age, family size, years of education, dummies for homeownership, marital status, two-earner status, defined benefit pension status, and individual retirement account participation status. To assess the reliability of this parametric PS model, we apply the specification test of SantAnna2018 with 1,000 bootstrap draws, and fail to reject the null of the propensity score model being correctly specified at the 10% level.

Panel A (Panel B) of Table (ref) shows the point estimates and standard errors (in parentheses) for the effect of 401(k) eligibility (participation among compliers) on net financial assets and total wealth. We present IPW estimators for the ATE, QTE$\left( 0.25\right) $, QTE$\left( 0.5\right) $ and QTE$\left( 0.75\right) $, and for the LATE, LQTE$\left( 0.25\right) $, LQTE$\left( 0.5\right) $ and LQTE$\left( 0.75\right) $ using the same PS estimation methods as in the simulation exercise in Section (ref), except the IPS and LIPS estimators based on the indicator weighting function, as they tend to be numerically unstable when the dimension of covariates is moderate.

The results in Panel A suggest that 401(k) eligibility has a positive and significant average impact on both net financial assets and total wealth and that the effect is more pronounced at the higher quantiles. When one compares the treatment effect measures across different PS estimation methods, we see that the results tend to be similar for net financial assets; for total wealth, we note that estimators based on the overidentified CBPS estimator suggest much larger effects of 401(k) eligibility at higher quantiles than those based on our proposed IPS estimators; see Figure (ref) for a more detailed comparison between the QTE estimates based on the IPS with projection weighting function, overidentified CBPS (the default in the CBPS R package), and those based on ML.

The results in Panel B paint a similar picture as those in Panel A: 401(k) participation tends to have a positive and significant average impact on both measures of wealth, and the effect is more pronounced at the right tail of the wealth measures. As we illustrated in Figure (ref), there are quantitative differences between the LQTE estimates based on different PS estimation methods, with those based on the overidentified instrument CBPS suggesting much larger effects than the other estimation methods, though the shape of the LQTE function is similar across specifications.

Conclusion

In this article, we proposed a framework to estimate propensity score parameters such that, instead of targeting to balance only some specific moments of covariates, it aims to balance all functions of covariates. The proposed estimator is of the minimum distance type, and is data-driven, $\sqrt{n}$-consistent, asymptotically normal, and admits an asymptotic linear representation that facilitates the study of inverse probability weighted estimators in a unified manner. Importantly, we have shown that our framework can accommodate the empirically relevant situation under which treatment allocation is endogenous. We derived the large sample properties of average, distributional and quantile treatment effect estimator based on the proposed integrated propensity scores, and illustrated its attractive properties via a Monte Carlo study and two empirical applications.

Although this paper devoted most of its attention to forming IPW-type treatment effect estimators, we note that sometimes researchers are willing to consider an outcome regression model, on top of the propensity score model. In such cases, we stress that one can easily combine our IPS estimation procedure with such outcome regression model to form doubly-robust, locally efficient treatment effect estimators, see, e.g., Sloczynski2018 and references therein. Perhaps even better, one can use the integrated moment approach adopted in this paper to estimate not only the propensity score, but also the outcome regression model. We leave the detailed discussion of such procedure for future research.

\onehalfspacing