Regularized Orthogonal Machine Learning for Nonlinear Semiparametric Models

Regularized Orthogonal Machine Learning for Nonlinear Semiparametric Models

Introduction

Conditional moment restrictions (CMRs) often emerge as natural restrictions summarizing conditional independence, exclusion, or structural assumptions in economic models. In such a setting, a major challenge is to allow for a data-driven selection among a (very) large number of conditioning covariates. The sparsity assumption, which requires the number of relevant covariates to be small, has appeared to be an interpretable and plausible alternative. If the target parameter minimizes some population loss, a natural way to impose sparsity is to add an $\ell_1$-penalty on the target parameter to the sample loss. As discussed in EHT, this Lasso approach has substantial computational and statistical advantages over its alternatives. However, in general, the loss may be difficult to find. Focusing on a class of single index CMRs (Ichimura:93, kleinspady), we make the Lasso approach feasible by deriving the loss for an arbitrary single index moment function and establish convergence rate for the Lasso estimator.

The starting point of the analysis is a single index CMR that, in addition to target parameter, may depend on a nuisance component, such as the propensity score, the conditional choice probability, the conditional density, or alike, that can be substantially more high-dimensional than the target parameter itself. A natural approach would be to plug a machine learning/regularized estimate of the nuisance parameter into the moment function, such as random forest or neural networks. We first adjust the moment function so that the gradient of the future loss function is insensitive (i.e., Neyman-orthogonal, chernozhukov2016double) with respect to the first-stage regularization bias, preserving the single index property. We then take the loss function to be an indefinite integral of the adjusted moment function with respect to the single-index. If the original moment function is monotone in the single-index, we report the global minimum of $\ell_1$-regularized $M$-estimator loss, following Negahban. Under mild conditions, the proposed estimator converges at the oracle rate, where the oracle knows the true value of the nuisance parameter and solves only the parametric problem.

We demonstrate the utility of our method with theoretical and empirical applications. First, we introduce a partially linear logistic model with heterogeneous treatment effects and derive an orthogonal loss. Second, we derive an orthogonal loss in conditional moment models with missing data, as studied in carroll:95, carroll:91, chen:08, lee:95, stepanski:93 and static games of incomplete information (e.g. see bajari:10 and bajari:13 among others). In all these settings, we give sufficient primitive conditions on the nuisance parameters to achieve oracle convergence. In the empirical application, we study the heterogeneous Jobs First effect on welfare participation decision via a partially linear logistic model. To detect treatment effect's heterogeneity, it is essential to use the orthogonal loss rather than the non-orthogonal one.

\paragraph{Literature review. }This paper is related to three lines of research: single index models, high-dimensional sparse models, and orthogonal/debiased inference based on machine learning methods. The first line of research concerns with the estimation of single index models (manski:75, powell, manski:85, Ichimura:93, kleinspady, ichimura:hall). This work focuses on a low-dimensional target parameter that can be treated as fixed. Focusing on smooth models, we allow the parameter's dimension to grow with sample size and even exceed it. We show that the single index property is sufficient to ensure the uniform convergence of the sample moments towards its population analog over an $\ell_1$-restricted ball, extending the generalization bounds in the machine learning literature (see e.g., Shalev2014) to single index CMRs.

The second line of research establishes the finite-sample bounds for a high-dimensional sparse parameter (belloni:11, Negahban, loh:13, geer, loh2017, zhu:17, Zhu2017). Our contribution is to allow the loss function to depend on a functional nuisance parameter. In the convex case, we establish global convergence of the $\ell_1$-regularized $M$-estimator, following Negahban. One could follow a similar path to establish local convergence, building on loh:13 and loh2017, in the non-convex case. Focusing on the double robustness property, tan2017regularized and tan2018modelassisted establish convergence guarantees in high-dimensional models that are potentially misspecified. After we released the working paper version of this article (arxiv.ID 1702.06240), many methods have proposed similar $M$-estimation approaches. Focusing on missing data, chakrabortty2019high develops an $M$-estimator, relying on a classical Robins orthogonal score. We derive an orthogonal loss function for an arbitrary single index moment restriction, including Robins as a special case. The follow-up paper by foster2019orthogonal extends our result to $M$-estimators with decomposable regularizers, including $\ell_1$ penalty as a special case. An alternative regularized minimum distance approach has been proposed in belloni2018highdimensional.

