Robust Orthogonal Machine Learning of Treatment Effects

Robust Orthogonal Machine Learning of Treatment Effects

\markboth{Robust Orthogonal Machine Learning of Treatment Effects} {Shell \MakeLowercase{et al.}: Bare Demo of IEEEtran.cls for Computer Society Journals}

\IEEEtitleabstractindextext{

}

\IEEEdisplaynontitleabstractindextext \IEEEpeerreviewmaketitle

\IEEEraisesectionheading{

Introduction

} \IEEEPARstart{H}{ow} to construct counterfactuals and estimate treatment effects correctly is the key to obtaining stable predictions and answering what if questions when one only has observational data. It is ubiquitous in many applications such as in healthcare and business when decision makers do not have the luxury to conduct randomized controlled trials (RCTs) because they are either costly, or time-consuming, or simply can not be done. There are three challenges in utilizing observational data to estimate the treatment-outcome relation. The first one is the selection bias. Selection bias arises during the selection process when subjects are differentially included or excluded from treatments. The control of selection bias requires data on all factors relevant to the selection mechanism for treatment. The second one is the confounding bias. Confounding bias results from inadequate adjustment of covariates that are simultaneously predictive of treatment and outcome. The control of confounding bias requires data on all covariates relevant to the outcome and the treatment assignment mechanism. The third one is the unknown and potentially nonlinear relationship between treatment and outcome, between treatment and covariates, and between outcome and covariates.

In healthcare applications, the selection mechanism for including or excluding a patient for treatment in RCTs is likely constrained by medical and ethical, and the outcomes are subject to covariates that are largely overlapping with those used in the selection process. For example, the Infant Health and Development Program (IHDP) hill2011bayesian in the United States conducted experiments to evaluate the efficacy of early intervention in reducing the developmental and health problems of low birth weight and premature infants. However, there is a wide range of variables (covariates) that can influence both infants' cognitive outcomes and whether or not they receive treatments, including pregnancy complications, child's gender, household composition, day care arrangements, source of health care, quality of the home environment, parents' race and ethnicity, and maternal age, education, IQ, and employment. They managed to obtain the study sample of infants which was stratified by birth weight (2,000 grams or less, 2,001-2,500 grams) and randomized to the Intervention Group or the Follow-Up Group. The goal is to study the treatment effect of the specialist visits (binary treatment) on the cognitive scores of children (continuous-valued outcome).

In an example of business applications, many e-commerce platforms provide shopping credit in the form of consumer loans to shoppers who frequent their marketplaces. In this new paradigm of consumer financing, platform lenders would like to know the precise impact of adjusting credit policies (the treatment), such as credit lines and loan interests, on shoppers' spending behavior. “It is an open question to what extent FinTech lenders are using randomization (or other techniques to recover causal effects) in order to improve the accuracy of their prediction models out-of-sample", as the authors of berg2021fintech put it. This is because conducting RCTs by boosting shoppers' credit lines is costly and waiting for changes in spending behavior to reveal is time-consuming. However, from observational data, it is difficult to assess whether the changes in spending amount (the outcome) come from covariates such as customer's age, gender, region, indebtedness, or seasonal effects, external shocks, or are indeed the results of credit policy changes (the treatment), not to mention that, in practice, similar covariates are also used by algorithms to select customers for adjusting credit lines and by how much.

In the recent literature, various machine learning-based methods have been widely adopted to estimate treatment effects. For example, researchers have suggested (i) classical linear models, e.g., Ordinary Least Square (OLS), and Least Absolute Shrinkage and Selection Operator (LASSO) belloni2014inference, bloniarz2016lasso; (ii) tree models and their ensembles, e.g., Regression Tree (RT) athey2016recursive and Causal Forests (CF) wager2018estimation, athey2019generalized; (iii) neural network models, e.g., TARNet shalit2017estimating and Dragonnet shi2019adapting. Other advanced approaches involve representation learning louizos2017causal, assaad2021counterfactual, yao2018representation, zhang2020learning, adversarial learning yoon2018ganite, du2021adversarial, and Bayesian learning hill2011bayesian, taddy2016nonparametric, alaa2017bayesian, ray2019debiased.

In general, causal machine learning methods concerning ATE estimations mainly include regression adjustment methods and re-weighting methods (see more details in 10.1145/3444944). Regression adjustment methods require an estimated feature-outcome relation (the outcome model) and directly average the predicted potential outcomes over the whole population to estimate ATE, so the associated estimator is called the direct regression (DR) estimator. The shortcoming of the DR estimator is that it overlooks the probabilistic impact of the covariates on the treatment assignment (i.e., the propensity score) and hence often results in biased estimations of ATE unless the outcome model is estimated accurately. Re-weighting methods mimic the principle of RCTs to make the re-weighted instances look like they receive alternative treatment. The Inverse Probability Weighting (IPW) is one of the prevalent re-weighting strategies. It involves the propensity scores rather than the outcome model. However, the IPW estimator is sensitive to the estimation of propensity scores, thus easily leading to high variance or even infinite estimates. Such occasion often occurs when the estimated propensity scores are close to 0 or 1. This is called the error-compounding issue.

The Debiased Machine Learning (DML) method, which is exploited by chernozhukov2018double based on neyman1979c, has become more and more popular and been applied in different fields since its appearance. Examples include the use of the DML estimators when the data is of panel type in semenova2018orthogonal, the application of DML for programme evaluation in knaus2020double, and the study of the action-response effects under the credit risk context in huang2021causal. Besides, researchers elaborate and investigate the DML estimators further. For example, DML has been extended for continuous treatment effects klosin2021automatic, dynamic treatment effects lewis2021double, and instrumental variable quantile regressions chen2021debiased. In addition, authors in oprescu2019orthogonal combine the DML estimators with the generalized random forest given in athey2019generalized to estimate the heterogeneous effects (or the conditional average treatment effects (CATEs)). Recently, mackey2018orthogonal firstly proposed orthogonal machine learning where the treatment effect estimator satisfying the higher-order orthogonality than in DML was constructed and nice theoretical properties were presented (our paper will follow the line of orthogonal machine learning works).

DML offsets the shortcomings of classical causal learning methods. The DML method combines the two classical approaches to ensure that the corresponding estimator is accurate as long as either the outcome model or the propensity score model, but not necessarily both, is correctly specified (see robins2017minimax,robins2008higher,mukherjee2017semiparametric,kang2007demystifying,van2014higher and the references therein). This notable merit is well known as the doubly robust property. However, the DML estimator can still suffer from the error compounding issue since the inverse term of the propensity score is still present. For example, assuming the true propensity score $\pi^i(\mathbf{z}_m)=0.01$ and the estimated one $\hat{\pi}^i(\mathbf{z}_m)=0.05$, they only differ with $|0.01-0.05|=0.04$. However, the error of their inverses can attain up to $|1/0.01-1/0.05|=80$. Inevitably, such situation occurs commonly in practice, especially when the distribution of the treated group is substantially different from that of the controlled group (see, for example, li2018balancing, busso2014new, knaus2020double, austin2015moving, dudik2011doubly). In practice, two variants of the DML estimator, Augmented Inverse Propensity Weighted (AIPW) robins1994estimation, rotnitzky1998semiparametric, bang2005doubly and DML-trim chernozhukov2018double, are often used to tackle this problem.

