Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Robust Causal Learning for the Estimation of Average Treatment Effects
abstractMany practical decision-making problems in economics and healthcare seek to estimate the average treatment effect (ATE) from observational data. The Double/Debiased Machine Learning (DML) is one of the prevalent methods to estimate ATE in the observational study. However, the DML estimators can suffer an error-compounding issue and even give an extreme estimate when the propensity scores are misspecified or very close to 0 or 1. Previous studies have overcome this issue through some empirical tricks such as propensity score trimming, yet none of the existing literature solves this problem from a theoretical standpoint. In this paper, we propose a Robust Causal Learning (RCL) method to offset the deficiencies of the DML estimators. Theoretically, the RCL estimators i) are as consistent and doubly robust as the DML estimators, and ii) can get rid of the error-compounding issue. Empirically, the comprehensive experiments show that i) the RCL estimators give more stable estimations of the causal parameters than the DML estimators, and ii) the RCL estimators outperform the traditional estimators and their variants when applying different machine learning models on both simulation and benchmark datasets.
IEEEkeywordstreatment effect estimation, causal inference, economics, healthcare
Introduction
Causal inference is ubiquitous for decision-making problems in various areas such as Healthcare glass2013causal, hill2013assessing, alaa2017bayesian and Economics belloni2014inference, farrell2015robust, chernozhukov2018double. At the core of causal machine learning, estimating the average treatment effect (ATE) from observational data is challenging because some features (covariates) can influence both treatment and outcome in most practical circumstances. For example, factors such as regions and races (covariates) can affect both the vaccination (treatment) assignment and the post-vaccination infection rate (outcome). To obtain a clean ATE, one can conduct the Randomized Controlled Trials (RCTs). RCTs are regarded as the gold standard to evaluate ATE, whereas conducting RCTs is often expensive and time-consuming. As a result, more and more researchers tend to estimate ATE from observational data.
In the observational study, classical causal 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 (aka 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 conundrum 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. Nevertheless, the IPW estimator is sensitive to the estimation of propensity scores and even leads to high variance estimates. Such occasion often occurs when the estimated propensity scores are close to 0 or 1. This is called an error-compounding issue.
The Debiased Machine Learning (DML) method, which is exploited by chernozhukov2018double based on neyman1979c, 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, 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 the error-compounding issue since the inverse term of the propensity score is still present. Inevitably, propensity scores usually fail to be correctly specified 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). This observation motivates us to go beyond the DML estimator and construct estimators that are more robust to the misspecification of the estimated propensity scores.
In this paper, we propose a Robust Causal Learning (RCL) method to establish the RCL estimators of the ATE. The contributions are summarized as follows:
enumerate• Our RCL estimators robustly ease the error-compounding issue exhibited by the DML estimator since the propensity scores in the RCL estimators are no longer in an inverse form.
• The RCL estimators inherit the consistency and doubly robust property of the DML estimator.
• The RCL methodology to construct an ATE estimator can also be applied to establish other prevalent treatment effect estimators.
• Extensive experiments show that the proposed RCL method achieves superior performance than DR, IPW, DML, and their variants across various combinations of machine learning regressors and classifiers.
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). Section (ref) reports extensive experimental results on simulation datasets and benchmark datasets. Due to the space limit, we defer all proof details and the code for reproducing our experiments to the full paper version.
Preliminaries
The Problem Setup
In this paper, we consider the potential outcome framework rubin1974estimating, rubin2005causal to study ATE. Let { $\mathbf{Z}$} be the covariates (aka confounders), and { $D$} be the treatment variable which can take values from { $\{d^{1},\dots, d^{n}\}$}. We denote { $Y$} as the outcome variable (aka 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 observed $N$ realizations of the 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
equation[equation omitted — 91 chars of source]
Identifying the 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 works.
assumption[Stable Unit Treatment Value Assumption (SUTVA)]
The potential outcomes for any individual do not vary regardless of the treatment status of other individuals.
assumption[Ignorability]
Given the covariates { $\mathbf{Z}$}, the potential outcome { $Y^{i}$} is independent to the treatment assignment { $D$}, i.e., { $(Y^{1},\cdots, Y^{n})\perp\!\!\!\perp D\mid \mathbf{Z} \;\; \forall i$}.
assumption[Positivity]
Treatment assignment is not deterministic regardless of the values covariates { $\mathbf{Z}$} take, i.e., { $0<\mathbb{P}\{D=d^{i}\mid \mathbf{Z}=\mathbf{z}\}<1$} { $\forall\;i$} and { $\forall\;\mathbf{z}$}.
assumption[Consistency]
If an individual receives treatment { $d^{i}$}, his factual outcome { $Y^{F}$} is equal to the potential outcome under treatment { $d^{i}$}, i.e., { $Y^{F}=Y^{i}$} if { $D=d^{i}$}.
These assumptions guarantee that the 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 (ref):
{
equation[equation omitted — 291 chars of source]
}
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 variables. 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=(\mathbf{Z}_m, D_m, Y_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).
definition[Moment Condition]
Let { $W=(\mathbf{Z}, D, Y)$} and { $\theta$} be the true causal parameter with { $\vartheta$} being a causal parameter that lies in the causal parameter set. Denoting the nuisance parameters as { $\varrho$} and the true nuisance parameters as { $\rho$}, we say a score { $\psi(W,\vartheta,\varrho)$} satisfies the moment condition if
\begin{equation}
\begin{aligned}
\mathbb{E}\left[\psi(W,\vartheta,\varrho)|_{\vartheta=\theta,\;\varrho=\rho}\right]=0.
\end{aligned}
\end{equation}
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.
example[The Score and Estimator for DR]
Let { $\varrho=\mathcal{g}^{i}$} and $\rho=g^{i}$. In the DR method, the score { $\psi_{DR}^{i}(W,\vartheta,\varrho)$} satisfying the moment condition and the associated estimator { $\hat{\theta}_{DR}^{i}$} are
{ \begin{equation*}
\begin{aligned}
\psi_{DR}^{i}(W,\vartheta,\varrho)=\vartheta-\mathcal{g}^{i}(\mathbf{Z}); \quad
\hat{\theta}_{DR}^{i}=\frac{1}{N}\underset{m=1}{\overset{N}{\sum}}\hat{g}^{i}(\mathbf{Z}_{m}).
\end{aligned}
\end{equation*}}
example[The Score and Estimator for IPW]
Let { $\varrho=a_{i}$} and { $\rho=\pi^{i}$}. In the IPW method, the score { $\psi_{IPW}^{i}(W,\vartheta,\varrho)$} satisfying the moment condition and the associated estimator { $\hat{\theta}_{IPW}^{i}$} are
{ \begin{equation*}
\begin{aligned}
\psi_{IPW}^{i}(W,\vartheta,\varrho)=\vartheta-\frac{Y\mathbf{1}_{\{D=d^{i}\}}}{a_{i}(\mathbf{Z})}; \
\hat{\theta}_{IPW}^{i}=\frac{1}{N}\underset{m=1}{\overset{N}{\sum}}\frac{Y_{m}\mathbf{1}_{\{D_{m}=d^{i}\}}}{\hat{\pi}^{i}(\mathbf{Z}_{m})}.
\end{aligned}
\end{equation*}}
Generally, the estimators established from the scores in Example 1 might be invalid unless { $\hat{g}^{i}$} and { $\hat{\pi}^{i}$} estimate { $g^{i}$} and { $\pi^{i}$} well. To obtain robust estimators, chernozhukov2018double suggest that we should construct scores which satisfy the Orthogonal Condition (Definition (ref)) apart from the moment condition.
definition[Orthogonal Condition]
Suppose that the nuisance parameters and the true nuisance parameters are { $\gamma$}-dimensional tuples, i.e., { $\varrho=(\mathcal{h}_{1},\cdots, \mathcal{h}_{\gamma})$} and { $\rho=(h_{1},\cdots, h_{\gamma})$}. Given { $S\subseteq \mathbb{Z}_{\geq 0}^{\gamma}$}, we say a score { $\psi(W,\vartheta,\varrho)$} satisfies the orthogonal condition if
\begin{equation}
{
\begin{aligned}
\mathbb{E}\left[\mathrm{\mathbf{D}}^{\alpha}\psi(W,\vartheta,\varrho)\mid_{\vartheta=\theta,\;\varrho=\rho}\mid\mathbf{Z}\right]=0\; \forall \alpha\in S.
\end{aligned}
}
\end{equation}
{ $S$} can be any subset of { $\mathbb{Z}_{\geq 0}^{\gamma}$}. Throughout the paper, for some positive integer { $k$}, we define { $S$} as
\begin{equation}
{
\begin{aligned}
S=\{\alpha\in\mathbb{Z}_{\geq 0}^{\gamma}\mid \left\|\alpha\right\|_{1}\leq k\}
\end{aligned}
}
\end{equation}
and { $\mathrm{\mathbf{D}}^{\alpha}\psi(W,\vartheta,\varrho)=\partial^{\alpha_{1}}_{\mathcal{h}_{1}}\partial^{\alpha_{2}}_{\mathcal{h}_{2}}\cdots\partial^{\alpha_{\gamma}}_{\mathcal{h}_{\gamma}}\psi(W,\vartheta,\mathcal{h}_{1},\cdots, \mathcal{h}_{\gamma})$}.
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 Examples (ref) and (ref) violate the orthogonal condition. Suppose $k$ in (ref) is $1$.
{
equation*[equation* omitted — 718 chars of source]
}
The above calculations show that { $\psi_{DR}^{i}(W,\vartheta,\varrho)$} and { $\psi_{IPW}^{i}(W,\vartheta,\varrho)$} do not satisfy the orthogonal condition. The 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.
example[The Score and Estimator for DML]
Let { $\varrho=(\mathcal{g}^{i}, a_{i})$} and { $\rho=(g^{i},\pi^{i})$}. In the DML method, the score { $\psi_{DML}^{i}(W,\vartheta,\varrho)$} that satisfies both the moment condition and orthogonal condition and the associated estimator { $\hat{\theta}_{DML}^{i}$} are
{
\begin{equation*}
\begin{aligned}
&\psi_{DML}^{i}(W,\vartheta,\varrho)=\vartheta-\mathcal{g}^{i}(\mathbf{Z})-\frac{\bm{1}_{\{D=d^{i}\}}}{a_{i}(\mathbf{Z})}(Y-\mathcal{g}^{i}(\mathbf{Z}));\\
&\hat{\theta}_{DML}^{i}=\frac{1}{N}\overset{N}{\underset{m=1}{\sum}}\hat{g}^{i}(\mathbf{Z}_{m})+\frac{1}{N}\overset{N}{\underset{m=1}{\sum}}\frac{\bm{1}_{\{D_{m}=d^{i}\}}(Y_{m}-\hat{g}^{i}(\mathbf{Z}_{m}))}{\hat{\pi}^{i}(\mathbf{Z}_{m})}.
\end{aligned}
\end{equation*}
}
We can prove that { $\psi_{DML}^{i}(W,\vartheta,\varrho)$} satisfies the orthogonal condition when { $k=1$} in (ref) (see chernozhukov2018double for detailed derivations) following similar calculation processes for DR and IPW. { $\psi_{DML}^{i}(W,\vartheta,\varrho)$} is therefore termed as the 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. Despite the doubly robust property, the DML estimator still suffers an error-compounding issue once the encompassed inverse propensity score is slightly misspecified for some data points. In real applications, one seldom encounters a situation that propensity scores are correctly estimated for all individuals. 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 misspecifications on propensity scores.
In the upcoming section, we will introduce a novel method, the Robust Causal Learning (RCL) method, to overcome the difficulties encountered by DR, IPW, and DML methods.
The Proposed Method
This section shows our main theoretical results. First, Section (ref) demonstrates the RCL scores. Then Section (ref) presents the detailed construction of the RCL estimators with an algorithm that describes how to obtain an estimate of { $\theta^{i}$} from observational data using the proposed RCL method.
Construction of The RCL Score
In this paper, we construct an orthogonal score, the RCL score, to derive an estimator of { $\theta^{i}$} along the lines of orthogonal machine learning works (e.g., mackey2018orthogonal, chernozhukov2018double). The relevant result is stated in Theorem (ref).
theorem[RCL score]
Suppose { $\varrho$} and { $\rho$} are { $2$}-dimensional tuples such that { $\varrho=(\mathcal{g}^{i},a_{i})$} and { $\rho=(g^{i},\pi^{i})$}. Let { $r$}, { $k$} be integers s.t. $1\leq k\leq r$. Assume the local moments { $\mathbb{E}\left[(\nu^{i})^{r}\mid\mathbf{Z}\right]\neq 0$} and { $\left|\mathbb{E}\left[(\nu^{i})^{q}\mid\mathbf{Z}\right]\right|<\infty\;a.s.\;\forall\;1\leq q\leq r$}. Under the assumptions on nuisance parameters and noise terms stated in chernozhukov2018double and mackey2018orthogonal, the RCL score { $\psi_{RCL}^{i}(W,\vartheta,\varrho)$} that satisfies the moment condition and the orthogonal condition is
{
\begin{align}
\psi_{RCL}^{i}(W,\vartheta,\varrho)&=\vartheta-\mathcal{g}^{i}(\mathbf{Z})-(Y^{i}-\mathcal{g}^{i}(\mathbf{Z}))A(D,\mathbf{Z};a_{i}).
\end{align}
}
Given an integer { $r$} and an integer { $k$}, we have
\begin{subequations}
\begin{equation}
{
\begin{aligned}
&A(D,\mathbf{Z};a_{i})=\bar{b}_{r}\left[\mathbf{1}_{\{D=d^{i}\}}-a_{i}(\mathbf{Z})\right]^{r}+\\
&\;\;\mathbf{1}_{\{k\neq 1\}}\left[\underset{q=1}{\overset{k-1}{\sum}}b_{q}\left(\left[\mathbf{1}_{\{D=d^{i}\}}-a_{i}(\mathbf{Z})\right]^{q}-\mathbb{E}\left[(\nu^{i})^{q}\mid\mathbf{Z}\right]\right)\right],
\end{aligned}
}
\end{equation}
where { $\bar{b}_{r}=\frac{1}{\mathbb{E}\left[(\nu^{i})^{r}\mid\mathbf{Z}\right]}$} and the coefficient { $b_{q}$} is computed by descending order for { $q \in \{k-1,\dots,1\}$}:
\begin{equation}
{
\begin{aligned}
b_{q}&=-\bar{b}_{r}\binom{r}{q}\mathbb{E}\left[(\nu^{i})^{r-q}\mid\mathbf{Z}\right]\\
&\quad-\overset{k-1-q}{\underset{u=1}{\sum}}b_{q+u}\binom{q+u}{q}\mathbb{E}\left[(\nu^{i})^{u}\mid\mathbf{Z}\right].
\end{aligned}
}
\end{equation}
\end{subequations}
From (ref), we can observe that $a_{i}(\cdot)$, the nuisance parameter of the propensity score, is no longer in an inverse form for the RCL score. As a consequence, the established RCL estimators from (ref) can avoid the error-compounding issue. Simultaneously, the RCL scores are orthogonal scores, so the RCL estimators are as doubly robust as the DML estimator.
Establishment of the RCL estimators
In this part, we will go into detail about the establishment of the RCL estimators. 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):
subequations{
\begin{align}
\tilde{\theta}^{i}=&\frac{1}{N}\underset{m=1}{\overset{N}{\sum}}g^{i}(\mathbf{Z}_{m}) \\
&+\frac{1}{N}\underset{m=1}{\overset{N}{\sum}}(Y^{i}_{m}-g^{i}(\mathbf{Z}_{m}))A(D_{m},\mathbf{Z}_{m};\pi^{i}).
\end{align}
}
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:
subequations{
\begin{align}
(ref) &=\frac{1}{N}\underset{m\in\mathscr{I}}{\overset{\sum}}(Y^{i}_{m}-g^{i}(\mathbf{Z}_{m}))A(D_{m},\mathbf{Z}_{m};\pi^{i}) \\
&+\frac{1}{N}\underset{m\in\mathscr{I}^{c}}{\overset{\sum}}(Y^{i}_{m}-g^{i}(\mathbf{Z}_{m}))A(D_{m},\mathbf{Z}_{m};\pi^{i}),
\end{align}
}
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 \}$},
{
equation*[equation* omitted — 370 chars of source]
}
holds for { $\forall j \in \{ 1,\dots,n \}$}. Thus, the sample mean of { $(Y^i-g^i(\mathbf{Z}))$} equals 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
equation[equation omitted — 183 chars of source]
Then, a replaced estimator of (ref) is obtained as follows:
enumerate• For the { $m^{\mathrm{th}}$} unit in the set { $\mathscr{I}^{c}$}, pick an element $\xi_{m}^{i}$ from { $\mathcal{A}$} and multiply it by { $A(D_{m},\mathbf{Z}_{m};\pi^{i})$}. Repeat the process until we go through all the individuals in the set { $\mathscr{I}^{c}$};
• Compute $\frac{1}{N}\underset{m\in\mathscr{I}^{c}}{\overset{}{\sum}}\xi_{m}^{i}A(D_{m},\mathbf{Z}_{m};\pi^{i})$;
• Repeat above steps { $R$} times to eliminate the randomness brought by the random picking procedure and return the substitute estimator { $\frac{1}{R}\overset{R}{\underset{u=1}{\sum}}\left[\frac{1}{N}\underset{m\in\mathscr{I}^{c}}{\overset{}{\sum}}\xi_{m,u}^{i} A(D_{m},\mathbf{Z}_{m};\pi^{i})\right]$}.
Consequently, (ref) can be inferred indirectly from observational data. With (ref) and (ref), the RCL estimator of { $\theta^{i}$} is finally established in Corollary (ref).
corollary[RCL estimator]
Let { $R\in\mathbb{Z}^{+}$}, { $(\hat{g}^{i},\hat{\pi}^{i})$} be the estimates of { $(g^{i},\pi^{i})$}, { $\hat{\mathcal{A}}$} be { $\mathcal{A}$} by replacing { $g^{i}$} with { $\hat{g}^{i}$} in (ref), and { $\hat{\xi}_{m,u}^{i}$} be the element that is randomly selected from the set { $\hat{\mathcal{A}}$} in the { $u^\mathrm{th}$} of $R$ repeated selections. The RCL estimator { $\hat{\theta}^{i}_{RCL}$} is given by
\begin{equation}
{
\begin{aligned}
\hat{\theta}^{i}_{RCL}=&\underbrace{\frac{1}{N}\underset{m=1}{\overset{N}{\sum}}\hat{g}^{i}(\mathbf{Z}_{m})}_{(a)}+\underbrace{\frac{1}{N}\underset{m\in\mathscr{I}}{\overset{\sum}}(Y^{i}_{m}-\hat{g}^{i}(\mathbf{Z}_{m}))A(D_{m},\mathbf{Z}_{m};\hat{\pi}^{i})}_{(b)}\\
&+\underbrace{\frac{1}{R}\overset{R}{\underset{u=1}{\sum}}\left[\frac{1}{N}\underset{m\in\mathscr{I}^{c}}{\overset{\sum}}\hat{\xi}_{m,u}^{i} A(D_{m},\mathbf{Z}_{m};\hat{\pi}^{i})\right]}_{(c)}.
\end{aligned}
}
\end{equation}
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. Due to the space limit, the proofs of Theorem (ref) and the consistency of { $\hat{\theta}^{i}_{RCL}$} can be seen in the full paper version. We also outline the procedures 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 will be only conducted on the training set, while Step 3 - Step 8 can be performed to obtain the estimates of { $\hat{\theta}^{i}_{RCL}$} on both the training set and the test set. The running complexity of our algorithm is at most { $O(NR)$}.
algorithm[algorithm omitted — 2,782 chars of source]
Numerical Studies
In this section, we compare the performances of our RCL estimators with the DML estimator and the DR estimator through simulation and empirical experiments. In both experiments, we consider three types of regressors: Lasso, Random Forests (RF), and Multi-layer Perceptron (MLP); and three types of classifiers: 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 inference: TARNet shalit2017estimating and Dragonnet shi2019adapting. All the experiments are run on Dell 3640 with Intel(R) Xeon(R) W-1290P CPU at 3.70GHz, and a set of NVIDIA GeForce RTX 2080Ti GPU.
For all the experiments throughout the paper, we use the following two metrics to evaluate the performance:
subequations\begin{equation}
{
\begin{gathered}
\epsilon_{ATE}=\frac{1}{M}\underset{m=1}{\overset{M}{\sum}}\epsilon_{ATE;m} \; ;
\end{gathered}
}
\end{equation}
\begin{equation}
{
\begin{gathered}
\sigma_{ATE}=\sqrt{\frac{1}{M-1}\underset{m=1}{\overset{M}{\sum}}\left[\epsilon_{ATE;m}-\epsilon_{ATE}\right]^{2}}.
\end{gathered}
}
\end{equation}
Here, { $\epsilon_{ATE;m}$} is the weighted relative error of the $m^{\mathrm{th}}$ experiment such that { $\epsilon_{ATE;m}=\frac{\underset{\underset{1 \leq i,j \leq n}{i \neq j}}{\overset{}{\sum}}\left|\hat{\theta}^{i,j;m}-\theta^{i,j;m}\right|}{\underset{\underset{1 \leq i,j \leq n}{i \neq j}}{\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$} of the { $m^{\mathrm{th}}$} experiment. { $n$} is the number of treatments and { $M$} is the number of experiments.
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
equation[equation omitted — 354 chars of source]
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 { $p\cdot r_c$} and { $\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
equation[equation omitted — 163 chars of source]
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 { $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.
\paragraph{Consistency of RCL estimators} In this part, we 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 estimator 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}$}.
figure[figure omitted — 274 chars of source]
figure[figure omitted — 290 chars of source]
figure[figure omitted — 293 chars of source]
figure[figure omitted — 260 chars of source]
figure[figure omitted — 263 chars of source]
\paragraph{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 errors of the ATE. Indeed, any ATE estimators that involve 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 classifiers. 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
table*[table* omitted — 4,830 chars of source]
table*[table* omitted — 5,047 chars of source]
\paragraph{Models} Similar to the simulation experiments, we choose Lasso, RF, and MLP as the regressors while 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.
\paragraph{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 the extreme values. 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 the 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 the higher-moment estimations. We use grid search to adjust the hyperparameters for those general 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.
\paragraph{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.
\paragraph{Twins} Twins dataset is introduced by louizos2017causal and it collects twin births in the USA between 1989 and 1991. The treatment { $D=1$} indicates the heavier twin while { $D=0$} indicates the lighter twin; 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 weigh 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 of the two twins as the observed sample based on the covariates of $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$} Twins datasets is split by the ratio of { $64\%/16\%/20\%$} as training/validation/test sets, which keeps consistent with yoon2018ganite.
\paragraph{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$} ({ $R_{DML}=\text{RCL}_{2,1}/\min(\text{IPW, AIPW, DML, DML-trim})-1$}) is used to evaluate the reduction ratio in { $\epsilon_{ATE}$} of RCL{ $_{2,2}$} relative to DR (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.
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 more minor { $\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 the 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, 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}$} 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, 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 benchmark datasets are typically generated by a parametric data generating process. Though the ground truth of treatment effects are accessable in this way, such semi-synthetic datasets fail to resemble the original real data sets.
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}$} ({ $\pi^{i}$}), the DML estimator can reduce the bias if { $\pi^{i}$} ({ $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 the ATE, reflecting that the DML estimator is volatile to the estimation 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 misspecification on 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 errors and slighter volatility to the estimated propensity scores.
Conclusion
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 estimators 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 estimation error due to the misspecification on propensity scores than the DML estimator and its variants. In the future research, we will i) investigate the optimal values of { $(r,k)$} in Theorem (ref); and ii) provide interpretability for deep learning models in causal inference using the RCL method.
Acknowledgements
Qi WU acknowledges the support from the Hong Kong Research Grants Council [General Research Fund 14206117, 11219420, and 11200219], CityU SRG-Fd fund 7005300, and the support from the CityU-JD Digits Laboratory in Financial Technology and Engineering, HK Institute of Data Science. The work described in this paper was partially supported by the InnoHK initiative, The Government of the HKSAR, and the Laboratory for AI-Powered Financial Technologies.
Shumin MA acknowledges the support from the Guangdong Provincial Key Laboratory of Interdisciplinary Research and Application for Data Science, BNU-HKBU United International College under project code 2022B1212010006, the support from Guangdong Higher Education Upgrading Plan (2021-2025) of “Rushing to the Top, Making Up Shortcomings and Strengthening Special Features" with UIC research grant R0400001-22, and the UIC grant UICR0700019-22.
Appendices
proofs
We present the theoretical proofs of Theorems and Corollaries given in the paper.
proof[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 estimators 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
\begin{equation}
{
\begin{aligned}
\psi^{i}(W,\vartheta,\varrho)&=\vartheta-\mathcal{g}^{i}(\mathbf{Z})-(Y^{i}-\mathcal{g}^{i}(\mathbf{Z}))A(D,\mathbf{Z};a_{i}),
\end{aligned}
}
\end{equation}
where
\begin{equation}
{
\begin{aligned}
A(D,\mathbf{Z};a_{i})&=\bar{b}_{r}\left[\mathbf{1}_{\{D=d^{i}\}}-a_{i}(\mathbf{Z})\right]^{r}\\
&+\underset{q=1}{\overset{k-1}{\sum}}b_{q}\big(\big[\mathbf{1}_{\{D=d^{i}\}}-a_{i}(\mathbf{Z})\big]^{q}-\mathbb{E}\big[(\nu^{i})^{q}\mid\mathbf{Z}\big]\big).
\end{aligned}
}
\end{equation}
Here, the coefficients $b_{1},\cdots,b_{k-1},\bar{b}_{r}$ depend on $\mathbf{Z}$ and $\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
\begin{equation*}
{
\begin{aligned}
&\mathbb{E}\left[\psi^{i}(W,\vartheta,\varrho)\mid_{\vartheta=\theta^{i},\;\varrho=\rho}\right]\\
=&\mathbb{E}\left[\theta^{i}-g^{i}(\mathbf{Z})-(Y^{i}-g^{i}(\mathbf{Z}))A(D,\mathbf{Z};\pi^{i})\right]\\
=&\mathbb{E}\left[\theta^{i}-g^{i}(\mathbf{Z})\right]-\mathbb{E}\left[(Y^{i}-g^{i}(\mathbf{Z}))A(D,\mathbf{Z};\pi^{i})\right]\\
=&-\mathbb{E}\left[\xi^{i}\times A(D,\mathbf{Z};\pi^{i})\right]\\
=&-\mathbb{E}\left[\mathbb{E}\left[\xi^{i}\times A(D,\mathbf{Z};\pi^{i})\mid D, \mathbf{Z}\right]\right]\\
=&-\mathbb{E}\left[A(D,\mathbf{Z};\pi^{i})\mathbb{E}\left[\xi^{i}\mid D,\mathbf{Z}\right]\right]=0.\\
\end{aligned}
}
\end{equation*}
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 of $b_{1},\cdots,b_{k-1},\bar{b}_{r}$ such that the score (ref) satisfies the $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},\cdots,b_{k-1},\bar{b}_{r}$ from
\begin{subequations}
{
\begin{empheq}[left=\empheqlbrace]{align}
0&=\mathbb{E}\left[\partial_{a_{i}}^{k}\psi^{i}(W,\vartheta,\varrho)\mid_{\vartheta=\theta^{i},\;\varrho=\rho}\mid\mathbf{Z}\right],\\
0&=\mathbb{E}\left[\partial_{\mathcal{g}^{i}}^{1}\partial_{a_{i}}^{q}\psi^{i}(W,\vartheta,\varrho)\mid_{\vartheta=\theta^{i},\;\varrho=\rho}\mid\mathbf{Z}\right]
\end{empheq}
}
\end{subequations}
$\forall q=0,\cdots,k-1$. However, (ref) always holds since
\begin{equation*}
{
\begin{aligned}
&\mathbb{E}\left[\partial_{a_{i}}^{k}\psi^{i}(W,\vartheta,\varrho)\mid_{\vartheta=\theta^{i},\;\varrho=\rho}\mid\mathbf{Z}\right]\\
=&\mathbb{E}\left[(Y^{i}-g^{i}(\mathbf{Z}))\times\partial_{a_{i}}^{k}A(D,\mathbf{Z};a_{i})\mid_{a_{i}=\pi^{i}}\mid\mathbf{Z}\right]\\
=&\mathbb{E}\left[\mathbb{E}\left[(Y^{i}-g^{i}(\mathbf{Z}))\times\partial_{a_{i}}^{k}A(D,\mathbf{Z};a_{i})\mid_{a_{i}=\pi^{i}}\mid D,\mathbf{Z}\right] \mid \mathbf{Z}\right]\\ =&\mathbb{E}\left[\partial_{a_{i}}^{k}A(D,\mathbf{Z};a_{i})\mid_{a_{i}=\pi^{i}}\mathbb{E}\left[(Y^{i}-g^{i}(\mathbf{Z}))\mid D,\mathbf{Z}\right] \mid \mathbf{Z}\right]\\
=&\mathbb{E}\left[\partial_{a_{i}}^{k}A(D,\mathbf{Z};a_{i})\mid_{a_{i}=\pi^{i}}\mathbb{E}\left[\xi^{i}\mid \mathbf{Z}\right] \mid D,\mathbf{Z}\right]=0.
\end{aligned}
}
\end{equation*}
Consequently, we need to find out the coefficients $b_{1}$, $b_{2}$, $\dots$, $b_{k-1}$, $\bar{b}_{r}$ from
{
\begin{align}\tag{(ref)}
\mathbb{E}\left[\partial_{\mathcal{g}^{i}}^{1}\partial_{a_{i}}^{q}\psi^{i}(W,\vartheta,\varrho)\mid_{\vartheta=\theta^{i},\;\varrho=\rho}\mid\mathbf{Z}\right]&=0
\end{align}
}
$\forall q=0,\cdots,k-1$. From (ref), there are $k$ equations and we need to solve the $k$ unknowns $b_{1},\cdots,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,\cdots,k-1$. Note that
\begin{equation*}
{
\begin{aligned}
\partial_{\mathcal{g}^{i}}^{1}\partial_{a_{i}}^{q}\psi^{i}(W,\vartheta,\varrho)&=-1+A(D,\mathbf{Z};a_{i})
\end{aligned}
}
\end{equation*}
when $q=0$ and
\begin{equation*}
{
\begin{aligned}
\partial_{\mathcal{g}^{i}}^{1}\partial_{a_{i}}^{q}\psi^{i}(W,\vartheta,\varrho)&=\bar{b}_{r}\frac{r!(-1)^{q}[\mathbf{1}_{\{D=d^{i}\}}-a_{i}(\mathbf{Z})]^{r-q}}{(r-q)!}\\
&\quad+\overset{k-1}{\underset{u=q}{\sum}}b_{u}\frac{u!(-1)^{q}[\mathbf{1}_{\{D=d^{i}\}}-a_{i}(\mathbf{Z})]^{u-q}}{(u-q)!}
\end{aligned}
}
\end{equation*}
when $1\leq q\leq k-1$. Consequently, we need to solve for $b_{1},\cdots,b_{k-1}$ and $\bar{b}_{r}$ simultaneously from
\begin{subequations}
\begin{equation}
{
\begin{aligned}
1&=\mathbb{E}\left[A(D,\mathbf{Z};\pi^{i})\mid\mathbf{Z}\right]
\end{aligned}
}
\end{equation}
and
\begin{equation}
{
\begin{aligned}
0&=\mathbb{E}\left[\bar{b}_{r}\frac{r!(-1)^{q}[\mathbf{1}_{\{D=d^{i}\}}-\pi^{i}(\mathbf{Z})]^{r-q}}{(r-q)!}\mid\mathbf{Z}\right]\\
&\quad+\mathbb{E}\left[\overset{k-1}{\underset{u=q}{\sum}}b_{u}\frac{u!(-1)^{q}[\mathbf{1}_{\{D=d^{i}\}}-\pi^{i}(\mathbf{Z})]^{u-q}}{(u-q)!}\mid\mathbf{Z}\right].
\end{aligned}
}
\end{equation}
\end{subequations}
From (ref), we have
\begin{equation}\tag{(ref)*}
{
\begin{aligned}
1&=\bar{b}_{r}\mathbb{E}\left[\left(\mathbf{1}_{\{D=d^{i}\}}-\pi^{i}(\mathbf{Z})\right)^{r}\mid\mathbf{Z}\right]\\
&\quad+\underset{q=1}{\overset{k-1}{\sum}}b_{q}\mathbb{E}\left[\left(\mathbf{1}_{\{D=d^{i}\}}-\pi^{i}(\mathbf{Z})\right)^{q}\mid\mathbf{Z}\right]\\
&\quad-\underset{q=1}{\overset{k-1}{\sum}}b_{q}\mathbb{E}\left[\mathbb{E}\left[(\nu^{i})^{q}\mid\mathbf{Z}\right]\mid\mathbf{Z}\right].
\end{aligned}
}
\end{equation}
Since
\begin{equation*}
\begin{aligned}
\mathbb{E}\left[\left(\mathbf{1}_{\{D=d^{i}\}}-\pi^{i}(\mathbf{Z})\right)^{q}\mid\mathbf{Z}\right]&=\mathbb{E}\left[(\nu^{i})^{q}\mid\mathbf{Z}\right]\quadand\\
\mathbb{E}\left[\mathbb{E}\left[(\nu^{i})^{q}\mid\mathbf{Z}\right]\mid\mathbf{Z}\right]&=\mathbb{E}\left[(\nu^{i})^{q}\mid\mathbf{Z}\right],
\end{aligned}
\end{equation*}
we understand that
\begin{equation*}
\begin{aligned}
\mathbb{E}\left[\left(\mathbf{1}_{\{D=d^{i}\}}-\pi^{i}(\mathbf{Z})\right)^{q}\mid\mathbf{Z}\right]-\mathbb{E}\left[\mathbb{E}\left[(\nu^{i})^{q}\mid\mathbf{Z}\right]\mid\mathbf{Z}\right]=0.
\end{aligned}
\end{equation*}
As such, (ref) can be reduced as
\begin{equation*}
{
\begin{aligned}
\bar{b}_{r}\mathbb{E}\left[\left(\mathbf{1}_{\{D=d^{i}\}}-\pi^{i}(\mathbf{Z})\right)^{r}\mid\mathbf{Z}\right]&=1\\
\Rightarrow \bar{b}_{r}\mathbb{E}\left[(\nu^{i})^{r}\mid\mathbf{Z}\right]&=1.
\end{aligned}
}
\end{equation*}
Hence, we can solve for $\bar{b}_{r}$ such that
\begin{equation*}
\begin{aligned}
\bar{b}_{r}=\frac{1}{\mathbb{E}\left[(\nu^{i})^{r}\mid\mathbf{Z}\right]}.
\end{aligned}
\end{equation*}
It remains to find out $b_{1},\cdots,b_{k-1}$ from (ref). Indeed, we can simplify (ref) as
\begin{equation}
{
\begin{aligned}
\bar{b}_{r}\mathbb{E}\left[\frac{r!(\nu^{i})^{r-q}}{(r-q)!}\mid\mathbf{Z}\right]+\overset{k-1}{\underset{u=q}{\sum}}b_{u}\mathbb{E}\left[\frac{u!(\nu^{i})^{u-q}}{(u-q)!}\mid\mathbf{Z}\right] &=0\\
\Rightarrow \bar{b}_{r}\binom{r}{q}\mathbb{E}\left[(\nu^{i})^{r-q}\mid\mathbf{Z}\right]+\overset{k-1}{\underset{u=q}{\sum}}b_{u}\binom{u}{q}\mathbb{E}\left[(\nu^{i})^{u-q}\mid\mathbf{Z}\right] &=0
\end{aligned}
}
\end{equation}
$\forall 1\leq q\leq k-1$. Now, we solve $b_{1},\cdots,b_{k-1}$. We start with finding out $b_{k-1}$, followed by $b_{k-2},\;b_{k-3},\cdots,b_{1}$ iteratively. When $q=k-1$, (ref) becomes
\begin{equation*}
{
\begin{aligned}
0&=\bar{b}_{r}\binom{r}{k-1}\mathbb{E}\left[(\nu^{i})^{r-k+1}\mid\mathbf{Z}\right]\\
&\qquad\qquad+b_{k-1}\binom{k-1}{k-1}\mathbb{E}\left[(\nu^{i})^{0}\mid\mathbf{Z}\right]\\
\Rightarrow b_{k-1}&=-\bar{b}_{r}\binom{r}{k-1}\mathbb{E}\left[(\nu^{i})^{r-k+1}\mid\mathbf{Z}\right].
\end{aligned}
}
\end{equation*}
Now, when $q=k-2$, (ref) becomes
\begin{equation*}
{
\begin{aligned}
0&=\bar{b}_{r}{r \choose k-2}\mathbb{E}\left[(\nu^{i})^{r-k+2}\mid\mathbf{Z}\right]\\
&\quad+b_{k-1}{k-1 \choose k-2}\mathbb{E}\left[(\nu^{i})^{(k-1)-(k-2)}\mid\mathbf{Z}\right]\\
&\quad+b_{k-2}\mathbb{E}\left[(\nu^{i})^{0}\mid\mathbf{Z}\right]\\
\Rightarrow b_{k-2}&=-b_{k-1}{k-1 \choose k-2}\mathbb{E}\left[(\nu^{i})^{1}\mid\mathbf{Z}\right]\\
&\quad-\bar{b}_{r}{r \choose k-2}\mathbb{E}\left[(\nu^{i})^{r-k+2}\mid\mathbf{Z}\right].
\end{aligned}
}
\end{equation*}
Now, suppose $b_{q+1},\cdots,b_{k-1}$ are known and we want to find out what $b_{q}$ is. We have to solve it from
\begin{equation*}
{
\begin{aligned}
0&=b_{q}\mathbb{E}\left[(\nu^{i})^{0}\mid\mathbf{Z}\right]\\
&\quad\quad+\bar{b}_{r}{r \choose q}\mathbb{E}\left[(\nu^{i})^{r-q}\mid\mathbf{Z}\right]+\overset{k-1}{\underset{u=q+1}{\sum}}b_{u}{u \choose q}\mathbb{E}\left[(\nu^{i})^{u-q}\mid\mathbf{Z}\right].
\end{aligned}
}
\end{equation*}
We can obtain $b_{q}$ from the above equation, which gives
\begin{equation*}
{
\begin{aligned}
b_{q}&=-\overset{k-1}{\underset{u=q+1}{\sum}}b_{u}{u \choose q}\mathbb{E}\left[(\nu^{i})^{u-q}\mid\mathbf{Z}\right]\\
&\quad\quad-\bar{b}_{r}{r \choose q}\mathbb{E}\left[(\nu^{i})^{r-q}\mid\mathbf{Z}\right]\\
\Rightarrow b_{q}&=-\overset{k-1-q}{\underset{u=1}{\sum}}b_{q+u}{q+u \choose q}\mathbb{E}\left[(\nu^{i})^{u}\mid\mathbf{Z}\right]\\
&\quad\quad-\bar{b}_{r}{r \choose q}\mathbb{E}\left[(\nu^{i})^{r-q}\mid\mathbf{Z}\right].
\end{aligned}
}
\end{equation*}
The proof is completed.
proof[Proof of Corollary (ref)]\
We have discussed the way to obtain the estimator in the main paper.
To facilitate the upcoming studies, we first introduce some notations. Recall that $Y^{i}$ is the potential outcome under the treatment $d^{i}$. We use $Y^{i;F}$ as the factual outcome if an individual receives $d^i$, and $Y^{i;CF}$ as the counterfactual outcome if an individual receives alternative treatments. Hence, we have
equation*[equation* omitted — 146 chars of source]
Based on the introduced notations, we can define two residual differences $\xi_{m}^{i;F}$ and $\xi_{m}^{i;CF}$ for the $m^{\text{th}}$ individual according to the sets $\mathscr{I}$ and $\mathscr{I}^{c}$ that the $m^{\text{th}}$ individual belongs to. Mathematically, $\xi_m^{i;F}:=Y_m^{i;F}-g^{i}(\mathbf{Z}_m)$ if $m \in \mathscr{I}$ and $\xi_m^{i;CF}:=Y_m^{i;CF}-g^{i}(\mathbf{Z}_m)$ if $m \in \mathscr{I}^c$.
We give the statistical properties between $\xi_{m}^{i;F}$ and $\xi_{\bar{m}}^{i;CF}$ for $m,\;\bar{m}$. First, $\xi_{m}^{i;F}\perp\!\!\!\perp\xi_{\bar{m}}^{i;CF}\mid \mathbf{Z}$ due to the SUTVA assumption for $m\neq\bar{m}$. Second. Regardless of the actual treatment the individual receives, the noises in terms of $d^{i}$ should be identical and independently distributed for different individuals. As a result, we should have
equation*[equation* omitted — 194 chars of source]
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$ and $m,\bar{m}$.
We give the 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}$. The properties are summarized in Proposition (ref).
propositionGiven the covariates $\mathbf{Z}$, the random variable $\xi_{m}^{i;F}$ and $\xi_{\bar{m}}^{i;CF}$ are independent and identically distributed, i.e., $\mathbb{E}\big[(\xi_{m}^{i;F})^{r}\mid D_{m},\mathbf{Z}\big]=\mathbb{E}\big[(\xi_{\bar{m}}^{i;CF})^{r} \mid D_{\bar{m}},\mathbf{Z}\big]$ and $\xi_{m}^{i;F} \perp\!\!\!\perp \xi_{\bar{m}}^{i;CF}\mid \mathbf{Z}$.
proofUsing the SUTVA assumption, we have $\xi_{m}^{i;F} \perp\!\!\!\perp \xi_{\bar{m}}^{i;CF}\mid\mathbf{Z}$. In addition, we have
\begin{equation*}
\resizebox{0.5\textwidth}{!}{$
\begin{aligned}
&\mathbb{E}\left[(\xi_{m}^{i;F})^{r} \mid D_{m}, \mathbf{Z}_{m}\right]=\mathbb{E}\left[(Y_{m}^{i;F}-g^{i}(\mathbf{Z}_{m}))^{r} \mid D_{m},\mathbf{Z}_{m}\right]\\
=&\underset{k=0}{\overset{r}{\sum}}\binom{r}{k}\mathbb{E}\left[(Y_{m}^{i;F})^{r} \mid D_{m}, \mathbf{Z}_{m}\right](-g^{i}(\mathbf{Z}_{m}))^{r-k}\\
\overset{\star}{=}&\underset{k=0}{\overset{r}{\sum}}\binom{r}{k}\mathbb{E}\left[(Y_{m}^{i;F})^{r} \mid \mathbf{Z}_{m}\right](-g^{i}(\mathbf{Z}_{m}))^{r-k}=\mathbb{E}\left[(\xi^i)^{r} \mid \mathbf{Z}\right].
\end{aligned}
$}
\end{equation*}
$\overset{\star}{=}$ is due to the ignorability assumption. Similarly, we also have
\begin{equation*}
\resizebox{0.5\textwidth}{!}{$
\begin{aligned}
\mathbb{E}\left[(\xi_{\bar{m}}^{i;CF})^{r} \mid D_{\bar{m}} \neq d^i, \mathbf{Z}_{\bar{m}}\right]=\mathbb{E}\left[(\xi_{\bar{m}}^i)^{r} \mid \mathbf{Z}_{\bar{m}}\right]=\mathbb{E}\left[(\xi^i)^{r} \mid \mathbf{Z}\right].
\end{aligned}
$}
\end{equation*}
In the remaining sequel, we investigate the consistency of our RCL estimators based on the basics of orthogonal machine learning theory. To start with, we give the assumptions on the nuisance parameters. Only the assumptions that are helpful in studying the consistency of our RCL estimators 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.
assumptionGiven that the nuisance parameters and the true nuisance parameters are $(\hat{g}^{i},\hat{\pi}^{i})$ and $(g^{i},\pi^{i})$, and $S=\{\bm{\alpha}=(\alpha_{1},\alpha_{2})\in\mathbb{Z}^{2}_{\geq 0}:\left\|\bm{\alpha}\right\|_{1}\leq k\}$, we have
\begin{enumerate}
• { $\mathbb{E}\big[\left|\hat{g}^{i}(\mathbf{Z})-g^{i}(\mathbf{Z})\right|^{4\alpha_{1}}\left|\hat{\pi}^{i}(\mathbf{Z})-\pi^{i}(\mathbf{Z})\right|^{4\alpha_{2}}\mid \hat{g}^{i},\hat{\pi}^{i}\big]\overset{p}{\rightarrow} 0$} { $\forall\;\bm{\alpha}\in S$}\\
• { $N^{\frac{1}{2}}\sqrt{\mathbb{E}\big[\left|\hat{g}^{i}(\mathbf{Z})-g^{i}(\mathbf{Z})\right|^{2\alpha_{1}}\left|\hat{\pi}^{i}(\mathbf{Z})-\pi^{i}(\mathbf{Z})\right|^{2\alpha_{2}}\mid \hat{g}^{i},\hat{\pi}^{i}\big]}\overset{p}{\rightarrow} 0$} { $\forall\;\bm{\alpha}\in \{\bm{\alpha}\in\mathbb{Z}^{2}_{\geq 0}:\left\|\bm{\alpha}\right\|_{1}\leq k+1\}\backslash S$}.
\end{enumerate}
For notational convenience, we rewrite our RCL estimator $\hat{\theta}^{i}_{RCL}$ and denote it as $\hat{\theta}^{i}_{N}$. Indeed, we have
equation[equation omitted — 479 chars of source]
Besides, we define two quantities $\hat{\tilde{\theta}}_{N}^{i}$ and $\hat{\bar{\theta}}_{N}^{i}$. They are
{
align[align omitted — 818 chars of source]
}
We also define
equation*[equation* omitted — 873 chars of source]
Then (ref) and (ref) can be rewritten as
{
align[align omitted — 681 chars of source]
}
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.
lemmaGiven two sequences of random variables $(X_{N})_{N=1}^{\infty}$ and $(Y_{N})_{N=1}^{\infty}$ such that $X_{N}\overset{d}{=}Y_{N}$. If $X_{N}\overset{p}{\rightarrow}c$ for some constant $c$, then $Y_{N}\overset{p}{\rightarrow}c$.
proofLet $f_{X_{N}}(\cdot)$ and $f_{Y_{N}}(\cdot)$ be the density functions of the random variables $X_{N}$ and $Y_{N}$ respectively. Since $X_{N}\overset{d}{=}Y_{N}$, $f_{X_{N}}(\cdot)=f_{Y_{N}}(\cdot)$. Hence, $\mathbb{P}\left\{\left|X_{N}-c\right|\geq \epsilon\right\}=\int_{|z-c|\geq\epsilon}f_{X_{N}}(z)dz=\int_{|z-c|\geq\epsilon}f_{Y_{N}}(z)dz=\mathbb{P}\left\{\left|Y_{N}-c\right|\geq \epsilon\right\}$. Consequently, $X_{N}\overset{p}{\rightarrow}c$ implies $Y_{N}\overset{p}{\rightarrow}c$.
lemmaGiven random variables $X$, $Y$, $E$, $Z$. If $(X \overset{d}{=} Y) \mid Z$, $(X \perp \!\!\! \perp E) \mid Z$, $(Y \perp \!\!\! \perp E) \mid Z$, then $Xh(E,Z) \overset{d}{=} Yh(E,Z)$ for any function $h$.
proofDefine $f_Z(z)$ as the density function of $Z$, $f_{X|Z}(x|z)$ is the conditional density function of $X|Z$, $f_{Y|Z}(y|z)$ is the conditional density function of $Y|Z$, $f_{E|Z}(e|z)$ is the conditional density function of $E|Z$, $f_{X,E|Z}(x,e|z)$ is the conditional joint density function of $X,E|Z$, and $f_{Y,E|Z}(y,e|z)$ is the conditional joint density function of $Y,E|Z$. For a measurable set $\mathcal{A}$, we have
\begin{equation*}
{
\begin{aligned}
&\mathbb{P}\{Xh(E,Z)\in \mathcal{A}\}\\
=&\int_{\Omega_{Z}}\left\{\iint_{\Omega_{X} \times \Omega_{E}} \mathbf{1}_{\{xh(e,z)\in \mathcal{A}\}}f_{X,E|Z}(x,e|z)dxde\right\}f_{Z}(z)dz\\
\overset{*}{=}&\int_{\Omega_{Z}}\left\{\iint_{\Omega_{X} \times \Omega_{E}} \mathbf{1}_{\{xh(e,z)\in \mathcal{A}\}}f_{X|Z}(x|z)f_{E|Z}(e|z)dxde\right\}f_{Z}(z)dz\\
\overset{\triangle}{=}&\int_{\Omega_{Z}}\left\{\iint_{\Omega_{Y} \times \Omega_{E}} \mathbf{1}_{\{yh(e,z)\in \mathcal{A}\}}f_{Y|Z}(y|z)f_{E|Z}(e|z)dyde\right\}f_{Z}(z)dz\\
\overset{\square}{=}&\int_{\Omega_{Z}}\left\{\iint_{\Omega_{Y} \times \Omega_{E}} \mathbf{1}_{\{yh(e,z)\in \mathcal{A}\}}f_{Y,E|Z}(y,e|z)dyde\right\}f_{Z}(z)dz\\
=&\mathbb{P}\{Yh(E,Z)\in \mathcal{A}\}
\end{aligned}
}
\end{equation*}
$*$ holds since $(X \perp \!\!\! \perp E) \mid Z$, $\triangle$ holds since $(X \overset{d}{=} Y) \mid Z$, and $\square$ holds since $(Y \perp \!\!\! \perp E) \mid Z$.
propositionSuppose { $\mathbb{E}\left[(\xi^{i;F}_{m})^{2}\mid\mathbf{Z}\right]$} and { $(A^{i}_m)^2$} exist such that $\mathbb{E}\left[(A_m^{i})^2\mathbb{E}\left[(\xi_m^{i;F})^{2}\mid\mathbf{Z}\right]\right]$ is finite for all $m$. We have { $\kappa_{N}^{i;CF}-\kappa_{N}^{i;F}\overset{p}{\rightarrow} 0$} when { $N\rightarrow \infty$}.
proof$\forall \epsilon>0$, we consider $\mathbb{P}\left\{\left|\kappa_{N}^{i;CF}-\kappa_{N}^{i;F}\right|\geq\epsilon\right\}$. Indeed, we have
\begin{equation*}
{
\begin{aligned}
&\mathbb{P}\left\{\left|\kappa_{N}^{i;CF}-\kappa_{N}^{i;F}\right|\geq\epsilon\right\}\leq \frac{\mathbb{E}\left[\left(\kappa_{N}^{i;CF}-\kappa_{N}^{i;F}\right)^{2}\right]}{\epsilon^{2}}\\
=&\frac{\frac{1}{N^{2}}\mathbb{E}\left[\left(\underset{m\in \mathscr{I}^{c}}{\sum}(\xi_{m}^{i;CF}-\xi_{m}^{i;F})A_{m}^{i}\right)^{2}\right]}{\epsilon^{2}}.\\
\end{aligned}
}
\end{equation*}
Denoting $\xi_{m}^{i;CF}-\xi_{m}^{i;F}$ as $\Xi_{m}^{i}$, we have
\begin{equation*}
{
\begin{aligned}
&\mathbb{E}\big[\big(\underset{m\in \mathscr{I}^{c}}{\sum}(\xi_{m}^{i;CF}-\xi_{m}^{i;F})A_{m}^{i}\big)^{2}\big]=\mathbb{E}\big[\big(\underset{m\in \mathscr{I}^{c}}{\sum}\Xi_{m}^{i}A_{m}^{i}\big)^{2}\big]\\
=&\mathbb{E}\big[\underset{m,\bar{m}\in \mathscr{I}^{c}}{\sum}\Xi_{m}^{i}A_{m}^{i}\Xi_{\bar{m}}^{i}A_{\bar{m}}^{i}\big]=\underset{m,\bar{m}\in \mathscr{I}^{c}}{\sum}\mathbb{E}\left[\Xi_{m}^{i}A_{m}^{i}\Xi_{\bar{m}}^{i}A_{\bar{m}}^{i}\right]\\
&\;+\underset{\substack{m,\bar{m}\in \mathscr{I}^{c} \\ m\neq \bar{m}}}{\sum}\mathbb{E}\left[A_{m}^{i}A_{\bar{m}}^{i}\mathbb{E}\left[\Xi_{m}^{i}\Xi_{\bar{m}}^{i}\mid D,\mathbf{Z}\right]\right]\\
=&\underset{m\in \mathscr{I}^{c}}{\sum}\mathbb{E}\left[(A_{m}^{i})^{2}\mathbb{E}\left[(\Xi_{m}^{i})^{2}\mid D,\mathbf{Z}\right]\right]\\
&\;+\underset{\substack{m,\bar{m}\in \mathscr{I}^{c} \\ m\neq \bar{m}}}{\sum}\mathbb{E}\left[A_{m}^{i}A_{\bar{m}}^{i}\mathbb{E}\left[\Xi_{m}^{i}\mid D,\mathbf{Z}\right]\mathbb{E}\left[\Xi_{\bar{m}}^{i}\mid D,\mathbf{Z}\right]\right]\\
=&\underset{m\in \mathscr{I}^{c}}{\sum}\mathbb{E}\left[(A_{m}^{i})^{2}\mathbb{E}\left[(\Xi_{m}^{i})^{2}\mid \mathbf{Z}\right]\right]=2\underset{m\in \mathscr{I}^{c}}{\sum}\mathbb{E}\left[(A_{m}^{i})^{2}\mathbb{E}\left[(\xi^{i;F}_{m})^{2}\mid \mathbf{Z}\right]\right]\\
\leq& 2N\mathbb{E}\left[(A^{i})^{2}\mathbb{E}\left[(\xi^{i;F})^{2}\mid \mathbf{Z}\right]\right].
\end{aligned}
}
\end{equation*}
The last equality follows from
\begin{equation*}
{
\begin{aligned}
&\mathbb{E}\big[(\Xi^{i}_{m})^{2}\mid\mathbf{Z}\big]=\mathbb{E}\big[(\xi^{i;CF}_{m}-\xi^{i;F}_{m})^{2}\mid\mathbf{Z}\big]\\
=&\mathbb{E}\big[(\xi^{i;CF}_{m})^{2}\mid\mathbf{Z}\big]-2\mathbb{E}\big[\xi^{i;F}_{m}\xi^{i;CF}_{m}\mid\mathbf{Z}\big]+\mathbb{E}\big[(\xi^{i;F}_{m})^{2}\mid\mathbf{Z}\big]\\
=&\mathbb{E}\big[(\xi^{i;CF}_{m})^{2}\mid\mathbf{Z}\big]-2\mathbb{E}\big[\xi^{i;F}_{m}\mid\mathbf{Z}\big]\mathbb{E}\big[\xi^{i;CF}_{m}\mid\mathbf{Z}\big]+\mathbb{E}\big[(\xi^{i;F}_{m})^{2}\mid\mathbf{Z}\big]\\
=&\mathbb{E}\big[(\xi^{i;CF}_{m})^{2}\mid\mathbf{Z}\big]+\mathbb{E}\big[(\xi^{i;F}_{m})^{2}\mid\mathbf{Z}\big]=2\mathbb{E}\big[(\xi^{i;F}_{m})^{2}\mid\mathbf{Z}\big].
\end{aligned}
}
\end{equation*}
As a consequence, we have
\begin{equation*}
{
\begin{aligned}
\mathbb{P}\left\{\left|\kappa_{N}^{i;CF}-\kappa_{N}^{i;F}\right|\geq\epsilon\right\}&\leq \frac{2N\mathbb{E}\left[(A^{i})^{2}\mathbb{E}\left[(\xi^{i;F})^{2}\mid\mathbf{Z}\right]\right]}{N^{2}\epsilon^{2}}\\
&=\frac{2\mathbb{E}\left[(A^{i})^{2}\mathbb{E}\left[(\xi^{i;F})^{2}\mid\mathbf{Z}\right]\right]}{N\epsilon^{2}}\rightarrow 0
\end{aligned}
}
\end{equation*}
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.
propositionSuppose 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}$} { $\forall u,\;\bar{u}\in\{1,2,\cdots,R\}$}. We have
\begin{equation*}
\kappa_{N}^{i;F}-\kappa_{R,N}^{i;F}\overset{p}{\rightarrow}0\quad when\quad N\rightarrow \infty.
\end{equation*}
proof\ \\
Write
\begin{equation*}
{
\begin{aligned}
\kappa_{R,N}^{i;F}&=\frac{1}{R}\overset{R}{\underset{u=1}{\sum}}\bigg[\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\overset{\sum}}\xi_{m,u}^{i;F} A_{m}^{i}\bigg]\\
&=\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\overset{\sum}}\left(\frac{1}{R}\overset{R}{\underset{u=1}{\sum}}\xi_{m,u}^{i;F}\right) A_{m}^{i}=\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\overset{\sum}}\mathscr{E}_{m}^{i} A_{m}^{i}.
\end{aligned}
}
\end{equation*}
$\forall\epsilon>0$, we have
\begin{equation*}
{
\begin{aligned}
\mathbb{P}\left\{\left|\kappa_{N}^{i;F}-\kappa_{R,N}^{i;F}\right|\geq\epsilon\right\}\leq \frac{\mathbb{E}\left[\left(\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\overset{\sum}}\left[\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right]A_{m}^{i}\right)^{2}\right]}{\epsilon^{2}}.
\end{aligned}
}
\end{equation*}
We simplify { $\mathbb{E}\left[\left(\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\overset{}{\sum}}\left[\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right]A_{m}^{i}\right)^{2}\right]$}. Note that
\begin{equation*}
\resizebox{0.5\textwidth}{!}{$
\begin{aligned}
&\mathbb{E}\left[\left(\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\overset{}{\sum}}\left[\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right]A_{m}^{i}\right)^{2}\right]\\
=&\frac{1}{N^{2}}\:\:\sum_{\mathclap{\substack{m,\bar{m}\in \mathscr{I}^{c}}}}\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)\left(\mathscr{E}_{\bar{m}}^{i}-\xi_{\bar{m}}^{i;F}\right)A_{\bar{m}}^{i}A_{m}^{i}\right]\\
=&\frac{1}{N^{2}}\:\:\sum_{\mathclap{\substack{m\in \mathscr{I}^{c}}}}\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)^{2}(A_{m}^{i})^{2}\right]\\
&\;+\frac{1}{N^{2}}\:\:\sum_{\mathclap{\substack{m,\bar{m}\in \mathscr{I}^{c}\\ m\neq\bar{m}}}}\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)\left(\mathscr{E}_{\bar{m}}^{i}-\xi_{\bar{m}}^{i;F}\right)A_{\bar{m}}^{i}A_{m}^{i}\right]\\
=&\frac{1}{N^{2}}\:\:\sum_{\mathclap{\substack{m\in \mathscr{I}^{c}}}}\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)^{2}(A_{m}^{i})^{2}\right]\\
&\;+\frac{1}{N^{2}}\:\:\sum_{\mathclap{\substack{m,\bar{m}\in \mathscr{I}^{c}\\ m\neq\bar{m}}}}\mathbb{E}\left[A_{\bar{m}}^{i}A_{m}^{i}\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)\left(\mathscr{E}_{\bar{m}}^{i}-\xi_{\bar{m}}^{i;F}\right)\mid D,\mathbf{Z}\right]\right]\\
=&\frac{1}{N^{2}}\:\:\sum_{\mathclap{\substack{m\in \mathscr{I}^{c}}}}\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)^{2}(A_{m}^{i})^{2}\right]\\
&\;+\frac{1}{N^{2}}\:\:\sum_{\mathclap{\substack{m,\bar{m}\in \mathscr{I}^{c}\\ m\neq\bar{m}}}}\mathbb{E}\left[A_{\bar{m}}^{i}A_{m}^{i}\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)\mid\mathbf{Z}\right]\mathbb{E}\left[\left(\mathscr{E}_{\bar{m}}^{i}-\xi_{\bar{m}}^{i;F}\right)\mid\mathbf{Z}\right]\right]\\
=&\frac{1}{N^{2}}\:\:\sum_{\mathclap{\substack{m\in \mathscr{I}^{c}}}}\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)^{2}(A_{m}^{i})^{2}\right].
\end{aligned}$
}
\end{equation*}
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_{u}^{i;F}\mid\mathbf{Z}\right]$ for any $m$ and $u\in\{1,2,\cdots,R\}$. Consequently, we have
\begin{equation*}
\resizebox{0.49\textwidth}{!}{$\begin{aligned}
&\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)\mid\mathbf{Z}\right]=\mathbb{E}\left[\left(\frac{1}{R}\overset{R}{\underset{u=1}{\sum}} \xi_{m,u}^{i;F}-\xi_{m}^{i;F}\right)\mid\mathbf{Z}\right]\\
=&\frac{1}{R}\overset{R}{\underset{u=1}{\sum}} \mathbb{E}\left[\xi_{m,u}^{i;F}\mid\mathbf{Z}\right]-\mathbb{E}\left[\xi_{m}^{i;F}\mid\mathbf{Z}\right]=\frac{1}{R}\overset{R}{\underset{u=1}{\sum}} \mathbb{E}\left[\xi_{m}^{i;F}\mid\mathbf{Z}\right]-\mathbb{E}\left[\xi_{m}^{i;F}\mid\mathbf{Z}\right]\\
=& \mathbb{E}\left[\xi_{m}^{i;F}\mid\mathbf{Z}\right]-\mathbb{E}\left[\xi_{m}^{i;F}\mid\mathbf{Z}\right]=0.
\end{aligned}
$}
\end{equation*}
In addition, we simplify the quantity { $\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)^{2}\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
\begin{equation*}
{
\begin{aligned}
&\mathbb{E}\left[\left(\overset{R}{\underset{u=1}{\sum}}\left[\xi_{m,u}^{i;F}-\xi_{m}^{i;F}\right]\right)^{2}\right]\\
=&\overset{R}{\underset{u,\bar{u}=1}{\sum}}\mathbb{E}\left[\left(\xi_{m,u}^{i;F}-\xi_{m}^{i;F}\right)\left(\xi_{m,\bar{u}}^{i;F}-\xi_{m}^{i;F}\right)\right]\\
=&\overset{R}{\underset{u,\bar{u}=1}{\sum}}\left\{\mathbb{E}\left[\xi_{m,u}^{i;F}\xi_{m,\bar{u}}^{i;F}\right]-\mathbb{E}\left[\xi_{m,u}^{i;F}\xi_{m}^{i;F}\right]\right.\\
&\quad\quad\quad\left.-\mathbb{E}\left[\xi_{m}^{i;F}\xi_{m,\bar{u}}^{i;F}\right]+\mathbb{E}\left[\xi_{m}^{i;F}\xi_{m}^{i;F}\right]\right\}\\
=&\overset{R}{\underset{u=1}{\sum}}\mathbb{E}\left[(\xi_{m,u}^{i;F})^{2}\right]-2R\overset{R}{\underset{u=1}{\sum}}\mathbb{E}\left[\xi_{m,u}^{i;F}\xi_{m}^{i;F}\right]+R^{2}\mathbb{E}\left[(\xi_{m}^{i;F})^{2}\right]\\
&\;+\overset{R}{\underset{\substack{u,\bar{u}=1 \\ u\neq\bar{u}}}{\sum}}\mathbb{E}\left[\xi_{m,u}^{i;F}\xi_{m,\bar{u}}^{i;F}\right]=\left[R^{2}+R\right]\mathbb{E}\left[\left(\xi_{m}^{i;F}\right)^{2}\right].
\end{aligned}
}
\end{equation*}
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
\begin{equation*}
\resizebox{0.5\textwidth}{!}{$
\begin{aligned}
&\mathbb{E}\left[\xi_{m,u}^{i;F}\xi_{m}^{i;F}\right]=\mathbb{E}\left[\mathbb{E}\left[\xi_{m,u}^{i;F}\xi_{m}^{i;F}\mid \mathbf{Z}\right]\right]=\mathbb{E}\left[\mathbb{E}\left[\xi_{m,u}^{i;F}\mid \mathbf{Z}\right]\mathbb{E}\left[\xi_{m}^{i;F}\mid \mathbf{Z}\right]\right]\\
=&\mathbb{E}\left[\mathbb{E}\left[\mathbb{E}\left[\xi_{m,u}^{i;F}\mid D, \mathbf{Z}\right]\mid\mathbf{Z}\right]\mathbb{E}\left[\mathbb{E}\left[\xi_{m}^{i;F}\mid D,\mathbf{Z}\right]\mid \mathbf{Z}\right]\right]=0
\end{aligned}
$}
\end{equation*}
and
\begin{equation*}
\resizebox{0.5\textwidth}{!}{$
\begin{aligned}
&\mathbb{E}\left[\xi_{m,u}^{i;F}\xi_{m,\bar{u}}^{i;F}\right]=\mathbb{E}\left[\mathbb{E}\left[\xi_{m,u}^{i;F}\xi_{m,\bar{u}}^{i;F}\mid \mathbf{Z}\right]\right]=\mathbb{E}\left[\mathbb{E}\left[\xi_{m,u}^{i;F}\mid \mathbf{Z}\right]\mathbb{E}\left[\xi_{m,\bar{u}}^{i;F}\mid \mathbf{Z}\right]\right]\\
=&\mathbb{E}\left[\mathbb{E}\left[\mathbb{E}\left[\xi_{m,u}^{i;F}\mid D,\mathbf{Z}\right]\mid \mathbf{Z}\right]\mathbb{E}\left[\mathbb{E}\left[\xi_{m,\bar{u}}^{i;F}\mid D,\mathbf{Z}\right]\mid \mathbf{Z}\right]\right]=0.
\end{aligned}
$}
\end{equation*}
Consequently, we have
\begin{equation*}
{
\begin{aligned}
\mathbb{E}\left[\left(\mathscr{E}_{m}^{i}-\xi_{m}^{i;F}\right)^{2}\right]=\left(1+\frac{1}{R}\right)\mathbb{E}\left[\left(\xi_{m}^{i;F}\right)^{2}\right].
\end{aligned}
}
\end{equation*}
Thus, we have
\begin{equation}
{
\begin{aligned}
\mathbb{P}\left\{\left|\kappa_{N}^{i;F}-\kappa_{R,N}^{i;F}\right|\geq\epsilon\right\}&\leq \frac{\frac{1}{N^{2}}\underset{m\in \mathscr{I}^{c}}{\overset{\sum}}\left(1+\frac{1}{R}\right)\mathbb{E}\left[\left(\xi_{m}^{i;F}\right)^{2}\right]}{\epsilon^{2}}\\
&\leq \frac{\left(1+\frac{1}{R}\right)\mathbb{E}\left[\left(\xi^{i;F}\right)^{2}\right]}{N\epsilon^{2}}.
\end{aligned}
}
\end{equation}
We notice that no matter we set $R\rightarrow \infty$ followed by $N\rightarrow \infty$ or vice versa, 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$}.
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$. Before presenting the proof, we notice that $\left(\hat{\xi}^{i;F}\overset{d}{=}\hat{\xi}^{i;CF}\right)\mid\mathbf{Z}$ when $\hat{g}^{i}$ and $g^{i}$ satisfie 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}$.
proof$\forall\epsilon>0$, we have
\begin{equation*}
{
\begin{aligned}
&\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\theta}_{N}^{i}-\theta^{i}\right| \geq \epsilon\right\}\\
=&\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\theta}_{N}^{i}-\hat{\bar{\theta}}_{N}^{i}+\hat{\bar{\theta}}_{N}^{i}-\theta^{i}\right| \geq \epsilon\right\}\\
\leq&\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\theta}_{N}^{i}-\hat{\bar{\theta}}_{N}^{i}\right| \geq \frac{\epsilon}{2}\right\}+\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\bar{\theta}}_{N}^{i}-\theta^{i}\right| \geq \frac{\epsilon}{2}\right\}.
\end{aligned}
}
\end{equation*}
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
\begin{equation}
{
\begin{aligned}
&\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\theta}_{N}^{i}-\hat{\bar{\theta}}_{N}^{i}\right| \geq \frac{\epsilon}{2}\right\}=\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\kappa}_{R,N}^{i;F}-\hat{\kappa}_{N}^{i;F}\right| \geq \frac{\epsilon}{2}\right\}\\
\leq&\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\kappa}_{R,N}^{i;F}-\kappa_{R,N}^{i;F}\right| \geq \frac{\epsilon}{8}\right\}}_{(a)}+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\left|\kappa_{R,N}^{i;F}-\kappa_{N}^{i;F}\right| \geq \frac{\epsilon}{8}\right\}}_{(b)}\\
&\;+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\left|\kappa_{N}^{i;F}-\kappa_{N}^{i;CF}\right| \geq \frac{\epsilon}{8}\right\}}_{(c)}+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\left|\kappa_{N}^{i;CF}-\hat{\kappa}_{N}^{i;F}\right| \geq \frac{\epsilon}{8}\right\}}_{(d)}.
\end{aligned}
}
\end{equation}
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 function { $\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
\begin{equation*}
{
\begin{aligned}
&\left|\kappa_{N}^{i;CF}-\hat{\kappa}_{N}^{i;F}\right|=\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\xi_{m}^{i;CF}A_{m}^{i}-\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}\hat{A}_{m}^{i}\right|\\
\leq&\underbrace{\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}(\xi_{m}^{i;CF}A_{m}^{i}-\xi_{m}^{i;F}A_{m}^{i})\right|}_{\Gamma_{1}}\\
&\;+\underbrace{\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}(\xi_{m}^{i;F}A_{m}^{i}-\hat{\xi}_{m}^{i;F}A_{m}^{i})\right|}_{\Gamma_{2}}+\underbrace{\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}(A_{m}^{i}-\hat{A}_{m}^{i})\right|}_{\Gamma_{3}},
\end{aligned}
}
\end{equation*}
((ref)d) is bounded above by
\begin{equation}
{
\begin{aligned}
&\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{1} \geq \frac{\epsilon}{24}\right\}}_{(a)}+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2} \geq \frac{\epsilon}{24}\right\}}_{(b)}+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{3} \geq \frac{\epsilon}{24}\right\}}_{(c)}.
\end{aligned}
}\\
\end{equation}
((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
\begin{equation*}
{
\begin{aligned}
\Gamma_{2}=&\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}(\xi_{m}^{i;F}A_{m}^{i}-\hat{\xi}_{m}^{i;F}A_{m}^{i})\right|\\
\leq&\underbrace{\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\xi_{m}^{i;F}A_{m}^{i}-\frac{N^{c}}{N}\mathbb{E}_{\hat{\rho}}\left[\xi^{i;F}A^{i}\right]\right|}_{\Gamma_{2;1}}\\
&\;+\underbrace{\left|\frac{N^{c}}{N}\mathbb{E}_{\hat{\rho}}\left[\xi^{i;F}A^{i}\right]-\frac{N^{c}}{N}\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]\right|}_{\Gamma_{2;2}}\\
&\;+\underbrace{\left|\frac{N^{c}}{N}\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]-\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}A_{m}^{i}\right|}_{\Gamma_{2;3}}.
\end{aligned}
}
\end{equation*}
We see that ((ref)b) can be further bounded by
\begin{equation}
{
\begin{aligned}
&\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2;1} \geq \frac{\epsilon}{72}\right\}}_{(a)}+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2;2} \geq \frac{\epsilon}{72}\right\}}_{(b)}+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2;3} \geq \frac{\epsilon}{72}\right\}}_{(c)}.
\end{aligned}
}
\end{equation}
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
\begin{equation*}
{
\begin{aligned}
&\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2;1}\geq\frac{\epsilon}{72}\right\}\\
=&\mathbb{P}_{\hat{\rho}}\left\{\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|\geq\frac{\epsilon}{72}\cdot\frac{N}{N^{c}}\right\}\\
\leq& \mathbb{P}_{\hat{\rho}}\left\{\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|\geq\frac{\epsilon}{72}\right\}\\
\leq&\frac{\mathbb{E}_{\hat{\rho}}\left[\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|^{2}\right]}{\left(\frac{\epsilon}{72}\right)^{2}}.
\end{aligned}
}
\end{equation*}
Consider { $\mathbb{E}_{\hat{\rho}}\left[\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|^{2}\right]$}.
Note that it equals
\begin{equation*}
\resizebox{0.5\textwidth}{!}{$
\begin{aligned}
&\frac{1}{(N^{c})^{2}}\underset{m\in \mathscr{I}^{c}}{\sum}\mathbb{E}_{\hat{\rho}}\left[\left|\xi_{m}^{i;F}A_{m}^{i}-\mathbb{E}_{\hat{\rho}}\left[\xi^{i;F}A^{i}\right]\right|^{2}\right]\\
&\;+\frac{1}{\left(N^{c}\right)^{2}}\underset{\substack{m,\bar{m}\in \mathscr{I}^{c}\\m\neq\bar{m}}}{\sum}\mathbb{E}_{\hat{\rho}}\left[(\xi_{m}^{i;F}A_{m}^{i}-\mathbb{E}_{\hat{\rho}}[\xi^{i;F}A^{i}])(\xi_{\bar{m}}^{i;F}A_{\bar{m}}^{i}-\mathbb{E}_{\hat{\rho}}[\xi^{i;F}A^{i}])\right]\\
=&\frac{1}{\left(N^{c}\right)^{2}}\underset{m\in \mathscr{I}^{c}}{\sum}\mathbb{E}_{\hat{\rho}}\left[\left(A_{m}^{i}\right)^{2}\mathbb{E}_{\hat{\rho}}\left[\left(\xi_{m}^{i;F}\right)^{2}\mid D,\mathbf{Z}\right]\right]\\
&\;+\frac{1}{\left(N^{c}\right)^{2}}\underset{\substack{m,\bar{m}\in \mathscr{I}^{c}\\m\neq\bar{m}}}{\sum}\mathbb{E}_{\hat{\rho}}\left[A_{m}^{i}A_{\bar{m}}^{i}\mathbb{E}_{\hat{\rho}}\left[\xi_{m}^{i;F}\mid D,\mathbf{Z}\right]\mathbb{E}_{\hat{\rho}}\left[\xi_{\bar{m}}^{i;F}\mid D,\mathbf{Z}\right]\right]\\
=&\frac{1}{\left(N^{c}\right)^{2}}\underset{m\in \mathscr{I}^{c}}{\sum}\mathbb{E}_{\hat{\rho}}\left[\left(A_{m}^{i}\right)^{2}\mathbb{E}_{\hat{\rho}}\left[\left(\xi_{m}^{i;F}\right)^{2}\mid\mathbf{Z}\right]\right]\\
=&\frac{1}{N^{c}}\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].
\end{aligned}
$}
\end{equation*}
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
\begin{equation*}
{
\begin{aligned}
\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2;1}\geq\frac{\epsilon}{72}\right\}\leq\frac{72^2\;\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]}{\epsilon^{2}N^{c}}\overset{p}{\rightarrow}0.
\end{aligned}
}
\end{equation*}
Now, we consider ((ref)b). Indeed, we have
\begin{equation*}
{
\begin{aligned}
&\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2;2}\geq\frac{\epsilon}{72}\right\}\\
=&\mathbb{P}_{\hat{\rho}}\left\{\left|\mathbb{E}_{\hat{\rho}}\left[\xi^{i;F}A^{i}\right]-\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]\right|\geq\frac{\epsilon}{72}\cdot\frac{N}{N^{c}}\right\}\\
\leq&\mathbb{P}_{\hat{\rho}}\left\{\left|\mathbb{E}_{\hat{\rho}}\left[(\xi^{i;F}-\hat{\xi}^{i;F})A^{i}\right]\right|\geq\frac{\epsilon}{72}\right\}\\
\leq&\frac{72^2\left\{\mathbb{E}_{\hat{\rho}}\left[(\xi^{i;F}-\hat{\xi}^{i;F})A^{i}\right]\right\}^{2}}{\epsilon^{2}}\\
\leq&\frac{72^2\;\left\{\mathbb{E}_{\hat{\rho}}\left[(\xi^{i;F}-\hat{\xi}^{i;F})^{4q}\right]\right\}^{\frac{1}{2q}}\left\{\mathbb{E}_{\hat{\rho}}\left[(A^{i})^{\frac{4q}{4q-1}}\right]\right\}^{2-\frac{1}{2q}}}{\epsilon^{2}}\overset{p}{\rightarrow}0.
\end{aligned}
}
\end{equation*}
Here, the last inequality follows from the H\"{o}lders inequality, while 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
\begin{equation*}
{
\begin{aligned}
\Gamma_{2;3}=\frac{N^{c}}{N}\left|\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]-\frac{1}{N^{c}}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}A_{m}^{i}\right|.
\end{aligned}
}
\end{equation*}
Now, we have
\begin{equation*}
{
\begin{aligned}
&\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{2;3}\geq\frac{\epsilon}{72}\right\}\\
\leq&\mathbb{P}_{\hat{\rho}}\left\{\left|\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]-\frac{1}{N^{c}}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}A_{m}^{i}\right|\geq\frac{\epsilon}{72}\right\}\\
\leq&\frac{72^2\;\mathbb{E}_{\hat{\rho}}\left[\left\{\underset{m\in \mathscr{I}^{c}}{\sum}\left(\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]-\hat{\xi}_{m}^{i;F}A_{m}^{i}\right)\right\}^{2}\right]}{\epsilon^{2}\left(N^{c}\right)^{2}}\\
=&\frac{72^2\;\underset{m\in \mathscr{I}^{c}}{\sum}\;\mathbb{E}_{\hat{\rho}}\left[\left(\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]-\hat{\xi}_{m}^{i;F}A_{m}^{i}\right)^{2}\right]}{\epsilon^{2}\left(N^{c}\right)^{2}}\\
&\;+\frac{72^2\;\underset{\substack{m,\bar{m}\in \mathscr{I}^{c}\\m\neq\bar{m}}}{\sum}\mathbb{E}_{\hat{\rho}}[(\hat{\xi}_{m}^{i;F}-\xi_{m}^{i;F})A_{m}^{i}]\mathbb{E}_{\hat{\rho}}[(\hat{\xi}_{\bar{m}}^{i;F}-\xi_{\bar{m}}^{i;F})A_{\bar{m}}^{i}]}{\epsilon^{2}\left(N^{c}\right)^{2}}\\
&\;-2\frac{72^2\;\underset{\substack{m,\bar{m}\in \mathscr{I}^{c}\\m<\bar{m}}}{\sum}\mathbb{E}_{\hat{\rho}}[(\hat{\xi}_{m}^{i;F}-\xi_{m}^{i;F})A_{m}^{i}]\mathbb{E}_{\hat{\rho}}[(\hat{\xi}^{i;F}-\xi^{i;F})A^{i}]}{\epsilon^{2}\left(N^{c}\right)^{2}}\\
&\;+\frac{72^2\;\underset{\substack{m,\bar{m}\in \mathscr{I}^{c}\\m\neq\bar{m}}}{\sum}\left\{\mathbb{E}_{\hat{\rho}}\left[\left(\hat{\xi}^{i;F}-\xi^{i;F}\right)A^{i}\right]\right\}^{2}}{\epsilon^{2}\left(N^{c}\right)^{2}}.
\end{aligned}
}
\end{equation*}
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
\begin{equation*}
{
\begin{aligned}
&\Gamma_{3}=\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}A_{m}^{i}-\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}\hat{A}_{m}^{i}\right|\\
\leq&\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}A_{m}^{i}-\frac{N^{c}}{N}\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]\right|\\
&\;+\left|\frac{N^{c}}{N}\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]-\frac{N^{c}}{N}\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}\hat{A}^{i}\right]\right|\\
&\;+\left|\frac{N^{c}}{N}\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}\hat{A}^{i}\right]-\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}\hat{A}_{m}^{i}\right|\\
=&\frac{N^{c}}{N}\underbrace{\left|\frac{1}{N^{c}}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}A_{m}^{i}-\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]\right|}_{\Gamma_{3;1}}\\
&\;+\frac{N^{c}}{N}\underbrace{\left|\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}A^{i}\right]-\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}\hat{A}^{i}\right]\right|}_{\Gamma_{3;2}}\\
&\;+\frac{N^{c}}{N}\underbrace{\left|\mathbb{E}_{\hat{\rho}}\left[\hat{\xi}^{i;F}\hat{A}^{i}\right]-\frac{1}{N^{c}}\underset{m\in \mathscr{I}^{c}}{\sum}\hat{\xi}_{m}^{i;F}\hat{A}_{m}^{i}\right|}_{\Gamma_{3;3}},
\end{aligned}
}
\end{equation*}
we see that ((ref)c) can be further bounded by
\begin{equation}
{
\begin{aligned}
&\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{3;1} \geq \frac{\epsilon}{72}\right\}}_{(a)}+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{3;2} \geq \frac{\epsilon}{72}\right\}}_{(b)}+\underbrace{\mathbb{P}_{\hat{\rho}}\left\{\Gamma_{3;3} \geq \frac{\epsilon}{72}\right\}}_{(c)}.
\end{aligned}
}
\end{equation}
Similarly, we can prove that ((ref)a) and ((ref)c) converge to $0$ in probability when $N\rightarrow \infty$ using the arguments in proving that ((ref)a) 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
\begin{subequations}
{
\begin{align}
&\mathbb{P}_{\hat{\rho}}\left\{\left|\hat{\kappa}_{R,N}^{i;F}-\kappa_{R,N}^{i;F}\right|\geq\frac{\epsilon}{8}\right\}\nonumber\\
\leq&\mathbb{P}_{\hat{\rho}}\left\{\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\overset{\sum}}\frac{1}{R}\underset{u=1}{\overset{R}{\sum}}\left(\hat{\xi}_{m,u}^{i;F}\right)(\hat{A}_{m}^{i}-A_{m}^{i})\right|\geq\frac{\epsilon}{16}\right\}\\
&\;+\mathbb{P}_{\hat{\rho}}\left\{\left|\frac{1}{N}\underset{m\in \mathscr{I}^{c}}{\overset{\sum}}A_{m}^{i}\frac{1}{R}\underset{u=1}{\overset{R}{\sum}}\left[(\hat{\xi}_{m,u}^{i;F}-\xi_{m,u}^{i;F})\right]\right|\geq\frac{\epsilon}{16}\right\}.
\end{align}
}
\end{subequations}
We can argue that (ref) converges to $0$ in probability as $N\rightarrow \infty$ using similar arguments when we prove that ((ref)b) 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)c) 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.