The third line of research obtains a $\sqrt{N}$-consistent and asymptotically normal estimator of a low-dimensional target parameter $\theta$ in the presence of a nonparametric nuisance parameter (Neyman:1959, Neyman:1979, HardleStoker1989,bickel:93, NeweyStoker, andrews1994,Newey1994, Robins, chen:03). A statistical procedure is called Neyman-orthogonal if it is locally insensitive with respect to the estimation error of the first-stage nuisance parameter. Combining orthogonality and sample splitting, the Double Machine Learning framework of LRSP and chernozhukov2016double has derived a root-N consistent asymptotically normal estimator of the target parameter based on the first-stage machine learning estimates. Extending this work, we establish convergence rates for $\ell_1$-regularized $M$-estimators whose loss function gradient is an orthogonal moment. Next, we also contribute to the literature that derives orthogonal moments starting from non-orthogonal ones. Specifically, we construct a bias correction term for a nuisance parameter that is identified by a general conditional exogeneity restriction, covering conditional expectation (Newey1994) and conditional quantile (Newey:2018) as leading special cases. An alternative approach based on automatic debiasing has been proposed in chernozhukov2021debiased and chernozhukov2021automatic. Finally, we also contribute to the growing literature on orthogonal/doubly robust estimation based on machine learning methods (sasaki2018estimation, chiang2019multiway, sasaki2020unconditional, Chiang), in particular, heterogeneous treatment effects estimation (nie2017, CGST, OprescuWu, Lieli, ZimLech, Colangelo, CherSem). In contrast to this work, our main example features partialling out inside the argument of nonlinear link function, which, to the best of our knowledge, is completely new.

This paper is organized as follows. Section (ref) introduces our main examples and gives a non-technical overview of the results. Section (ref) formally states our results. Section (ref) derives the concrete conditions for our applications. Section (ref) gives an empirical application. Appendix (ref) generalizes Theorem (ref) to the case of extremum estimators beyond $M$-estimators. Appendix (ref) proves Theorem (ref) and Theorem (ref). Appendix (ref) verifies the conditions of Theorem (ref) for each of the three applications.

Set-Up

We start with the description of a single index conditional moment restriction (CMR) framework. The CMR takes the form

where the first argument $W \in \mathcal W \subseteq \mathrm{R}^{\text{dim} W}$ is the data vector, the second argument $t \in \mathrm{R}$ is the single index, and the third one $\gamma \in \mathrm{R}^d$ is the output of the functional nuisance parameter. The parameter of interest $\theta \in \mathrm{R}^p$ enters the moment function only via its inner product

where the index function $\Lambda(z,\gamma): \mathcal{Z} \bigtimes \mathrm{R}^d \rightarrow \mathrm{R}^p$ is known up to $\gamma$. In many cases (e.g., Ichimura:93, kleinspady), the index function $\Lambda(z,\gamma)$ reduces to $\Lambda(z,\gamma)=z$. Examples of the nuisance vector-function $$ g_0 = g_0(z) $$ include the propensity score, the conditional choice probability, and the regression function, or a combination of these functions. Given the CMR (ref), our goal is to find a loss function $Q(\theta, \gamma): \mathrm{R}^p \bigtimes \mathrm{R}^d \rightarrow \mathrm{R}$ so that the true parameter value $\theta_0$ obeys

We find the loss function in two steps. We first solve the ordinary differential equation (ODE)

We then plug $t=\Lambda(Z,\gamma)'\theta$ and $\gamma = g(Z)$ into the sketch of the loss function $\ell(w, t, \gamma)$ to obtain a population loss

whose gradient at $g=g_0$ and $\theta = \theta_0$ is

In what follows, we refer to $\ell (w, t, \gamma)$ as the loss sketch, to $Q(\theta, g)$ and $\nabla_{\theta} Q(\theta, g)$ as the population loss and the population gradient, to

and to $\nabla_{\theta} \widehat{Q}(\theta, \widehat{g})$ as the sample loss and the sample gradient. The conditional independence assumption (ref) ensures that (ref) is a valid moment equation for $\theta_0$. If $m(w, t, \gamma)$ is non-decreasing (non-increasing) in $t$, $\theta_0$ is the unique minimizer of $Q(\theta, g_0)$ ($-Q(\theta, g_0)$).

