Stable Probability Weighting: Large-Sample and Finite-Sample Estimation and Inference Methods for Heterogeneous Causal Effects of Multivalued Treatments Under Limited Overlap

Stable Probability Weighting Large-Sample and Finite-Sample Estimation and Inference Methods for Heterogeneous Causal Effects of Multivalued Treatments Under Limited Overlap

\onehalfspacing

Introduction

Inverse Probability Weighting (IPW), or its augmented doubly robust version, is now a widely used technique for estimating causal effects and for dealing with missing data. The technique is based on the influential work of horvitz1952generalization, robins1994estimation, robins1995semiparametric, hahn1998role, hirano2003efficient, and several others. However, even in its early days, the method was not without its critics. basu1971essay illustrates some practical issues with it using an amusing tale, in which a circus statistician gets sacked by the circus owner after proposing an unsound IPW estimator of the total mass of circus elephants (given by the mass of a randomly chosen elephant multiplied by the reciprocal of an associated extreme selection probability). Unfortunately, “Basu's elephants” are not confined to fictional circuses; they manifest themselves in many policy-relevant empirical settings, and extreme weights based on IPW can result in questionable and unstable estimates of causal effects.\footnote{See, e.g., kang2007demystifying, robins2007comment, tsiatis2007comment, kang2007rejoinder, frolich2004finite, khan2010irregular, ma2020robust, heiler2021valid, sasaki2022estimation, ma2022testing, d2021overlap, crump2006moving, li2019addressing, petersen2012diagnosing, and the references therein.} This problem with IPW shows up in many important datasets across various fields, including the health sciences\footnote{See, e.g., crump2009dealing, scharfstein1999adjusting, heiler2021valid.} and the social sciences.\footnote{See, e.g., ma2020robust, huber2013performance, busso2014new.} In fact, the pioneers of the modern form of (augmented) IPW have themselves repeatedly issued warnings regarding highly variable, extreme inverse probability weights.\footnote{See, e.g., robins1995semiparametric, scharfstein1999adjusting, robins2000inference, robins2007comment.}

In this paper, I provide new large-sample and finite-sample methods, which serve as alternatives to the (augmented) IPW framework, for estimating and inferring heterogeneous causal effects under unconfoundedness, i.e., the standard assumption of strongly ignorable treatment assignment rosenbaum1983central. Even though I primarily work with binary treatment variables to avoid cumbersome notation, I show that my framework readily extends to the setting of multivalued treatments with generalized propensity scores imbens2000role. While studying the existing literature on limited overlap (extreme selection on observables), I have come across only papers that deal with the average treatment effect (ATE) parameter in this context (to various extents). To the best of my knowledge, this paper is the first to propose estimation and inference methods for higher-dimensional parameters that characterize heterogeneity in the (possibly multivalued) treatment effects when limited overlap is empirically important and non-negligible.\footnote{As far as I am aware, the existing literature (from a frequentist or machine learning or policy learning perspective) on heterogeneous treatment effect estimation and inference relies on strong overlap. See wang2020debiased, zhao2019selective, zhao2022selective, athey2019generalized, powers2018some, nie2021quasi, kunzel2019metalearners, wager2018estimation, athey2016recursive, kennedy2022towards, kennedy2022minimax, semenova2021debiased, imai2013estimating,belloni2017program,singh2020kernel,chernozhukov2018generic,knaus2021machine, oprescu2019orthogonal,semenova2022estimation,crump2008nonparametric,sant2021nonparametric, athey2021policy, kitagawa2018should, kitagawa2022stochastic, sun2021empirical,abrevaya2015estimating,fan2022estimation,lee2017doubly. There is some discussion on limited overlap in the Bayesian literature on heterogeneous treatment effect estimation li2022bayesian. hill2013assessing propose some heuristic strategies based on Bayesian additive regression trees to deal with poor overlap. hahn2020bayesian propose finding the regions with limited overlap and then trimming them or to use the spline-based extrapolation procedures suggested by nethery2019estimating. The existing literature on causal estimation and inference in the case of multi-valued treatments also relies on the assumption of strong overlap. See, ai2021unified, ai2022estimation, bugni2019inference, cattaneo2010efficient, kennedy2017non,farrell2015robust,su2019non,colangelo2020double,yang2016propensity, zhang2022towards. A very recent exception is the work of dias2021inference, who allow the proportion of treated units to diminish asymptotically. They provide some interesting results in a specialized setting with discrete covariates and some shape restrictions, especially a “conditional rank invariance restriction [that] is analogous to the control function framework for quantile regression.” They also provide valid inferential methods for a narrow class of null hypotheses. Since dias2021inference allow for dynamic treatment effects, their methods are very useful in some panel data settings, but their model setup is specialized and different from the standard version of unconfoundedness that I use.}