AIPW has the same mathematical form as the DML estimator, but can avoid an infinite estimate. When IPW is infinite, the two terms in AIPW that contain IPW will cancel each other out linden2016estimating. DML-trim removes the individuals with very high or low estimated propensity scores crump2009dealing or re-adjusts the propensity scores based on a specific range (e.g., in experiments of this paper, following chernozhukov2018double, we trim the estimated propensity scores with $0.01$ and $0.99$). Nevertheless, such incomplete solutions are trick-based, not from a theoretical or methodological perspective. Although we can deal with the error-compounding issue with them from a practical perspective, stabilizing the compounded error caused by IPW from a theoretical standpoint still remains a challenge. This observation motivates us to go beyond the DML estimator and construct estimators that are more robust to the extremes of the propensity scores.

In this paper, we propose a Robust Causal Learning (RCL) method to establish the RCL estimators of ATE. The contributions are summarized as follows:

The rest of the paper is organized as follows. Section (ref) introduces the problem setup and the background of orthogonal scores. Section (ref) presents the main theoretical results, including the RCL score in Theorem (ref) and the RCL estimator in Corollary (ref). In Section (ref), we theoretically prove the consistency of the RCL estimator. Section (ref) reports extensive experimental results on simulation datasets, benchmark datasets, and a mimic consumer credit dataset generated by WGAN. We defer some proof details in the appendix.

Problem Formulation and Assumptions

The Problem Setup

In this paper, we consider the potential outcome framework rubin1974estimating, rubin2005causal to study the ATE estimation. Let {$\mathbf{Z}$} be the covariates (confounders), and {$D$} be the treatment variable which can take values from {$\{d^{1},\dots, d^{n}\}$}. We denote {$Y$} as the outcome variable (the response), and {$Y^{i}$} represents the potential outcome under the treatment {$d^{i}$}. If the observed treatment is {$d^{i}$}, then the factual outcome {$Y^{F}$} equals {$Y^{i}$}. We denote {$\{w_m=(y_m, d_m, \mathbf{z}_m)\}_{m=1}^{N}$} as the observed $N$ realizations of i.i.d. random variables {$\{W_m=(Y_m, D_m, \mathbf{Z}_m)\}_{m=1}^{N}$}. Given the true causal parameter {$\theta^{i}:=\mathbb{E}[Y^{i}]$}, the target quantity ATE between treatment $d^{i}$ and treatment $d^{j}$ is defined as

Identifying ATE {$\theta^{i,j}$} under the potential outcome framework requires some fundamental assumptions to ensure that {$\theta^{i}$} and {$\theta^{j}$} are identifiable. Thus, we impose the following assumptions as stated in existing causal inference literature.

These assumptions guarantee that ATE can be inferred if we specify the relation {$\mathbb{E}\left[Y \mid D, \mathbf{Z}\right]$}, which is equivalent to estimating {$g^{i}(\mathbf{Z})$} for each {$i \in \{1, \dots, n\}$} in the generalized propensity score model setting\footnote{This model setting allows {$D$} to be multi-valued. It can be reduced to the iteractive model in chernozhukov2018double once the treatment {$D$} takes binary values.} imbens2000role,tu2013using:

Here, {$g^i(\cdot)$} and {$\pi^i(\cdot)$} are true nuisance parameters. {$\xi^{i}$} and {$\nu^{i}$} are the noise terms. {$\pi^i(\mathbf{Z})=\mathbb{E}\left[\mathbf{1}_{\{D=d^{i}\}} \mid \mathbf{Z}\right]$} is known as the generalized propensity score (GPS) with multi-valued treatment variable. Notice that $\xi^{i} \perp\!\!\!\perp D\mid \mathbf{Z}$, so $\mathbb{E}\left[\xi^{i}\mid \mathbf{Z}\right]=0.$ Finally, the true causal parameter {$\theta^i$} for {$i \in \{1, \dots, n\}$} can be computed by {$\theta^i:=\mathbb{E}\left[Y^{i}\right]=\mathbb{E}\left[g^i(\mathbf{Z})\right]$} and the true ATE can be computed by {$\theta^{i,j}=\theta^{i}-\theta^{j}$}.

Non-Orthogonal Scores and Orthogonal Scores

We aim to estimate the true causal parameters $\theta^{i}$ given $N$ i.i.d. samples $\{W_m=(Y_m, D_m, \mathbf{Z}_m)\}^{N}_{m=1}$. According to chernozhukov2018double, the standard procedure to obtain the estimated causal parameter $\hat{\theta}^i$ is: 1) getting the estimated nuisance parameters $\hat{\rho}$, e.g., $\hat{\rho}=(\hat{g}^i, \hat{\pi}^i)$; 2) constructing a score that satisfies the moment condition (Definition (ref)); 3) establishing the estimator of $\theta^i$, which is solved from the moment condition (ref).

The moment condition guarantees that the estimator derived from the score is unbiased if the nuisance parameters equal the true ones. Here, we give the scores which satisfy the moment condition of two classical causal learning methods (DR and IPW) introduced before.

Generally, the estimators established from the scores in Example (ref) and (ref) might be invalid unless {$\hat{g}^{i}$} and {$\hat{\pi}^{i}$} estimate {$g^{i}$} and {$\pi^{i}$} well. To obtain robust estimators, chernozhukov2018double,mackey2018orthogonal suggest that we should construct scores which satisfy the (higher-order) Orthogonal Condition (Definition (ref)) apart from the moment condition.

The orthogonal condition ensures that the established estimators can still be valid even though some nuisance parameters are misspecified (see chernozhukov2018double, robins2008higher, mackey2018orthogonal for more details). Below we demonstrate how to utilize {$\mathbf{D}^{\alpha}\psi(W,\vartheta,\varrho)$} to justify that the scores in Example (ref) and (ref) violate the orthogonal condition. Suppose $k$ in (ref) is $1$, then {

} The above calculations show that {$\psi_{DR}^{i}(W,\vartheta,\varrho)$} and {$\psi_{IPW}^{i}(W,\vartheta,\varrho)$} do not satisfy the orthogonal condition. Such scores are usually termed as the non-orthogonal scores. As a consequence, their associated estimators are not doubly robust. To obtain a doubly robust estimator, chernozhukov2018double propose the DML method to construct the DML score as we summarize below.