The proposed estimator $\widehat{\theta}$ has two stages. First, on the auxiliary sample, we construct an estimate $\widehat{g}$ of the nuisance parameter $g_0$, using a machine learning estimator capable of dealing with the high-dimensional covariate vector $Z$. Second, on the main sample, the target parameter's estimate is taken to be the minimizer of $\ell_1$-regularized sample loss. For a fixed vector, its $\ell_2$ norm is denoted by $\| \cdot \|_2$, the $\ell_1$ norm is denoted by $\| \cdot \|_1$, the $\ell_{\infty}$ norm is denoted by $\| \cdot \|_{\infty}$, and $\ell_0$ norm is denoted by $\| \cdot \|_{0}$.

As discussed in Negahban, the sample gradient $\nabla_{\theta} \widehat{Q} (\theta_0, \widehat{g})$ at $\theta = \theta_0$ summarizes the noise of the problem. If the population gradient (ref) possesses the orthogonality property (Neyman:1959)

the biased first-stage estimation error $\widehat{g}(Z) - g_0(Z)$ has no first-order effect on the sample gradient. As a result, there exists a moderate penalty choice $\lambda=\lambda_{\text{mod}}$ obeying

that is sufficiently large to dominate the noise

If (ref) does not hold, the event (ref) requires an aggressive choice $\lambda=\lambda_{\text{agg}}$ obeying

To sum up, for $C$ sufficiently large, there exists an admissible penalty level $\lambda = \lambda_{\text{adm}} $ obeying (ref)

which reduces to $\lambda_{\text{mod}} $ if (ref) holds (i.e., $B_0 = 0$) and to $\lambda_{\text{agg}} $ otherwise. For an admissible penalty choice, Theorem (ref) establishes the following bounds

In particular, if $B_0 = 0$ and $g_n^2= o( \log p/n)^{1/2}$ and $\lambda =\lambda_{\text{adm}}$ obeying (ref), $\widehat{\theta}$ converges at the oracle rate, where the oracle knows the nuisance parameter $g_0$ and estimates only $\theta_0$. In what follows, if the condition (ref) holds, we refer to the $m(w, t, \gamma)$ and $\ell(w, t, \gamma)$ as the orthogonal moment and orthogonal loss, respectively.

We conclude this section by studying the partially linear logistic model. The model takes the form

where $D \in \mathbb{R}$ is a one-dimensional base treatment, $X \in \mathbb{R}^p$ is a vector of controls, $Y$ is a binary outcome, $W=(D,X,Y)$ is the data vector, and $G(t)$ is the logistic link function. The treatment variable $D$ affects the outcome $Y$ via its interactions with the controls $(1,X)$. In addition, $X$ affects $Y$ via the confounding function $f_0(X)$ that enters (ref) in an additively separable way. To make progress, most papers (see e.g., belloni2016postselection) require this function to be linear so that $\theta_0$ and the nuisance parameter can be estimated under a joint sparsity assumption. We describe below how to circumvent this bottleneck.

Inspired by robinson:88, we propose to partial out the controls inside the link function's argument

where

is the conditional expectation of the treatment and

is the conditional expectation of the link function's argument. In contrast to (ref), the nuisance parameters $p_0(x)$ and $q_0(x)$ are identified separately from $\theta_0$ and permit a wider class of approaches to estimate them. Thus, equation (ref) is a special case of (ref) with the conditioning vector $Z=(D,X)$, the moment function

the single index $t= (d - \gamma_1) \cdot (1,x)'\theta$ and the nuisance parameter $g_0(x) = \{ p_0(x),q_0(x) \}$ whose output is denoted by $\gamma = (\gamma_1, \gamma_2)$.

One may be tempted to proceed with the moment function (ref). Solving the ODE (ref) gives

The sample loss $\widehat{Q}(\theta, (\widehat{p}, \widehat{q}))$ in (ref) coincides with the negative logistic likelihood. However, the population gradient (ref) does not obey the orthogonality condition (ref)

As a result, the biased estimation error $\widehat{q}(X)-q_0(X)$ has a first-order effect on the sample gradient. Thus, the admissible penalty choice reduces to $\lambda_{\text{adm}}=\lambda_{\text{agg}}$ in (ref), which makes the estimator's rate in (ref) a slow one.