In the setting of rosenbaum1983central, if the assumption of strong overlap (also known as strict overlap, i.e., propensity scores being bounded away from zero and one) does not hold, then inverse weight estimation is problematic; specifically, khan2010irregular show that IPW may not lead to an asymptotically normal regular $\sqrt{n}$-consistent estimator of even just the ATE, which is a simpler parameter than the heterogeneous treatment effect parameters. Depending on how the propensity scores are distributed, the IPW estimator may be non-Gaussian (e.g., asymmetric L\'evy stable) asymptotically ma2020robust. Although strong overlap is often invoked in the literature hirano2003efficient, chernozhukov2022automatic or is implicit in other assumptions d2021overlap, it rules out many simple, plausible models (e.g., some probit models) for the treatment assignment variable and also contradicts many empirically observed propensity score distribution shapes heiler2021valid, ma2022testing,lei2021distribution.

In the presence of limited overlap, some researchers suggest focusing attention instead on alternative estimands and parameters that can be estimated efficiently using heavy trimming or winsorization of the propensity scores crump2009dealing, zhou2020propensity or using other weighting procedures.\footnote{See li2018balancing, graham2012inverse, zubizarreta2015stable, wang2020minimal, athey2018approximate, hainmueller2012entropy, zhao2019covariate, imai2014covariate, imai2015robust, robins2000marginal, ai2021unified, yang2018asymptotic,wong2018kernel,ning2020robust,hirshberg2021augmented,ben2022using,wang2020debiased,matsouaka2022overlap,khan2021adaptive, chen2021robust,li2019propensity. Some methods, such as those of chen2008semiparametric and hirshberg2021augmented, make weaker assumptions than strong overlap but still restrict limited overlap. For example, hirshberg2021augmented require the first moment of the inverse propensity score to be bounded, but this testable assumption may lack empirical support ma2022testing.} lee2021bounding propose partial identification methods. ma2020robust develop bias-correction and robust inference procedures for trimmed IPW estimators.\footnote{chaudhuri2014heavy, yang2018asymptotic, sasaki2022estimation, khan2022uniform,khan2022heteroscedasticity also propose related trimming approaches for ATE estimation and inference under limited overlap.} heiler2021valid propose an empirically appealing method (a modified $m$-out-of-$n$ bootstrap procedure that does not rely on trimming) for robust inference on the (augmented) IPW estimators of the ATE in the presence of limited overlap. However, heiler2021valid restrict “the occurrence of extreme inverse probability weighted potential outcomes or conditional mean errors.” Although the focus of my paper is on the heterogeneity of causal effects, my proposed methods also offer an alternative to the aforementioned existing large-sample methods for inference on the ATE.

In Section (ref), I formalize the observational setup and the parameters of interest that are used throughout the paper. I do not impose the empirically restrictive assumption of strict overlap, but the setting is otherwise quite standard in the literature on causal inference under unconfoundedness. In Section (ref), I review issues with the popular (augmented) inverse probability weighting techniques under weak overlap. I illustrate these issues by considering a very simple setting where the propensity scores are known and a linear model describes heterogeneity in average causal effects. This particular setting is designed to help us focus on the core issues related to limited overlap, paving the way for a more thorough consideration of the problem from semi/non-parametric perspectives in later sections. Even in the simple setting of Section (ref), IPW approaches can have undesirable statistical properties when overlap is not strong. However, I show that there exist simple alternatives: “Non-Inverse Probability Weighting” (NPW) estimators, which are consistent and asymptotically normal while also being robust to limited overlap. I define a general estimator class that nests IPW and NPW estimators and show that they can be expressed as weighted versions of one another. I also discuss NPW-based estimation of the ATE.

My estimators in Section (ref) are based on a general principle called “Stable Probability Weighting” (SPW), which I develop more fully in the next section, for learning about the conditional average treatment effect (CATE) function that characterizes the heterogeneity in the average causal effects. In Section (ref), I present the relevant conditional moment restrictions that enable estimation and inference procedures that are robust to limited overlap and provide several examples. (If strict overlap holds, my general framework nests IPW as a special case.) I also show how to augment the basic conditional moment model in order to achieve double robustness (i.e., consistent estimation despite misspecification of one of the nuisance parameters) and compare my approach with some moment-based methods that are widely used in the existing literature. I then generalize the SPW approach for analyzing multivalued treatments and distributional or quantile treatment effects. I also discuss SPW-based estimation and inference from parametric, semiparametric, nonparametric, machine learning, and policy learning perspectives. Thanks to the currently rich econometric literature on conditional moment models, I do not have to reinvent the wheel for large-sample estimation and inference. Next, I proceed to develop finite-sample methods in Section (ref).

Even under strict overlap, the standard IPW estimators may have poor finite-sample behavior. But under limited overlap, the problem is much worse, and those estimators typically have very undesirable small-sample properties armstrong2021finite, rothe2017robust, hong2020inference. hong2020inference provide restrictions on limited overlap that allow the use of standard methods for valid asymptotic inference on the finite-population version of the ATE when propensity scores degenerate to zero asymptotically. In a setting with finite strata, rothe2017robust assumes normally distributed potential outcomes and proposes valid inference for the empirical ATE (EATE), which is the sample average of the expected treatment effect conditional on the covariates, by turning the problem into a general version of the Behrens--Fisher problem.\footnote{An alternative to the approach of rothe2017robust is that of dias2021inference, who also use discrete covariates. Rather than assuming normally distributed potential outcomes like rothe2017robust does, dias2021inference impose some shape restrictions, especially a “conditional rank invariance restriction [that] is analogous to the control function framework for quantile regression.” My finite-sample methods in this paper do not impose either set of assumptions.} armstrong2021finite also develop some finite-sample inference methods for the EATE under the assumption of normally distributed errors with known variances for potential outcomes whose conditional means lie within a known convex function class.

Both rothe2017robust and armstrong2021finite acknowledge the restrictiveness of the normality assumption for finite-sample inference.\footnote{rothe2017robust says, “We work with normality since without some restriction of this type, it would seem impossible to obtain meaningful theoretical statements about the distribution of (studentized) average outcomes in covariate-treatment cells with very few observations,” in addition to saying that the normality assumption “is clearly restrictive; but without imposing some additional structure it would seem impossible to conduct valid inference in the presence of small groups.”} They also use fixed designs not only for covariates but also treatment status. While it is standard in the literature on finite-sample inference to condition on covariates for testing purposes,\footnote{See, e.g., lehmann1993fisher, lehmann2005testing,zhang2022randomization.} it is not standard to also condition on the treatment status. Conditioning on the treatment status may be even more problematic in the limited overlap setting, because such an approach may understate actual uncertainty in the resulting estimates.\footnote{In the limited overlap scenario, very few treated (or untreated) observations are used to estimate the conditional means of the relevant potential outcomes. However, the number of those few treated (or untreated) observations is itself a random variable. Conditioning on the treatment statuses of the observations ignores this source of randomness, potentially leading to underreported uncertainty in the estimates of the conditional means.}

The finite-sample framework of my paper goes beyond the EATE parameter, which is the focus of rothe2017robust, armstrong2021finite, and hong2020inference.\footnote{To be more precise, rothe2017robust and armstrong2021finite focus on the EATE, and hong2020inference focus on the finite-population ATE, which is related to but conceptually different from EATE. Another important note is that rothe2017robust and armstrong2021finite use different terms for the same EATE parameter. rothe2017robust calls it the “sample average treatment effect” (SATE), but this term is used for a different parameter in much of the literature. armstrong2021finite use the term “conditional average treatment effect” (CATE) to refer to the EATE parameter, but their usage is very nonstandard in the literature. They themselves say, “We note that the terminology varies in the literature. Some papers call this object the sample average treatment effect (SATE); other papers use the terms CATE and SATE for different objects entirely.” To avoid all this confusion, I ensure that the names and abbreviations of some common parameters used in this paper, as defined in Section (ref), are the same as those used in the vast majority of the causal inference literature.} I propose finite-sample estimation and inference methods for each term comprising the EATE, which is the sample average of heterogeneous treatment effect means, rather than just the aggregate EATE object. Unlike armstrong2021finite, I do not impose assumptions on the smoothness or the shape of the conditional means of the outcomes (for finite-sample results), and so my focus is on finite-sample unbiased estimation and valid finite-sample inference rather than “optimal” procedures. To allow for fine partitions of the covariate space (e.g., using unsupervised feature learning), I use the strata setup of rothe2017robust, hong2020inference, and dias2021inference but with two generalizations. Unlike their setups, I allow the number of strata to potentially grow with the overall sample size, while also allowing overlap to be arbitrarily low within each stratum. I do not impose the restrictive normality assumptions used by rothe2017robust and armstrong2021finite. In addition, unlike these papers, I do not condition on the treatment status variables, which are the source of the fundamental missing (counterfactual) data problem. I instead use a design-based finite-sample approach\footnote{See, e.g., lehmann2005testing, wu2021randomization, abadie2020sampling, imbens2021causal, ding2016randomization, young2019channeling, fisher1925statistical, fisher1935design,zhang2022randomization,xie2013confidence.} that I believe is simpler, broader, and more transparent.

Even though the heterogeneous average treatment effects are point-identified in the above setting, it is typically not straightforward to obtain unbiased point-estimates of the objects of interest. This difficulty is a result of some practical statistical issues with the reciprocal of the estimated propensity score. Thus, in Section (ref), I develop a “Finite-Sample Stable Probability Weighting” (FPW) set-estimator for analyzing heterogeneous average effects of multivalued treatments. I define a notion of finite-sample unbiasedness for set-estimators (a generalization of the usual concept for point-estimators), and I show that the property holds for the FPW set-estimator. An interesting feature of the FPW set-estimator is that it actually reduces to a point-estimator for many practical purposes and thus serves as a simpler alternative to other available set-estimators that are wider lee2021bounding in the (high-dimensional) strata setting.

After discussing finite-sample-unbiased estimation of average effects within strata, I propose new finite-sample inference methods for testing a general class of “weak” null hypotheses, which specify a hypothesized value for average effects. Except for the case of binary outcomes,\footnote{See rigdon2015randomization and li2016exact. See caughey2021randomization for another exception.} currently there only exist asymptotically robust tests\footnote{See wu2021randomization,chung2013exact,chung2016multivariate.} of general weak null hypotheses that have finite-sample validity under some “sharp” null hypotheses, which restrictively specify hypothesized values of counterfactual outcomes for all the observations. I thus tackle the question of how to conduct finite-sample tests of some general weak null hypotheses. I also use partial identification for this purpose. The finite-sample inference methods that I propose are computationally convenient, and they can be used to construct valid confidence sets and to bound the finite-sample confidence distribution.\footnote{Since sharp null hypotheses are a subset of weak null hypotheses, my proposals also offer a computationally fast way to compute the finite-sample confidence distributions when sharp null hypotheses are of interest.} These contributions (to the literature on finite-sample inference) are of independent interest. An appealing aspect of my finite-sample methods is that they are valid under weak assumptions and are yet simple conceptually and computationally.

Standard Observational Setup and Parameters of Interest

In this section, I formalize the observational setup and parameters of interest that are used throughout the paper. The parameters and the associated empirical setting that I consider are quite standard in the literature on causal inference under unconfoundedness imbens2000role,hirano2004propensity, rosenbaum1983central,imbens2015causal. However, I do not impose the restrictive assumption of strict overlap (boundedness of the reciprocals of propensity scores), which is crucial for many of the results in the existing literature. Thus, the following setup is more general than the standard one that is usually used hirano2003efficient, kennedy2022towards, kennedy2022minimax, ai2021unified, kennedy2017non,semenova2021debiased.

Let $\mathbb{W} \subset \mathbb{R}$ be a bounded set of possible treatment states, and let $Y^*_w$ be the random variable representing the potential outcome under the treatment state $w \in \mathbb{W}$. However, the fundamental problem of causal inference is that the random element $(Y^*_w)_{w\in \mathbb{W}}$ is not completely observable; only one component of it is not latent. The observable component is $Y = Y^*_W$, where $W$ is itself a random variable on $\mathbb{W}$ so that $W$ represents the treatment status. When $\mathbb{W} = \{0,1\}$ and $W$ is the associated binary random variable, $Y = W\,Y^*_1+(1-W)\,Y^*_0$. More generally, when $\mathbb{W}$ is discrete, $Y = \sum_{w \in \mathbb{W}} \mathbb{I}\{W = w\}\,Y^*_w$. When $\mathbb{W}$ is an open interval, $Y = \int_{\mathbb{W}}\, Y^*_{\tilde{w}}\,\delta(W - \tilde{w})\,d\tilde{w} = Y^*_W$, where $\delta(\cdot)$ is the Dirac delta function. The random elements $(Y^*_w)_{w \in \mathbb{W}}$ and $W$ may all depend on characteristics or covariates $X$, where $X$ is a random vector on a subset $\mathbb{X}$ of some Euclidean space. The only observable random elements in this setting are $(Y, X, W)$, but a lot can be learned from them. The relationship between $W$ and $X$ is characterized by an object called the propensity score, which plays a crucial role in causal analysis rosenbaum1983central. Building on imbens2000role and hirano2004propensity, the propensity score can be generally defined as follows.

Note that Definition (ref) is applicable even if $W$ is not a purely discrete or a purely continuous random variable. For example, if $W$ is a mixed continuous and discrete random variable, then $\tilde{p}(\cdot \,|\,x)$ in the above definition is the Radon--Nikodym derivative of the probability measure (induced by $W \,|\,X = x$) with respect to a combined measure, which equals the Lebesgue measure plus another measure on $\mathbb{R}$ for which the measure of any Borel set is equal to the number of integers on that Borel set. Definition (ref) can be applied even in a setting where $W$ is a random vector.

It is relatively less difficult to learn about the marginal distributions of the potential outcomes (under some assumptions) than about the joint distribution of the potential outcomes. Thus, it is useful to focus on parameters such as the “Conditional Average Response” (CAR) function and the `Conditional Average Contrast” (CAC), which are defined as follows.

In addition to CAR and CAC, one may also be interested in parameters that depend on the conditional marginal distributions of the potential outcomes. For this purpose, it is useful to define the “Conditional Response Distribution” (CRD) as follows.

It is difficult to identify the above parameters (CAR, CAC, and CRD) without further assumptions. To make parameter identification feasible, a large strand of literature on causal inference---as reviewed by imbens2004nonparametric, imbens2015causal, and wager2020causal---assumes unconfoundedness (also known as selection on observables, conditional exogeneity, or strongly ignorable treatment assignment), population overlap (i.e., at least weak overlap), and SUTVA (stable unit treatment value assumption). These notions are formalized in Assumptions (ref)--(ref) below.

The above assumptions are sufficient to identify CAR as follows for all $(w, x)\in \mathbb{W} \times \mathbb{X}$: $$\mu_w(x) = \mu(w, x) = \mathbb{E}[Y^*_w \mid X = x] = \mathbb{E}[Y^*_W \mid W = w, X = x] = \mathbb{E}[Y \mid W = w,\, X = x].$$ Note that the conditioning event in the above expression is $\{W = w,\, X = x\}$. Assumption (ref), i.e., $p(w, x) > 0$, ensures that $\{W = w,\, X = x\}$ is not a null event. It is also possible to identify CAR using the conditioning event $\{X = x\}$. For example, when $\mathbb{W}$ is discrete, CAR can be identified as follows for all $(w, x)\in \mathbb{W} \times \mathbb{X}$: $$\frac{\mathbb{E}[\,\mathbb{I}\{W = w\}\,Y \mid X = x]}{\mathbb{E}[\,\mathbb{I}\{W = w\}\mid X = x]} = \frac{\mathbb{E}[\,\mathbb{I}\{W = w\}\,Y^*_w \mid X = x]}{p(w, x)} = \frac{p(w, x)\,\mathbb{E}[Y^*_w \mid X = x]}{p(w, x)} = \mu_w(x).$$ When $\mathbb{W}$ is an interval, if $K_h(W - w) \equiv h^{-1}\,K\big(\frac{W - w}{h}\big)$, where $K(\cdot)$ is a kernel and $h > 0$, then $$\mathrm{lim}_{\,h \,\downarrow\, 0}\frac{\mathbb{E}[\,K_h(W-w)\,Y \mid X = x]}{\mathbb{E}\,[K_h(W-w) \mid X = x]} = \frac{\int_\mathbb{W}\,\delta(\tilde{w} - w)\,\mu(\tilde{w}, x)\,p(\tilde{w}, x)\,d\tilde{w}}{\mathbb{E}[\delta(W-w) \mid X = x]} = \frac{\mu(w, x)\,p(w, x)}{p(w, x)} = \mu_w(x).$$ Thus, the positivity of the propensity score plays a key role in the identification of CAR. However, note that Assumption (ref) is not the same as the more restrictive condition of strict (or strong) overlap defined below. Although a vast number of results in the existing literature depend on strict overlap, it is a testable assumption that may lack empirical support ma2022testing,lei2021distribution. In this paper, I only assume population overlap (Assumption (ref)), which allows for arbitrarily low propensity scores, i.e., “weak\,/\,limited overlap,” rather than the following condition.

Since CAR is identified under Assumptions (ref)--(ref), CAC and CRD are also identified because they both can be expressed using different forms of CAR. Although these three assumptions are adequate for identifying the parameters of interest, estimating them is still challenging without more restrictions (on the outcomes, covariates, and observations), such as the following assumptions used by, e.g., kitagawa2018should,kennedy2022towards,kennedy2022minimax.

Several authors kitagawa2018should,kennedy2022towards,kennedy2022minimax use Assumption (ref), which restricts the support of the potential outcomes, because it is convenient analytically; it also rules out cases that are practically unimportant but theoretically pathological. It implies that CAR and CAC are also bounded in magnitude. Assumption (ref) may seem restrictive, but it is indeed satisfied in practice for most economic or health outcomes of interest and also, of course, for binary outcomes. Nevertheless, Assumption (ref) is not strictly needed; it is possible to do the analysis in this paper by dropping it and using a weaker version of it. For example, one could instead assume that the conditional marginal distributions of the potential outcomes (conditional on covariates) satisfy standardized uniform integrability romano2012uniform. Alternatively, one could use an even weaker assumption that $\mathbb{E}[(Y^*_w)^2 \,|\,X = x] < \infty$ for all $(w, x) \in \mathbb{W} \times \mathbb{X}$. Assumptions (ref) and (ref) are also quite standard in the literature hirano2003efficient.

The majority of the literature on causal inference under unconfoundedness focuses on the empirically important setting where observations can be classified as “treated” and “untreated” (or “control”) units. In this setting, $\mathbb{W} = \{0, 1\}$, and $W$ is a binary indicator of treatment status; the event $\{W = 0\}$ refers to the control state, and the event $\{W = 1\}$ refers to the treatment state. In this case, the parameters CAR, CAC, and the related parameters all have special names as follows.

Assumptions (ref), (ref), and (ref) identify all of the above parameters except for SATE, which is not nonparametrically identified because it directly involves the unobserved counterfactual outcomes rather than expectations. However, note that PATE equals the expected value of both EATE and SATE. Other parameters, such as quantile treatment effects, are briefly considered in a later subsection, but the rest of the paper mostly focuses on aspects of the CATE function. Although it would be ideal to know CATE at every point on $\mathbb{X}$, an interpretable low-level summary of CATE would be much more useful in practice. For this purpose, an object called “Best Summary of CATE” (BATE) is defined as follows and includes the “Group Average Treatment Effect” (GATE) vector (defined below) and the PATE as special cases. GATEs are common parameters of interest in the literature because they are easily interpretable. The finite-dimensional BATE parameter is of interest for another obvious reason: estimates of the higher-dimensional (or infinite-dimensional) CATE function may have too much statistical uncertainty to be practically useful.

Since CATE is a specific version of CAR, which is identified, it follows that BATE is also identified. The propensity of being treated is useful for expressing BATE in terms of observables.

Specifically, since $e_i = \mathbb{E}[W_i \mid X_i]$ and the expected value of $\mathcal{I}_i \equiv \frac{(W_i - e_i)\,Y_i}{e_i\,(1-e_i)}= \frac{W_i\,Y_i}{e_i} - \frac{(1-W_i)\,Y_i}{1 - e_i}$ given $X_i$ is $\mathbb{E}[\,\mathcal{I}_i \mid X_i] = \mathbb{E}\Big[\frac{W_i\,Y_i}{e_i} \Big\lvert \,X_i\Big] - \mathbb{E}\Big[\frac{(1-W_i)\,Y_i}{1 - e_i} \Big\lvert \,X_i\Big] = \frac{e_i\,\mathbb{E}[Y^*_{1,i} \mid X_i]}{e_i} - \frac{(1-e_i)\,\mathbb{E}[Y^*_{0,i} \mid X_i]}{1 - e_i}= \tau(X_i)$, $$\mathbb{E}[Z_i\,Z'_i]^{-1}\mathbb{E}[Z_i\,\mathcal{I}_i] = \mathbb{E}[Z_i\,Z'_i]^{-1}\mathbb{E}[Z_i\,\mathbb{E}[\,\mathcal{I}_i \mid X_i]] =\mathbb{E}[Z_i\,Z'_i]^{-1}\mathbb{E}[Z_i\,\tau(X_i)] = \beta$$ is identified under Assumptions (ref)--(ref). However, these conditions alone do not guarantee that the sample analog of $\mathbb{E}[Z_i\,Z'_i]^{-1}\mathbb{E}[Z_i\,\mathcal{I}_i]$ is $\sqrt{n}$-consistent and asymptotically normal, even in the simple case where $\widetilde{Z}(\cdot) \equiv 1$ khan2010irregular. This motivates the next section.

Non-Inverse Alternatives to Inverse Probability Weighting

If strict overlap (in Definition (ref)) holds in addition to Assumptions (ref)--(ref), then it is quite straightforward to prove that $\widehat{\mathbb{E}}_n[Z_i\,Z'_i]^{-1}\widehat{\mathbb{E}}_n[Z_i\,\mathcal{I}_i]$, where $\widehat{\mathbb{E}}_n[\,\cdot\,] \equiv \frac{1}{n}\sum_{i = 1}^n [\,\cdot\,]$, is a $\sqrt{n}$-consistent and asymptotically normal estimator of $\beta$ in the case where $W$ is binary. This estimator is based on a simple regression of $\mathcal{I}_i$ on $Z_i$; strict overlap (together with the other assumptions) ensures that $0 \prec \mathbb{V}[Z_i(\mathcal{I}_i - Z'_i\beta)] \prec \infty$, where $\mathbb{V}[\,\cdot\,]$ denotes the covariance matrix. The regressand $\mathcal{I}_i$ is called a “pseudo-outcome” kennedy2022towards based on “Inverse Probability Weighting” (IPW), which is a widely used tool for causal inference. There are two versions of it hirano2003efficient, hahn1998role that are practically different but are based on the same underlying idea.

In the above expressions, the propensity score is a nuisance function but plays a central role in the identification of CAR and CAC. However, it is possible to identify CAC (and thus also CAR) in another manner, for which the propensity score is relatively less central. A widely used strategy in the literature kennedy2022towards,semenova2021debiased,su2019non,chernozhukov2022locally, based on the proposal of robins1994estimation, augments the IPW pseudo-outcome to make identification of CAC robust to misspecification of either the propensity score function or the conditional means of the potential outcomes. The introduction of these additional nuisance functions limits the importance of the propensity score for identifying CAC.

To understand the implications of limited overlap (i.e., scenario where strict overlap does not necessarily hold) for IPW or AIPW pseudo-outcome regressions kennedy2022towards, it is useful to consider the case with binary treatment rather than the multivalued treatment setting. This helps in keeping the focus on key ideas rather than cumbersome notation. Hence, I deliberately do not return to the (conceptually similar) multivalued treatment setting until later in Section (ref).

For flexible estimation of CATE in the binary treatment setting, the latest methods in the literature kennedy2022towards involve nonparametric regressions of $\mathcal{A}[\{0, 1\}, (-1, 1); \widehat{\gamma}, \widehat{\varphi}, X_i]$ on $X_i$, where $\widehat{\gamma}$ and $\widehat{\varphi}$ are nonparametric estimates of the respective nuisance parameters that are constructed using a separate auxiliary sample. Convergence rates of such estimators crucially depend on the risk of the “oracle estimator,” which nonparametrically regresses $\mathcal{A}[\{0, 1\}, (-1, 1); \gamma, \varphi, X_i]$ on $X_i$. This “oracle risk” further depends on the conditional variance of $\mathcal{A}[\{0, 1\}, (-1, 1); \gamma, \varphi, X]$, which is the AIPW pseudo-outcome. Similarly, the conditional variance of the IPW pseudo-outcome $\mathcal{I}[\{0, 1\}, (-1, 1); \gamma, \varphi, X]$ determines the risk of the oracle estimator that nonparametrically regresses $\mathcal{I}[\{0, 1\}, (-1, 1); p, X_i]$ on $X_i$. However, without strict overlap, the conditional variances of the IPW and AIPW pseudo-outcomes may not necessarily be bounded, as discussed below. Thus, the associated oracle estimators may themselves be unstable under limited overlap.

The above arguments, related to the points made by khan2010irregular and others, show that the (augmented) IPW technique is a double-edged sword, because even simple propensity score models heiler2021valid, such as $e(X) = X = \Phi(2\,X^*)$ with $X^* \sim \mathcal{N}(0,1)$, violate strict overlap. The main component of the (A)IPW pseudo-outcome (i.e., the denominator equalling the conditional variance of the binary treatment indicator) is powerful and empirically useful (under strong overlap) because it allows for nonparametric estimation of PATE without strong restrictions on the heterogeneity in the unit-level treatment effects hirano2003efficient. However, the inverse weighting component is also the reason for the potential instability of the (A)IPW technique in many empirically relevant scenarios (involving limited overlap). One need not look beyond a very simple setting to understand the issue: a situation where the propensity score function is known and the CATE has a linear form. My approach to the problem in this setting, as demonstrated in this section, illustrates and motivates my general framework in Section (ref) that tackles the challenge from semi/non-parametric perspectives. Hence, throughout this section, I make the following assumption and also treat the propensity score function $e(\cdot)$ as known.

I now define a class of “General Probability Weighting” (GPW) estimators that may be used to estimate the above parameter of interest $\beta$. This estimator class includes the usual IPW estimator but also “Non-Inverse Probability Weighting” (NPW) estimators that are defined below. As shown and demonstrated later, NPW estimators are asymptotically normal under Assumptions (ref)--(ref) and Assumption (ref) even if the IPW estimator may not be $\sqrt{n}$-consistent and asymptotically normal.

Even in the simple case where $Z_i = 1$ “a.s.” (almost surely), the IPW estimator $\widehat{\beta}_{-1}$ may not be $\sqrt{n}$-estimable khan2010irregular under Assumptions (ref)--(ref) and Assumption (ref). Specifically, ma2020robust show that the distribution of $\widehat{\beta}_{-1}$ depends on the distribution of the propensity scores and may thus be non-Gaussian (e.g., asymmetric L\'evy stable) asymptotically. More generally, the “usual” (multivariate) central limit theorems may not necessarily be applicable to the sample mean of $(e_i(1-e_i))^{\nu}Z_i(W_i-e_i)Y_i$ when $\nu < 0$, in which case the random variable may potentially have infinite variance (if overlap is weak\,/\,limited). In this case, generalized central limit theorems may need to be applied instead, but they typically imply asymptotic non-Gaussianity of $\widehat{\mathbb{E}}_n[Z_i(\nu)Y_i(\nu)]$ (and thus of $\widehat{\beta}_\nu$) when $\nu < 0$ under limited overlap. However, I show in this section that the NPW estimators are asymptotically normal unlike the IPW estimators (or, more generally, GPW estimators with negative indices) using the following proposition.

Figures (ref) through (ref) in the Appendix illustrate Proposition (ref) by contrasting the sampling distributions of the IPW estimator $\widehat{\beta}_{-1}$ and the NPW estimator $\widehat{\beta}_1$ when $n = 10^5$ (sample size) in the following simple setting with limited overlap: $X \sim \mathrm{Uniform}(0,1)$, $W \sim \mathrm{Bernoulli}(X^4)$, $(\upsilon_1, \upsilon_2) \sim \mathrm{Uniform}[(-2, 2)^2]$, $Y = 10\,(1-X^4) + X^4\,\upsilon_1 + W\,(\tau(X) + 2\,\upsilon_2)$, where $\tau(X) = \beta'Z$ with $\beta = (3, -2)$ and $Z = (1,X)$ so that $\tau^* = \mathbb{E}[\tau(X)] = 3 - 2\,\mathbb{E}[X] = 2$. As the figures show, the NPW estimator in this setting is approximately Gaussian but the IPW estimator is not.

The theoretical underpinning of GPW (and thus NPW) estimators is the moment condition $$\mathbb{E}[Z_i\{e_i(1-e_i)\}^{\nu}\{(W_i - e_i)Y_i - e_i(1-e_i)Z'_i\beta\}] = 0$$ for any $\nu \in \mathbb{R}$, since $\mathbb{E}[Z_i\{e_i(1-e_i)\}^{\nu}\,\mathbb{E}[(W_i - e_i)Y_i - e_i(1-e_i)Z'_i\beta \mid X_i]] = 0$ by Proposition (ref). However, an alternative moment condition is $\mathbb{E}[Z_i\{(W_i - e_i)Y_i - (W_i - e_i)^2Z'_i\beta\}] = 0$, which holds because $\mathbb{E}[(W_i - e_i)^2Z'_i\beta \mid X_i] = e_i(1-e_i)Z'_i\beta$ almost surely, resulting in the estimator $\widehat{\mathbb{E}}_n[(W_i - e_i)^2Z_iZ'_i]^{-1}\,\widehat{\mathbb{E}}_n[Z_i(W_i-e_i)Y_i],$ which can be obtained by regressing $Y_i$ on $(W_i - e_i)Z_i$. (The connection between this estimator and the so-called “robinson1988root transformation” is discussed in the next section.) The equality $\mathbb{E}[Z_i\{(W_i - e_i)Y_i - \frac{1}{2}[W_i(1-e_i)+e_i(1-W_i)]Z'_i\beta\}] = 0$ could also be exploited because $\mathbb{E}[\frac{1}{2}[W_i(1-e_i)+e_i(1-W_i)]Z'_i\beta \mid X_i] = e_i(1-e_i)Z'_i\beta$ almost surely, and so $\widehat{\mathbb{E}}_n[\frac{1}{2}[W_i(1-e_i)+e_i(1-W_i)]Z_iZ'_i]^{-1}\,\widehat{\mathbb{E}}_n[Z_i(W_i-e_i)Y_i]$ is another possible estimator. Similar reasoning as in Proposition (ref) can be used to establish the asymptotic normality of both of these alternative estimators even under limited overlap when Assumptions (ref)--(ref) and (ref) hold. In addition, when one-sided overlap holds, e.g., if $1 - e_i$ is bounded away from zero almost surely, the estimator $\widehat{\mathbb{E}}_n[W_i\,Z_iZ'_i]^{-1}\widehat{\mathbb{E}}_n[Z_i(W_i - e_i)Y_i/(1-e_i)]$ is asymptotically normal under the same assumptions (even if may have bad properties in finite samples or under weak identification asymptotics) because $\mathbb{E}[W_i\tau(X_i) - (W_i - e_i)Y_i/(1-e_i) \mid X_i] = e_i\tau(X_i) - e_i(1-e_i)\tau(X_i)/(1-e_i) = 0$ almost surely. All of these estimators, including the NPW estimator(s), are based on a more general principle called “Stable Probability Weighting” (SPW), which I develop in the next section.

Stable Probability Weighting for Estimation and Inference

Under Assumptions (ref)--(ref), the Conditional Average Treatment Effect (CATE) function $\tau(\cdot)$ is nonparametrically identified without assuming strict overlap and without using propensity scores because $\mathbb{E}[Y \mid W = 1,\, X = x] - \mathbb{E}[Y \mid W = 0,\, X = x] = \mathbb{E}[Y^*_1 \mid X = x] - \mathbb{E}[Y^*_0 \mid X = x] = \mu_1(x) - \mu_0(x) = \mu(1,x) - \mu(0,x) = \tau(x)$ for all $x \in \mathbb{X}$. This suggests that one could potentially learn about CATE by simply estimating $\mu$ through a nonparametric regression of $Y$ on $(W, X)$ and then use a plug-in version of $\mu(1,x) - \mu(0,x)$ to estimate CATE. However, the resulting naive plug-in estimator may be too noisy even if the CATE function itself is smooth and even when strict overlap holds; see Section 2.2 of kennedy2022towards. Such an approach may be especially problematic from an inferential perspective under weak overlap; see Section 3.2 in the Supporting Information of kunzel2019metalearners. An alternative is to use a parametric (or semiparametric) model for $\mu$, but this approach may be highly sensitive to model specification. A more fruitful approach is to treat $\mu$ as a nuisance parameter and to focus on the target parameter $\tau$ in addition to exploiting the data structure (i.e., unconfoundedness) more fully. kennedy2022towards and others (cited in Section (ref)) take such an approach under the assumption of strict overlap. However, to the best of my knowledge, the empirically relevant setting of limited overlap is largely ignored in the existing literature on heterogeneous causal effects. Hence, this section develops the relevant large-sample theory.

In Subsection (ref), I develop a basic form of a general principle called “Stable Probability Weighting” (SPW) for learning about CATE under limited overlap. I then further develop a more robust version of it called Augmented SPW (ASPW) in Subsection (ref). In Subsection (ref), I generalize the (A)SPW framework for learning about the Conditional Average Contrast (CAC) function as well as distributional or quantile treatment effects in settings with multivalued treatments. In Subsection (ref), I discuss methods for estimating and inferring parameters of interest from parametric, semiparametric, nonparametric, machine learning, and policy learning perspectives.

Main Conditional Moment Restrictions Robust to Limited Overlap

Without further ado, I define the main Stable Probability Weighting (SPW) principle for CATE.

It would be desirable for the generalized residual $\Psi(Y, X, W; \tau, e, \eta)$ to satisfy an important property called “Neyman orthogonality” chernozhukov2017double, chernozhukov2018double, mackey2018orthogonal, which enables application of some recent (but increasingly popular) nonparametric or machine learning methods, in addition to having zero conditional mean and finite conditional variance. The residual $\Psi(Y, X, W; \tau, e, \eta)$ satisfies Neyman orthogonality if the G\^ateaux derivative of its conditional expectation is zero almost surely with respect the last two (nuisance) arguments (evaluated at their true values $e$ and $\eta$) so that the conditional moment equality $\mathbb{E}[\Psi(Y, X, W; \tau, e, \eta) \mid X] = 0$ (a.s.) is insensitive to local deviations in the values of the nuisance parameters, as defined below.

The discussion so far has focused on the SPW principle, which relies on the propensity score function $e(\cdot)$, but it is also possible to construct generalized residuals that are robust to limited overlap without using $e(\cdot)$. For this purpose, it is useful to define a broader “Stable Residualization Principle” (SRP), which nests the SPW principle.

The generalized residuals discussed so far treat the CATE function $\tau(\cdot)$ as the target parameter and the functions $(\mu_0(\cdot), \mu_1(\cdot))$ as nuisance parameters (in addition to $e(\cdot)$, the propensity score function). However, when $(\mu_0(\cdot), \mu_1(\cdot))$ is itself of interest, another SPW principle can be used.

Augmented Conditional Moment Restrictions for Double Robustness

I now return to case where CATE is of interest and CTM\,/\,CCM are treated as nuisance components. Robustifying the conditional moment restrictions (for CATE) to the nuisance parameters is an important consideration in the context of semi/non-parametric (as well as parametric) statistical methods. While the Neyman orthogonality property discussed in the previous subsection provides a certain form of “local” robustness (i.e., low sensitivity to marginally minute errors in the nuisance components around their true values), it does not cover the case where there may be systematic mistakes in one of the nuisance functions. It would be desirable to construct conditional moment restrictions that have “double robustness” to misspecification of one of the nuisance components. There is a basic version of such property (assessed at the true value of CATE) and also a more global version (assessed in the space where CATE lies). Both versions are formalized below.

The BDR property in Definition (ref) is simply a conditional version of the unconditional version usually defined in the literature on “doubly robust” methods robins2001comment,rothe2019properties,chaudhuri2019review. Similarly, the GDR property is a generalization of the unconditional version defined by chernozhukov2022locally and is much stronger than BDR.

The above remark shows that the Neyman orthogonal generalized residual based on the Robinson transformation does not satisfy the ASPW principle, even under strict overlap. When overlap is strong, the AIPW residual (defined as CATE minus the AIPW pseudo-outcome) can be used instead; see Definition (ref) and the remark below it. However, even under limited overlap, there is a simple SPW-based way to fix the usual form of Robinson transformation in order to achieve BDR. There are also several other ways to construct generalized residuals that satisfy the ASPW principle under one-sided or two-sided limited overlap, as discussed below.

The above proof not only shows that GNPW residual has basic double robustness (BDR) but also implicitly shows that it does not have global double robustness (GDR) in general. Nevertheless, it does satisfy the ASPW principle. The above theorem also provides a way to fix the Robinson transformation-based residual (RT residual) to make it doubly robust. The generalized residual $(W - e(X))^2\tau(X) - (W - e(X))(Y - e(X)\mu_1(X) - (1-e(X))\mu_0(X))$ is the GNPW residual with $(\vartheta_1, \vartheta_2, \vartheta_3, \vartheta_4) = (1, 0, -2, 1)$ and has a vector-valued nuisance function $(\mu_0(\cdot), \mu_1(\cdot))$ in addition to $e(\cdot)$. Since this residual satisfies the ASPW principle, it can be used instead of the usual RT residual given by $(W - e(X))^2\tau(X) - (W - e(X))(Y - \eta(X))$, which only uses one nuisance function $\eta(X) \equiv \mathbb{E}[Y \mid X]$ in addition to $e(\cdot)$. As argued in Remark (ref), the RT residual does not have double robustness because correctly specifying $\eta(\cdot)$, which implicitly involves the true value of the propensity score function, while misspecifying $e(\cdot)$ leads to a systematic imbalance that leads to a nonzero conditional expectation of the RT residual in that case. Thus, surprisingly, when $e(\cdot)$ is misspecified as $\tilde{e}(\cdot)$, it is better to use the residual $(W - \tilde{e}(X))^2\tau(X) - (W - \tilde{e}(X))(Y - \tilde{e}(X)\mu_1(X) - (1-\tilde{e}(X))\mu_0(X))$ than to use $(W - \tilde{e}(X))^2\tau(X) - (W - \tilde{e}(X))(Y - \eta(X))$, where $\eta(X) \equiv e(X)\mu_1(X) + (1-e(X))\mu_0(X)$ is the correct conditional expectation rather than the incorrect version $\tilde{e}(X)\mu_1(X) + (1-\tilde{e}(X))\mu_0(X)$; the latter is nevertheless deliberately exploited by the GNPW residual to achieve double robustness. Another generalized residual that satisfies the ASPW principle (but not GDR) is given below.

The residuals discussed in Theorems (ref) and (ref) are valid under two-sided limited overlap and satisfy the BDR condition but not the GDR condition in general. The reason is that they multiply the main parameter of interest $\tau(\cdot)$ by a function of $e(\cdot)$. However, there are two alternatives that do satisfy the GDR principle and the ASPW principle in addition to satisfying Neyman orthogonality.

Augmented Stable Probability Weighting for Multivalued Treatments

In this section, I show that the SPW principle readily extends to the case where general conditional structural parameters (defined using conditional outcome means or quantiles or distribution functions) are of interest in settings with multivalued treatments. However, to keep the discussion concise, the version of (Augmented) Stable Probability Weighting that I formulate for this purpose is less broad than the version developed previously for CATE.

Recall that the Conditional Average Response (CAR) and the Conditional Response Distribution (CRD) functions are given by $\mu(w, x) = \mathbb{E}[Y^*_w \mid X = x]$ and $\varrho(u; w, x) = \mathbb{P}\{Y^*_w \leq u \mid X = x\}$, respectively, for all $(u, w, x) \in \mathbb{R} \times \mathbb{W} \times \mathbb{X}$. The Conditional Average Contrast (CAC) function is given by $\textstyle \theta[\widetilde{\mathbb{W}}, \kappa] = \sum_{w \,\in\, \widetilde{\mathbb{W}}} \,\,\kappa_w\,\mu_w$ when $\widetilde{\mathbb{W}}$ is a finite subset of $\mathbb{W}$ and $\kappa \equiv (\kappa_w)_{w\,\in\,\mathbb{W}}$ is a $|\widetilde{\mathbb{W}}|$-dimensional vector of real constants. In addition, let $q_v(w,x) = \mathrm{inf}\{u \in \mathbb{R}: \varrho(u; w, x) \geq v\}$ be the Conditional Quantile Response (CQR) function for all $v \in (0, 1)$ and $(w, x) \in \mathbb{W} \times \mathbb{X}$. In the below discussion, $\widetilde{K}(W - w)$ denotes either $\mathbb{I}\{W - w = 0\}$ or $\mathrm{lim}_{\,h\,\downarrow\,0} \,K_h(W-w)$ for all $w \in \mathbb{W}$, depending on whether $\mathbb{W}$ is a discrete set or a continuum, respectively.

Large-Sample Statistical Methods Using Stable Probability Weighting

In this subsection, I discuss how the (A)SPW principle can be empirically implemented using several existing methods for conditional moment models. Thanks to this rich literature, this paper does not attempt to reinvent the wheel for large-sample estimation and inference.\footnote{Of course, future research can try to improve methods for efficient identification/singularity-robust semiparametric or nonparametric estimation and inference using infinite conditional moment equalities based on exogeneity.}

Let $\rho(Y, X, W; g)$, where $g$ includes the main parameters of interest and any nuisance parameters, be a (finite-dimensional) vector containing the (A)SPW-based generalized residuals for the parameters of interest as well as generalized residuals that identify the nuisance functions such that $\mathbb{E}[\rho(Y, X, W; g) \mid X] = 0$ almost surely. For example, when $\mathbb{W} = \{0, 1\}$ and CATE is of interest, the required conditional moment restriction is satisfied by setting $g = (\tau, e, \mu_0, \mu_1)$ and $\rho(Y, X, W; g) = (W - e(X),\,(1-W)(Y - \mu_0(X)),\,W(Y - \mu_1(X)),\,\Psi(Y, X, W; \tau, e, (\mu_0, \mu_1)))$, where $\Psi(Y, X, W; \tau, e, (\mu_0, \mu_1))$ is an (A)SPW generalized residual for CATE.\footnote{Sometimes a few nuisance functions, such as the $r(\cdot)$ functions in Theorems (ref) and (ref) need to be excluded from $g$. When $r(\cdot)$ is unknown in such cases, data-driven assumptions on them can be made using auxiliary sample that is not used for the main analysis. However, the main analysis will then be conditional on the auxiliary sample.} BATE and PATE satisfy $\mathbb{E}[\tau(X) - Z'\beta \mid Z] = 0$ (a.s.) and $\mathbb{E}[\tau(X) - \tau^*] = 0$. chen2014local provide a sufficient condition for (local) identification of $g$.\footnote{For example, invertibility of $\mathbb{E}[(e_i(1-e_i))^{\nu + 1}Z_iZ'_i]$, i.e., no perfect collinearity in $\sqrt{(e_i(1-e_i))^{\nu + 1}}Z_i$, where $\nu \geq 0$, is crucial for NPW-based identification of $\beta$ in Section (ref) when CATE has linearity and $e(\cdot)$ is known.} Typically, further restrictions on $g$ (or further assumptions on the setting) are needed to satisfy the sufficient identification condition of chen2014local. I cannot list them out with high specificity because they depend on the empirical context.

Nonparametric estimation of $g$ is feasible when $\mathbb{E}[\rho(Y, X, W; g) \mid X] = 0$ (a.s.) and other regularity conditions hold. One such important condition is that $\mathbb{E}[\lVert\rho(Y, X, W; g)\rVert^2 \mid X]$ is bounded, as required by ai2003efficient and newey2003instrumental for smooth generalized residuals $\rho(\cdot)$. For example, when $g = (\tau, e, \mu_0, \mu_1)$, boundedness of $\mathbb{E}[\lVert\rho(Y, X, W; g)\rVert^2 \mid X]$ under Assumptions (ref)--(ref) crucially depends on the boundedness of $\mathbb{E}[\Psi(Y, X, W; \tau, e, (\mu_0, \mu_1))^2 \mid X]$, which holds when the (A)SPW principle is used. The sieve minimum distance (SMD) estimation procedure of ai2003efficient or other alternatives can be implemented for consistent estimation of $\tau$ and to obtain $\sqrt{n}$-consistent and asymptotically normal estimators of $(\tau^*, \beta)$. When the main parameters of interest are finite-dimensional, two-step semiparametric procedures may be used.\footnote{See, e.g., bravo2020two,chernozhukov2022locally,cattaneo2010efficient,chenliao2015sieve,ai2012semiparametric,ackerberg2012practical,ackerberg2014asymptotic,otsu2007penalized,hahn2018nonparametric,andrews1994asymptotics,newey1994asymptotic,chernozhukov2018double.} When all the parameters are finite-dimensional, there exist simpler methods.\footnote{See, e.g., dominguez2004consistent,kitamura2004empirical,antoine2014conditional,andrews2013inference,jun2012testing,jun2009semiparametric,andrews2017inference. If unconditional moment equalities identify the parameters, one can also use the methods of andrews2019identification,andrews2012estimation,andrews2016conditional.}

In the general case where $\rho(\cdot)$ is possibly nonsmooth, penalized SMD estimators chen2012estimation can be used, and sieve-based tests chen2015sieve,hong2017inference help with inference on functionals of $g$. chen2016methods provide a useful survey of other available methods. chernozhukov2022constrained propose valid procedures when constraints, such as shape restrictions, are imposed on $g$. When $\tau$ (or the main structural parameter) is weakly identified but a strongly identified functional of $\tau$, such as $(\tau^*, \beta)$, is of interest, the penalized minimax (adversarial) approach of bennett2022inference can be used for inference. Some recent high-dimensional methods dong2021high,chang2018new,chang2021high,nekipelov022regularized may also be used.\footnote{High-dimensional statistical methods can be used by turning the conditional moment restrictions into unconditional moment equalities using high-dimensional instruments and parameters. For example, when PATE, BATE, and CATE are of interest so that $g = (\tau^*, \beta, \tau, e, \mu_0, \mu_1)$, high-dimensional logit and linear regression models may be specified for $e$ and $(\mu_0, \mu_1, \tau)$, respectively, and the corresponding high-dimensional regressors may be used as instruments. Then, under some restrictions on the high-dimensionality of the model (e.g., parameter sparsity), estimation and inference are possible using (penalized) sieve generalized method of moments dong2021high or using high-dimensional penalized empirical likelihood methods chang2018new,chang2021high. Regularized machine learning methods nekipelov022regularized may also be used based on (A)SPW generalized residuals.} ASPW-based nonparametric estimators of CATE may be used for policy learning using, e.g., kitagawa2018should's (kitagawa2018should) hybrid Empirical Welfare Maximization (EWM) approach.

When $g$ only has infinite-dimensional components, e.g., $g = (\tau, e, \mu_0, \mu_1)$, and the main parameter evaluated at some $x \in \mathbb{X}$, e.g., $\tau(x)$, is of interest, then local estimation and inference methods are available. Kernel and local polynomial methods can be used for estimation.\footnote{See, e.g., lewbel2007local,carroll1998local,han2020identification,zhang2003local.} In this context of localized moment restrictions, andrews2014nonparametric and xu2020inference provide useful inference methods that are robust to identification strength. Machine learning-based methods, such as orthogonal random forests oprescu2019orthogonal, may also be used based on ASPW residuals.

Finite-Sample Stable Probability Weighting and Inference

In this section, I develop finite-sample estimation and inference methods for heterogeneous average treatment effects of finite multivalued treatments $\mathbb{W}$. The setup for finite-sample results uses a “fixed design” that treats the covariates as given.\footnote{This type of conditioning on covariates is quite common in finite-sample statistical theory. For the results on design-based inference, I also condition on the (realized or latent) potential outcomes, following common practice.} I assume that the dataset $\mathcal{D}_n \equiv (Y_i, X_i, W_i)_{i = 1}^n$ has independent observations and follows the model specified in Assumptions (ref)--(ref) below.

The above assumptions essentially specify a (potentially) high-dimensional model (conditional on covariates) with fine strata $\mathbb{X}_n$ (and $\mathbb{Z}_n$) that could depend on the sample size.\footnote{Note that the strata (i.e., the covariate set $\mathbb{X}_n$) in Assumption (ref) is different from the set $\mathbb{X}$ used for large-sample theory. Since the finite-sample analysis is conditional on covariates, the researcher may use unsupervised feature learning to form the strata $\mathbb{X}_n$ based on a richer set of background variables. For example, if a scalar background variable is available, it can be partitioned (say, based on quantiles) into disjoint bins, which then are the strata.} Henceforth, I assume without loss of generality that $|\mathbb{Z}_n | = 1$ so that $\mathbb{E}[Y^*_{w,i} \mid X_i] = \mu_w$ for some real number $\mu_w \in [\underline{C}_w, \overline{C}_w]$ for all $w\in \mathbb{W}$ and $i \in \{1, \dots, n\}$.\footnote{Assumption (ref) uses a partitioning-based high-dimensional model (similar to a a zero-degree piecewise polynomial or spline) for the propensity score with potentially $(|\mathbb{W}| - 1) \times K_n$ unknown nuisance parameters. Assumption (ref) uses (potentially) coarser strata for the Conditional Average Response (CAR) function but still could be high-dimensional. If we are interested in estimating CAR at a particular value $z \in \mathbb{Z}_n$, then only the units $i$ with $X_i \in \tilde{Z}^{-1}(z)$ need to be used, and those units would then comprise the dataset for finite-sample analysis. Since the dataset $\mathcal{D}_n$ could be redefined to contain only those units, we may assume without loss of generality that $\tilde{Z} \equiv 1$ so that $|\mathbb{Z}_n | = 1$.} I am interested in “unbiased” estimation of the Conditional Average Contrast (CAC) parameter $\theta = \sum_{w \in \mathbb{W}} \kappa_w\,\mu_w$ for a given $(\kappa_w)_{w \in \mathbb{W}} \in \mathbb{R}^{|\mathbb{W}|}$.

When $K_n = 1$, an unbiased estimator of $\mu_w$ ($w \in \mathbb{W}$) is $\frac{1}{n}\sum_{i = 1}^n \mathbb{I}\{W_i = w\}\,Y_i/\lambda_{w,1}$, but this is not computable because $\lambda_{w,1}$ is unknown. A conventionally used alternative estimator of $\mu_w$ is the subsample mean\footnote{I do not call the subsample mean the IPW estimator for reasons that will become clear later.} $\big(\sum_{i = 1}^n \mathbb{I}\{W_i = w\}\,Y_i\big)/\big(\sum_{i = 1}^n \mathbb{I}\{W_i = w\}\big) = \frac{1}{n}\sum_{i = 1}^n \mathbb{I}\{W_i = w\}\,Y_i/\hat{\lambda}_{w,1}$, where $\hat{\lambda}_{w,1} = \frac{1}{n}\sum_{i = 1}^n \mathbb{I}\{W_i = w\}$. However, the expectation of this subsample mean does not exist because the denominator $\hat{\lambda}_{w,1}$ can be zero with positive probability, especially under limited overlap. Since unbiased point-estimation may not be generally feasible in this setting, I focus on set-estimators, for which I define a generalized notion of “unbiasedness” as follows.

To inform unbiased set-estimation, it is useful to compute the finite-sample bias of a (biased) estimator whose expectation does exist (unlike the subsample mean mentioned before). For this purpose, some sensible estimator (with a finite expectation) of the reciprocal of the propensity score is needed. However, in general, it is impossible to construct an unbiased estimator of the reciprocal of a Bernoulli parameter voinov1993unbiased. Nevertheless, a (biased) shrinkage estimator fattorini2006applying is sometimes used for practical purposes. I use a leave-one-out version of it below.

An interesting feature of the estimator $\widetilde{\mu}_{w,k}$ is that it is numerically equivalent to a modified subsample mean, i.e., $\widetilde{\mu}_{w,k} \equiv \big(\sum_{i \in \mathcal{X}_k} \mathbb{I}\{W_i = w\}\,Y_i\big)/\mathrm{max}\{1, \sum_{i \in \mathcal{X}_k} \mathbb{I}\{W_i = w\}\}$, which does have a finite expectation, for all $(w, k) \in \mathbb{W} \times \mathbb{X}_n$. If it is possible to construct an unbiased set-estimator of the bias of $\widetilde{\mu}_{w,k}$, then an unbiased set-estimator of $\mu_w$ can also be constructed for any $w \in \mathbb{W}$. I provide such a construction below but without pooling information across strata.

Although the set-estimators provided above are unbiased, it may be more efficient to pool information across strata, especially when there is strong overlap in some strata. I do this using the “Finite-Sample Stable Probability Weighting” (FPW) set-estimator below and show that it is unbiased. It can be seen as a counterpart of ASPW estimators in the finite-sample setting.

The contribution $\bar{\mu}_{w,i}(t)$ of each observation $i \in \{1, \dots, n\}$ to the FPW set-estimator $\widehat{\mathcal{U}}_w$ (when $\theta = \mu_w$) for some $w \in \mathbb{W}$ is very intuitive. When the stratum of $i$ has information about $\mu_w$, then $\bar{\mu}_{w,i}(t)$ uses a type of stabilized probability weighting $\widehat{R}_{w,i}\,\mathbb{I}\{W_i = w\}\,Y_i$ using the stable reciprocal weight $\widehat{R}_{w,i}$. Otherwise, $\bar{\mu}_{w,i}(t)$ uses a certain type of imputation based on information (about $\mu_w$) in other strata. Theorem (ref) shows that this sort of case-specific imputation is unbiased in a sense. The FPW approach can be interpreted as a finite-sample counterpart of the large-sample ASPW approach, which is conceptually similar.\footnote{My FPW approach is similar in spirit to that of lee2021bounding, but they consider a broader setting where the parameters of interest may not even be point-identified, so the set-estimators they propose are generally wider. Under my assumptions, the parameters are indeed point-identified. Thus, in this setting, the FPW set-estimator, which actually reduces to a point-estimator for many practical purposes, is much simpler conceptually and computationally.} The FPW estimator can also be seen as a bias-corrected version of the Weighted Modified Difference (WMD) estimator $\big[\sum_{w\in \mathbb{W}} \kappa_w \big(\sum_{k \in \mathbb{X}_n} \frac{N_k}{n}\,\widetilde{\mu}_{w,k}\big)\big] \equiv \big[\sum_{w\in \mathbb{W}} \kappa_w \big(\sum_{k \in \mathbb{X}_n} \frac{N_k}{n}\,\big(\sum_{i \in \mathcal{X}_k} \mathbb{I}\{W_i = w\}\,Y_i\big)/\mathrm{max}\{1, \sum_{i \in \mathcal{X}_k} \mathbb{I}\{W_i = w\}\}\big)\big]$. Both the FPW and WMD estimators can be contrasted with the IPW estimator $\big[\sum_{w\in \mathbb{W}} \kappa_w \big(\sum_{k \in \mathbb{X}_n} \frac{N_k}{n}\,\breve{\mu}_{w,k}\big)\big]$, where $\breve{\mu}_{w,k} \equiv \frac{1}{N_k}\sum_{i \in \mathcal{X}_k} \mathbb{I}\{W_i = w\}\,Y_i/\mathrm{max}\{\widehat{P}_{w,i}, 1/(2N_k - 2)\}$ that uses the leave-one-out propensity score estimate $\widehat{P}_{w,i} \equiv (\sum_{j \neq i\,|\, X_j = X_i} \mathbb{I}\{W_j = w\})/(N_{X_i} - 1)$. It follows from Jensen's inequality\footnote{Note that $\mathbb{E}[1/\mathrm{max}\{\widehat{P}_{w,i}, [2(N_k - 1)]^{-1}\}] \geq 1/\mathbb{E}[\widehat{P}_{w,i}] = 1/\lambda_{w,X_i}$ if $\lambda_{w,X_i} \geq 1/[2(N_k - 1)]$, because the following function is convex for all $k \in \mathbb{X}_n$: $1/P$ if $P \geq 1/(N_k - 1)$, and $2(N_k-1) - (N_k - 1)^2P$ otherwise.} that $\breve{\mu}_{w,k}$ is a biased estimator of $\mu_w$ for all $(w, k) \in \mathbb{W} \times \mathbb{X}_n$. Figures (ref) through (ref) in the Appendix illustrate the unbiasedness property of the FPW estimator (compared with the IPW and WMD estimators that can have very large biases and variances) in a simple setting.

The leave-one-out propensity score estimate $\widehat{P}_{w,i} \equiv (\sum_{j \neq i\,|\, X_j = X_i} \mathbb{I}\{W_j = w\})/(N_{X_i} - 1)$ is used by the FPW and IPW estimators (in different ways) for finite-sample estimation of the CAC parameter $\theta$. The IPW estimator is a point-estimator but it is biased. The FPW is unbiased but it is a set-estimator. A natural question is whether it is possible to obtain an unbiased point-estimator (potentially based on $\widehat{P}_{w,i}$) for any particular forms of CAC $\theta = \sum_{w \in \mathbb{W}} \kappa_w\,\mu_w$, where the constants $(\kappa_w)_{w \in \mathbb{W}} \in \mathbb{R}^{|\mathbb{W}|}$ may themselves potentially depend on unknown (nuisance) parameters. The below theorem addresses this question and also paves the way for practical finite-sample inference.

Note that the moment equality in the above theorem strongly resembles the conditional moment equality for SPW generalized residuals; this is not surprising because they are based on the same underlying idea. Clearly, Theorem (ref) implies that $\widehat{T}_{a,b}/\mathcal{G}_{a,b}$ is an (infeasible) unbiased point-estimator of the average treatment effect $\mu_a - \mu_b$. Despite the infeasibility of $\widehat{T}_{a,b}/\mathcal{G}_{a,b}$ for finite-sample estimation, $\widehat{T}_{a,b}/\mathcal{G}_{a,b}$ is highly useful for finite-sample inference. Suppose $T^{**}_{a,b}$ is a random draw from the finite-sample distribution corresponding to $\widehat{T}_{a,b}$. Then, despite the fact that $\widehat{T}_{a,b}/\mathcal{G}_{a,b}$ is latent and not computable, note that $\mathbb{P}\{T^{**}_{a,b}/\mathcal{G}_{a,b} \geq \widehat{T}_{a,b}/\mathcal{G}_{a,b}\} = \mathbb{P}\{T^{**}_{a,b} \geq \widehat{T}_{a,b}\}$. Thus, the $p$-values associated with the test statistics $\widehat{T}_{a,b}$ and $\widehat{T}_{a,b}/\mathcal{G}_{a,b}$ are indeed equivalent.

I now propose a conceptually and computationally simple finite-sample inference method based on the above statistic to test null hypotheses regarding the average effects of multivalued treatments. Specifically, I am interested in “weak” null hypotheses that concern the average treatment effect vector $T^\circ = (T^\circ_w)_{w = 1}^L$, where $\mathbb{W} = \{0,1,\dots,L\}$ without loss of generality and $T^\circ_w \equiv \mu_w - \mu_0$, so that $T^\circ \equiv (\mu_1 - \mu_0, \dots, \mu_L - \mu_0)$. In contrast, “sharp” null hypotheses hypothesize a specific set of restricted values of all the missing (counterfactual) potential outcomes. Currently, there only exist asymptotically valid tests\footnote{See chung2013exact,chung2016multivariate,wu2021randomization.} of weak null hypotheses in my general setting.\footnote{There exist some exceptions, such as the case of binary outcomes, for which rigdon2015randomization and li2016exact propose some methods. caughey2021randomization provide useful methods for testing “bounded” null hypotheses with one-sided treatment effect heterogeneity, but the weak null hypotheses I consider are more general.} I use partial identification and a simple simulation-based method to bound the entire “$p$-value function” luo2021leveraging, which may be called a finite-sample version of the “confidence distribution” xie2013confidence, using only a single Monte Carlo sample.

As is standard in the literature on design-based finite-sample inference, I treat the following as fixed for inferential purposes: the (partially observed) data $\mathcal{D}^*_n \equiv ((Y^*_{w,i})_{w \in \mathbb{W}}, X_i)_{i = 1}^n$, and a set $\Lambda$ that partially identifies the nuisance parameters $\{\lambda_{w,k}: (w, k) \in \mathbb{W} \times \mathbb{X}_n\}$ by assuming bounds for them, implying the assumption that $\{\lambda_{w,k}: (w, k) \in \mathbb{W} \times \mathbb{X}_n\} \in \Lambda$. Thus, in the following results, the notation $\mathbb{P}_l\{\cdot\}$ is used to represent probability when the only source of randomness is the vector of treatment statuses that are generated according to the model $l \in \Lambda$.

For purposes of inference on the parameter $\mu_a - \mu_b$, I consider “linear” test statistics that can be expressed as $\breve{T}_{a,b} = \frac{1}{n}\sum_{i = 1}^n Q_i^{a,b}(\mathcal{W}_n)\,Y_i$, where the weight $Q_i^{a,b}(\mathcal{W}_n)$ may depend on the entire set $\mathcal{W}_n \equiv (X_i, W_i)_{i = 1}^n$ for all $i \in \{1, \dots, n\}$, for a distinct pair $(a, b) \in \mathbb{W}^2$. I prefer to choose $\breve{T}_{a,b} \equiv \widehat{T}_{a,b}$, but $\breve{T}_{a,b}$ may also be chosen to be the appropriate IPW estimator or WMD estimator. For testing weak null hypotheses, I use a flexible restriction on unit-level effect heterogeneity.

To avoid cumbersome notation, I henceforth concentrate on the case where $\mathbb{W} = \{0, 1\}$ in order to construct bounds for the $p$-value function in Theorem (ref) below. However, a conceptually (and mechanically) similar construction is possible for the general case where $\mathbb{W} = \{0, \dots, L\}$.

Stochastic approximation may be used to compute $\mathbb{P}_l\{\cdot\}$ for all $l \in \Lambda$ without affecting the finite-sample validity of the test lehmann2005testing. Only a single set of Monte Carlo draws are needed to bound the entire $p$-value function based on the above result. Valid confidence sets can also be computed by simply “inverting” luo2021leveraging all possible $p$-value functions.

\setstretch{1} {2pt}