We can prove that {$\psi_{DML}^{i}(W,\vartheta,\varrho)$} satisfies the orthogonal condition when {$k=1$} in Equation (ref) (see chernozhukov2018double for detailed derivations) following similar calculation processes for DR and IPW. {$\psi_{DML}^{i}(W,\vartheta,\varrho)$} is therefore termed as an orthogonal score. The orthogonal condition assures that the DML estimator is doubly robust, i.e., the estimator is locally unbiased and consistent as long as either {$g^{i}$} or {$\pi^{i}$} is correctly specified. However, in spite of the doubly robust property, the DML estimator still suffers from the error-compounding issue once some encompassed inverse propensity scores are extreme. In real applications, one seldom encounters a situation that propensity scores are correctly estimated for all individuals and with no extreme values.

This dilemma motivates us to construct scores such that: 1) the scores are orthogonal scores, i.e., they satisfy the moment condition (Definition (ref)) and the orthogonal condition (Definition (ref)); 2) the estimators established from the scores can stabilize the estimation error due to the extremes of propensity scores. In the upcoming section, we will introduce a novel method, the Robust Causal Learning (RCL) method, to achieve this goal and overcome the difficulties encountered by DR, IPW, and DML methods.

The Proposed Method

This section introduces our main contribution, the RCL method proposed by us. First, Section (ref) shows the constructed RCL score. Then Section (ref) presents the detailed construction of the RCL estimator with an algorithm that describes how to obtain such an estimator of {$\theta^{i}$} from observational data.

Construction of The RCL Score

In this subsection, we construct an orthogonal score, the RCL score, to derive an estimator of {$\theta^{i}$} along the line of orthogonal machine learning works (e.g., chernozhukov2018double,mackey2018orthogonal). The main result is stated in Theorem (ref).

We might denote the consequent estimator constructed later using the above score as RCL{$_{r,k}$}. From (ref), we can observe that $a_{i}(\cdot)$, the nuisance parameter of the propensity score, is no longer in an inverse form in the RCL score. As a consequence, the established RCL estimator from (ref) can avoid the error-compounding issue. Simultaneously, the RCL score is an orthogonal score, so the RCL estimator is as doubly robust as the DML one. The proof of Theorem (ref) can be found in Appendix (ref). From the proof, one can find some intuition behind the construction of such a score.

Establishment of The RCL Estimator

In this part, we will go into detail about the establishment of the RCL estimator. To begin with, we can solve the estimator $\tilde{\theta}^{i}$ from (ref) using the emprical version of the moment condition for the RCL score (ref):

Equation (ref) is referred to as the DR estimator when the true nuisance parameter {$g^i$} is replaced by the estimated one {$\hat{g}^i$}. Equation (ref) can then be divided into two parts:

where {$\mathscr{I}$} is the sample set in which the units are all treated with {$d^{i}$} while {$\mathscr{I}^c$} is the sample set in which the units are not treated with {$d^{i}$}. It is obvious that (ref) and (ref) can be directly calculated from observational data, whereas (ref) that contains the counterfactual outcomes is unavailable to compute in a direct manner. Instead of pursuing the unobservable counterfactuals, we realize that given {$i \in \{ 1,\dots,n \}$}, {

} holds for {$\forall j \in \{ 1,\dots,n \}$}. Thus, the sample means of {$Y^i-g^i(\mathbf{Z})$} and {$(Y^i-g^i(\mathbf{Z}))A(D,\mathbf{Z};\pi^{i})$} equal zero regardless the samples come from {$\mathscr{I}$} or {$\mathscr{I}^c$}. This observation allows us to replace the sample mean of the counterfactuals in (ref) with that of the factual ones. To be specific, we first define the set {$\mathcal{A}$} such that

Then, a replaced estimator of (ref) is obtained as follows:

Consequently, (ref) can be inferred indirectly from the observational data. With (ref) and (ref) as well, the RCL estimator of {$\theta^{i}$} is finally established in Corollary (ref).

The proposed RCL estimator {$\hat{\theta}^{i}_{RCL}$} is a consistent estimator of {$\theta^{i}$} if {$(\hat{g}^{i},\hat{\pi}^{i})$} satisfy the assumptions stated in mackey2018orthogonal and chernozhukov2018double. The detailed consistency results of {$\hat{\theta}^{i}_{RCL}$} are in the next section and some pivotal proofs are left in the appendix due to the space limit. We also outline the procedure of estimating {$\theta^{i}$} from observational data using the proposed RCL method in Algorithm (ref). Note that if the whole dataset is split into the training set and the test set, Step 2 in Algorithm (ref) will be only conducted on the training set, while Step 3-8 can be performed to obtain the estimate {$\hat{\theta}^{i}_{RCL}$} on either the training set or the test set. The running complexity of our algorithm is obviously at most {$O(NR)$}.

Theoretical Analysis

To facilitate the upcoming studies, we first introduce some notations. These will assist us in the following consistency analysis of the RCL estimator. Recall that $Y^{i}$ is the potential outcome under the treatment $d^{i}$, and {$\mathscr{I}$} is the sample set in which the units are all treated with {$d^{i}$} while {$\mathscr{I}^c$} is the sample set in which the units are not treated with {$d^{i}$}. We use $Y^{i;F}$ as the factual outcome if an individual from {$\mathscr{I}$} actually receives $d^i$, and $Y^{i;CF}$ as the counterfactual outcome if an individual from {$\mathscr{I}^c$} receives $d^i$. Based on the introduced notations, we can define two residual differences $\xi_{m}^{i;F}$ and $\xi_{\bar{m}}^{i;CF}$ for $m\in\mathscr{I}$ and $\bar{m}\in\mathscr{I}^{c}$. Mathematically, $\xi_m^{i;F}:=Y_m^{i;F}-g^{i}(\mathbf{Z}_m)$ if $m \in \mathscr{I}$ and $\xi_{\bar{m}}^{i;CF}:=Y_{\bar{m}}^{i;CF}-g^{i}(\mathbf{Z}_{\bar{m}})$ if $\bar{m} \in \mathscr{I}^c$.

We give some assumptions on statistical properties between $\xi_{m}^{i;F}$ and $\xi_{\bar{m}}^{i;CF}$. To be precise, we study $\xi_{m}^{i;F}$ and $\xi_{\bar{m}}^{i;CF}$ conditioning on $\mathbf{Z}$. It is an abbreviation of conditioning on $(\mathbf{Z_m},\mathbf{Z_{\bar{m}}})$ or on all $\mathbf{Z_m}$ and $\mathbf{Z_{\bar{m}}}$, and this applies throughout in the rest of the paper. First, $\xi_{m}^{i;F}$ and $\xi_{\bar{m}}^{i;CF}$ are independent conditioning on $\mathbf{Z}$ due to the SUTVA assumption and have the same distribution, i.e.,

From the above assumptions, we have $\mathbb{E}\big[(\xi_{m}^{i;F})^{r}\mid\mathbf{Z}\big]=\mathbb{E}\big[(\xi_{\bar{m}}^{i;CF})^{r} \mid\mathbf{Z}\big]$, $\forall\;r,\;m$, and $\bar{m}$.