To restore orthogonality, we reweigh the moment function as in belloni2016postselection. The sketch of the new moment function is

where $\gamma=(\gamma_1,\gamma_2,\gamma_3)$ corresponds to the output of $g_0(z) = \{ p_0(x),q_0(x),V_0(d,x)\}$ and the weighting function $V_0(d,x)$ is the conditional variance

Solving the ODE (ref) gives

The sample loss coincides with the negative weighted logistic likelihood (belloni2016postselection). In contrast to (ref), the population gradient obeys the orthogonality condition (ref)

As a result, the admissible penalty choice (ref) reduces to $\lambda_{\text{adm}}=\lambda_{\text{mod}}$ in (ref), which makes the estimator's rate in (ref) a fast one.

Examples

Remarks

Following the sparsity bounds established in belloni2016postselection, we conjecture that the final estimator $\widehat{\theta}$ is sparse. Therefore, one can interpret the Lasso estimator $\widehat{\theta}$ as a model selector. We expect the post-Lasso-logistic based on single-selection procedure

to be too sensitive to moderate model selection mistakes, occurring when the non-zero coefficients of $\theta_0$ are statistically indistinguishable from zero. According to LeebPotcher, such inference procedures do not provide a Gaussian approximation that is uniform over the space of $\theta$ and $g_0$ and are not honest. Instead, we discuss the following the debiasing procedure of geer.

Theoretical results

\paragraph{Notation. } We will use the following notation. For two sequences of random variables $a_n, b_n, n \geqslant 1: a_n \lesssim_{P} b_n$ means $a_n = O_{P} (b_n)$. For two sequences of numbers $a_n, b_n, n \geqslant 1$, $a_n \lesssim b_n$ means $a_n = O (b_n)$. Let $T = \{ j: \quad \theta_{0, j} \neq 0 \}$ be the set of coordinates of $\theta_0$ that are not equal to zero, and let $T^c$ be the complement of $T$. For a vector $\delta \in \mathbb{R}^p$, let $(\delta_T)_j = \delta_j$ for each $j \in T$ and $(\delta_T)_j = 0$ for $j \in T^c$. For a vector-valued function $g(z):\mathcal{Z} \rightarrow \mathbb{R}^{d}$, denote its $\ell_2$ and $\ell_{\infty}$ norms as

Furthermore, assume that the index function $\Lambda(z,\gamma): {\mathcal Z} \bigtimes \mathbb{R}^d \rightarrow \mathbb{R}^p$ is sufficiently smooth with respect to $\gamma$, so that the gradient $ \nabla_{\gamma} \Lambda_j(z, \gamma)$ and the Hessian $ \nabla_{\gamma \gamma} \Lambda_j(z, \gamma)$ of each coordinate $j \in \{1,2,\dots, p\}$ are well-defined. Finally, we will use the empirical process notation $${\mathbb{E}_n} f(W_i) := \dfrac{1}{n} \sum_{i=1}^n f(W_i).$$

Assumption (ref) ensures that the loss sketch $\ell(w,t,\gamma)$ is a convex function of $t$ for any $w$ and $\gamma$. As a result, the sample loss $\theta \rightarrow \widehat{Q}(\theta, \widehat{g})$ defined in (ref) is a convex function of $\theta$. By convexity, on the event (ref), the vector of errors $ \nu = \widehat{\theta} - \theta_0 $ belongs to the restricted cone

as shown in Negahban (see Lemma (ref)). Define the restricted set as

Define the sample curvature as

and let the population curvature be the analog of (ref) based on the population Hessian $\nabla_{\theta \theta} Q(\theta, g_0) $ instead of $\nabla_{\theta \theta} \widehat{Q}(\theta, \widehat{g})$.

Assumptions (ref) with $C_{\text{min}}$ and (ref) with $B_{\text{min}}$ imply that the population curvature is bounded from below by $\bar{\gamma} = B_{\text{min}} \cdot C_{\text{min}}$ (see Lemma (ref)). On the event

the sample curvature (ref) is bounded from below by $\bar{\gamma}/2$ (see Lemma (ref)). Since

it suffices to show that each summand above is $o_P(1)$. We bound the first and the second summand in Lemmas (ref) and (ref), respectively.

