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.
97,119 characters · 19 sections · 57 citation commands
Orthogonal Series Estimation for the Ratio of Conditional Expectation Functions
In various fields of data science, researchers often face problems of estimating the ratio of conditional expectation functions (CEFR). Although CEFR appears not only in causal inference studies, there are several important examples of CEFR in the treatment effect estimation literature. When an outcome of interest is the relapse rate of a specific disease, it is more natural to consider the ratio of the conditional expectation of potential outcomes $E[Y_1|X]/E[Y_0|X]$, rather than the difference $E[Y_1|X]-E[Y_0|X]$, as a measure of causal effects of a medical treatment, where $Y_1$ and $Y_0$ are the potential outcome with and without a treatment, respectively, and $X$ is a vector of baseline covariates. Likewise, ratio-based treatment effects such as the odds ratio and hazard ratio have been widely used especially in clinical settings. Furthermore, in a data combination setting where an outcome and a treatment status are only separately observed, both conditional average treatment effect (CATE) and local average treatment effect (LATE) are identified as in the form of CEFR yamane,shinoda2022estimation, while these effects are defined as the difference of the potential outcomes.
In this article, we start by developing a novel series estimator for CEFR in a very simple setting without selection bias in observed data. This series estimator is itself useful in estimating treatment effects when data can be collected completely at random from the population of interest, but such randomized data are often not available in practice. Technically, when there is selection bias in collected data, we need to estimate potentially infinite-dimensional nuisance parameters to adjust for the bias, but these parameters may be hard to estimate with a “sufficiently high quality” in observational studies on complex systems since they can be very high-dimensional and/or highly nonlinear. The highly complex nuisance parameters do not satisfy the traditional assumptions that limit the complexity of a function class, and therefore the resulting semiparametric estimator fails to be $\sqrt{N}$-consistent. We employ debiased machine learning (DML), a set of techniques to enable the use of flexible machine learning (ML) methods for nuisance estimation, to develop a simple and general framework for constructing a high-quality estimator for CEFR even in the presence of selection bias in observed data.
The major contribution of this study is the development of a novel inference framework for CEFR with theoretical guarantees. We derive the general asymptotic results for estimation and uniform inference on the best linear approximation to the target CEFR under the proposed framework, including the validity of the Gaussian bootstrap. It is worth noting that we do not have to add stronger regularity conditions than the assumptions previously known in the literature to establish the theoretical results in this article. In addition to the general results, we provide a set of low-level sufficient conditions to apply the proposed framework to several specific settings. Besides the asymptotic analysis, we conduct numerical simulations to evaluate the performance of the proposed method on finite samples. We also illustrate the use of the framework in an empirical example.
This study builds upon three important bodies of research within the semiparametric literature: DML, series estimation and CEFR estimation. DML chernozhukov2018,chernozhukov2022locally enables inference on the finite-dimensional parameter in the presence of infinite-dimensional nuisance parameters by using the Neyman orthogonal moment conditions. When the moment condition satisfies the Neyman orthogonality neyman1959optimal, bias in the nuisance estimates has no first-order asymptotic effect on the estimator of the target parameter. DML is powerful enough to allow the use of a broad class of ML estimators ---such as $\ell_1$-penalized methods in sparse models bickel2009simultaneous,buhlmann2011statistics,belloni2011square,belloni2012sparse,belloni2011ell1,belloni2013least, $L_2$-boosting methods in sparse linear models luo2016high, neural nets schmidt2020nonparametric,farrell2021deep,kohler2021on, trees and random forest wager2015adaptive,syrgkanis2020estimation--- under a broad range of data generating processes and for various causal parameters. Moreover, the extensions of DML have been proposed to estimate a function, e.g. CATE or continuous treatment effect, with the existence of complex infinite-dimensional nuisance parameters jacob2019group,zimmert2019nonparametric,colangelo2020double,semenova2021debiased,fan2022estimation. The present study provides a novel method for the orthogonal estimation under a variety of realistic settings that cannot be covered by the previous works, such as ratio-based treatment effects, and treatment effects in the data combination settings.
The second foundation of this study is the least squares series estimation. Series estimation is a type of nonparametric estimation method that approximates a function of interest by a linear combination of multiple basis functions. It is especially useful when the exact functional form of the target function is unknown. Its asymptotic properties have been investigated intensively in the literature van1990estimating,andrews1991asymptotic,eastwood1991adaptive,gallant1991asymptotic,newey1997convergence,van2002m,huang2003local,chen2007large,cattaneo2013optimal,belloni2015some,chen2015optimal, among which we will mainly rely on the results in belloni2015some. This study contributes to the series estimation literature by establishing the asymptotic theory for the CEFR problems rather than the standard regression setting.
This study is also closely related to CEFR estimation in causal inference. A simple yet important example is ratio-based treatment effects such as the odds and hazard ratio. To the best of our knowledge, the previous works on the estimation of the ratio-based treatment effects all impose assumptions on the functional form somewhere in the model dukes2018note,liang2020relative,yadlowsky2021estimation,lee2022survival. For example, liang2020relative supposes that the ratio-based treatment effects can be expressed as the monotone single index model, and yadlowsky2021estimation imposes a stronger condition that the monotone link function is an exponential function. On the other hand, this study considers a fully nonparametric model for treatment effects, imposing no assumptions on the functional form of CEFR.
Furthermore, yamane and shinoda2022estimation show that difference-based CATE and LATE are identified in the form of CEFR in the data combination setting where we cannot observe an outcome and a treatment status simultaneously in a single dataset, and develop estimation methods for CEFR. This study accommodates the generalized version of their problem settings, where each separate dataset contains selection bias. Their estimators cannot handle the selection bias without introducing additional nuisance parameters, but their methods are not orthogonal to the nuisance parameters.
The rest of this article is organized as follows. In Section (ref), we develop a simple series estimator for CEFR. We propose the general framework for estimation and inference of CEFR in Section (ref). The main theoretical results for the asymptotic properties of the series estimator under the proposed framework are presented in Section (ref), and the application of the framework to some specific treatment effects problems is shown in Section (ref). In Section (ref), we illustrate the finite sample performance of the proposed method by numerical simulations. We also apply the proposed method to the empirical example of estimating a causal effect of participation in 401(k) on net financial assets in Section (ref). Finally, we provide further discussion on the limitation and future direction of the present study in Section (ref).
Consider the following model:
where $X\in\mathcal{X}\subset\mathbb{R}^q$ is a $q$-dimensional vector of covariates, and $\zeta_0(x)\neq0$ for all $x\in\mathcal{X}$ is necessary for $\theta_0(x)$ to be well-defined. Now we are interested in estimating $\theta_0(x)$ from separately observed iid samples: $\{u_i,x_i\}_{i=1}^{N_U}$ of $(U,X)$ and $\{t_i,x_i\}_{i=1}^{N_T}$ of $(T,X)$. Suppose $N=N_U=N_T$ without loss of generality. Also, assume that the samples $\{u_i,x_i\}_{i=1}^N$ and $\{t_i,x_i\}_{i=1}^N$ are complete random draws from the population of interest, i.e. there is no selection bias in these separate samples. Then, regarding $U$ and $T$ as outcome variables with and without a treatment, respectively, this problem coincides with the estimation of ratio-based treatment effects ogburn2015,dukes2018note,yadlowsky2021estimation in randomized control trials. Moreover, estimation of difference-based CATE and LATE from separately observed samples yamane,shinoda2022estimation is included in the model ((ref)). Although the estimator developed below primarily aims at estimating $\theta_0(x)$ from the separate samples, it can also be applicable to the estimation from joint samples $\{u_i,t_i,x_i\}_{i=1}^N$ of $(U,T,X)$. Thereafter, we omit arguments of functions if they are obvious from the context.
It follows from the model ((ref)) that $E[U-\zeta_0(X)\theta_0(X)|X]=0$. Then, using this moment condition, we can consider the linear approximation to the target function $\theta_0(x)$ by a vector of $k$ basis functions $p(x):=(p_1(x),\ldots,p_k(x))'$ as $\theta(x)=p(x)'\beta$, where
However, as the true denominator function $\zeta_0(x)$ is unknown, the above equation cannot be directly applied for estimation. To deal with this issue, shinoda2022estimation proposed to plug-in the estimator $\hat{\zeta}(x)$ of $\zeta_0(x)$, while yamane introduced an auxiliary function and reformulate the problem into minimax optimization. Both of these approaches have problems: the former lacks the theoretical justification to use flexible nonparametric methods for $\hat{\zeta}(x)$ since simply plugging-in nonparametric estimates generally does not ensure $\sqrt{N}$-consistency of the target estimator, and the latter suffers from unreliable hyperparameter selection caused by the minimax objective. In what follows, we develop a simpler one-step estimator for CEFR in the spirit of vapnik1999nature's principle: When solving a given problem, try to avoid solving a more general problem as an intermediate step.
Define $f(x):=\zeta_0(x)p(x)$ and its series estimator
Substituting $f$ for ((ref)) gives
Then, by replacing $f(x)$ with its estimator $\hat{f}(x)$, we obtain
and its sample analog estimator $\hat{\beta}=\hat{Q}^{-1}E_N[p(x_i)u_i]$, where $E_N[g(x_i)]:=N^{-1}\sum_{i=1}^N f(x_i)$ for some function $g$ and $\hat{Q}:=E_N[p(x_i)p(x_i)'t_i]$. The key observation here is that we can compute $\hat{\beta}$ even when $\{u_i,x_i\}_{i=1}^N$ and $\{t_i,x_i\}_{i=1}^N$ are only separately observed. Henceforth, we refer to this estimator as the Direct Series Ratio (DSR) estimator.
In this section, we propose a general framework for the CEFR problems by using DSR developed in the previous section as one of the building blocks. The proposed framework is based on the generalized version of the model ((ref)), and therefore it can accommodate a variety of real-world applications.
In observational studies of complex systems, possibly high-dimensional covariates $X$ may be necessary for conditional independence to hold, which is one of the critical conditions in many causal inference problems. However, researchers are often not interested in estimating treatment effects as a function of all covariates, but they need to focus on the heterogeneous relationships between treatment effects and only a few important factors. In such a case, estimating the target function on the full vector $X$ is unnecessarily hard due to the curse of dimensionality, and directly estimating a function of a subvector can stabilize the estimation. Therefore, suppose that our parameter of interest $\theta_0$ is now a function of $V\in\mathcal{V}$, a subvector of covariates $X$. Also, suppose that $\nu_0$ and $\zeta_0$ in ((ref)) are defined as
where the random variables $U$ and $T$ depend on a vector of observed random variables $O$ and a set of infinite-dimensional nuisance parameters $\eta_0(x)$. Note that $\eta_0$ is a function of $X$, reflecting that $X$ is indispensable for adjusting for selection biases. We refer to $U$ and $T$ as signals, following semenova2021debiased, and use the notations $U_0:=U(O,\eta_0)$ and $T_0:=T(O,\eta_0)$.
Among many possible choices of the signals $U$ and $T$ that satisfy ((ref)), we focus on the signals with the Neyman orthogonal property, which is defined later. The Neyman orthogonality is the key property to deliver high-quality inference on the target function $\theta_0(v)$ even when modern ML estimators that do not satisfy Donsker conditions are used for the nuisance functions $\eta_0$. The classic semiparametric theories ensure that the target estimator achieves $\sqrt{N}$-consistency by bounding the complexity of the nuisance space, but estimators in such a restricted class are not suitable for fitting high-dimensional and/or highly nonlinear nuisance functions in data-rich environments. On the other hand, modern ML estimators are so flexible to fit any function that there is no need to specify the functional form of the nuisance parameters in advance. ML estimators also perform well in high-dimensional settings by employing regularization to reduce variance at the cost of regularization bias. These properties are especially important in observational studies in which the distribution of an outcome is a complex mixture of treated and non-treated populations, and there are many confounding factors.
Although ML estimators are effective in predicting the values of nuisance functions, their prediction is inherently biased by model selection and regularization. Therefore, the naive plug-in estimator of the target function fails to be $\sqrt{N}$-consistent. To ensure the desirable theoretical properties of the target estimator while using flexible ML estimators for nuisance functions, we need signals that are locally insensitive to the bias of the nuisance estimators. Formally, the Neyman orthogonality is defined in terms of the pathwise (Gateaux) derivative as follows:
for all $v\in\mathcal{V}$, where $\partial$ is the partial derivative operator and $s\in[0,1]$ is a constant.
We describe some empirically relevant examples for the generalized setting. The detailed discussions on conditions for identification and inference of estimands in each example are given in Section (ref).
We develop a two-stage estimator using the orthogonal signals and DSR proposed in Section (ref). We refer to this two-stage estimator as the Orthogonal Series Ratio (OSR) estimator. The first stage consists of constructing the orthogonal signals by cross-fitting. Cross-fitting is another important technique to eliminate biases in the orthogonal signals constructed from finite samples. The cross-fitting of the orthogonal signals is implemented as follows:
Cross-fitting uses different samples for estimating the nuisance functions and constructing the orthogonal signals. Such sample-splitting allows the nuisance estimators to be treated as non-random when constructing the orthogonal signals $\hat{U}$ and $\hat{T}$, which helps the debiasing of the target estimator.
In the second stage, the estimated orthogonal signals $\hat{U}$ and $\hat{T}$ are plugged-in to the DSR estimator:
where we redefine $Q$ and $\hat{Q}$ as $Q=E[pp'T_0]$ and $\hat{Q}=E_N[p_ip_i'\hat{T}_i]$, respectively. The asymptotic covariance matrix of the OSR estimator is
where $\varepsilon_U:=U_0-\nu_0(V)$ and $\varepsilon_T:=T_0-\zeta_0(V)$ are the stochastic errors. The sample analogue of the asymptotic variance is
where $\hat{\theta}(v)=p(v)'\hat{\beta}$ is the OSR estimator of the target function $\theta_0$.
For statistical inference, denote the standard deviation of the OSR estimator at $v$ and its sample analogue as $\sigma(v)=\sqrt{p(v)'\Omega p(v)}$ and $\hat{\sigma}(v)=\sqrt{p(v)'\hat{\Omega}p(v)}$, respectively. Then, we can write $t$-statistic as
and the bootstrapped $t$-statistic as
where $\mathcal{N}_k^b$ is a bootstrap draw from $N(0,I_k)$. We will show that the Gaussian bootstrap is valid for the OSR estimator in the next section. We can calculate the confidence bands for $\theta_0(v)$ as
where the critical value $c_N(1-\delta)$ is the $(1-\delta)$-quantile of $N(0,1)$ for the pointwise bands, and the $(1-\delta)$-quantile of $\sup_{v\in\mathcal{V}}|\hat{\tau}_N^b(v)|$ for the uniform bands.
Here, we present the main theoretical results on the asymptotic properties of the proposed OSR estimator. These results heavily rely on the theoretical analyses in belloni2015some and semenova2021debiased.
First, we set up some additional notations. Define a random variable that takes the value of the stochastic errors with the larger absolute value:
and the lower and upper bounds on its second moments:
We denote the best linear approximation to the target function $\theta_0(v)$ by $\theta_k(v):=p(v)'\beta_k$, where
and the approximation error by $r(v):=\theta_0(v)-\theta_k(v)$. We use the same notation $\|\cdot\|$ for the $\ell_2$-norm of a vector and the operator norm of a matrix. The notation $a\lesssim b$ is used when $a\leq cb$ for some positive constant $c$ which does not depend on $N$, and $a\lesssim_P b$ is used when $a=O_P(b)$. $a\land b$ and $a\lor b$ mean $\min\{a,b\}$ and $\max\{a,b\}$, respectively. The scaled and demeaned sample average for some function $f$ is denoted by
The following assumptions are the collection of regularity conditions on the covariates distribution, basis functions, error terms, nuisance estimators, and the target function. We use these assumptions to establish the pointwise asymptotic theory for the orthogonal estimator.
The uniform boundedness of the functions $\nu_0$ and $\zeta_0$ is anyway included in the low-level conditions stated in Section (ref). Assumption (ref) therefore imposes virtually no extra restriction on these functions.
The following is the result on the pointwise convergence rate and linearization, which extends the results obtained in the standard regression setting belloni2015some,semenova2021debiased.
Lemma (ref) states that the OSR converges to the pseudo-target function at the same rate as in the standard regression setting under mild conditions. However, the linearization result is slightly different from one obtained in semenova2021debiased, where the linearization is possible even when the term $p_ir_i$ is included because $E[pr]=0$. In the CEFR problems, we have the term $p_i\zeta_{0i}r_i$ instead of $p_ir_i$, and $E[p\zeta_0r]\neq0$ in general.
The following theorem establishes the pointwise normality of the OSR estimator with additional conditions to satisfy Lindeberg's condition for the central limit theorem.
Stronger conditions than needed for the pointwise results are required to establish the uniform asymptotic theory for the OSR estimator. The following conditions control the behaviour of the stochastic errors, basis functions and nuisance errors more strictly.
Denote by $\alpha(v):=p(v)/\|p(v)\|$ the normalized value of the basis $p(v)$. Define the Lipschitz constant for $\alpha(v)$ as:
The following lemma is about the uniform convergence rate and linearization of the OSR estimator.
Likewise in Lemma (ref), we cannot include the term $p_i\zeta_{0i}r_i$ in the linearization result. However, its impact is negligible when the basis is sufficiently rich so that $l_kr_k=o(\sqrt{\log N})$.
The following theorem establishes a strong approximation of the OSR estimator's series process.
Theorem (ref) derives the convergence rate of the covariance matrix estimator $\hat{\Omega}$.
Theorem (ref) establishes the validity of Gaussian bootstrap.
Theorem (ref) is on the validity of the uniform confidence bands and their width.
We apply the general theoretical results presented in the previous section to the examples in Section (ref).
Consider the setting of Example (ref), and recall that the parameter of interest is $\theta_0(v)=E[Y_1-Y_0|D_1>D_0,V=v]$. Below, we provide the identification assumptions for LATE, and the orthogonal signals in this setting.
Assumption (ref).1 states that $Z$ is randomly assigned within a subpopulation sharing the same level of the covariates. Assumption (ref).2 restricts the treatment assignment probability from taking extreme values. Assumption (ref).3 states that the treatment assignment has non-zero effects on the actual treatment status. Assumption (ref).4 excludes defiers, those who never follow the given treatment assignment, from our analysis. Assumption (ref).5 relates the potential variables to their realized counterparts.
Consider the following doubly robust signals:
where $\rho_0(x)=P(Z=1|X=x)=E[Z|X=x]$ is a propensity score. The following theorem establishes the identification of LATE and the Neyman orthogonality of the signals.
Then, we give a set of sufficient conditions the nuisance estimators must satisfy so that the general results in Section (ref) hold. Given the true nuisance functions $\eta_0=\{\mu_0,\pi_0,\rho_0\}$ and sequences of shrinking neighborhoods $\mathcal{S}_N^\mu$ of $\mu_0$, $\mathcal{S}_N^\pi$ of $\pi_0$, $\mathcal{S}_N^\rho$ of $\rho_0$, define the following rates:
where $q\geq2$ or $q=\infty$.
Consider the setting of Example (ref), and recall that the parameter of interest is $\theta_0(v)=E[Y_1|V=v]/E[Y_0|V=v]$. Below, we provide the identification assumptions for ratio-based CATE, and the orthogonal signals in this setting.
Assumption (ref).1 suppose there are no unobserved confounders. Assumption (ref).2 restricts the treatment probability from taking extreme values. Assumption (ref).3 is a necessary condition to ensure that ratio-based CATE is well-defined uniformly over $\mathcal{V}$. Assumption (ref).4 relates the potential outcomes to the observed outcome. Note that when we are interested in the odds ratio or hazard ratio, CV based on the criterion ((ref)) is possible because the outcome $Y$ is positive in these cases.
Consider the following doubly robust signals:
where $\pi_0(x)=E[D|X=x]$ is a propensity score. The following theorem establishes the identification of ratio-based CATE and the Neyman orthogonality of the signals.
Then, we give a set of sufficient conditions the nuisance estimators must satisfy so that the general results in Section (ref) hold. Given the true nuisance functions $\eta_0=\{\mu_0,\pi_0\}$ and sequences of shrinking neighborhoods $\mathcal{S}_N^\mu$ of $\mu_0$, $\mathcal{S}_N^\pi$ of $\pi_0$, define the following rates:
where $q\geq2$ or $q=\infty$.
Consider the setting of Example (ref), and recall that the parameter of interest is $\theta_0(v)=E[Y_1-Y_0|V=v]$. Below, we provide the identification assumptions for CATE in the IDID setting, and the orthogonal signals in this setting.
Assumption (ref).1 relates the potential variables to their realized counterparts. Assumption (ref).2 restricts the propensity score from taking extreme values. Assumption (ref).3 states that the distributions of the potential variables remain the same over time when conditioned on $Z$ and $X$. Assumption (ref).4, 5 and 6 are weaker than the usual IV assumptions, where $Z$ is correlated with $D$ but independent of every potential variable. Instead, IDID only assumes $Z$ is a valid IV for the difference of the potential variables. Assumption (ref) states that the magnitude of a treatment effect does not change over time.
Consider the following signals
where $\rho_0(w,z,x)=P(W=w,Z=z|X=x)$ and $1_A$ is an indicator function that takes 1 if $A$ is true. Note that estimating $\rho$ can be done by constructing a four-class classifier and then obtaining posterior probabilities. The following theorem establishes the identification of CATE in IDID and the Neyman orthogonality of the signals.
Then, we give a set of sufficient conditions the nuisance estimators must satisfy so that the general results in Section (ref) hold. Given the true nuisance functions $\eta_0=\{\mu_0,\pi_0,\rho_0\}$ and sequences of shrinking neighborhoods $\mathcal{S}_N^\mu$ of $\mu_0$, $\mathcal{S}_N^\pi$ of $\pi_0$, $\mathcal{S}_N^\rho$ of $\rho_0$, define the following rates:
where $q\geq2$ or $q=\infty$.
Consider the setting of Example (ref), and recall that the parameter of interest is CATE $\theta_0(v)=E[Y_1-Y_0|V=v]$ or LATE $\theta_0(v)=E[Y_1-Y_0|D_1>D_0,V=v]$. Below, we provide the identification assumptions for CATE and LATE in the data combination setting, and the orthogonal signals in this setting.
Assumption (ref).1 states that $H$ and $W$ are randomly assigned within the subpopulation sharing the same level of covariates. Likewise, Assumption (ref).2 states that $Z$ is randomly assigned within strata defined by $X$. Assumption (ref).3 can be automatically satisfied as long as we have four different datasets. Assumption (ref).4 is the key assumption in this setting, which states that the two treatment regimes are different, and the treatment assignment has some impact on the individual's decision to receive a treatment or not. Assumption (ref).5 relates the potential variables to their realized counterparts. We need $P(D_1\geq D_0|V=v)=1$ for $v\in\mathcal{V}$ when the parameter of interest is LATE, and $P(D_1>D_0|V=v)=1$ for $v\in\mathcal{V}$ when the parameter of interest is CATE.
Consider the following signals
where $\rho_0(h,w,x)=P(H=h,W=w|X=x)$. The same comment for $\rho_0$ in Example (ref) also applies here. The following theorem establishes the identification of CATE and LATE in the data combination setting and the Neyman orthogonality of the signals.
Then, we give a set of sufficient conditions the nuisance estimators must satisfy so that the general results in Section (ref) hold. Given the true nuisance functions $\eta_0=\{\mu_0,\pi_0,\rho_0\}$ and sequences of shrinking neighborhoods $\mathcal{S}_N^\mu$ of $\mu_0$, $\mathcal{S}_N^\pi$ of $\pi_0$, $\mathcal{S}_N^\rho$ of $\rho_0$, define the following rates:
where $q\geq2$ or $q=\infty$.
In this section, we conduct numerical simulations to evaluate the finite-sample performance of the proposed method. We first compare DSR developed in Section (ref) to the existing methods, and then illustrate the performance of OSR and its uniform confidence band.
Setting. In this simulation study, we focus on CATE estimation by data combination explained in Example (ref). We generate random samples of $O=(HY+(1-H)D,W,H,X)$ with a sample size $N=500,1000,2000$ using the following two data generating processes (DGPs):
where $X_1,X_2,\epsilon$ are mutually independent random variables drawn from $N(0,1)$, and $\varsigma(x)=(1+e^{-x})^{-1}$ is the logistic function. Note that values of $H$ and $W$ are determined completely at random independent of $X$ so that DSR is applicable, which is the original setting considered in yamane. Also, we consider the case where no intervention takes place in regime 0 ($P(D_0=1|X)=0$) since it enables the CV method in Remark (ref) by ensuring the denominator is strictly positive. See Remark (ref) for this point. The number of replication is set $1000$.
We consider the following candidate basis vectors:
and regularization parameters $\lambda\in\{0.001,0.01,0.1,1\}$. These hyperparameters are selected based on 5-fold CV in each replication using the criterion ((ref)). For comparison with DSR, we consider a naive separate estimation (SEP), the Direct Least Squares (DLS) proposed in yamane and the Directly Weighted Least Squares (DWLS) proposed in shinoda2022estimation. The estimation of $\pi_0$, when necessary, is implemented by the $\ell_2$-regularized logistic regression with a regularization parameter $\lambda$.
\noindentResults. Table (ref) summarises the hyperparameters most frequently selected by 5-fold CV. SEP1 and SEP2 in the table represent the estimation of the numerator and denominator for the SEP estimator, respectively, while DLS1 and DLS2 represent the estimation of CATE and the auxiliary function, respectively. It can be seen that CV for DSR is working well as a smaller $\lambda$ is chosen as the sample size grows in DGP-L, and the optimal basis is selected in DGP-Q. However, for other estimators, CV sometimes fails to select the small regularization and optimal basis even with the large sample size.
Table (ref) and (ref) are the simulation results for DGP-L and DGP-Q, respectively. Bias and SD in the tables represent the bias and standard deviation of the estimates at $(X_1,X_2)=(1,1)$, respectively, and MSE is calculated on the test samples. In most cases, DSR outperforms the other estimators, while it has a slightly larger bias than DWLS in DGP-Q. However, DSR has a significantly smaller MSE than DWLS as the sample size grows, which implies that the former is more efficient. DSR is also found to have an advantage over other estimators in terms of computational efficiency. As shown in Table (ref), DSR requires the shortest computation time for estimation (without CV), reflecting its simple estimation procedures. Moreover, the computational advantage of DSR becomes even clearer by taking CV into account. For example, since there are 3 possible basis vectors and 4 possible regularization parameters for each function in this simulation, we must repeat CV 24 times for SEP and DWLS and 144 times for DLS, while DSR requires only 12 times as it has no nuisance estimation.
Then, we conduct the sensitivity analysis to evaluate the impact of the hyperparameters $k,\lambda$ on the performance of the estimators. The detailed results for all estimators are summarised in Table (ref) and (ref) in Section (ref) of Appendix. Figure (ref) and (ref) illustrate the change in MSE of DSR with different hyperparameter combinations. It can be seen that smaller $\lambda$ gives smaller MSE as the number of samples increases in both DGPs. However, in DGP-L, MSE changes significantly with the choice of $\lambda$ as the sample size grows, whereas in DGP-Q, MSE does not respond to $\lambda$ to the same extent at a large sample size. Therefore, when the number of samples is more than 1000 in DGP-L, MSE can be small depending on the choice of $\lambda$ even if the model is not optimal ($k=3$), but in DGP-Q, if the model is too simple ($k=3$) or too complex ($k=10$), MSE is considerably larger and the importance of the choice of $k$ is relatively high. Also, by comparing the results of the sensitivity analysis and Table (ref), we can see that the CV method for DSR successfully selected the optimal hyperparameters in all cases except when $N=500$ in DGP-Q.
Setting. As in the previous simulation, we consider CATE estimation by data combination, but now $H$ and $W$ depend on covariates. Therefore, simply applying DSR must result in biased estimation. We generate random samples of $O=(HY+(1-H)D,W,H,X)$ with a sample size $N=1000,2000,3000$ from the following DGP:
where $X=(X_1,\ldots,X_5)'$ is a five-dimensional covariate vector with all elements having zero mean, and $\gamma$ is a five-dimensional vector with one in all elements. The covariance matrix of $X$ takes one for diagonal elements, and the absolute value of non-diagonal elements is randomly chosen from $[0.1,0.3]$ except that $X_1$ is independent of all the other covariates. The number of replication is set $1000$.
The estimand is CATE as a function of only the first covariate: $\theta_0(v)=E[Y_1-Y_0|X_1=v]=0.4v$. We use the polynomial basis up to the third order to construct OSR, and the order is selected based on 5-fold CV in each replication. ML estimators used for the nuisance estimation are random forest (RF), gradient boosting trees (GBT) and multi-layer perceptron (MLP). ML estimators are implemented using scikit-learn 1.0.2, and we use the default value for all options and hyperparameters. We use 5-fold cross-fitting to construct the orthogonal signals $\hat{U}$ and $\hat{T}$. Unlike the previous one, regularization is not employed in this simulation.
\noindentResults. Table (ref), (ref) and (ref) present the simulation results when using RF, GBT and MLP for the nuisance estimation, respectively. Bias in the tables denotes the estimation bias evaluated at $v=1$, Width is the largest width of the 95% uniform confidence band, and CVR is the empirical coverage of the 95% uniform confidence band, namely, the proportion of the times out of 1000 replications that the true target function $\theta_0(v)=0.4v$ is included in the estimated confidence band. In addition to mean and standard deviation, we report the $q$-quantiles of MSE, Bias and Width for $q=0.2,0.4,0.6,0.8$ because OSR produces rare but extremely inaccurate estimates. This instability may be due to the structure of the orthogonal signals, where the inverse of the estimated probability is included, and it has been pointed out in the literature singh2019automatic,chernozhukov2022automatic,chernozhukov2022riesz. Further discussion on this point can be found in Section (ref).
The results diverge for the different ML methods when the sample size is small, but it seems that the methods for nuisance estimates have less impact on the performance of OSR as the sample size grows. The mean of MSE and the width of the confidence band get smaller as $N$ gets larger for all ML methods. Although the mean of bias remains almost the same across sample sizes, the standard deviation decreases, which also indicates that an increase in sample size stabilises the estimation. The empirical coverage is reasonably close to the nominal rate of 95% for all the cases with RF and when GBT is used with $N=1000,2000$. However, when using GBT with $N=3000$ and MLP with all the sample sizes, the empirical coverage deviates from the nominal coverage, possibly reflecting the slightly large bias and small width in these cases. The large bias in GBT and MLP may be due to the fixed hyperparameters. Although we used fixed hyperparameters to reduce the computational burden in this simulation, we could select hyperparameters based on CV in practice to lower the bias and obtain the confidence band with the correct coverage.
In this section, we apply OSR to estimate a causal effect of participation in 401(k) on household assets.
\noindentSetting. In the US, 401(k) is an employer-sponsored personal pension program first implemented in 1978. It is tax-deferred to encourage household savings for their retirement. Data from the US Census Bureau's 1991 Survey of Income and Program Participation (SIPP) has been used in poterba1994401,poterba1995401 and many subsequent studies to examine the effect of 401(k) participation on household savings. The key technical challenge in estimating the causal effect of 401(k) participation is there are not enough covariates available in the SIPP data to explain the self-selected participation among those who are eligible to participate in 401(k). poterba1994401,poterba1995401 argue that eligibility for the 401(k) program can be regarded as exogenous after conditioning on some important variables because whether an employer offered 401(k) would not affect people's job selection at least at the time the program just started, but they would instead make a decision based on other aspects of the job such as salary. Adopting this argument, we can estimate LATE of 401(k) participation in household savings with 401(k) eligibility as IV and other variables related to job choice to control selection bias.
In this empirical example, we use the data analyzed in chernozhukov2004effects and chernozhukov2018, consisting of samples of 9915 households with the reference person of 25-64 years old, and at least one member is employed but no one is self-employed. We use net financial assets ---defined as the sum of IRA (Individual Retirement Account) balances, 401(k) balances, checking accounts, US saving bonds, other interest-earning accounts in financial institutions, other interest-earning assets, stocks, and mutual funds minus non-mortgage debt--- as the outcome $Y$, a binary indicator for 401(k) eligibility as IV $Z$, and a binary indicator for 401(k) participation as the treatment $D$. The covariate vector $X$ used in this analysis includes age, income, family size, years of education, marital status, two-earner status, defined benefit pension status, IRA participation status, and home ownership status.
We conduct analyses based on the usual one-sample LATE estimation as in Example (ref) and two-sample LATE estimation explained in Section (ref). The parameter of interest is LATE as a function of income. For two-sample LATE estimation, we generate a dataset indicator $H$ such that $P(H=1|X=x)=\varsigma(0.1\gamma'\tilde{x})$, where $\tilde{x}$ is a vector of the covariates scaled so that the values fall within $[0,1]$. We use the polynomial basis and the order is selected based on CV from $k=1,2,3$. We use a relatively large number of partitions $G=20$ to mitigate the impact of random sample-splitting on the performance of OSR. GBT is used for nuisance estimation.
\noindentResults. Figure (ref)(a) and (b) illustrate the estimated LATE function of income and its 95% uniform confidence band constructed by one-sample estimation and two-sample estimation, respectively. As can be seen, $k=1$ is selected in both one-sample and two-sample estimations. The estimated LATE is linearly increasing in household income and the slope is approximately 0.36 in one-sample estimation and 0.45 in two-sample estimation. Furthermore, the point estimate of LATE for households with annual income \$50000 is \$16225 in one-sample estimation and \$17141 in two-sample estimation, which are consistent with the analysis in ogburn2015 whose estimate is \$14910. Statistically significant positive LATE is indicated for households with annual income \$24000-\$68000 in one-sample estimation, while LATE is significantly positive for households with income \$26000-61000 in two-sample estimation, reflecting the wider confidence band in two-sample estimation.
This section discusses the limitation and future direction of the present study. The first topic is model selection in the proposed framework. As explained in Remark (ref), we can perform model selection based on CV using the criterion ((ref)) only when $\zeta_0$ is strictly positive. We explained the several examples where $\zeta_0$ is necessarily positive, and it is shown in the simulation study in Section (ref) that CV using the criterion ((ref)) works well. However, there are also situations where $\zeta_0$ can take both positive and negative values, such as LATE estimation and IDID. Therefore, a model selection method for the general situation is key to increasing the practicality of the proposed framework. Despite its practical importance, little attention has been paid to model selection in the treatment effects estimation schuler2018comparison,caron2020estimating. To the best of our knowledge, there exist only a few attempts in the literature to develop a flexible method for selecting the treatment effect model brookhart2006semiparametric,rolling2014model,saito2020. Although we may be able to extend the ideas of these previous studies for the CEFR problems, further research on model selection is essential to enhance the feasibility of causal inference methods.
One of the drawbacks of the proposed framework is the large variability found in the simulation of Section (ref). It may be due to the structure of the orthogonal signals, in which the inverse of the estimated propensity score is used. Inverse probability weighting (IPW) is known to suffer from unstable estimates especially when the propensity score is close to zero wooldridge2002,Wooldridge2007-vn,robins2007comment,Seaman2013-mw. This problem is especially acute in the data combination settings including two-sample estimation of LATE and IDID, where we have to perform four-class or eight-class classification to estimate propensity scores. Although we can increase stability by trimming small probabilities, determining the optimal threshold is nontrivial lee2011, and trimming can cause additional bias. Recently developed automatic debiased machine learning (Auto-DML) chernozhukov2022automatic can be an effective solution to the problem because it avoids the estimation of propensity scores. Auto-DML directly estimates the Riesz representer of the orthogonal signals rather than constructing it with the inverse of the estimated propensity score. A much smaller variance of Auto-DML compared to the original DML has been empirically verified in numerical experiments singh2019automatic,chernozhukov2022riesz. Thus, extending the procedures and theory of the proposed framework to accommodate signals obtained by Auto-DML is a promising future direction.
The present study proposed a general and flexible framework for the CEFR problems, but more efficient estimation and inference may be possible in some specific settings. For example, the orthogonal moment condition for LATE using the interaction term of $Y$ and $D$ has been proposed in singh2019automatic, while OSR uses $Y$ and $D$ only separately. Comparison of the efficiency of the method in singh2019automatic and OSR is beyond the scope of this study, but intuitively, leveraging information expressed in the form of interaction of $Y$ and $D$ can improve efficiency. However, the contribution of this study for offering the flexible inference framework in the data combination settings is significant, as the method in singh2019automatic is not applicable to situations where $Y$ and $D$ are separately observed.