From the construction of our RCL estimator, it is easy to find that for $\bar{m}\in\mathscr{I}^c$, we replace $\xi_{\bar{m}}^{i;CF}$ with a $\xi_{m}^{i;F}$ for a random $m$ picked from $\mathscr{I}$. For the convenience of mathematical derivations, we denote an i.i.d. copy of this $\xi_{m}^{i;F}$ as $\xi_{\bar{m}}^{i;F}$. So for any $\bar{m}\in\mathscr{I}^c$, both $\xi_{\bar{m}}^{i;CF}$ and $\xi_{\bar{m}}^{i;F}$ are well-defined now. Moreover, $\xi_{\bar{m},u}^{i;F}$ is the $u^{\text{th}}$ i.i.d. copy of $\xi_{m}^{i;F}$ similarly.

In the remaining sequel, we investigate the consistency of our RCL estimator based on some basics of orthogonal machine learning theory. To start with, we state the assumptions on the consistency rates for the nuisance parameters in Assumption (ref). Only the assumptions that are helpful in studying the consistency of our RCL estimator are stated. Other assumptions that concentrate on the conditions of the scores, including orthogonality, identifiability, non-degeneracy, smoothness, and the regularity of moments can be found in mackey2018orthogonal and references therein.

Consistency Results

For notational convenience, we rewrite our RCL estimator $\hat{\theta}^{i}_{RCL}$ and denote it as $\hat{\theta}^{i}_{N}$. Indeed, we have

Here, $A^i_m=A(D_m,\mathbf{Z}_m;\pi^i)$ and $\hat{A}^i_m$ is the estimate of $A^i_m$. Besides, we define two quantities $\hat{\tilde{\theta}}_{N}^{i}$ and $\hat{\bar{\theta}}_{N}^{i}$. They are {

} We also define

Then (ref), (ref), and (ref) can be rewritten as {

} for simplicity. In addition, we have to use two lemmas and two propositions to study the consistency of $\hat{\theta}_{N}^{i}$. We state them with the proofs of the two lemmas as follows. The full proofs of the remaining Proposition (ref) and Proposition (ref) are left in Appendix (ref) and Appendix (ref) respectively.

Now, we are ready to investigate if the estimator $\hat{\theta}_{N}^{i}$ is a consistent estimator of $\theta^{i}$. Our goal is to show that $\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\theta}_{N}^{i}-\theta^{i}\right| \geq \epsilon\right\}\overset{p}{\rightarrow} 0$, where $\hat{\rho}=(\hat{g}^i, \hat{\pi}^i)$ are the estimated nuisance parameters. Before presenting the proof, we notice that for $m\in \mathscr{I}^{c}$, $\left(\hat{\xi}^{i;F}\overset{d}{=}\hat{\xi}^{i;CF}\right)\mid\mathbf{Z}$ when $\hat{g}^{i}$ and $g^{i}$ satisfy Assumption (ref). This is because that $\hat{\xi}^{i;F}=Y^{i;F}-\hat{g}^{i}(\mathbf{Z})=\xi^{i;F} + g^{i}(\mathbf{Z})-\hat{g}^{i}(\mathbf{Z})$, $\hat{\xi}^{i;CF}=Y^{i;CF}-\hat{g}^{i}(\mathbf{Z})=\xi^{i;CF} + g^{i}(\mathbf{Z})-\hat{g}^{i}(\mathbf{Z})$, {$(\xi^{i;F} \perp \!\!\! \perp D) \mid \mathbf{Z}$}, {$(\xi^{i;CF} \perp \!\!\! \perp D) \mid \mathbf{Z}$}, and $\left(\xi^{i;F}\overset{d}{=}\xi^{i;CF}\right)\mid\mathbf{Z}$. Our concluding result is that under the assumptions stated before, the estimator $\hat{\theta}_{N}^{i}$ is indeed a consistent estimator of $\theta^{i}$. Due to the space limit, we defer the full proof of consistency in Appendix (ref).

Numerical Studies

In this section, we compare the performances of our RCL estimators with the DR, IPW, DML estimators and their variants through simulation and empirical experiments. In both experiments, we consider three types of regressors for $g^i$: Lasso, Random Forests (RF), and Multi-layer Perceptron (MLP); and three types of classifiers for $\pi^i$: Logistic Regression (LR), RF, and MLP. We combine the regression model A and the classification model B to estimate {$g^{i}$} and {$\pi^i$} respectively, and denote the combination as A+B, e.g., Lasso+LR. In the empirical experiments, we consider two additional state-of-the-art neural network models in causal learning: TARNet shalit2017estimating and Dragonnet shi2019adapting. All the experiments are run on Dell 3640 with an Intel(R) Xeon(R) W-1290P CPU at 3.70GHz, and a NVIDIA GeForce RTX 2080Ti GPU.

The RCL estimator with the values of {$r$} and {$k$} (see Theorem (ref)) is denoted as RCL{$_{r,k}$}. For all the experiments throughout the paper, we use the following two metrics to evaluate the performance:

Here, {$\epsilon_{ATE;m}$} is the relative error in the $m^{\mathrm{th}}$ experiment or simulation trial where $$\epsilon_{ATE;m}=\frac{\underset{1 \leq i<j \leq n}{\overset{}{\sum}}\left|\hat{\theta}^{i,j;m}-\theta^{i,j;m}\right|}{\underset{1 \leq i<j \leq n}{\overset{}{\sum}}\left|\theta^{i,j;m}\right|}$$ with {$\theta^{i,j;m}$} and {$\hat{\theta}^{i,j;m}$} being the true ATE and the estimated ATE between the treatment {$d^i$} and the treatment {$d^j$} in the {$m^{\mathrm{th}}$} experiment. {$n$} is the number of treatments and {$M$} is the number of experiments or simulation trials. {$\epsilon_{ATE}$} and {$\sigma_{ATE}$} are the average and standard deviation of these relative errors.

Numerical Studies on Simulation Datasets

We first introduce the data generating process (DGP) for the simulation experiments. Given the covariates {$\mathbf{Z}=(Z_{1},...,Z_{p})^{T}$} which follow a standard multivariate Gaussian distribution, the treatment variable {$D$} has the treatment space {$\{d^{1}, d^{2}, d^{3}\}$} with the corresponding probability

where the values of coefficients {$\beta_{iu}$} are randomly picked from the uniform distribution {$\mathcal{U}(-0.1, 0.1)$}. {$r_c$} is the confounding ratio ranging from {$0$} to {$1$}, and the number of covariates in {$\mathbf{Z}$} used to generate {$D$} is {$\lfloor p\cdot r_{c}\rfloor\in\mathbb{N}$}. For example, if {$p=10$} and {$r_{c} = 0.56$}, then {$p\cdot r_c = 5.6$} and {$\lfloor p\cdot r_{c}\rfloor= 5$}. We generate the potential outcome {$Y^{i}$} for treatment indices {$i \in \{1,2,3\}$} as