Assumption (ref) ensures that the moment function is sufficiently smooth in its second and third arguments. Let $\mathcal W$ be an open bounded set containing the support of the data vector $W$. Likewise, let $\mathcal{T}$ be an open set containing the support of $\Lambda(Z,g_0(Z))'\theta_0$. Finally, let $\Gamma$ be an open bounded set containing the support of vector $g(Z)$, when $g \in {\mathcal G}_n$. In what follows, we assume that the sets $\mathcal W, \mathcal{T}, \Gamma$ do not change with $n$.

Lemma (ref) establishes a convergence rate for the sample Hessian uniformly over the restricted set (ref). A gradient version of Lemma (ref) is available in the concurrent work by belloni2018highdimensional.

Assumption (ref) formalizes the convergence of the nuisance parameter's estimator. It introduces a sequence of nuisance realization sets ${\mathcal G}_{n} \subseteq \mathcal{G}$ that contain the true value $g_0$ and the estimator $\widehat{g}$ with probability $1-\delta_n$. As the sample size $n$ increases, the sets ${\mathcal G}_{n}$ shrink. The shrinkage speed is measured by the rate $g_{n}$ and is referred to as the first-stage rate.

Assumption (ref) is satisfied by many machine learning estimators under structural assumptions on the model in $\ell_2$ and/or $\ell_{\infty}$ norm. For example, it holds for Lasso (belloni:11, belloni2013) in linear and generalized linear models, for $L_2$-boosting in sparse models (Luo), for neural network (ChenWhite1991), and for random forest in low-dimensional (wagerwalther) and high-dimensional sparse (fastRF) models. For a broad overview of low-level primitive conditions that covers neural networks and random forest see, e.g., jeong2020robust, Appendix 1.

Theorem (ref) is our main result. It establishes the convergence rate for the Regularized $M$-estimator, covering the orthogonal (i.e., $B_0 = 0$) and the non-orthogonal (i.e., $B_0 \neq 0$) cases. In the former case, the admissible penalty choice $ \lambda_{\text{adm}}$ reduces to $ \lambda_{\text{mod}}$ in (ref), which makes the estimator's rate (ref) fast. Otherwise, $ \lambda_{\text{adm}}$ reduces to $ \lambda_{\text{agg}}$ in (ref), which makes the estimator's rate (ref) slow. In all our examples, we invoke Theorem (ref) twice: first, for the preliminary estimate $\check{\theta}$ based on a non-orthogonal CMR (ref) and, second, for the final estimate $\widehat{\theta}$ based on the orthogonal CMR.

Suppose that the moment function $m_{\text{pre}} (w, t, \gamma_1): \mathcal W \bigtimes \mathbb{R} \bigtimes \mathbb{R} \rightarrow \mathbb{R}$ and/or the index function $\Lambda(z, \gamma_1): \mathcal W \bigtimes \mathbb{R} \rightarrow \mathbb{R}^p$ depend on a one-dimensional functional nuisance parameter $p_0(z)$:

Furthermore, suppose the nuisance parameter is identified by a conditional exogeneity restriction

For example, $ R(W, p(X)) = D - p(X) $ defines the conditional expectation function as in (ref), and $ R(W,p(X)) = 1_{ \{ V \leqslant p(X) \}} - \tau$ defines the conditional $\tau$-quantile function. Starting from an arbitrary CMR (ref), Theorem (ref) derives a loss whose gradient obeys the orthogonality condition (ref).

Starting from an arbitrary CMR (ref), Theorem (ref) adjusts the moment function so that population gradient $\nabla_{\theta} Q(\theta_0, g_0)$ obeys orthogonality condition (ref). Special cases of equation (ref), such as the conditional expectation function and conditional quantile function, are available in Newey1994, LRSP and Newey:2018 for unconditional moment problems. We extend the results above to allow for an arbitrary conditional exogeneity restriction.

Applications

In all the applications below, our starting point is the CMR (ref) which does not obey the orthogonality condition (ref). Invoking either the weighting idea (Example (ref)) or Theorem (ref) (Examples (ref)-(ref)), we derive an orthogonal CMR obeying (ref). However, this CMR requires an additional preliminary estimate of $\theta_0$ on top of the nuisance parameters involved in (ref). Definition (ref) describes the suggested two-step procedure, spelling out the choices of main/auxiliary samples and the penalty parameters in each step.

Nonlinear Treatment Effects

Consider Example (ref). Define the covariance matrix

Missing Data

Consider Example (ref). Define the covariance matrix

and the function $q(t, x)$ as

Static Games of Incomplete Information

Consider Example (ref). Define the covariance matrix

Empirical Application

In this section, we study the short-term impact of Connecticut's Jobs First welfare reform experiment on women's labor supply and welfare participation decisions. Jobs First, a welfare-to-work assistance program, was introduced in Connecticut in the late 1990s as an alternative to the federal Aid to Families with Dependent Children (AFDC) program. Imposing revealed preference restrictions, KT's nonparametric bounds show that Jobs First induced many women to work but led some others to reduce their earnings in order to receive assistance. To gain more insight into the question, we postulate a partially linear logistic model for women's welfare participation decision and estimate heterogeneous Jobs First effects.

The data for our analysis are the same as in KT. They come from Manpower Demonstration Research Corporation (MDRC). In 1996, MDRC conducted a randomized trial that randomly selected a set of eligible female applicants into Jobs First (the treatment group), leaving the remaining females eligible for AFDC (the control group). The outcome of interest is the binary indicator that is equal to one if a woman receives any type of welfare (e.g., AFDC or food stamps) at the fourth quarter after random assignment. Baseline characteristics include constant, age, education level, quarterly history of employment, earnings, AFDC, and food stamps during $8$ quarters before random assignment (RA). We postulate a logistic specification with a partially linear index

where $D=1$ if a woman is assigned to Jobs First and $Y=1$ if a woman is on welfare. We assume that the treatment $D$ and the controls $X$ affect the outcome $Y$ via an additively separable index

entering the logistic link function $G(\cdot): \mathbb{R} \rightarrow \mathbb{R}$. The base treatment $D$ enters the index (ref) through its interactions with the control vector $(1,X)$. In addition, the control vector $X$ enters the index (ref) through an unknown function $f_0(x)$, which summarizes the confounding effect of the controls on $Y$. Our main object of interest is the heterogeneous treatment effects vector $$\theta_0= (\theta_{0,1}, \theta_{0,-1})',$$ where $\theta_{0,1}$ is the baseline treatment effect and $\theta_{0,-1}$ is the treatment interaction effect. We assume that, out of $20$ treatment effect interactions, only few have non-zero value, but do not know their identities.

A standard approach to this problem is to require the confounding function $f_0(x)$ to be a linear sparse function of the controls

where $B(x)$ is a vector of basis functions of $x$ and $n$ is the sample size. We take $B(x)$ to be a vector of $p_{\alpha} = 1, 600$ pairwise interactions of the controls. The joint estimate of $\theta_0$ and $\alpha_0$ is

where the penalty parameter $\lambda_{\text{direct}}$ recommended by belloni2016postselection is $$ \lambda_{\text{direct}} = \dfrac{1.1}{2\sqrt{n}} \Phi^{-1}(1-0.05/ (p+p_{\alpha}) \log (n) ) = 0.036. $$ Following tan2018modelassisted, we do not penalize the first coordinate $\theta$, imposing the sparsity assumption on the heterogeneous modification effects but not the level of the treatment effect for the baseline category.

We compare the direct estimator to the Regularized $M$-Estimator $\widehat{\theta}$, described in Section (ref). This estimator no longer requires the control function $f_0(x)$ to be sparse with respect to the chosen basis $B(x)$. Instead, we assume that the conditional probability of welfare $$ G_0(d,x)= \mathbb{E}[ Y | D=d, X=x] $$ can be well-approximated by trees. While this assumption is less interpretable, it allows $G_0(d,x)$ to have a nonlinear argument. We estimate $G_0(d,x)$ by the probability random forest as in Malley with $100$ trees and default size of leaf node. As for $q_0(x)$, we take $$ \widehat{q}(x) = \widehat{\mathbb{E}} \bigg[ \log \dfrac{\widehat{G}(d,x)}{1-\widehat{G}(d,x)} \bigg| X=x \bigg], $$ where the outer expectation function is estimated by regular random forest of Breiman. In addition, we assume that the propensity score $p_0(x)$ in equation (ref) is a sufficiently smooth function of $x$ and estimate it by simple logistic regression. The weighting function $\widehat{V}(d,x)$ is estimated as in (ref). In Definition (ref), we set $\lambda$ to be $$ \lambda=\lambda_{\text{ortho}} = \dfrac{1.1}{2\sqrt{n}} \Phi^{-1}(1-0.05/n ) = 0.031, $$ standardizing each covariate after interacting it with treatment. The Lasso estimate $\widehat{\theta}$ and its post-penalized analog post-Lasso-logistic $\widehat{\theta}_{\text{PL}}$ are reported in Table (ref), Columns (2)-(3). Finally, we report the unpenalized version of the orthogonal estimator defined as