where {$\bm{a}_{i}$} is a {$p\times 1$} constant vector whose elements are randomly chosen from {$\mathcal{U}(0.1, 0.5)$}. We also set {$d^{1}=0.1$}, {$d^{2}=0.5$}, {$d^{3}=1$}, {$\xi^{1} \sim \mathcal{N}(0,9)$}, {$\xi^{2} \sim \mathcal{N}(0,4)$}, and {$\xi^{3} \sim \mathcal{N}(0,1)$}. Next, we generate {$N$} i.i.d. observations based on the DGP. Suppose the realized covariates of the {$m^{\mathrm{th}}$} individual are {$\mathbf{z}_{m}$}, then the actual treatment {$d_{m}$} will be $d^{k}$, where $k$ is determined by {$k=\underset{u \in \{1,2,3\}}{\arg\max} \; \pi^{u}(\mathbf{z}_{m})$}. Under the actual treatment $d^{k}$, the observed factual outcome {$y_{m}$} will correspondingly be a realization of {$Y^{k}$}.

For the simulation experiments, we compute the DR, DML, and our RCL estimators with different values of {$r$} and {$k$} (see Theorem (ref)), which is denoted by RCL{$_{r,k}$}. We then use {$\epsilon_{ATE}$} in (ref) with {$n=3$} and {$M=100$} to evaluate the performance of different estimators for each combination of the regressor {$g^{i}$} and the classifier {$\pi^{i}$} (denoted as regressor+classifier). We split every dataset by the ratio {$56\%/14\%/30\%$} as training/validation/test sets.

Consistency of RCL estimators. In this part, we first set {$r_c=1$}, {$p=5$}, and let the number of observations {$N$} vary in {$\{1, 2, 4, 8, 16\}\times 10000$}. We check the consistency of RCL estimators through simulations and report {$\epsilon_{ATE}$} in Fig. (ref). The result indicates that the error reduces when the sample size increases for our RCL estimators. Besides, we also find that when {$g^{i}$} is fitted well, {$\forall i$} (e.g., when the regressor is chosen as Lasso or RF), RCL{$_{2,2}$} performs better than DR, DML, and other RCL estimators. On the other hand, when {$g^{i}$} is not fitted well for some {$i$}, e.g., when the regressor is chosen as MLP, the DML and the RCL estimators with $k=1$ can significantly correct the bias thanks to the doubly robust property. In this case, despite similar performances produced by RCL{$_{2,1}$} estimator and the DML estimator, RCL{$_{2,1}$} still has a smaller {$\epsilon_{ATE}$}.

Varying {$r_{c}$} and {$p$}. In the following experiments, we mainly compare RCL{$_{2,2}$} with DR, and RCL{$_{2,1}$} with DML since the above simulation experiments indicate that RCL{$_{2,2}$} and DR perform similarly, and RCL{$_{2,1}$} has similar trends to DML. We set {$N=10000$} and plot {$\epsilon_{ATE}$} produced by each model combination A+B versus i) different {$r_{c}$} with {$p = 100$} in Fig. (ref) and Fig. (ref); ii) different {$p$} with {$r_{c}=1$} in Fig. (ref) and Fig. (ref). From all the four figures, we observe that the DML estimator is sensitive to the change of {$r_{c}$} and {$p$}, especially when the classifier is MLP. As analyzed before, if the estimation error of {$\pi^{i}(\mathbf{z})$} is non-negligible for some $\mathbf{z}$, the term {$\frac{1}{\hat{\pi}^i(\cdot)}$} of the DML estimator often gives extreme values especially when {$\hat{\pi}^{i}(\mathbf{z})$} is small, leading to pronounced estimation error of ATE. Indeed, any ATE estimator that involves the inverse propensity score term might face this error-compounding issue. By contrast, our RCL estimators are less volatile to the variation of {$r_{c}$} and {$p$} regardless of the choice of the classifier. For example, in Fig. (ref), we notice that when the classifier is MLP, the error of DML rises dramatically as {$p$} increases, while RCL{$_{2,1}$} performs more steadily. In addition, our RCL{$_{2,2}$} estimator overall has a smaller {$\epsilon_{ATE}$} than the DR estimator no matter how {$r_{c}$} or {$p$} varies.

Numerical Studies on Benchmark Datasets

Models. Similar to the simulation experiments, we choose Lasso, RF, and MLP as the regressors and choose LR, RF, and MLP as the classifiers. Additionally, two prevalent neural network models, TARNet and Dragonnet, are also considered for learning the nuisance parameters. According to shi2019adapting, these two neural network structures can incorporate the estimations of both {$g^{i}$} and {$\pi^{i}$} using the representation learning technique.

Settings. We implement the above methods on two widely adopted benchmark datasets for causal inference, i.e., IHDP and Twins, and then compare RCL estimators with DR, IPW, DML, and their variants AIPW and DML-trim estimators. Mathematically, both the AIPW estimator and the DML-trim estimator are the same as the DML estimator. However, empirically, AIPW and DML-trim are less prone to suffer from the extreme values of inverse propensity scores. To be precise, AIPW decomposes the estimator into two parts that both contain the IPW term (see linden2016estimating), while DML-trim trims estimated propensity scores at the cutoff points of $0.01$ and $0.99$ (see chernozhukov2018double).

We take RCL{$_{2,1}$} and RCL{$_{2,2}$} as the representatives of the general RCL estimators because for the real datasets with a relatively small sample size and a large dimension of features, the second moment estimation of {$\nu^i$} is more reliable compared to higher moment estimations. We use grid search to adjust the hyperparameters on the validation set for those machine learning models. For TARNet and Dragonnet, we use the same network structures (layers, units, regularization, batch size, learning rate, and stopping criterion) as suggested in shalit2017estimating and shi2019adapting.

IHDP. It is a widely used benchmark dataset for causal inference introduced by hill2011bayesian. IHDP dataset is constructed based on the randomized controlled experiment conducted by Infant Health and Development Program. The collected 25-dimensional confounders from the 747 samples are associated with the properties of infants and their mothers, such as birth weight and mother's age. Our aim is to study the treatment effect of the specialist visits (binary treatment) on the cognitive scores (continuous-valued outcome). By removing a subset of the treated group, the selection bias in the IHDP dataset occurs. There are 1000 IHDP datasets given in hill2011bayesian. Each dataset is split by the ratio of {$63\%/27\%/10\%$} as training/validation/test sets, which keeps consistent with shalit2017estimating.

Twins. Twins dataset is introduced by louizos2017causal and it collects twin births in USA between 1989 and 1991. The treatment {$D=1$} indicates the heavier twin while {$D=0$} indicates the lighter one; the outcome {$Y$} is a binary variable defined as the mortality in the first year; the covariates {$\mathbf{Z}$} include 30 features relevant to the parents, the pregnancy, and the birth. Similar to yoon2018ganite, we only select twins that have the same gender and both weights less than 2kg. Finally, we have {$11440$} pairs of twins whose mortality rates are {$17.7\%$} for lighter twin and {$16.1\%$} for heavier twin. To simulate an observational dataset with selection bias, we selectively choose one from the twins as the observed sample based on the covariates of the $m^{\text{th}}$ individual: {$D_m|\mathbf{Z}_m \sim$} Bernoulli(Sigmoid({$\mathbf{w}^T\mathbf{Z}_m+b$})), where {$\mathbf{w}^T \sim \mathcal{U}_{30}((-0.01,0.01)^{30})$} and {$b \sim \mathcal{N}(0,0.01)$}. We repeat this process {$100$} times, and each of the generated {$100$} datasets is split by the ratio of {$64\%/16\%/20\%$} as training/validation/test sets, which keeps consistent with yoon2018ganite.

Analysis. In Table (ref) and Table (ref), we report the performance of every model combination, measured by {$\epsilon_{ATE} \; (\pm \sigma_{ATE})$}, for IHDP and Twins experiments respectively. The smaller {$\epsilon_{ATE}$}, the better. The metric {$R_{DR}=\text{RCL}_{2,2}/\text{DR}-1$} is used to evaluate the reduction ratio in {$\epsilon_{ATE}$} of RCL{$_{2,2}$} relative to DR, and {$R_{DML}=\text{RCL}_{2,1}/\min(\text{IPW, AIPW, DML, DML-trim})-1$} is to evaluate the reduction ratio in RCL{$_{2,1}$} relative to the best estimator among IPW, AIPW, DML, and DML-trim. The negative {$R_{DR}$} ({$R_{DML}$}) indicates that the RCL estimator has a smaller {$\epsilon_{ATE}$} than the DR (IPW, AIPW, DML, and DML-trim) estimator. The smaller {$R_{DR}$} ({$R_{DML}$}), the better for our estimators.

Table (ref) reports the experimental results on IHDP datasets. It illustrates that although the DR estimator produces reasonable estimates, the RCL{$_{2,2}$} estimator has a smaller {$\epsilon_{ATE}$} than the DR estimator, with the error reduced relatively by {$0.9\%-4.8\%$}. Simultaneously, RCL{$_{2,2}$} achieves the best performance among all the estimators across all the model combinations. We also notice that even though AIPW and DML-trim avoid extreme values encountered by DML (e.g., when the classifier is chosen as RF, the inverse propensity score is estimated with an infinity value for some data points), the RCL{$_{2,1}$} estimator is still at most $67\%$ better than the best of IPW, AIPW, DML, and DML-trim estimators. More importantly, when the variance of inverse propensity scores is large (e.g., when the classifier is MLP), the improvement of RCL{$_{2,1}$} relative to DML becomes more substantial.

Table (ref) presents the experimental results on Twins datasets. It can be observed that the RCL{$_{2,2}$} estimator has a significantly smaller {$\epsilon_{ATE}$} compared with other estimators for all model combinations, and it can reduce the estimation error relatively by {$2.2\%-6.3\%$} compared with the DR method. Besides, the RCL{$_{2,1}$} estimator can reduce the estimation error by {$0.3\%-10\%$} relative to the best of IPW, AIPW, DML, and DML-trim estimators. It is also noticeable that when {$\pi^{i}$} is well specified (e.g., the case on using Dragonnet in Table (ref)), our RCL{$_{2,1}$} estimator still outperforms the DML estimator even though the error {$\epsilon_{ATE}$} produced by the DML estimator is small enough.

Numerical Studies on Credit dataset

Causal inference methods are popularly assessed with simulation datasets which are typically generated by a parametric data generating process. Though the ground truth of treatment effects is accessible in this way, such synthetic datasets fail to resemble the original real datasets. Some benchmark real datasets are used for evaluating methods by hiding some part of the full dataset with a selecting formula for some specific purpose, like IHDP and Twins. The resulting semi-synthetic datasets are used in many existing works. But such original datasets are difficult and expensive to generate. Recently, athey2021using suggest using Wasserstein Generative Adversarial Network (WGAN) to generate datasets to evaluate various treatment effect estimators. Such artificial data can closely mimic the original real-world data, and the non-parametric generating way provides a fairer evaluation to alleviate concerns of researchers choosing some specific parametric generating processes to support the proposed methods.

Credit. Following the data generating strategy in athey2021using, we generate the Credit dataset based on a consumer loan dataset collected from an online consumer credit lending platform. In the Credit dataset, the treatment is whether to increase the credit limit for an individual in June or not, and the outcome is the monthly consumption in June. The covariates of the Credit dataset are age, using days, credit status in the last month; and the average credit line, order amount, borrowing amount, repayment amount, number of orders, number of borrowings, number of repayments, handling fee, and loan term in the last three months. We randomly shuffle the samples to obtain $100$ different Credit datasets with $2000$ samples generated by WGAN, and then split by the ratio of $64\%/16\%/20\%$ as training/validation/test sets. The comparisons between the generated data and the original data are shown in Table (ref).

Result analysis. We keep the same setting as in the IHDP experiments and report results on test sets in Table (ref). It is obvious that RCL$_{2,2}$ and RCL$_{2,1}$ are overwhelmingly superior compared with DR and variants of DML, respectively. The improvements are especially pronounced for DML variants, with the estimation error reduced by $94.5\%$ at most. Moreover, our estimators protect the estimates from high variance or infinity values caused by the error-compounding issue. Table (ref) reveals that our estimators can significantly improve the ATE estimations compared to traditional estimators on datasets generated in a nonparametric way suggested by athey2021using.

In summary, DML is recognized as a better method than DR and IPW because when the DR (IPW) estimator has a notable bias due to the misspecification on {$g^{i}$} (or {$\pi^{i}$}), the DML estimator can reduce the bias if {$\pi^{i}$} (or {$g^{i}$}) is well estimated. However, the advantages of DML are not easy to achieve in practice. First, the DML estimator, which incorporates the inverse propensity score term, may give a very large estimation or even infinite value of ATE, reflecting that the DML estimator is volatile to the estimation error of propensity scores. Second, if {$g^{i}$} is approximated well enough, DML will not assuredly perform better than DR due to the high variance of the IPW term. By contrast, our RCL estimators are more practical since i) they can stabilize the error caused by the extremes of propensity scores; ii) if {$g^{i}$} is well approximated, the RCL{$_{2,2}$} estimator will outperform the DR estimator owing to the RCL scores are orthogonal scores; iii) if {$g^{i}$} is not well approximated, but {$\pi^{i}$} is correctly specified, the RCL estimator with {$k=1$} performs better than the DML estimator with smaller estimation error and slighter volatility to the estimation error of propensity scores.

Conclusions