Since the number of treatment interactions $p=20$ is less than the sample size, this estimator is well-defined. We report this estimate and the standard errors in Table (ref), Columns (4)-(5). Our method is straightforward to implement using the \url{glmnet} function in the \url{glmnet} $R$ package, setting the \url{offset} argument to $\widehat{q}(x)$, the \url{weights} argument to $\widehat{V}(d,x)$, and the \url{penalty.factor}=$(0,1,1,\dots, 1)$ to accommodate the penalty in (ref)

.

Our empirical findings are as follows. The direct estimator (Column (1)) implies that Jobs First unambiguously increased the fraction of women on welfare for all women. In contrast, the orthogonal estimator (Column (2)) implies that women with long history of prior AFDC receipt (at least 5 months) were pushed out of the assistance. This pattern makes sense. The Jobs First (treatment group) faced a time limit of 21 months on welfare while the AFDC (control group) had no time limit. The orthogonal estimator picks up this pattern, while the direct one fails to do so. Finally, the unpenalized model in Columns (4) and (5) finds evidence in favor of treatment effect heterogeneity, but most interaction effects are not significant at $\alpha=0.05$. Therefore, orthogonal Lasso captures a more nuanced pattern of treatment effect that is left out by direct Lasso due to a substantially larger complexity of the nuisance parameter relative to the target one, failure of the sparsity assumption (ref), or both.

\setcounter{equation}{0} \setcounter{theorem}{0} \setcounter{lemma}{0} \setcounter{assumption}{0}

General Theory of Extremum Estimators

In this section, we introduce Regularized Extremum Estimator and establish its properties. Suppose there exists a population loss $Q(\theta, g)$ and its sample analog $\widehat{Q}(\theta, \widehat{g})$. Define the target parameter $\theta_0$ as the unique minimizer of $Q(\theta, g_0$) (i.e., (ref) holds), and the Regularized Extremum Estimator as in (ref). In contrast to (ref), $\widehat{Q}(\theta, \widehat{g})$ may not be a sample average. Finally, let the set $\mathbb{B}$ be as in (ref) and the event $\mathcal{V}$ be as in (ref).

We describe the influence of the estimation error $\widehat{g}-g_0$ on the gradient $\nabla_{\theta} Q(\theta_0, \cdot)$ using pathwise derivatives w.r.t the nuisance parameter $g_0$. The first-order derivative is

and the second-order derivative is

Theorem (ref) establishes the convergence rate for Regularized Extremum Estimator. Its building blocks are described in the following lemmas.

Proofs of Section 3

In this section, we prove Theorem (ref) as a special case of Theorem (ref). We also prove Theorem (ref). Define an event $\mathcal{E}_n := \{ \widehat{g}_{k} \in {\mathcal G}_n \quad \forall k \in [K] \}$, such that the nuisance parameter estimate $\widehat{g}_k$ belongs to the realization set ${\mathcal G}_n$ for each fold $k \in [K]$. By union bound, this event holds w.h.p. $${\mathrm{P}} (\mathcal{E}_n ) \geqslant 1- K \epsilon_n = 1-o(1).$$ For a given partition $k$ in $\{1,2, \dots, K\}$, define the partition-specific averages $${\mathbb{E}_{n,k}} f(W_i) := \dfrac{1}{n_k} \sum_{i \in J_k} f(W_i), \quad {\mathbb{G}_{n,k}} f(W_i) := \dfrac{1}{\sqrt{n_k}} \sum_{i \in J_k} [f(W_i) - \int f(w) dP (w)] .$$

Proofs of Section 4

Let $\check{\theta}$ be a preliminary estimator of $\theta_0$ converging at rate $\theta_n$. Lemma (ref) establishes the shrinkage rates of the nuisance realization sets introduced below.