This paper constructs the RCL scores and establishes the RCL estimators for the ATE estimation. Theoretically, we prove that the RCL scores are orthogonal scores and the RCL estimators are consistent. Numerically, the comprehensive experiments have shown that our estimators outperform the commonly used methods such as DR, IPW, AIPW, DML, and DML-trim estimators. In addition, the proposed RCL estimators have the same merit, i.e., the doubly robust property, as the DML estimator. However, unlike the DML estimator, the RCL estimators are more stable to the extremes of propensity scores than the DML estimator and its variants.

To be more specific, RCL{$_{2,2}$} always outperforms DR and RCL{$_{2,1}$} always works better than the IPW-based and DML-based methods. The comparison result between RCL{$_{2,2}$} and RCL{$_{2,1}$} varies across datasets and different combinations of machine learning regressor and classifier. Notice that there is also no consistent superiority between DR and DML (it may depend on the severity of the confounding effect or the difference between the distributions of the treated group and controlled group). We suggest choosing RCL{$_{2,2}$} or RCL{$_{2,1}$} based on the performance difference between DR and DML. For RCL{$_{r,k}$} with {$r>2$}, in the current stage it is not superior than the case with {$r=2$}. In the future research, we will investigate the optimal values of {$(r,k)$} for the RCL estimator.

Future research also includes providing interpretability for deep learning models in causal inference context using the RCL method, since deep learning has become the most attractive methodology and technique in the big data era and the interpretability is critical for researchers. In addition, the applications of RCL method in decision makings are worthwhile future works, because better causality models can answer well more {what if} problems and thus are inherently linked and helpful to decision makings.

\appendices

Proof of Theorem (ref)

Given the nuisance parameters $\varrho=(\mathcal{g}^{i},a_{i})$ and the true nuisance parameters $\rho=(g^{i},\pi^{i})$, we find out the RCL score $\psi^{i}(W,\vartheta,\varrho)$ w.r.t. the nuisance parameters $\varrho=(\mathcal{g}^{i},a_{i})$ which can be used to construct the estimator of the causal parameter $\theta^{i}:=\mathbb{E}\left[g^{i}(\mathbf{Z})\right]$. We try an ansatz of $\psi^{i}(W,\vartheta,\varrho)$ such that

where

Here, the coefficients $b_{1},\dots,b_{k-1},\bar{b}_{r}$ depend on $\mathbf{Z}$ and the moments of $\nu^{i}$ only. Using the ansatz, we notice that $\psi^{i}(W,\vartheta,\varrho)$ satisfies the moment condition, i.e., $\mathbb{E}\left[\psi^{i}(W,\vartheta,\varrho)\mid_{\vartheta=\theta^{i},\;\varrho=\rho}\right]=0$. Indeed, we have

The second last equality comes from the fact that $A(D,\mathbf{Z};\pi^{i})$ is a function of $(D,\mathbf{Z})$. The last equality comes from the fact that $(\xi^{i} \perp \!\!\! \perp D) \mid \mathbf{Z}$. Now, we aim to find out the coefficients $b_{1},\dots,b_{k-1},\bar{b}_{r}$ such that the score (ref) satisfies the orthogonal score condition in Definition (ref) for $\alpha\in S_k$ (abbreviated as a $k^{\mathrm{th}}$ score). Indeed, we need to have $\mathbb{E}\left[\partial_{\mathcal{g}^{i}}^{\alpha_{1}}\partial_{a_{i}}^{\alpha_{2}}\psi^{i}(W,\vartheta,\varrho)\mid_{\vartheta=\theta^{i},\;\varrho=\rho}\mid\mathbf{Z}\right]=0$ for all $\alpha_{1}$ and $\alpha_{2}$ which are non-negative integers such that $1\leq \alpha_{1}+\alpha_{2}\leq k$. Since $\partial_{\mathcal{g}^{i}}^{\alpha_{1}}\partial_{a_{i}}^{\alpha_{2}}\psi^{i}(W,\vartheta,\varrho)=0$ when $\alpha_{1}\geq 2$, we only need to solve the coefficients $b_{1},\dots,b_{k-1},\bar{b}_{r}$ from

$\forall q=0,\dots,k-1$. However, (ref) always holds since

Consequently, we need to find out the coefficients $b_{1}$, $b_{2}$, $\dots$, $b_{k-1}$, $\bar{b}_{r}$ from {

} $\forall q=0,\dots,k-1$. From (ref), there are $k$ equations and we need to solve the $k$ unknowns $b_{1},\dots,b_{k-1},\bar{b}_{r}$ from the $k$ equations. Generally, the $k$ unknowns could be solved uniquely.

To start with, we compute $\partial_{\mathcal{g}^{i}}^{1}\partial_{a_{i}}^{q}\psi^{i}(W,\vartheta,\varrho)$ for $q=0,\dots,k-1$. Note that

when $q=0$ and

when $1\leq q\leq k-1$. Consequently, we need to solve for $b_{1},\dots,b_{k-1}$ and $\bar{b}_{r}$ simultaneously from

From (ref), we have

Since

we understand that

As such, (ref) can be reduced as

Hence, we can solve for $\bar{b}_{r}$ such that

It remains to find out $b_{1},\dots,b_{k-1}$ from (ref). Indeed, we can simplify (ref) as

$\forall 1\leq q\leq k-1$. Now, we solve $b_{1},\dots,b_{k-1}$. We start with finding out $b_{k-1}$, followed by $b_{k-2},\;b_{k-3},\dots,b_{1}$ iteratively. When $q=k-1$, (ref) becomes

Now, when $q=k-2$, (ref) becomes

Now, suppose $b_{q+1},\dots,b_{k-1}$ are known and we want to find out what $b_{q}$ is. We have to solve it from

We can obtain $b_{q}$ from the above equation, which gives

The proof is completed.

\ifCLASSOPTIONcaptionsoff \fi

\onecolumn Now we present the theoretical proofs of some theorems, propositions, and the consistency results given in the paper. Appendix (ref) is in the main body of the paper.

Proof of Proposition (ref)

$\forall \epsilon>0$, we consider $\mathbb{P}\left\{\left|\kappa_{N}^{i;CF}-\kappa_{N}^{i;F}\right|\geq\epsilon\right\}$. Indeed, we have

Denoting $\xi_{m}^{i;CF}-\xi_{m}^{i;F}$ as $\Xi_{m}^{i}$, we have

The last equality follows from

As a consequence, we have

when $N\rightarrow \infty$. As a result, we have {$\kappa_{N}^{i;CF}-\kappa_{N}^{i;F}\overset{p}{\rightarrow} 0$}. The proof is completed.

Proof of Proposition (ref)

First, we write

$\forall\epsilon>0$, we have

We simplify the numerator term. Note that

The last equality in the above derivation follows from the fact that, conditioning on {$\mathbf{Z}$}, {$\xi_{m,u}^{i;F}$} are i.i.d. of {$\xi_{m}^{i;F}$} for any {$u\in\{1,2,\cdots,R\}$}. Indeed, we have $\mathbb{E}\left[\xi_{m,u}^{i;F}\mid\mathbf{Z}\right]=\mathbb{E}\left[\xi_{m}^{i;F}\mid\mathbf{Z}\right]$ for any $m$ and $u\in\{1,2,\cdots,R\}$. Consequently, we have $\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)\mid\mathbf{Z}\right]=0$.

In addition, we simplify the quantity {$\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)^{2}\mid\mathbf{Z}\right]$}. Note that {$\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}=\frac{1}{R}\underset{u=1}{\overset{R}{\sum}}[\xi_{m,u}^{i;F}-\xi_{m}^{i;F}]$}. We therefore have

We justify the last equality. The last equality follows from the fact that, conditioning on {$\mathbf{Z}$}, {$\xi_{m,u}^{i;F}$} are i.i.d. of {$\xi_{m}^{i;F}$} and {$\xi_{m,u}^{i;F}$} are i.i.d. of {$\xi_{m,\bar{u}}^{i;F}$} for any {$u,\;\bar{u}\in\{1,2,\cdots,R\}$}. Indeed, under the given fact, we have

and

Consequently, we have

Thus, we have

We notice that no matter we set $R\rightarrow \infty$ followed by $N\rightarrow \infty$, or we fix $R$ but let $N\rightarrow \infty$, we see that {$\mathbb{P}\left\{\left|\kappa_{N}^{i;F}-\kappa_{R,N}^{i;F}\right|\geq\epsilon\right\} \rightarrow 0$}. The proof is completed.

Proof of The Consistency of $\hat{\theta}_{N}^{i}$

$\forall\epsilon>0$, we have

Since {$\left(\xi^{i;F}\overset{d}{=}\xi^{i;CF}\right)\mid\mathbf{Z}$} and {$\left(\hat{\xi}^{i;F}\overset{d}{=}\hat{\xi}^{i;CF}\right)\mid\mathbf{Z}$}, we have {$\hat{\bar{\theta}}_{N}^{i}\overset{d}{=}\hat{\tilde{\theta}}_{N}^{i}$} by Lemma (ref). Moreover, we know that {$\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\tilde{\theta}}_{N}^{i}-\theta^{i}\right| \geq \frac{\epsilon}{2}\right\}\overset{p}{\rightarrow}0$} under the assumptions given in mackey2018orthogonal. Together with the fact that {$\hat{\bar{\theta}}_{N}^{i}\overset{d}{=}\hat{\tilde{\theta}}_{N}^{i}$}, we have {$\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\bar{\theta}}_{N}^{i}-\theta^{i}\right| \geq \frac{\epsilon}{2}\right\}\overset{p}{\rightarrow}0$} by Lemma (ref). We turn to consider the quantity {$\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\theta}_{N}^{i}-\hat{\bar{\theta}}_{N}^{i}\right| \geq \frac{\epsilon}{2}\right\}$}, and we aim to show that {$\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\theta}_{N}^{i}-\hat{\bar{\theta}}_{N}^{i}\right| \geq \frac{\epsilon}{2}\right\}\overset{p}{\rightarrow} 0$}. Notice that

Note that {$\left|\kappa_{R,N}^{i;F}-\kappa_{N}^{i;F}\right|$} and {$\left|\kappa_{N}^{i;F}-\kappa_{N}^{i;CF}\right|$} do not incorporate any terms related to the estimated {$\hat{\rho}$}. From Proposition (ref) and Proposition (ref), we conclude that ((ref)b) and ((ref)c) converge to {$0$} in probability respectively. It remains to show the convergence of ((ref)a) and ((ref)d). Consider ((ref)d) first. Since

((ref)d) is bounded above by

((ref)a) converges to $0$ in probability due to Proposition (ref). We study the quantities ((ref)b) and ((ref)c).

((ref)b) can be further bounded. If {$N^{c}$} is the size of {$\mathscr{I}^{c}$}, then we have

We see that ((ref)b) can be further bounded by

We investigate if ((ref)a), ((ref)b), and ((ref)c) converge to $0$ in probability. We consider ((ref)a) first. Recall the assumptions that {$(\xi^{i;F} \perp \!\!\! \perp D) \mid \mathbf{Z}$}, {$\left(\xi^{i;F}\overset{d}{=}\xi^{i;CF}\right) \mid \mathbf{Z}$}, and {$(\xi^{i;CF} \perp \!\!\! \perp D) \mid \mathbf{Z}$}, we have {$\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\xi_{m}^{i;F}A_{m}^{i}\overset{d}{=}\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\xi_{m}^{i;CF}A_{m}^{i}$} by Lemma (ref). Since {$\Gamma_{2;1}=\frac{N^{c}}{N}\left|\frac{1}{N^{c}}\underset{m\in \mathscr{I}^{c}}{\sum}\xi_{m}^{i;F}A_{m}^{i}-\mathbb{E}_{\hat{\rho}}\left[\xi^{i;F}A^{i}\right]\right|$}, we have

Consider the numerator term, note that it equals

Since {$A^{i}$} and {$\xi^{i;F}$} do not include the estimated nuisance parameters, {$\mathbb{E}_{\hat{\rho}}\left[\left(A^{i}\right)^{2}\mathbb{E}_{\hat{\rho}}\left[\left(\xi^{i;F}\right)^{2}\mid\mathbf{Z}\right]\right]$} is a constant. Moreover, note that {$N^{c}\rightarrow \infty$} when {$N\rightarrow \infty$}, we have

Now, we consider ((ref)b). Indeed, we have

Here, the last inequality follows from the H\"{o}lders inequality, and thereby the convergence holds $\forall q \in \{1,2,\dots,k\}$ according to Assumption 1.5 of mackey2018orthogonal. Finally, we consider ((ref)c). We can rewrite {$\Gamma_{2;3}$} as

Now, we have

Using Assumption 1.5 of mackey2018orthogonal, we can conclude that {$\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2;3}\geq\frac{\epsilon}{72}\right\}\overset{p}{\rightarrow}0$}.

Next, we come to bound ((ref)c). Since {$N^{c}$} is the size of {$\mathscr{I}^{c}$} and

we see that ((ref)c) can be further bounded by

Similarly, we can prove that ((ref)a), ((ref)b), and ((ref)c) converge to $0$ in probability when $N\rightarrow \infty$ using the arguments in proving that ((ref)a), ((ref)b), and ((ref)c) converge to $0$. As a result, the quantity ((ref)d) converges to $0$ in probability when $N\rightarrow \infty$.

Lastly, we turn to consider the quantity ((ref)a). In fact, we have

We can argue that (ref) converges to $0$ in probability as $N\rightarrow \infty$ using similar arguments when we prove that ((ref)c) converges to $0$ in probability. Simultaneously, we can argue (ref) converges to $0$ in probability as $N\rightarrow \infty$ using similar arguments when we prove that ((ref)b) converges to $0$ in probability. Consequently, we have $\hat{\kappa}_{R,N}^{i;F}-\kappa_{R,N}^{i;F}$ converges to $0$ in probability.

The proof is completed.