⁰⁰footnotetext: This paper is a complete and more advanced version of a recent technical note by the author. See the note from https://arxiv.org/abs/2008.12467.¹¹footnotetext: Department of Biostatistics, Harvard Chan School of Public Health.²²footnotetext: Department of Statistics, Harvard University³³footnotetext: Department of Statistics, University of California, Davis.

Double/Debiased Machine Learning for Logistic Partially
Linear Model

Molei Liu¹, Yi Zhang², Doudou Zhou³

Abstract

We propose double/debiased machine learning approaches to infer (at the parametric rate $n^{-1/2}$ ) the parametric component of a logistic partially linear model with the binary response following a conditional logistic model of a low dimensional linear parametric function of some key (exposure) covariates and a nonparametric function adjusting for the confounding effect of other covariates. We consider a Neyman orthogonal (doubly robust) score equation consisting of two nuisance functions: nonparametric component in the logistic model and conditional mean of the exposure on the other covariates and with the response fixed. To estimate the nuisance models, we separately consider the use of high dimensional (HD) sparse parametric models and more general (typically nonparametric) machine learning (ML) methods. In the HD case, we derive certain moment equations to calibrate the first order bias of the nuisance models and grant our method a model double robustness property in the sense that our estimator achieves the desirable $O(n^{-1/2})$ rate when at least one of the nuisance models is correctly specified and both of them are ultra-sparse. In the ML case, the non-linearity of the logit link makes it substantially harder than the partially linear setting to use an arbitrary conditional mean learning algorithm to estimate nuisance component of the logistic model. We handle this obstacle through a novel “full model refitting” procedure that is easy-to-implement and facilitates the use of nonparametric ML algorithms in our framework. Our ML estimator is rate doubly robust in the same sense as Chernozhukov et al., 2018a . We evaluate our methods through simulation studies and apply them in assessing the effect of emergency contraceptive (EC) pill on early gestation foetal with a policy reform in Chile in 2008 (Bentancor and Clarke,, 2017).

1 Introduction

Consider a logistic partially linear model. Let $\{(Y_{i},A_{i},\boldsymbol{X}_{i})\mathrel{\mathop{\mathchar 58\relax}}i=1,2,\ldots,n\}$ be independent and identically distributed samples of $Y\in\{0,1\}$ , $A\in\mathbb{R}$ and $\boldsymbol{X}\in\mathbb{R}^{p}$ . Assume that

\mathbb{P}(Y=1\mid A,\boldsymbol{X})={\rm expit}\{\beta_{0}A+r_{0}(\boldsymbol{X})\},

(1)

where ${\rm expit}(\cdot)={\rm logit}^{-1}(\cdot)$ , ${\rm logit}(a)=\log\{a/(1-a)\}$ and $r_{0}(\cdot)$ is an unknown nuisance function of $\boldsymbol{X}$ . In an experimental or observational study with $A$ taken as the exposure variable, $Y$ being the binary response of interests and $\boldsymbol{X}$ representing the observed confounding variables, parameter $\beta_{0}$ is of particular interests as it measures the conditional effect of $A$ on $Y$ in the scale of logarithmic odds ratio. And as the most common and natural way to characterize the conditional model of a binary outcome against some exposure, model (1) has been extensively used in economics and policy science studies.

Our goal is to estimate and infer $\beta_{0}$ asymptotic normally at the rate $n^{-1/2}$ . When $\boldsymbol{X}$ is a scalar and $r_{0}(\cdot)$ is smooth, classic semiparametric kernel or sieve regression works well for this purpose (Severini and Staniswalis,, 1994; Lin and Carroll,, 2006). When $\boldsymbol{X}$ is of high dimensionality, these approaches can have poor performance due to curse of dimensionality and one would rather estimate $r_{0}(\cdot)$ with modern HD (parametric) or ML (nonparametric)¹¹1Our HD setting refers to HD parametric (linear or generalized linear) model and the ML setting refers to ML models of conditional mean estimation (prediction/classification) that is blackbox and usually nonparametric. methods that are much more resistant to the growing dimensionality and complexity of $\boldsymbol{X}$ . However, unlike the partially linear model scenario (Chernozhukov et al., 2018c, ; Dukes and Vansteelandt,, 2020), robust and efficient inference of $\beta_{0}$ in (1) with HD or ML nuisance models is yet to be well-studied.

In recent, Tan, 2019a proposed a simple and flexible doubly robust estimator to enhance the robustness to the potential misspecification of $r(\boldsymbol{x})$ specified as a fixed-dimensional parametric function: $r(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}$ . They introduced a parametric model $m(\boldsymbol{x})=g(\boldsymbol{x}^{\intercal}\boldsymbol{\alpha})$ with a (known) link function $g(\cdot)$ for the conditional mean model $m_{0}(\boldsymbol{x})=\mathbb{E}(A\mid Y=0,\boldsymbol{X}=\boldsymbol{x})$ and proposed a doubly robust estimating equation:

\frac{1}{n}\sum_{i=1}^{n}\widehat{\phi}(\boldsymbol{X}_{i})\left\{Y_{i}e^{-\beta A_{i}-\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\gamma}}}-(1-Y_{i})\right\}\left\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\right\}=0,

(2)

where $\widehat{\phi}(\boldsymbol{x})$ is an estimation of some scalar nuisance function $\phi(\boldsymbol{x})$ affecting the asymptotic efficiency of the estimator, and $\widehat{\boldsymbol{\alpha}}$ and $\widehat{\boldsymbol{\gamma}}$ are two fixed dimensional nuisance model estimators. Estimator $\widehat{\beta}$ solved from (2) is doubly robust that it is valid when either $r(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}$ is correctly specified for the logistic model nonparametric component, or $m(\boldsymbol{x})=g(\boldsymbol{x}^{\intercal}\boldsymbol{\alpha})$ is correctly specified for the conditional mean model $m_{0}(\boldsymbol{x})=\mathbb{E}(A\mid Y=0,\boldsymbol{X}=\boldsymbol{x})$ . Prior to this, the doubly robust semiparametric estimation of odds ratio was built upon $p(A\mid\boldsymbol{X},Y=0)$ , the conditional density of $A$ given $\boldsymbol{X}$ and $Y=0$ (Chen,, 2007; Tchetgen Tchetgen et al.,, 2010, e.g.), requiring a stronger model assumption²²2Chen, (2007) and Tchetgen Tchetgen et al., (2010) estimate the conditional distribution of $A$ given $\boldsymbol{X}$ while Tan, 2019a only needs to specify a conditional mean model of $A$ . than Tan, 2019a for continuous $A$ .

Nevertheless, Tan, 2019a focuses on fixed dimensional parametric nuisance models that are still prone to model misspecification. And their proposed approach is not readily applicable to the recently developed and exploited HD ( $p\gg n$ and the two nuisance components are specified as parametric models with sparse coefficients) or ML (the nuisance functions are estimated by arbitrary blackbox learning algorithm of condition mean) approaches. For the HD realization, this is because simply using some regularized nuisance estimators in (2) would typically incur excessive bias and not guarantee the parametric rate of convergence. We realize bias reduction with respect to the regularization errors by constructing certain dantzig moment equations to estimate the nuisance parameters. With the ultra-sparse nuisance parameters, i.e. their sparsity level is $o(n^{1/2}/\log p)$ , our proposed estimator preserves the model double robustness property that it approaches $\beta_{0}$ at the rate $O_{p}(n^{-1/2})$ when either $r(\cdot)$ or $m(\cdot)$ is correctly specified. Under the ML framework, the non-linearity and unextractablity of the logit link makes it impossible to naturally estimate $r_{0}(\cdot)$ with a learning algorithm of conditional mean³³3Without purposed modification, the natural form of most contemporary supervised learning methods, e.g. random forest, support vector machine and neural network, can be conceptualized as a blackbox algorithm of conditional mean estimation because their task is prediction for continuous responses and classification for categorical responses. as trivially done under the partially linear model (Chernozhukov et al., 2018a, ). We handle this challenge through an easy-to-implement “full model refitting” (FMR) procedure that facilitates flexible implementation of arbitrary conditional mean learning algorithms in our framework. And our DML estimator for $\beta_{0}$ is rate doubly robust in the same sense as Chernozhukov et al., 2018a , i.e. being asymptotically normal at rate $n^{-1/2}$ when the two nuisance ML estimators are consistent for the true models and their mean squared errors are controlled by $o_{p}(n^{-1/2})$ .

In recent years, there has been a large body of literature developed for semiparametric inference (of a low-dimensional parameter) with HD and ML nuisance models, which gets arising attention and application in economics and policy sciences (Athey and Imbens,, 2017; Knaus,, 2018, 2020; Yang et al.,, 2020). Similar to the structure of our paper, there are two different frameworks among the literature, HD setting and ML setting, though sometimes both of them referred as “machine learning” approaches. As the main difference between them, the HD setting imposes parametric assumptions on the nuisance models and may allow for their potential misspecification while ML setting uses nonparametric ML estimation supposed to approach the true nuisance models at certain geometric rate.

To estimate low dimensional parameters like average treatment effect (ATE) and conditionl treatment effect in linear or log-linear model under the HD setting with potentially misspecified nuisance models, recent work including Smucler et al., (2019); Tan, 2020a ; Tan, 2020b ; Ning et al., (2020); Dukes and Vansteelandt, (2020) constructed $\ell_{1}$ -regularized estimating equations with certain $\ell_{\infty}$ -constraints to simultaneously estimate the nuisance parameters and calibrate their first order bias. In comparison, Bradic et al., (2019) proposed a more sparsity-rate robust ATE estimator that requires substantially weaker sparsity assumptions but needs both HD parametric nuisance models to be correctly specified. Among these literature, our work is the first to consider the logistic partially linear model under a similar regime. In parallel to these but closely relevant to our target, existing approaches including debiased (desparsified) LASSO (Van de Geer et al.,, 2014; Janková and Van De Geer,, 2016) and regularized Riesz representer (Chernozhukov et al., 2018c, ; Belloni et al.,, 2018) used the empirical inverse of the information matrix obtained with $\ell_{1}$ -regularized regression to correct for the bias of logistic LASSO estimator. They imposed a technical ultra-sparsity condition on inverse of the information matrix, which has been criticized as unreasonable and unverifiable (Xia et al.,, 2020). Compared with them, our sparsity assumption is model-specific and thus more reasonable and explainable as will be shown later.

We note that near the finishing of our paper, a parallel paper, Ghosh and Tan, (2020) is published on arxiv⁴⁴4As a preliminary version of our paper, our previous technical note https://arxiv.org/abs/2008.12467 appears slightly earlier than Ghosh and Tan, (2020).. Compared with their work, our HD part is studying the same problem but using different doubly robust estimating equation and calibrated procedures for the nuisance models. While the two proposals have similar theoretical properties and numerical implementation strategy (see our Sections 3.1, 4.1 and Appendix D). Also, they should have similar numerical performance. Certainly, our method under the ML setting introduced below has no overlapping with their work.

For the nonparametric ML realization, Chernozhukov et al., 2018a established a double machine learning (DML) framework utilizing Neyman orthogonal scores and cross-fitting to construct parametrically efficient ML based casual estimator. Their framework has been playing a central role in semiparametric inference with ML. As complements or extensions of it, recent work including Chernozhukov et al., 2018b ; Zimmert and Lechner, (2019); Colangelo and Lee, (2020) localized the orthogonal score function to estimate conditional average treatment effect; Semenova and Chernozhukov, (2020) constructed the best linear approximation of a structural function with ML; Farrell et al., (2018) used deep neural networks for DML estimation; Wager and Athey, (2018); Oprescu et al., (2019) proposed and studied the tree-based ML approaches for causal inference; and Cui and Tchetgen, (2019) proposed a minimax data-driven model selection approach to choose the ML nuisance models with the lowest bias on the DML estimator. The above mentioned work elaborated on specific inference problems including partially linear model, ATE, and heterogeneous treatment effect with their nuisance models directly estimable with arbitrary (supervised) ML algorithms. While as summarized above, $r_{0}(\cdot)$ cannot be estimated likewise with general ML algorithms, due to the non-linear structure of (1). To the best of our knowledge, this paper is the first one to solve this non-trivial technical problem through the proposed FMR procedure.

The rest of this paper is organized as follows. In Section 2, we define the Neyman orthogonal (doubly robust) score equation for logistic partially linear model. In Section 3, we introduce the realization of debiased and DML inference for $\beta_{0}$ under the HD and ML settings separately. In Section 4, we present and justify the asymptotic property of our HD and ML estimator separately. In Sections 5 and 6, we conduct simulation studies on empirical performance of our method and apply it to assess the effect of EC pill on early gestation foetal.

2 Neyman orthogonal score

Before coming to the specific approaches in Section 3, we introduce a Neyman orthogonal (doubly robust) score function for logistic partially linear model and derive its first order bias, which plays a central role in motivating and guiding our method construction and theoretical analysis. Let observation $\boldsymbol{D}_{i}=\{Y_{i},A_{i},\boldsymbol{X}_{i}\}$ for $i=1,\ldots,n$ and $\boldsymbol{D}=\{Y,A,\boldsymbol{X}\}$ be a realization of $\boldsymbol{D}_{i}$ . Motivated by equation (2) proposed and studied in Tan, 2019a , we define the Neyman orthogonal score as

h(\boldsymbol{D};\beta,\eta)=\psi(\boldsymbol{X})\{Ye^{-\beta A}-(1-Y)e^{r(\boldsymbol{X})}\}\{A-m(\boldsymbol{X})\},

where $\eta=\{r(\cdot),m(\cdot),\psi(\cdot)\}$ represents the whole set of nuisance functions. In analog to (2), $r(\cdot)$ and $m(\cdot)$ corresponds to the nonparametric component $r_{0}(\cdot)$ defined in (1) and the nuisance model $m_{0}(\boldsymbol{x})=\mathbb{E}(A\mid Y=0,\boldsymbol{X}=\boldsymbol{x})$ , respectively. And $\psi(\boldsymbol{x})$ is a nuisance function affecting the asymptotic variance of the estimator that may depend on $r(\boldsymbol{x})$ and $m(\boldsymbol{x})$ and actually corresponds to $\phi(\boldsymbol{x})e^{-r(\boldsymbol{x})}$ with $\phi(\boldsymbol{x})$ defined by (2)⁵⁵5We rewrite (2) with $\psi(\boldsymbol{x})=\phi(\boldsymbol{x})e^{-r(\boldsymbol{x})}$ to induce the form of score function with both its partial derivatives on $r$ and $\psi$ are Neyman orthogonal..

Remark 1.

According to Tan, 2019a , the score function $h(\boldsymbol{D};\beta,\eta)$ is doubly robust in the sense that when $r(\cdot)=r_{0}(\cdot)$ or $m(\cdot)=m_{0}(\cdot)$ , $\beta_{0}$ solves $\mathbb{E}h(\boldsymbol{D};\beta,\eta)=0$ . We shortly demonstrate this as follows. When either $r(\cdot)=r_{0}(\cdot)$ or $m(\cdot)=m_{0}(\cdot)$ holds, we have

		$\displaystyle\mathbb{E}\psi(\boldsymbol{X})(1-Y)\{e^{r(\boldsymbol{X})}-e^{r_{0}(\boldsymbol{X})}\}\{A-m(\boldsymbol{X})\}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left[\psi(\boldsymbol{X})\{e^{r(\boldsymbol{X})}-e^{r_{0}(\boldsymbol{X})}\}\{A-m(\boldsymbol{X})\}\Big{\|}Y=0,\boldsymbol{X}\right]=0,$

which combined with (1) lead to that

	$\displaystyle\mathbb{E}h(\boldsymbol{D};\beta_{0},\eta)=$	$\displaystyle\mathbb{E}\psi(\boldsymbol{X})\{Ye^{-\beta_{0}A}-(1-Y)e^{r(\boldsymbol{X})}\}\{A-m(\boldsymbol{X})\}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\psi(\boldsymbol{X})e^{r_{0}(\boldsymbol{X})}\{Ye^{-\beta_{0}A-r_{0}(\boldsymbol{X})}-(1-Y)\}\{A-m(\boldsymbol{X})\}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\psi(\boldsymbol{X})e^{r_{0}(\boldsymbol{X})}\{A-m(\boldsymbol{X})\}\frac{Y-\mathbb{P}(Y=1\mid A,\boldsymbol{X})}{\mathbb{P}(Y=1\mid A,\boldsymbol{X})}=0.$

Suppose the nuisance models $r_{0}(\boldsymbol{x})$ and $m_{0}(\boldsymbol{x})$ are estimated by $\widehat{r}(\boldsymbol{x})$ and $\widehat{m}(\boldsymbol{x})$ converging to $\bar{r}(\boldsymbol{x})$ and $\bar{m}(\boldsymbol{x})$ and $\widehat{\psi}(\boldsymbol{x})$ represents the estimator for $\psi(\boldsymbol{x})$ approaching $\bar{\psi}(\boldsymbol{x})$ . Denote by $\bar{\eta}=\{\bar{r}(\cdot),\bar{m}(\cdot),\bar{\psi}(\cdot)\}$ and $\widehat{\eta}=\{\widehat{r}(\cdot),\widehat{m}(\cdot),\widehat{\psi}(\cdot)\}$ . We then write the Gateaux (partial) derivative of the score function $h(\boldsymbol{D};\beta_{0},\bar{\eta})$ as

\begin{split}&\partial_{\eta}h(\boldsymbol{D};\beta_{0},\bar{\eta})[\eta-\bar{\eta}]\\ =&\partial_{\psi}h(\boldsymbol{D};\beta_{0},\bar{\eta})[\psi-\bar{\psi}]+\partial_{r}h(\boldsymbol{D};\beta_{0},\bar{\eta})[r-\bar{r}]+\partial_{m}h(\boldsymbol{D};\beta_{0},\bar{\eta})[m-\bar{m}]\\ =\mathrel{\mathop{\mathchar 58\relax}}&\{Ye^{-\beta_{0}A}-(1-Y)e^{\bar{r}(\boldsymbol{X})}\}\{A-\bar{m}(\boldsymbol{X})\}\{\psi(\boldsymbol{X})-\bar{\psi}(\boldsymbol{X})\}\\ &-(1-Y)\bar{\psi}(\boldsymbol{X})e^{\bar{r}(\boldsymbol{X})}\{A-\bar{m}(\boldsymbol{X})\}\{r(\boldsymbol{X})-\bar{r}(\boldsymbol{X})\}\\ &-\bar{\psi}(\boldsymbol{X})\{Ye^{-\beta_{0}A}-(1-Y)e^{\bar{r}(\boldsymbol{X})}\}\{m(\boldsymbol{X})-\bar{m}(\boldsymbol{X})\}.\end{split}

(3)

Remark 2.

We evaluate the Neyman orthogonal score on some limiting parameters $\bar{r}(\cdot)$ and $\bar{m}(\cdot)$ instead of on $r_{0}(\cdot)$ and $m_{0}(\cdot)$ as in Chernozhukov et al., 2018a . This is because different from their ML framework (and ours) assuming both nuisance estimators converge to the true models, i.e. $\bar{r}(\cdot)=r_{0}(\cdot)$ and $\bar{m}(\cdot)=m_{0}(\cdot)$ , our HD realization allows at most one nuisance model to be wrongly specified, and thus the score function to be analyzed with may not be evaluated at the true models.

Inspired by our deduction in Remark 1, $\mathbb{E}\partial_{\psi}h(\boldsymbol{D};\beta_{0},\bar{\eta})[\psi-\bar{\psi}]=0$ for any $\psi$ whenever $\bar{r}(\cdot)=r_{0}(\cdot)$ or $\bar{m}(\cdot)=m_{0}(\cdot)$ . Also, $\mathbb{E}\partial_{r}h(\boldsymbol{D};\beta_{0},\bar{\eta})[r-\bar{r}]=0$ when $\bar{m}(\cdot)=m_{0}(\cdot)$ and $\mathbb{E}\partial_{m}h(\boldsymbol{D};\beta_{0},\bar{\eta})[m-\bar{m}]=0$ when $\bar{r}(\cdot)=r_{0}(\cdot)$ . Thus, under the ML setting, $h(\boldsymbol{D};\beta_{0},\bar{\eta})$ satisfies Neyman orthogonality defined in Chernozhukov et al., 2018a and the first order (over-fitting) bias of the estimating equation: $n^{-1}\sum_{i=1}^{n}h(\boldsymbol{D}_{i};\beta,\widehat{\eta})=0$ can be removed through cross-fitting (introduced in Section 3.2) and concentration. While under the HD parametric setting with $\bar{r}(\cdot)\neq r_{0}(\cdot)$ or $\bar{m}(\cdot)\neq m_{0}(\cdot)$ , one needs to carefully construct the moment equations for $\bar{r}(\cdot)$ and $\bar{m}(\cdot)$ to ensure the orthogonality conditions. While similar to existing literature (Chernozhukov et al., 2018a, ; Tan, 2020a, ), removal of the second order (and beyond) bias relies on the assumptions that quality of the nuisance estimators $\widehat{r}(\cdot)$ and $\widehat{m}(\cdot)$ are good enough (see Section 4).

3 Method

In this section, we separately present our specific construction procedures for HD parametric and ML nonparametric realization of the debiased/DML estimator for $\beta_{0}$ , based on the Neyman orthogonal score derived in Section 2.

3.1 High dimensional parametric model realization

Consider the setting with $p\gg n$ , each $\boldsymbol{X}_{i}$ has its first element being $1$ , $r(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}$ and $m(\boldsymbol{x})=g(\boldsymbol{x}^{\intercal}\boldsymbol{\alpha})$ where $g(\cdot)$ is a monotone and smooth link function with derivative $g^{\prime}(\cdot)$ . Inspired by existing work including Smucler et al., (2019); Tan, 2020a ; Dukes and Vansteelandt, (2020), we construct dantzig moment equations to ensure the Neyman orthogonality empirically: $\partial_{r}h(\boldsymbol{D};\beta_{0},\bar{\eta})[r-\bar{r}]=0$ and $\partial_{m}h(\boldsymbol{D};\beta_{0},\bar{\eta})[m-\bar{m}]=0$ , under potential misspecification of (at most one) the nuisance models.

First, we obtain $\widetilde{\boldsymbol{\gamma}}$ as some initial estimator for $\boldsymbol{\gamma}$ that converges to some limiting parameter $\boldsymbol{\gamma}^{*}$ equaling to the true model parameter $\boldsymbol{\gamma}_{0}$ when $r(\boldsymbol{x})$ is correctly specified: $r(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}_{0}$ . And let $\widehat{\psi}(\boldsymbol{x})$ be some estimator of the nuisance function $\psi(\boldsymbol{x})$ depending on $\widetilde{\boldsymbol{\gamma}}$ with its limiting function being $\bar{\psi}(\boldsymbol{x})$ , whose choice will be discussed in Section 3.3 later. According to (3), we obtain $\widehat{\boldsymbol{\alpha}}$ through the dantzig moment equation:

{\rm min}_{\boldsymbol{\alpha}\in\mathbb{R}^{p}}\|\boldsymbol{\alpha}\|_{1}\quad{\rm s.t}\quad\left\|n^{-1}\sum_{i=1}^{n}(1-Y_{i})\widehat{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}}}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\alpha})\}\boldsymbol{X}_{i}\right\|_{\infty}\leq\lambda_{\alpha},

(4)

where $\lambda_{\alpha}$ is a tuning parameter controlling the regularization bias. Finally, we obtain the nuisance estimator $\widehat{\boldsymbol{\gamma}}$ and the targeted HD estimator $\widehat{\beta}_{\sf\scriptscriptstyle HD}$ simultaneously from:

\begin{split}{\rm min}_{\beta\in\mathbb{R},\boldsymbol{\gamma}\in\mathbb{R}^{p}}\|\boldsymbol{\gamma}\|_{1}\quad{\rm s.t}\quad\left\|n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\beta A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}}\}g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\boldsymbol{X}_{i}\right\|_{\infty}&\leq\lambda_{\gamma};\\ n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\beta A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}&=0.\end{split}

(5)

Let the limits of $\{\widehat{\boldsymbol{\alpha}},\widehat{\boldsymbol{\gamma}}\}$ be $\{\bar{\boldsymbol{\alpha}},\bar{\boldsymbol{\gamma}}\}$ , and $\bar{\eta}=\{\bar{r}(\cdot),\bar{m}(\cdot),\bar{\psi}(\cdot)\}$ where $\bar{r}(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\bar{\boldsymbol{\gamma}}$ , $\bar{m}(\boldsymbol{x})=g(\boldsymbol{x}^{\intercal}\bar{\boldsymbol{\alpha}})$ and $\bar{\psi}(\boldsymbol{x})$ as given in Section 3.3. We shall comment on the orthogonality (moment) conditions of our proposal in Remark 3, compare our method with Dukes and Vansteelandt, (2020) in Remark 4, and discuss its numerical implementation with a weighted LASSO formation in Remark 5.

Remark 3.

Neglect the second (and above) order error terms for now. When $r(\boldsymbol{x})$ is correct (see Assumption HD1), i.e. $r_{0}(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}_{0}$ , it naturally holds that $\mathbb{E}\partial_{m}h(\boldsymbol{D};\beta_{0},\bar{\eta})[\widehat{m}-m_{0}]=0$ and $\boldsymbol{\gamma}^{*}=\bar{\boldsymbol{\gamma}}=\boldsymbol{\gamma}_{0}$ . Then the $\ell_{\infty}$ -constraint in (4) imposes that

\mathbb{E}\partial_{r}h(\boldsymbol{D};\beta_{0},\bar{\eta})[\widehat{r}-r_{0}]\approx\mathbb{E}(1-Y)\bar{\psi}(\boldsymbol{X})e^{\boldsymbol{X}^{\intercal}\boldsymbol{\gamma}_{0}}\{A-g(\boldsymbol{X}^{\intercal}\bar{\boldsymbol{\alpha}})\}\boldsymbol{X}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{0})=\mathbf{0}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}_{0}).

When $\bar{m}(\boldsymbol{x})=m_{0}(\boldsymbol{x})=g(\boldsymbol{x}^{\intercal}\boldsymbol{\alpha}_{0})$ , we have $\mathbb{E}\partial_{r}h(\boldsymbol{D};\beta_{0},\bar{\eta})[\widehat{r}-r_{0}]$ and $\bar{\boldsymbol{\alpha}}=\boldsymbol{\alpha}_{0}$ in turn. And the $\ell_{\infty}$ -constraint of (5) corresponds to

\mathbb{E}\partial_{m}h(\boldsymbol{D};\beta_{0},\bar{\eta})[\widehat{m}-m_{0}]\approx\mathbb{E}\bar{\psi}(\boldsymbol{X})\{Ye^{-\beta_{0}A}-(1-Y)e^{\boldsymbol{X}^{\intercal}\bar{\boldsymbol{\gamma}}}\}g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\alpha}_{0})\boldsymbol{X}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}_{0})=\mathbf{0}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}_{0}).

Thus, the Neyman orthogonality condition $\partial_{\eta}h(\boldsymbol{D};\beta_{0},\bar{\eta})[\eta-\bar{\eta}]=0$ as introduced in Section 2 is satisfied under our construction when either $r(\cdot)$ or $m(\cdot)$ is correctly specified.

Remark 4.

Similar the HD partially linear (or log-linear) setting studied in Dukes and Vansteelandt, (2020), estimating equation for the nuisance parameter $\boldsymbol{\gamma}$ in our framework involves the unknown $\beta$ . Unlike their construction procedure that plug-in $\beta_{0}$ as every $\beta\in\mathbb{R}$ to estimate $\boldsymbol{\gamma}$ and invert the resulted score-test $p$ -values for interval estimation of $\beta_{0}$ , we solve for $\widehat{\beta}$ and $\widehat{\boldsymbol{\gamma}}$ jointly from (5) with the moment equation for $\widehat{\boldsymbol{\beta}}$ being doubly robust, as demonstrated in Remark 1. Compared with their method, ours is more friendly in computation and implementation, additionally provides $n^{-1/2}$ -consistent point estimator of $\beta_{0}$ and preserves similar theoretical guarantee (see Section 4.1).

Remark 5.

As is detailed in Appendix D, one could also construct LASSO problems with the same Karus–Kuhn–Tucker (KKT) conditions as the $\ell_{\infty}$ -norm constraints in (4) and (5), to obtain the estimators $\widehat{\boldsymbol{\alpha}}$ and $\widehat{\boldsymbol{\gamma}}$ , which has equivalent theoretical properties as the dantzig equations.⁶⁶6Here we present the dantzig equation version because it is more intuitive and directly connected with the Neyman orthogonal conditions. Due to the non-convexity⁷⁷7Its partial derivative on $\beta$ , $n^{-1}\sum_{i=1}^{n}-\widehat{\psi}(\boldsymbol{X}_{i})Y_{i}e^{-\beta A_{i}}A_{i}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}$ is not always positive or negative definite, through empirically and theoretically, the partial derivative should be negative with very high chances. of the equation in the second row of (5), numerical solution of the LASSO counterpart of (5) cannot be obtained with existing software like R package “glmnet” (Friedman et al.,, 2010) and “RCAL” (Tan, 2019b, ). A direct solution to this is programming an optimization procedure such as the Fisher scoring descent algorithm used by Tan, 2020b . While we also find a convenient way that moderately modifies the construction procedure to make the regularized estimating equations solvable with R package “RCAL”, and use it for the numerical implementation of our method. In Appendix D, we outline this modification and demonstrate its theoretical guarantee.

3.2 Machine learning realization

We turn to a (nonparametric) ML setting under which any learning algorithms of conditional mean could potentially be applied to estimate the nuisance functions. Similar to Chernozhukov et al., 2018a , we randomly split the $n$ samples into $K$ folds: $\mathcal{I}_{1},\mathcal{I}_{2},\ldots,\mathcal{I}_{K}$ of equal size, to assist removing the first order (over-fitting) bias through concentration. Then the cross-fitted estimating equation for $\beta$ is constructed as

n^{-1}\sum_{k=1}^{K}\sum_{i\in\mathcal{I}_{k}}h(\boldsymbol{D}_{i};\beta,\widehat{\eta}^{[\text{-}k]})=0,

(6)

where $\widehat{\eta}^{[\text{-}k]}=\{\widehat{r}^{[\text{-}k]}(\cdot),m^{[\text{-}k]}(\cdot),\psi^{[\text{-}k]}(\cdot)\}$ , representing ML estimators converging to $\bar{r}(\cdot)=r_{0}(\cdot)$ , $\bar{m}(\cdot)=m_{0}(\cdot)$ and $\bar{\psi}(\cdot)$ , obtained with the samples in $\mathcal{I}_{\text{-}k}=\{1,\ldots,n\}\setminus\mathcal{I}_{k}$ and independent of the samples in $\mathcal{I}_{k}$ . Now we present the specific construction procedures for $\widehat{r}^{[\text{-}k]}(\cdot)$ and $m^{[\text{-}k]}(\cdot)$ with the choice of $\psi^{[\text{-}k]}(\cdot)$ again discussed in Section 3.3.

Suppose there is a blackbox learning algorithm $\mathscr{L}(R_{i},\boldsymbol{C}_{i};\mathcal{I})$ that inputs samples $\mathcal{I}\subseteq\{1,2,\ldots,n\}$ with some response $R_{i}$ and covariates $\boldsymbol{C}_{i}$ , and outputs an estimation of $\mathbb{E}[R_{i}\mid\boldsymbol{C}_{i},i\in\mathcal{I}]$ . We outline as follows our approach utilizing $\mathscr{L}$ to estimate the nuisance functions. Corresponding to the definition of $m_{0}(\cdot)$ , it can be estimated by: $\widehat{m}^{[\text{-}k]}(\cdot)=\mathscr{L}(A_{i},\boldsymbol{X}_{i};\mathcal{I}_{\text{-}k}\cap\{i\mathrel{\mathop{\mathchar 58\relax}}Y_{i}=0\})$ . Compared to the partially linear setting in Chernozhukov et al., 2018a , estimation of $r_{0}(\cdot)$ with $\mathscr{L}$ is more sophisticated since it is defined through a non-linear form: $\mathbb{P}(Y=1\mid A,\boldsymbol{X})={\rm expit}\{\beta_{0}A+r_{0}(\boldsymbol{X})\}$ . One could modify some ML approaches, e.g. neural network⁸⁸8By setting the last layer of the neural network to be the combination of a complex network of $\boldsymbol{X}$ and a linear function of $A$ and linking it with the outcome through an expit link. to accommodate this form. However, such modification is not readily available in general, and typically requires additional human efforts on its implementing and validating if it exists. So alternatively, we propose a “full model refitting” (FMR) procedure using an arbitrary $\mathscr{L}$ to estimate $r_{0}(\cdot)$ . Our method is motivated by a simple proposition:

Proposition 1.

Let $M_{0}(A,\boldsymbol{X})=\mathbb{P}(Y=1\mid A,\boldsymbol{X})={\rm expit}\{\beta_{0}A+r_{0}(\boldsymbol{X})\}$ . We have:

\beta_{0}={\rm argmin}_{\beta\in\mathbb{R}}\mathbb{E}\left[{\rm logit}\{M_{0}(A,\boldsymbol{X})\}-\beta(A-\mathbb{E}[A|\boldsymbol{X}])\right]^{2}.

Proof.

For any $\beta\in\mathbb{R}$ , we have

		$\displaystyle\mathbb{E}\left[{\rm logit}\{M_{0}(A,\boldsymbol{X})\}-\beta(A-\mathbb{E}[A\|\boldsymbol{X}])\right]^{2}=\mathbb{E}\left\{\beta_{0}A+r_{0}(\boldsymbol{X})-\beta(A-\mathbb{E}[A\|\boldsymbol{X}])\right\}^{2}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left\{(\beta_{0}-\beta)(A-\mathbb{E}[A\|\boldsymbol{X}])+\eta(\boldsymbol{X})\right\}^{2}=(\beta_{0}-\beta)^{2}\mathbb{E}(A-\mathbb{E}[A\|\boldsymbol{X}])^{2}+\mathbb{E}\{\eta(\boldsymbol{X})\}^{2},$

where $\eta(\boldsymbol{X})=r_{0}(\boldsymbol{X})+\beta_{0}\mathbb{E}[A|\boldsymbol{X}]$ . Thus, $\beta_{0}$ minimizes $\mathbb{E}\left[{\rm logit}\{M_{0}(A,\boldsymbol{X})\}-\beta(A-\mathbb{E}[A|\boldsymbol{X}])\right]^{2}$ . ∎

Further randomly each split $\mathcal{I}_{\text{-}k}$ into $K$ folds $\mathcal{I}_{\text{-}k,1},\ldots\mathcal{I}_{\text{-}k,K},$ of equal size and denote by $\mathcal{I}_{\text{-}k,\text{-}j}=\mathcal{I}_{\text{-}k}\setminus\mathcal{I}_{\text{-}k,j}$ . Motivated by Proposition 1, we first estimate the “full” model $M_{0}(A,\boldsymbol{X})$ with $\mathcal{I}_{\text{-}k,\text{-}j}$ as:

\widehat{M}^{[\text{-}k,\text{-}j]}(\cdot)=\mathscr{L}(Y_{i},(A_{i},\boldsymbol{X}_{i}^{\intercal})^{\intercal};\mathcal{I}_{\text{-}k,\text{-}j}),

and learn $a_{0}(\boldsymbol{x})=\mathbb{E}[A|\boldsymbol{X}=\boldsymbol{x}]$ as $\widehat{a}^{[\text{-}k,\text{-}j]}(\cdot)=\mathscr{L}(A_{i},\boldsymbol{X}_{i};\mathcal{I}_{\text{-}k,\text{-}j})$ . Then we fit the (cross-fitted) least square regression to obtain:

\breve{\beta}^{[\text{-}k]}={\rm argmin}_{\beta\in\mathbb{R}}\frac{1}{|\mathcal{I}_{\text{-}k}|}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}\left[{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(A_{i},\boldsymbol{X}_{i})\}-\beta\{A_{i}-\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{X}_{i})\}\right]^{2},

(7)

as an estimator approaching $\beta_{0}$ at certain rate typically larger than $n^{-1/2}$ . Then $r_{0}(\cdot)$ could be identified through $r_{0}(\boldsymbol{X}_{i})={\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}-\beta_{0}A_{i}$ . Note that the empirically estimated version of ${\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}-\beta_{0}A_{i}$ typically involves $A_{i}$ due to the discrepancy of $\beta_{0}$ and $M(\cdot)$ from their empirical estimation. This can essentially impede removal of the over-fitting bias since $\partial_{r}h(\boldsymbol{D};\beta_{0},\bar{\eta})$ is not orthogonal to the error functions dependending on $A$ . So we further estimate the conditional mean of ${\rm logit}\{M_{0}(A,\boldsymbol{X})\}-\beta_{0}A$ on $\boldsymbol{X}$ to estimate $r_{0}(\cdot)$ . Denote by $W_{i}={\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(A_{i},\boldsymbol{X}_{i})\}$ for each $i\in\mathcal{I}_{\text{-}k,j}$ and get $\widehat{t}^{[\text{-}k]}(\cdot)=\mathscr{L}(W_{i},\boldsymbol{X}_{i};\mathcal{I}_{\text{-}k})$ to estimate $t_{0}(\boldsymbol{x})=\mathrel{\mathop{\mathchar 58\relax}}\mathbb{E}[{\rm logit}\{M_{0}(A,\boldsymbol{X})\}|\boldsymbol{X}=\boldsymbol{x}]$ . Then the estimator of $r_{0}(\cdot)$ is given by:

\widehat{r}^{[\text{-}k]}(\boldsymbol{x})=\widehat{t}^{[\text{-}k]}(\boldsymbol{x})-\breve{\beta}^{[\text{-}k]}\widehat{a}^{[\text{-}k]}(\boldsymbol{x}),\quad\mbox{where}\quad\widehat{a}^{[\text{-}k]}(\boldsymbol{x})=\frac{1}{K}\sum_{j=1}^{K}\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{x}).

(8)

Alternatively, one can estimate $r_{0}(\cdot)$ through

\widehat{r}^{[\text{-}k]}(\cdot)=\log\left(\frac{\mathscr{L}(e^{-\breve{\beta}^{[\text{-}k]}A_{i}},\boldsymbol{X}_{i};\mathcal{I}_{\text{-}k}\cap\{i\mathrel{\mathop{\mathchar 58\relax}}Y_{i}=1\})}{\mathscr{L}(1-Y_{i},\boldsymbol{X}_{i};\mathcal{I}_{\text{-}k})}\right),

(9)

inspired by the moment condition that is sufficient to identify $r_{0}(\cdot)$ :

\mathbb{E}\left[Ye^{-\beta_{0}A}-(1-Y)e^{r_{0}(\boldsymbol{X})}\Big{|}\boldsymbol{X}\right]=\mathbb{E}\left[e^{-\beta_{0}A}\Big{|}\boldsymbol{X},Y=1\right]-e^{r_{0}(\boldsymbol{X})}\mathbb{E}\left[(1-Y)\big{|}\boldsymbol{X}\right]=0.

We refer the estimation step for $\breve{\beta}^{[\text{-}k]}$ and $\widehat{r}^{[\text{-}k]}(\cdot)$ introduced above as “refitting”, and the whole procedure as FMR since we “refit” the least square problem (7) and ML models $\mathscr{L}$ to estimate $r_{0}(\cdot)$ with the initially estimated full model ${\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}$ as a psuedo-outcome. Finally, we solve (6) based on $\widehat{\eta}^{[\text{-}k]}$ to obtain the DML estimator $\widehat{\beta}_{\sf\scriptscriptstyle ML}$ .

Remark 6.

We further use cross-fitting in FMR to avoid over-fitting of the models $\widehat{M}^{[\text{-}k,\text{-}j]}(\cdot)$ and $a^{[\text{-}k,\text{-}j]}(\cdot)$ when they are used to obtain the estimators $\breve{\beta}^{[\text{-}k]}$ , $\widehat{t}^{[\text{-}k]}(\boldsymbol{x})$ and $\widehat{r}^{[\text{-}k]}(\boldsymbol{x})$ .

Remark 7.

The FMR implicitly assumes that $\mathscr{L}$ should perform similarly well on different learning objects with the covariates set as either $\boldsymbol{X}$ or $(A,\boldsymbol{X}^{\intercal})^{\intercal}$ . Classic nonparametric approaches like kernel smoothing or sieve may not satisfy this since including one more covariate $A$ in addition to the very low dimensional $\boldsymbol{X}$ can have substantial impact on estimation performance. Thus, we recommend using more dimensionality-robust modern ML approaches, such as random forest and neural network, in our ML framework. While the classic “plug-in” sieve or kernel method has been well-studied in existing literature (Severini and Staniswalis,, 1994; Lin and Carroll,, 2006).

3.3 Efficiency consideration

The nuisance function $\psi(\cdot)$ in our framework is included and chosen in consideration of estimation efficiency. Tan, 2019a proposed and studied two options for $\phi(\cdot)$ used and defined in (2), with the corresponding function $\psi(\cdot)$ taken as:

\psi_{\rm opt}(\boldsymbol{x})=\frac{e^{-r(\boldsymbol{x})}\mathbb{E}[\{A-m(\boldsymbol{X})\}^{2}|\boldsymbol{X}=\boldsymbol{x},Y=0]}{\mathbb{E}[\{A-m(\boldsymbol{X})\}^{2}/{\rm expit}\{\beta_{0}A+r(\boldsymbol{X})\}|\boldsymbol{X}=\boldsymbol{x},Y=0]};\quad\psi_{\rm simp}(\boldsymbol{x})={\rm expit}\{-r(\boldsymbol{x})\}.

It was shown that when both nuisance models are correctly specified, the estimator solved with the weight $\psi_{\rm opt}(\cdot)$ achieves the minimum asymptotic variance. However, computation of $\psi_{\rm opt}(\cdot)$ involves numerical integration with respect to $\boldsymbol{X}$ given $Y=0$ , making it sometimes inconvenient to implement. So Tan, 2019a proposed a simplified but reasonable choice $\psi_{\rm simp}(\boldsymbol{x})$ defined as above that is obtained by evaluating $\psi_{\rm opt}(\boldsymbol{x})$ at $\beta_{0}=0$ . In the following theoretical and numerical studies, we stick to $\psi(\boldsymbol{x})=\psi_{\rm simp}(\boldsymbol{x})$ , $\widehat{\psi}(\boldsymbol{x})={\rm expit}(-\boldsymbol{x}^{\intercal}\widetilde{\boldsymbol{\gamma}})$ and correspondingly $\bar{\psi}(\boldsymbol{x})={\rm expit}(-\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}^{*})$ under the HD setting, and $\widehat{\psi}^{[\text{-}k]}(\boldsymbol{x})={\rm expit}\{-\widehat{r}^{[\text{-}k]}(\boldsymbol{x})\}$ and $\bar{\psi}(\boldsymbol{x})={\rm expit}\{-r_{0}(\boldsymbol{x})\}$ under the ML setting.

4 Asymptotic analysis

Let $o(\alpha_{n})$ , $O(\alpha_{n})$ , $\omega(\alpha_{n})$ , $\Omega(\alpha_{n})$ and $\Theta(\alpha_{n})$ represent the sequences growing at a smaller, equal/smaller, larger, equal/larger and equal rate of $\alpha_{n}$ , respectively. And let $o_{\mathbb{P}}$ , $O_{\mathbb{P}}$ , $\omega_{\mathbb{P}}$ , $\Omega_{\mathbb{P}}$ and $\Theta_{\mathbb{P}}$ be the corresponding rates with probability approaching $1$ as $n\rightarrow\infty$ . Let $\mathcal{X}\subseteq\mathbb{R}^{p}$ be the domain of $\boldsymbol{X}$ . First, we introduce the regularity condition for $\beta$ and its estimating equation used under both HD and ML settings as Assumption REG, which is standard and can be commonly found in literature of the asymptotic analysis of $M$ -estimator (Van der Vaart,, 2000, Chapter 5). And we shall then study the asymptotic properties of $\widehat{\beta}_{\sf\scriptscriptstyle HD}$ and $\widehat{\beta}_{\sf\scriptscriptstyle ML}$ in Sections 4.1 and 4.2 respectively.

Assumption REG (Regularity of estimating equation for $\beta$ ).

Parameter $\beta$ belongs to a compact set $\mathcal{B}\subseteq\mathbb{R}$ and there exists $\delta_{n}=\Omega(n^{-1/2}\log n)$ such that $(\beta_{0}-\delta_{n},\beta_{0}+\delta_{n})\subseteq\mathcal{B}$ . Exposure $A$ belongs to a compact set $\mathcal{A}$ and $\sup_{\boldsymbol{x}\in\mathcal{X}}|\mathbb{E}[A|\boldsymbol{X}=\boldsymbol{x},Y=y]|=O(1)$ for $y=0,1$ . Also, it is satisfied that

\mathbb{E}\bar{\psi}(\boldsymbol{X})Ye^{-\beta_{0}A}A\{A-\bar{m}(\boldsymbol{X})\}=\Theta(1)\quad\mbox{and}\quad\mathbb{E}h^{2}(\boldsymbol{D};\beta_{0},\bar{\eta})=\Theta(1).

4.1 High dimensional (parametric) setting

Let ${\rm expit}^{\prime}(\cdot)$ be the derivative function of ${\rm expit}(\cdot)$ , $\|\cdot\|_{0}$ represents the number of non-zero elements in a vector and $s=\max\{\|\boldsymbol{\gamma}^{*}\|_{0},\|\bar{\boldsymbol{\gamma}}\|_{0},\|\bar{\boldsymbol{\alpha}}\|_{0}\}$ . We introduce following assumptions to regularize the covariates and nuisance estimators.

Assumption HD1 (Model double robustness).

At least one of the following conditions hold: (a) there exists $\boldsymbol{\gamma}_{0}\in\mathbb{R}^{p}$ such that $r_{0}(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}_{0}$ and $\boldsymbol{\gamma}^{*}=\bar{\boldsymbol{\gamma}}=\boldsymbol{\gamma}_{0}$ ; (b) there exists $\boldsymbol{\alpha}_{0}\in\mathbb{R}^{p}$ such that $m_{0}(\boldsymbol{x})=g(\boldsymbol{x}^{\intercal}\boldsymbol{\alpha}_{0})$ and $\bar{\boldsymbol{\alpha}}=\boldsymbol{\alpha}_{0}$ .

Assumption HD2 (Concentration rate).

It holds that

\displaystyle\left\|n^{-1}\sum_{i=1}^{n}(1-Y_{i}){\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{*})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}\boldsymbol{X}_{i}\right\|_{\infty}=O_{\mathbb{P}}\{(\log p/n)^{1/2}\};

	$\displaystyle\left\\|n^{-1}\sum_{i=1}^{n}{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{*})\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\boldsymbol{X}_{i}\right\\|_{\infty}=O_{\mathbb{P}}\{(\log p/n)^{1/2}\};$
	$\displaystyle\left\\|n^{-1}\sum_{i=1}^{n}{\rm expit}^{\prime}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{*})\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}\boldsymbol{X}_{i}\right\\|_{\infty}=O_{\mathbb{P}}\{(\log p/n)^{1/2}\}.$

Assumption HD3 (Smooth link function).

There exists $L=\Theta(1)$ that for any $u,v\in\mathbb{R}$ ,

|g^{\prime}(u)-g^{\prime}(v)|\leq L|u-v|.

Assumption HD4 (Risk of the $\ell_{1}$ -regularized estimators).

There exists tuning parameters $\lambda_{\alpha},\lambda_{\gamma}=\Theta\{(\log p/n)^{1/2}\}$ such that (4) and (5) have feasible solutions with probability approaching $1$ and

	$\displaystyle\sup_{i\in\{1,\ldots,n\}}\|g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\|=O_{\mathbb{P}}(1);\quad\\|\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*}\\|_{1}+\\|\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}}\\|_{1}+\\|\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}}\\|_{1}=O_{\mathbb{P}}\{s(\log p/n)^{1/2}\};$
	$\displaystyle n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\left[\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}^{2}+\{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\}^{2}\right]+(\widehat{\beta}_{\sf\scriptscriptstyle HD}-\beta_{0})^{2}=O_{\mathbb{P}}(s\log p/n);$
	$\displaystyle n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\left[\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})\}^{2}+\{g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}^{2}\right]=O_{\mathbb{P}}(s\log p/n).$

Assumption HD5 (Ultra-sparsity).

It holds that $s=o(n^{1/2}/\log p)$ .

Remark 8.

Under Assumption HD1 and our constructions (4) and (5) (or the one introduced in Appendix D), the expectations of the terms to be concentrated in Assumption HD2 are $\mathbf{0}$ by Remark 3. Then their maximum norms can be controlled by $O_{\mathbb{P}}\{(\log p/n)^{1/2}\}$ as assumed in HD2, when the covariates $\boldsymbol{X}_{i}$ are bounded, subgaussian or beyond (Kuchibhotla and Chakrabortty,, 2018), using the concentration results derived in existing literature (Giné and Nickl,, 2016).

Remark 9.

Rates of the prediction and estimation risk of the nuisance estimators in Assumption HD4 can be derived following the general theoretical framework for $\ell_{1}$ -regularized estimation introduced in Candes et al., (2007); Bickel et al., (2009); Bühlmann and Van De Geer, (2011); Negahban et al., (2012). The same rate properties has been used for analyzing doubly robust estimators with HD nuisance models in existing literature (Smucler et al.,, 2019; Tan, 2020a, ; Dukes and Vansteelandt,, 2020). Note that (4) and (5) involves the estimators $\widetilde{\boldsymbol{\gamma}}$ or $\widehat{\boldsymbol{\alpha}}$ obtained beforehand. This will require some additional effort on removing the “plug-in” errors of $\widetilde{\boldsymbol{\gamma}}$ or $\widehat{\boldsymbol{\alpha}}$ when deriving the risk rates for $\widehat{\boldsymbol{\alpha}}$ and $\widehat{\boldsymbol{\gamma}}$ , compared to the standard analysis procedures. One could see Tan, 2020a for a similar technical issue and relevant technical details being used to handle it.

In addition, $\sup_{i\in\{1,\ldots,n\}}|g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})|=O_{\mathbb{P}}(1)$ imposed in HD4 is not a standard assumption but is mild and very likely to hold: $\sup_{\boldsymbol{x}\in\mathcal{X}}|g(\boldsymbol{x}^{\intercal}\bar{\boldsymbol{\alpha}})|=O(1)$ according to Assumption REG so we only need $g(\boldsymbol{x}^{\intercal}\widehat{\boldsymbol{\alpha}})-g(\boldsymbol{x}^{\intercal}\bar{\boldsymbol{\alpha}})$ to be $O_{\mathbb{P}}(1)$ uniformly.

Remark 10.

The ultra-sparsity assumption HD5 was also imposed in existing literature including Tan, 2020a and Dukes and Vansteelandt, (2020), to control the rate of bias incurred by the HD estimators: $s\log p/n$ below the parametric rate. For linear nuisance model, existing work like Zhu et al., (2018) and Dukes and Vansteelandt, (2020) suggested to add additional moment (KKT) constraints to relax the ultra-sparsity assumption. However, their approach has not been shown to be feasible for the case with non-linear models yet, so it may be promising but still remains unclear for our framework.

We present the asymptotic property of $\widehat{\beta}_{\sf\scriptscriptstyle HD}$ in Theorem 2 and its proof in Appendix A.

Theorem 1.

Denote by $\bar{I}=\mathbb{E}\bar{\psi}(\boldsymbol{X})Ye^{-\beta_{0}A}A\{A-\bar{m}(\boldsymbol{X})\}$ and $\bar{\sigma}^{2}=\bar{I}^{-2}\mathbb{E}h^{2}(\boldsymbol{D};\beta_{0},\bar{\eta})$ . Under Assumptions REG and HD1–HD5, we have

\sqrt{n}\bar{\sigma}^{-1}(\widehat{\beta}_{\sf\scriptscriptstyle HD}-\beta_{0})=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\bar{\sigma}\bar{I})^{-1}h(\boldsymbol{D}_{i};\beta_{0},\bar{\eta})+o_{\mathbb{P}}(1),

which weakly converge to ${\rm N}(0,1)$ .

Recently, logistic debiased LASSO (Van de Geer et al.,, 2014; Janková and Van De Geer,, 2016) has been criticized on that its sparse inverse information matrix condition is not explainable and justifiable, leading to a subpar performance theoretically and numerically (Xia et al.,, 2020). Interestingly, we find the model sparsity assumption of our method is more reasonable than debiased LASSO and present a simple comparison of these two approaches in Remark 11.

Remark 11.

Assume the logistic model $\mathbb{P}(Y=1\mid A,\boldsymbol{X})={\rm expit}\{\beta_{0}A+\boldsymbol{X}^{\intercal}\boldsymbol{\gamma}_{0}\}$ is correctly specified. As is argued by Xia et al., (2020), assuming the information matrix of the logistic model has an ultra-sparse¹⁰¹⁰10Or approximately sparse (see recent work like Belloni et al., (2018); Ma et al., (2020); Liu et al., (2020)). inverse is crucial to ensure the desirable parametric rate of the debiased LASSO estimator for $\beta_{0}$ . However, this assumption is not explainable or convincing for the common gaussian design with sparse precision matrix, due to the presence of logistic canonical link. In comparison, we require that $\mathbb{E}(A\mid Y=0,\boldsymbol{X}=\boldsymbol{x})=g(\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\alpha}_{0})$ with $\|\boldsymbol{\alpha}_{0}\|_{0}=o(n^{1/2}/\log p)$ . This assumption has two advantages. First, it accommodates nonlinear link function $g(\cdot)$ and can be more reasonable for a categorical $A$ . Second, it is imposed on a conditional model directly and thus more explainable. For example, consider a conditional gaussian model: $(A,\boldsymbol{X}^{\intercal})^{\intercal}\mid\{Y=j\}\sim\mathcal{N}(\boldsymbol{\mu}_{j},\boldsymbol{\Sigma})$ for $j=0,1$ . Then we have $r_{0}(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}_{0}$ where $(\beta_{0},\boldsymbol{\gamma}_{0}^{\intercal})^{\intercal}=\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0})$ and $A\mid\boldsymbol{X},Y=0$ follows a gaussian linear model with the coefficient $\boldsymbol{\alpha}_{0}$ determined by $\boldsymbol{\Sigma}^{-1}$ . Therefore, our sparsity assumptions on $\boldsymbol{\alpha}_{0}$ and $\boldsymbol{\gamma}_{0}$ actually assumes the data generation parameters $\boldsymbol{\Sigma}^{-1}$ and $\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0}$ to be sparse, which are more explainable and verifiable in practice.

4.2 Machine learning (nonparametric) setting

Define that $\|f(\cdot)\|_{Q,q}=\mathrel{\mathop{\mathchar 58\relax}}\|f(\boldsymbol{U})\|_{Q,q}=\mathrel{\mathop{\mathchar 58\relax}}\{\int|f(u)|^{q}dQ(u)\}^{1/q}$ for any real number $q>0$ , function $f(\cdot)$ , random variables $\boldsymbol{U}$ and probability measure $Q$ . Let $P$ denote the probability measure of the observed $\boldsymbol{D}$ . We assume that $K=\Theta(1)$ and introduce the following assumption.

Assumption ML1 (Quality of the ML nuisance estimators).

For each $k\in\{1,2,\ldots,K\}$ ,

	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X}}\|\widehat{r}^{[\text{-}k]}(\boldsymbol{x})-r_{0}(\boldsymbol{x})\|+\|\widehat{m}^{[\text{-}k]}(\boldsymbol{x})-m_{0}(\boldsymbol{x})\|$	$\displaystyle=o_{\mathbb{P}}(1);$
	$\displaystyle\\|\widehat{r}^{[\text{-}k]}(\cdot)-r_{0}(\cdot)\\|_{P,2}+\\|\widehat{m}^{[\text{-}k]}(\cdot)-m_{0}(\cdot)\\|_{P,2}$	$\displaystyle=o_{\mathbb{P}}(n^{-1/4}).$

Remark 12.

Similar to Assumptions 3.2 and 3.4 of Chernozhukov et al., 2018a , our Assumption ML1 requires that the ML estimators for $r_{0}(\cdot)$ and $m_{0}(\cdot)$ are uniformly consistent and their mean squared errors (MSE) achieve the rate $o_{\mathbb{P}}(n^{-1/4})$ . This assumption is also referred as rate doubly robust property (Smucler et al.,, 2019) in that it requires production of the MSEs of $\widehat{r}^{[\text{-}k]}(\cdot)$ and $\widehat{m}^{[\text{-}k]}(\cdot)$ to be $o_{\mathbb{P}}(n^{-1/2})$ . In Appendix C, we provide justification for our proposed FMR procedure to derive that the resulted $\widehat{r}^{[\text{-}k]}(\cdot)$ satisfies Assumption ML1 as long as the learning algorithm $\mathscr{L}$ satisfies the same strong convergence properties as Assumption ML1, on all the learning tasks in FMR. Thus, the use of FMR procedure does not actually clip the wings of the ML algorithm being used in our framework.

We present the asymptotic property of $\widehat{\beta}_{\sf\scriptscriptstyle ML}$ in Theorem 2 with its proof found in Appendix B.

Theorem 2.

Denote by $I_{0}=\mathbb{E}\bar{\psi}(\boldsymbol{X})Ye^{-\beta_{0}A}A\{A-m_{0}(\boldsymbol{X})\}$ and $\sigma_{0}^{2}=I_{0}^{-2}\mathbb{E}h^{2}(\boldsymbol{D};\beta_{0},\eta_{0})$ . Under Assumptions REG and ML1, we have

\sqrt{n}\sigma_{0}^{-1}(\widehat{\beta}_{\sf\scriptscriptstyle ML}-\beta_{0})=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\sigma_{0}I_{0})^{-1}h(\boldsymbol{D}_{i};\beta_{0},\eta_{0})+o_{\mathbb{P}}(1),

which weakly converge to ${\rm N}(0,1)$ .

Remark 13.

Since $\bar{\psi}(\boldsymbol{x})={\rm expit}\{-r_{0}(\boldsymbol{x})\}$ and $\widehat{\psi}^{[\text{-}k]}(\boldsymbol{x})={\rm expit}\{-\widehat{r}^{[\text{-}k]}(\boldsymbol{x})\}$ in our ML case, one could show that $\widehat{\psi}^{[\text{-}k]}(\boldsymbol{x})$ achieves the same strong convergence and rate properties as $\widehat{r}^{[\text{-}k]}(\cdot)$ under Assumption ML1. While generally speaking, uniform consistency of $\widehat{\psi}^{[\text{-}k]}(\cdot)$ is sufficient for the desirable conclusion in Theorem 2 so our framework accommodates more flexible choices on $\psi(\boldsymbol{x})$ , for example, $\psi_{\rm opt}(\boldsymbol{x})$ as introduced in Section 3.3. We demonstrate this point during the proof in Appendix B.

5 Simulation study

We conduct simulation studies for HD and ML settings separately in Sections 5.1 and 5.2, to study the point and interval estimation performance of our method.

5.1 High dimensional (parametric) setting

For the HD parametric setting, we design three data generation configurations introduced as follows to simulate different scenarios of model specification:

(i)

First, generate $Y$ following $P(Y=1)=1/2$ . Then generate $(A,\boldsymbol{X}^{\intercal})^{\intercal}\mid\{Y=j\}\sim\mathcal{N}(\boldsymbol{\mu}_{j},\boldsymbol{\Sigma})$ for $j=0,1$ . Specification of $\boldsymbol{\mu}_{j}$ and $\boldsymbol{\Sigma}$ are presented in Appendix E such that $\beta_{0}=0.5$ , $r_{0}(\boldsymbol{X})=-0.22(X_{1}+X_{2})+0.08(X_{3}+X_{4})$ , and $m_{0}(\boldsymbol{X})=-0.1\dot{3}(X_{1}+X_{2}+X_{3}+X_{4})$ .
(ii)

Generate $Y$ following $P(Y=1)=1/2$ and $(A,\boldsymbol{X}^{\intercal})^{\intercal}\mid\{Y=j\}\sim\mathcal{N}(\boldsymbol{\mu}_{j},\boldsymbol{\Sigma})$ for $j=0,1$ . Specification of $\boldsymbol{\mu}_{j}$ and $\boldsymbol{\Sigma}_{j}$ are presented in Appendix E such that $\beta_{0}=0.5$ , $r_{0}(\boldsymbol{X})=-0.22(X_{1}+X_{2})+0.08(X_{3}+X_{4})-0.15(X_{1}X_{2}+X_{1}X_{3}+X_{2}X_{3})$ , and $m_{0}(\boldsymbol{X})=-0.1\dot{3}(X_{1}+X_{2}+X_{3}+X_{4})$ .

(iii)

First generate $\boldsymbol{X}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma})$ with $\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$ given in Appendix E. Then we generate $A$ given $\boldsymbol{X}$ from a gaussian linear model with unit variance and conditional mean

E(A\mid\boldsymbol{X})=0.15X_{1}+0.15X_{2}+0.15X_{3}+0.15X_{4}+0.075X_{1}X_{2}+0.075X_{1}X_{3}+0.075X_{2}X_{3}.

Finally, generate $Y$ by $\mathbb{P}(Y=1\mid A,\boldsymbol{X})={\rm expit}(0.5A+0.25X_{1}+0.25X_{2}+0.1X_{3}+0.1X_{4})$ .

We realize configurations (i)–(iii) with the sample size $n=1000$ , $1500$ or $2000$ separately and the dimension of $\boldsymbol{X}$ fixed as $p=200$ . Under all these settings, we specify the nuisance models as: $r(\boldsymbol{X})=\boldsymbol{X}^{\intercal}\boldsymbol{\gamma}$ and $m(\boldsymbol{X})=\boldsymbol{X}^{\intercal}\boldsymbol{\alpha}$ . Then both nuisance models are correctly specified under (i), only $m(\boldsymbol{X})$ is correctly specified under (ii), and only $r(\boldsymbol{X})$ is correct under (iii). Note that we could not extract the explicit form of $m_{0}(\boldsymbol{X})$ under (iii) since $A$ is generated conditional on $\boldsymbol{X}$ without fixed $Y=0$ . While we still expect the linear model $m(\boldsymbol{X})=\boldsymbol{X}^{\intercal}\boldsymbol{\alpha}$ is misspecified under (iii) as there is non-linear terms included in $E(A|\boldsymbol{X})$ . Implementing details of our HD approach are presented in Appendix D. Specifically, all the tuning parameters in $\ell_{1}$ -regularized regression are selected using cross-validation among $[0.2(\log p/n)^{1/2},2(\log p/n)^{1/2}]$ . We conducted each setting with $300$ repeated simulations.

Table 1 evaluates the performance of our estimator $\widehat{\beta}_{\sf\scriptscriptstyle HD}$ under all settings on its mean square error (MSE), absolute bias and coverage probability (CP) of the 95% confidence interval (CI) estimated using Gaussian bootstrap multiplier. Under all the settings, our method outputs low root-MSE and bias respectively occupying at most $18\%$ and $7\%$ of the magnitude of the true $\beta_{0}(=0.5)$ when $n=1000$ , and at most $12\%$ and $4\%$ of the $\beta_{0}$ when $n=2000$ . As the sample size $n$ grows, one could see a trend of decaying on the MSEs and bias of our estimator. In addition, under all the settings, our interval estimation has proper CP locating in $\pm 0.03$ range of the nominal level $0.95$ . Thus, our HD estimator performs steadily well under different model specification scenarios as long as at least one nuisance models are correctly specified.

Table 1: Average mean square error (MSE), average absolute bias (Bias), and average coverage probability (CP) of 95% CI of our HD estimator with the sample size set as

1000

1500

and

2000

, under configurations (i)–(iii) described in Section 5.1. Number of repetition for each setting is

300

	Configuration (i)			Configuration (ii)			Configuration (iii)
$n$	$1000$	$1500$	$2000$	$1000$	$1500$	$2000$	$1000$	$1500$	$2000$
MSE	$0.007$	$0.006$	$0.004$	$0.007$	$0.006$	$0.004$	$0.008$	$0.004$	$0.003$
Bias	$0.024$	$0.016$	$0.017$	$0.032$	$0.014$	$0.020$	$0.034$	$0.016$	$0.012$
CP	$0.95$	$0.94$	$0.93$	$0.96$	$0.92$	$0.93$	$0.93$	$0.95$	$0.95$

5.2 Machine learning (nonparametric) setting

To study our proposed method under the ML setting, we let $\boldsymbol{\Sigma}\in R^{p\times p}$ with $\Sigma_{ii}=1$ ; $\Sigma_{ij}=0.2$ for $i\neq j$ , and generate $\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma})$ random vectors and truncate their entries by $(-2,2)$ to obtain $\boldsymbol{X}$ . We then generate $A$ from the gaussian model given $\boldsymbol{X}$ with unit variance and conditional mean $a_{0}(\boldsymbol{X})=\boldsymbol{\zeta}_{a}^{\intercal}f_{a}(\boldsymbol{X})$ where $f_{a}(\boldsymbol{x})$ is a non-linear basis function of $\boldsymbol{x}$ including various types of effects (interaction, indicator and trigonometric function, etc.), defined in Appendix E and $\boldsymbol{\zeta}_{a}$ represents its loading coefficients also given in Appendix E. Finally, we set $\beta_{0}=1$ , $r_{0}(\boldsymbol{X})=\boldsymbol{\zeta}_{r}^{\intercal}f_{r}(\boldsymbol{X})$ (see Appendix E), and generate $Y$ with $\mathbb{P}(Y=1\mid A,\boldsymbol{X})={\rm expit}\{\beta_{0}A+r_{0}(\boldsymbol{X})\}$ . We fix $p=20$ and set $n=1000$ or $n=2000$ separately as two settings.

To estimate the nuisance function $r_{0}(\boldsymbol{x})$ , we use the FMR procedure with its last step being (8). And the number of fold for cross-fitting is set as $K=5$ . For choice of the learning algorithms $\mathscr{L}$ , we consider four ML methods and a hybrid method of the ML models introduced as follows.

(a)

Gradient boosted machines (GBM): an ensemble approach of classification and regression tree (CART) using gradient boosting. Implemented by R-package “gbm” (Greenwell et al.,, 2020).
(b)

Random forest (RF): ensemble of CART with bagging. Implemented with R package “RandomForest” (Liaw and Wiener,, 2002).
(c)

Support vector machine (SVM): specified with linear kernel and implemented using R package “e1071” (Dimitriadou et al.,, 2004).
(d)

Neural network (NN): single hidden layer neural network implemented with R package “nnet” (Ripley et al.,, 2016).
(e)

Best nuisance models (Best): similar to Chernozhukov et al., 2018a , we use a simple hybrid method choosing the ML estimator for each nuisance component with the best prediction performance evaluated by cross-validated sum-squared loss.

All the above mentioned ML algorithms are popular in the field of ML and have been considered in the literature of DML (Chernozhukov et al., 2018a, ; Cui and Tchetgen,, 2019). Tuning parameters of the ML models including the number of trees of GBM and RF, the margin of SVM, and the number of units and the weight decay of NN, are selected using the resampling approach of R package “caret” (Kuhn et al.,, 2020). We conducted each setting of $n$ again with $300$ repeated simulations.

Table 2 presents the resulted average MSE, bias and CP of 95% CI of the estimator $\widehat{\beta}_{\sf\scriptscriptstyle ML}$ obtained with the five ML modelling strategies for $n=1000$ and $n=2000$ separately. The five approaches have relatively consistent performance in terms of MSE, bias and CP under both settings, with the variation of their MSEs smaller than $0.002$ . This demonstrates that performance of our framework is robust to the choice of ML algorithms. While to certain degree, NN has the best performance (with the lowest bias and MSE) when $n=1000$ , and GBM and Best have the best performance when $n=2000$ . Also, interval estimators of all the approaches achieve proper coverage rates.

Table 2: Average mean square error (MSE), average absolute bias (Bias), and average coverage probability (CP) of 95% CI of our ML estimator with sample sizes set as

1000

and

2000

, and the nuisance models estimated by the four ML algorithms and the “Best” approach described in Section 5.2. Number of repetition for each setting is

300

	$n=1000$					$n=2000$
	GBM	RF	SVM	NN	Best	GBM	RF	SVM	NN	Best
MSE	$0.013$	$0.013$	$0.012$	$0.011$	$0.013$	$0.006$	$0.007$	$0.007$	$0.007$	$0.006$
Bias	$0.036$	$0.048$	$0.045$	$0.015$	$0.037$	$0.035$	$0.049$	$0.046$	$0.041$	$0.035$
CP	$0.93$	$0.94$	$0.92$	$0.95$	$0.93$	$0.94$	$0.94$	$0.92$	$0.94$	$0.94$

6 A real example: effect of EC pill on early gestation foetal

In this section, we implement our proposed HD and ML methods to study the effect of emergency contraceptive (EC) pill on the rate of new birth and early gestation foetal death (abortion), by revisiting and exploring the data of a quasi-experimental study based on the policy reform of EC pill’s legislation in Chile (Bentancor and Clarke,, 2017). In their original study, the authors collected all records of birth and foetal death in Chile, as well as a number of municipality level features (education, salary and healthcare, etc) of woman at reproductive age (15–34), in the years around 2008, during which the country was experiencing a reform on the legislation of EC pills. As the consequence of this reform, there is about half of the municipalities in Chile started to provide EC pill freely in 2009, while in the remaining half, EC pill is still not available or restricted in use during that period. This policy was mostly dependent on the political, economic and public health factors characterized by totally $16$ features (denoted as $\boldsymbol{Z}$ ) such as education spending, public health spending, condom use, and political conservativeness. Thus, the treatment of EC pill ( $A=1$ for EC pill accessible; $A=0$ for EC pill not accessible) can be regarded as exogenous for the individuals.

Let $Y^{\scriptscriptstyle(1)}$ denote the indicator for the status of early gestation foetal death in each individual record and $Y^{\scriptscriptstyle(2)}$ indicate giving new birth (pregnant and did not incur foetal death). Assume that

	$\displaystyle\mathbb{P}(Y^{\scriptscriptstyle(1)}=1\mid A,\boldsymbol{Z})={\rm expit}\{\beta_{0}^{\scriptscriptstyle(1)}A+r_{0}^{\scriptscriptstyle(1)}(\boldsymbol{Z})\};$
	$\displaystyle\mathbb{P}(Y^{\scriptscriptstyle(2)}=1\mid A,\boldsymbol{Z})={\rm expit}\{\beta_{0}^{\scriptscriptstyle(2)}A+r_{0}^{\scriptscriptstyle(2)}(\boldsymbol{Z})\},$

where $r_{0}^{\scriptscriptstyle(1)}(\cdot)$ and $r_{0}^{\scriptscriptstyle(2)}(\cdot)$ are two unknown functions. We are interested in inferring the two parameters $\beta_{0}^{\scriptscriptstyle(1)}$ and $\beta_{0}^{\scriptscriptstyle(2)}$ characterizing the log odd ratios (log-OR) of abortion (among the pregnant individuals) and birth (among all individuals) to the treatment of EC pill respectively. To investigate $\beta_{0}^{\scriptscriptstyle(1)}$ , we follow a similar strategy as Bentancor and Clarke, (2017) that focuses on the individual records at age between 15 and 25, on which early gestation foetal death can be viewed as a reasonable proxy for illegal abortion. Note that the prevalence of $Y^{\scriptscriptstyle(1)}$ and $Y^{\scriptscriptstyle(2)}$ in their corresponding populations are both less than $5\%$ , which could cause logistics model unstable to fit. So we randomly downsample the $0$ ’s in both analysis to make the prevalence of $Y^{\scriptscriptstyle(1)}$ and $Y^{\scriptscriptstyle(2)}$ being $1/4$ . This procedure only changes intercepts of the logistic models and does not affect the target parameters. The resulted data set for analyzing $\beta_{0}^{\scriptscriptstyle(1)}$ (abortion) has $n^{\scriptscriptstyle(1)}=5,824$ samples. While we only take a subset with $n^{\scriptscriptstyle(2)}=10,000$ samples for $\beta_{0}^{\scriptscriptstyle(2)}$ so that our algorithms will not require excessive computation times.

For our HD approach, we let $\boldsymbol{X}$ be the $p=175$ dimensional bases joining $\boldsymbol{Z}$ , all the interaction terms of $\boldsymbol{Z}$ and the three-dimensional natural splines of all the continuous variable in $\boldsymbol{Z}$ . And we specify the nuisance functions as $m(\boldsymbol{x})={\rm expit}(\boldsymbol{x}^{\intercal}\boldsymbol{\alpha}^{\scriptscriptstyle(\ell)})$ , and $r(\boldsymbol{x})=\boldsymbol{x}^{\intercal}\boldsymbol{\gamma}^{\scriptscriptstyle(\ell)}$ for $j=1,2$ . We take $\boldsymbol{X}=\boldsymbol{Z}$ as the input covariates of our ML method. The choice and implementation of the ML algorithms are basically the same as in Section 5.2, except that we additionally introduce dropout¹¹¹¹11A common and flexible technique in ML research used for regularization and avoiding over-fitting. Here we randomly and independently set each entry of the training covariates matrix as ${\rm N}(0,1)$ variable with probability $0.4$ . for the two tree-based ML algorithms, GBM and RF, to avoid their over-fitting due to that most covariates are of municipality level but the records are of individual level.

Tables 3 and 4 present the point estimation, 95% CI estimation and $p$ -values of our HD and ML approaches for $\beta^{\scriptscriptstyle(1)}_{0}$ and $\beta^{\scriptscriptstyle(2)}_{0}$ , respectively. For $\beta^{\scriptscriptstyle(1)}_{0}$ , log-OR of early gestation foetal death to EC pill, point estimations of all methods are positive and around $-0.175\pm 0.03$ . Their interval estimations are also internally consistent except that SVM outputs a relatively narrow CI and NN includes $0$ near its CI lower bound. And all methods reject the null: $\beta^{\scriptscriptstyle(1)}_{0}=0$ at level $0.05$ . This is because NN outputs slightly worse prediction model for $A\mid\boldsymbol{X},Y=0$ . Note that the result of our hybrid method “Best” is consistent enough with HD, indicating that our methods under both settings lead to basically the same conclusion. Similar situation occurs to the estimation of $\beta^{\scriptscriptstyle(2)}_{0}$ as well. For $\beta^{\scriptscriptstyle(2)}_{0}$ , all methods reject the null at level $0.05$ and their point and CI estimations are internally consistent (SVM shows a moderate variation from other methods).

Our analysis results reveal that distribution of EC pill could significantly reduce the rate of illegal abortion (in the age group 15–25) and new birth. This is consistent with the results of Bentancor and Clarke, (2017) obtained through their municipality level analysis. Although the estimated effect sizes are of different scales¹²¹²12Their effect is defined in a partially linear model of the abortion/birth rate against the treatment and control variables. While we are measuring the effect of EC pill in a logistic model at the individual level. and thus incomparable across the two studies, our $p$ -values appear to show more significance in that nearly all of them are below $0.05$ while their estimated $p$ -values are between $0.05$ and $0.1$ . This is because that we use more complex and robust nuisance models to adjust for the confounding effects of $\boldsymbol{Z}$ and perform our analysis from the individual level granting us to have larger sample sizes.

Table 3: Point estimations, estimated 95% CI lower/upper bounds (LB/UB) and (two-sided)

p

-values for

\beta^{\scriptscriptstyle(1)}_{0}

(log odds ratio of early gestation foetal death to the treatment of EC pill) of our HD and ML (with the five different realization described in Section 5.2) approaches.

Method	HD	GBM	SVM	RF	NN	Best
$\hat{\beta_{0}}$	-0.175	-0.207	-0.144	-0.181	-0.154	-0.181
CI LB	-0.343	-0.348	-0.267	-0.342	-0.356	-0.342
CI UB	-0.004	-0.067	-0.020	-0.020	0.047	-0.020
$p$ -value	0.046	0.004	0.023	0.028	0.133	0.028

Table 4: Point estimations, estimated 95% CI lower/upper bounds (LB/UB) and (two-sided)

p

-values for

\beta^{\scriptscriptstyle(2)}_{0}

(log odds ratio of new birth to the treatment of EC pill) of our HD and ML (with the five different realization described in Section 5.2) approaches.

Method	HD	GBM	SVM	RF	NN	Best
$\hat{\beta_{0}}$	-0.181	-0.132	-0.099	-0.150	-0.131	-0.150
CI LB	-0.272	-0.235	-0.194	-0.279	-0.255	-0.279
CI UB	-0.091	-0.030	-0.005	-0.021	-0.008	-0.021
$p$ -value	0.000	0.012	0.039	0.023	0.037	0.023

Acknowledgements

The authors thanks their advisor, Tianxi Cai for helpful discussion and comments on this paper.

References

Athey and Imbens, (2017) Athey, S. and Imbens, G. W. (2017). The state of applied econometrics: Causality and policy evaluation. Journal of Economic Perspectives, 31(2):3–32.
Belloni et al., (2018) Belloni, A., Chernozhukov, V., Chetverikov, D., Hansen, C., and Kato, K. (2018). High-dimensional econometrics and regularized GMM. arXiv preprint arXiv:1806.01888.
Bentancor and Clarke, (2017) Bentancor, A. and Clarke, D. (2017). Assessing plan B: The effect of the morning after pill on children and women. The Economic Journal, 127(607):2525–2552.
Bickel et al., (2009) Bickel, P. J., Ritov, Y., Tsybakov, A. B., et al. (2009). Simultaneous analysis of lasso and dantzig selector. The Annals of statistics, 37(4):1705–1732.
Bradic et al., (2019) Bradic, J., Wager, S., and Zhu, Y. (2019). Sparsity double robust inference of average treatment effects. arXiv preprint arXiv:1905.00744.
Bühlmann and Van De Geer, (2011) Bühlmann, P. and Van De Geer, S. (2011). Statistics for high-dimensional data: methods, theory and applications. Springer Science & Business Media.
Candes et al., (2007) Candes, E., Tao, T., et al. (2007). The dantzig selector: Statistical estimation when $p$ is much larger than $n$ . The annals of Statistics, 35(6):2313–2351.
Chen, (2007) Chen, H. Y. (2007). A semiparametric odds ratio model for measuring association. Biometrics, 63(2):413–421.
(9) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018a). Double/debiased machine learning for treatment and structural parameters.
(10) Chernozhukov, V., Newey, W., Robins, J., and Singh, R. (2018b). Double/debiased machine learning of global and local parameters using regularized riesz representers. arXiv preprint arXiv:1802.08667.
(11) Chernozhukov, V., Newey, W. K., and Robins, J. (2018c). Double/debiased machine learning using regularized riesz representers. Technical report, cemmap working paper.
Colangelo and Lee, (2020) Colangelo, K. and Lee, Y.-Y. (2020). Double debiased machine learning nonparametric inference with continuous treatments. arXiv preprint arXiv:2004.03036.
Cui and Tchetgen, (2019) Cui, Y. and Tchetgen, E. T. (2019). Bias-aware model selection for machine learning of doubly robust functionals. arXiv preprint arXiv:1911.02029.
Dimitriadou et al., (2004) Dimitriadou, E., Hornik, K., Leisch, F., Meyer, D., and Weingessel, A. (2004). R package: e1071: Misc functions of the department of statistics.
Dukes and Vansteelandt, (2020) Dukes, O. and Vansteelandt, S. (2020). Inference on treatment effect parameters in potentially misspecified high-dimensional models. Biometrika.
Farrell et al., (2018) Farrell, M. H., Liang, T., and Misra, S. (2018). Deep neural networks for estimation and inference: Application to causal effects and other semiparametric estimands. arXiv preprint arXiv:1809.09953.
Friedman et al., (2010) Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1):1–22.
Ghosh and Tan, (2020) Ghosh, S. and Tan, Z. (2020). Doubly robust semiparametric inference using regularized calibrated estimation with high-dimensional data. arXiv preprint arXiv:2009.12033.
Giné and Nickl, (2016) Giné, E. and Nickl, R. (2016). Mathematical foundations of infinite-dimensional statistical models, volume 40. Cambridge University Press.
Greenwell et al., (2020) Greenwell, B., Boehmke, B., Cunningham, J., Developers, G., and Greenwell, M. B. (2020). Package ‘gbm’.
Janková and Van De Geer, (2016) Janková, J. and Van De Geer, S. (2016). Confidence regions for high-dimensional generalized linear models under sparsity. arXiv preprint arXiv:1610.01353.
Knaus, (2018) Knaus, M. C. (2018). A double machine learning approach to estimate the effects of musical practice on student’s skills. arXiv preprint arXiv:1805.10300.
Knaus, (2020) Knaus, M. C. (2020). Double machine learning based program evaluation under unconfoundedness. arXiv preprint arXiv:2003.03191.
Kuchibhotla and Chakrabortty, (2018) Kuchibhotla, A. K. and Chakrabortty, A. (2018). Moving beyond sub-gaussianity in high-dimensional statistics: Applications in covariance estimation and linear regression. arXiv preprint arXiv:1804.02605.
Kuhn et al., (2020) Kuhn, M., Wing, J., Weston, S., Williams, A., Keefer, C., Engelhardt, A., Cooper, T., Mayer, Z., Kenkel, B., Team, R. C., et al. (2020). Package ‘caret’. The R Journal.
Liaw and Wiener, (2002) Liaw, A. and Wiener, M. (2002). Classification and regression by randomforest. R News, 2(3):18–22.
Lin and Carroll, (2006) Lin, X. and Carroll, R. J. (2006). Semiparametric estimation in general repeated measures problems. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(1):69–88.
Liu et al., (2020) Liu, M., Xia, Y., Cho, K., and Cai, T. (2020). Integrative high dimensional multiple testing with heterogeneity under data sharing constraints. arXiv preprint arXiv:2004.00816.
Ma et al., (2020) Ma, R., Tony Cai, T., and Li, H. (2020). Global and simultaneous hypothesis testing for high-dimensional logistic regression models. Journal of the American Statistical Association, pages 1–15.
Negahban et al., (2012) Negahban, S. N., Ravikumar, P., Wainwright, M. J., Yu, B., et al. (2012). A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers. Statistical Science, 27(4):538–557.
Ning et al., (2020) Ning, Y., Sida, P., and Imai, K. (2020). Robust estimation of causal effects via a high-dimensional covariate balancing propensity score. Biometrika, 107(3):533–554.
Oprescu et al., (2019) Oprescu, M., Syrgkanis, V., and Wu, Z. S. (2019). Orthogonal random forest for causal inference. In International Conference on Machine Learning, pages 4932–4941. PMLR.
Ripley et al., (2016) Ripley, B., Venables, W., and Ripley, M. B. (2016). Package ‘nnet’. R package version, 7:3–12.
Semenova and Chernozhukov, (2020) Semenova, V. and Chernozhukov, V. (2020). Debiased machine learning of conditional average treatment effects and other causal functions. The Econometrics Journal.
Severini and Staniswalis, (1994) Severini, T. A. and Staniswalis, J. G. (1994). Quasi-likelihood estimation in semiparametric models. Journal of the American statistical Association, 89(426):501–511.
Smucler et al., (2019) Smucler, E., Rotnitzky, A., and Robins, J. M. (2019). A unifying approach for doubly-robust $\ell_{1}$ -regularized estimation of causal contrasts. arXiv preprint arXiv:1904.03737.
(37) Tan, Z. (2019a). On doubly robust estimation for logistic partially linear models. Statistics $\&$ Probability Letters, 155:108577.
(38) Tan, Z. (2019b). RCAL: Regularized calibrated estimation. R package version 1.0.
(39) Tan, Z. (2020a). Model-assisted inference for treatment effects using regularized calibrated estimation with high-dimensional data. Annals of Statistics, 48(2):811–837.
(40) Tan, Z. (2020b). Regularized calibrated estimation of propensity scores with model misspecification and high-dimensional data. Biometrika, 107(1):137–158.
Tchetgen Tchetgen et al., (2010) Tchetgen Tchetgen, E. J., Robins, J. M., and Rotnitzky, A. (2010). On doubly robust estimation in a semiparametric odds ratio model. Biometrika, 97(1):171–180.
Van de Geer et al., (2014) Van de Geer, S., Bühlmann, P., Ritov, Y., Dezeure, R., et al. (2014). On asymptotically optimal confidence regions and tests for high-dimensional models. The Annals of Statistics, 42(3):1166–1202.
Van der Vaart, (2000) Van der Vaart, A. W. (2000). Asymptotic statistics, volume 3. Cambridge University Press.
Wager and Athey, (2018) Wager, S. and Athey, S. (2018). Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association, 113(523):1228–1242.
Xia et al., (2020) Xia, L., Nan, B., and Li, Y. (2020). A revisit to debiased lasso for generalized linear models. arXiv preprint arXiv:2006.12778.
Yang et al., (2020) Yang, J.-C., Chuang, H.-C., and Kuan, C.-M. (2020). Double machine learning with gradient boosting and its application to the big $n$ audit quality effect. Journal of Econometrics.
Zhu et al., (2018) Zhu, Y., Bradic, J., et al. (2018). Significance testing in non-sparse high-dimensional linear models. Electronic Journal of Statistics, 12(2):3312–3364.
Zimmert and Lechner, (2019) Zimmert, M. and Lechner, M. (2019). Nonparametric estimation of causal heterogeneity under high-dimensional confounding. arXiv preprint arXiv:1908.08779.

Appendix

Appendix A Proof of Theorem 1

Proof.

By (5), we have

n^{-1}\sum_{i=1}^{n}h(\boldsymbol{D}_{i};\widehat{\beta}_{\sf\scriptscriptstyle HD},\widehat{\eta})=n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\widehat{\beta}_{\sf\scriptscriptstyle HD}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}=0.

Our main involvement is to remove the approximation error $n^{-1}\sum_{i=1}^{n}h(\boldsymbol{D}_{i};\widehat{\beta}_{\sf\scriptscriptstyle HD},\widehat{\eta})-h(\boldsymbol{D}_{i};\widehat{\beta}_{\sf\scriptscriptstyle HD},\bar{\eta})$ , asymptotically. Note that

		$\displaystyle n^{-1}\sum_{i=1}^{n}h(\boldsymbol{D}_{i};\widehat{\beta}_{\sf\scriptscriptstyle HD},\widehat{\eta})-h(\boldsymbol{D}_{i};\widehat{\beta}_{\sf\scriptscriptstyle HD},\bar{\eta})$
	$\displaystyle=$	$\displaystyle n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})(1-Y_{i})\{e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}-e^{\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\widehat{\beta}_{\sf\scriptscriptstyle HD}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}\{{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}})-{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{*})\}\{Y_{i}e^{-\widehat{\beta}_{\sf\scriptscriptstyle HD}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}$
	$\displaystyle=\mathrel{\mathop{\mathchar 58\relax}}$	$\displaystyle\Delta_{1}+\Delta_{2}+\Delta_{3}.$

We handle the terms $\Delta_{1}$ , $\Delta_{2}$ and $\Delta_{3}$ separately as follows. First, we have

	$\displaystyle\Delta_{1}=$	$\displaystyle n^{-1}\sum_{i=1}^{n}(1-Y_{i})\{\widehat{\psi}(\boldsymbol{X}_{i})-\bar{\psi}(\boldsymbol{X}_{i})\}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{1-e^{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}(1-Y_{i})\widehat{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{1-e^{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})}\}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}(1-Y_{i})\bar{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{1-e^{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})}-\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}(1-Y_{i})\bar{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}(1-Y_{i})\bar{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})$
	$\displaystyle=\mathrel{\mathop{\mathchar 58\relax}}$	$\displaystyle\Delta_{11}+\Delta_{12}+\Delta_{13}+\Delta_{14}+\Delta_{15}.$

As $\sup_{i\in\{1,\ldots,n\}}|\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})|=O_{\mathbb{P}}(1)$ by Assumption HD4, there exists $M_{1}=O(1)$ such that with probability approaching $1$ ,

|1-e^{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})}|\leq M_{1}|\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})|,\quad|1-e^{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})}-\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})|\leq M_{1}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}^{2};

(A1)

|\widehat{\psi}(\boldsymbol{X}_{i})-\bar{\psi}(\boldsymbol{X}_{i})|=\frac{e^{\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}*}|1-e^{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}*)}|}{(1+e^{\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}*})(1+e^{\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}}})}\leq M_{1}|\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})|.

(A2)

And by Assumptions REG and HD4, there exists $M_{2}=\Theta(1)$ that $\sup_{i\in\{1,\ldots,n\}}|A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})|+|A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})|\leq M_{2}$ . Consequently, by Assumptions HD4 and boundness of $\psi(\cdot)$ , we have

	$\displaystyle\|\Delta_{11}\|\leq$	$\displaystyle n^{-1}\sum_{i=1}^{n}M_{1}^{2}M_{2}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\|\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\|\|\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\|$
	$\displaystyle\leq$	$\displaystyle M_{1}^{2}M_{2}\left[n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\}^{2}\cdot n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}^{2}\right]^{1/2}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{12}\|\leq$	$\displaystyle n^{-1}\sum_{i=1}^{n}M_{1}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\|\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\|\|g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\|$
	$\displaystyle\leq$	$\displaystyle M_{1}\left[n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}^{2}\cdot n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}^{2}\right]^{1/2}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{13}\|\leq$	$\displaystyle n^{-1}\sum_{i=1}^{n}M_{1}M_{2}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}^{2}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{14}\|\leq$	$\displaystyle n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\|\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\|\|g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\|=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right),~{}\mbox{similar to }\|\Delta_{12}\|.$

By Assumptions HD2 and HD4,

|\Delta_{15}|\leq\left\|n^{-1}\sum_{i=1}^{n}(1-Y_{i})\bar{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}\boldsymbol{X}_{i}\right\|_{\infty}\cdot\|\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}}\|_{1}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right).

Thus, we have $|\Delta_{1}|=O_{\mathbb{P}}(s\log p/n)$ . For $\Delta_{2}$ , we have

	$\displaystyle\Delta_{2}=$	$\displaystyle n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})Y_{i}(e^{-\widehat{\beta}_{\sf\scriptscriptstyle HD}A_{i}}-e^{-\beta_{0}A_{i}})\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}\{\widehat{\psi}(\boldsymbol{X}_{i})-\bar{\psi}(\boldsymbol{X}_{i})\}\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}\bar{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})-g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}\bar{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})$
	$\displaystyle=\mathrel{\mathop{\mathchar 58\relax}}$	$\displaystyle\Delta_{21}+\Delta_{22}+\Delta_{23}+\Delta_{24}.$

Again using Assumptions REG and HD4, there exists $M_{3}=\Theta(1)$ such that

|e^{-\widehat{\beta}_{\sf\scriptscriptstyle HD}A_{i}}-e^{-\beta_{0}A_{i}}|\leq M_{3}|\widehat{\beta}_{\sf\scriptscriptstyle HD}-\beta_{0}|.

(A3)

And by Assumption HD3 and the mean value theorem, for each $i$ , there exits $t_{i}$ lying between $\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}}$ and $\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}}$ such that

|g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})-g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})|\leq|g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g^{\prime}(t_{i})||\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})|\leq L\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})\}^{2}.

(A4)

These combined with (A1), (A2) and Assumptions REG and HD4 lead to that

	$\displaystyle\|\Delta_{21}\|$	$\displaystyle=O\left(\|\widehat{\beta}_{\sf\scriptscriptstyle HD}-\beta_{0}\|\left[n^{-1}\sum_{i=1}^{n}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}^{2}\right]^{1/2}\right)=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{22}\|$	$\displaystyle=O\left(\left[n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}\bar{\boldsymbol{\gamma}}}\}\{\widehat{\psi}(\boldsymbol{X}_{i})-\bar{\psi}(\boldsymbol{X}_{i})\}^{2}\cdot n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}\bar{\boldsymbol{\gamma}}}\}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}^{2}\right]^{1/2}\right)$
		$\displaystyle=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{23}\|$	$\displaystyle=O\left(n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}\bar{\boldsymbol{\gamma}}}\}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})\}^{2}\right)=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right).$

And by Assumptions HD2 and HD4,

|\Delta_{24}|\leq\left\|n^{-1}\sum_{i=1}^{n}\bar{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\boldsymbol{X}_{i}\right\|_{\infty}\cdot\|\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}}\|_{1}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right).

So we also have $|\Delta_{2}|=O_{\mathbb{P}}(s\log p/n)$ . For $\Delta_{3}$ , we have

	$\displaystyle\Delta_{3}=$	$\displaystyle n^{-1}\sum_{i=1}^{n}\{{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}})-{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{*})\}Y_{i}(e^{-\widehat{\beta}_{\sf\scriptscriptstyle HD}A_{i}}-e^{-\beta_{0}A_{i}})\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}\Big{[}\{{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}})-{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{})-{\rm expit}^{\prime}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{})\boldsymbol{X}_{i}^{\intercal}(\boldsymbol{\gamma}^{*}-\widetilde{\boldsymbol{\gamma}})\}$
		$\displaystyle\quad\quad\quad\quad~{}\cdot\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}\Big{]}$
		$\displaystyle+n^{-1}\sum_{i=1}^{n}{\rm expit}^{\prime}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{})\boldsymbol{X}_{i}^{\intercal}(\boldsymbol{\gamma}^{}-\widetilde{\boldsymbol{\gamma}})\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}$
	$\displaystyle=\mathrel{\mathop{\mathchar 58\relax}}$	$\displaystyle\Delta_{31}+\Delta_{32}+\Delta_{33}.$

Again using the mean value theorem and the fact that $|{\rm expit}^{\prime}(u)-{\rm expit}^{\prime}(v)|\leq|u-v|$ for any $u,v\in\mathbb{R}$ (by $|{\rm expit}^{\prime\prime}(\cdot)|\leq 1$ ), we have

|{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}})-{\rm expit}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{*})-{\rm expit}^{\prime}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{*})\boldsymbol{X}_{i}^{\intercal}(\boldsymbol{\gamma}^{*}-\widetilde{\boldsymbol{\gamma}})|\leq\{\boldsymbol{X}_{i}^{\intercal}(\boldsymbol{\gamma}^{*}-\widetilde{\boldsymbol{\gamma}})\}^{2}.

Then by (A2), (A3), and Assumptions REG and HD4, we have

	$\displaystyle\|\Delta_{31}\|=$	$\displaystyle O\left(\|\widehat{\beta}_{\sf\scriptscriptstyle HD}-\beta_{0}\|\left[n^{-1}\sum_{i=1}^{n}\{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\}^{2}\right]^{1/2}\right)=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{32}\|=$	$\displaystyle O\left(n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}\bar{\boldsymbol{\gamma}}}\}\{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\}^{2}\right)=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right).$

And by Assumptions HD2 and HD4,

|\Delta_{33}|\leq\left\|n^{-1}\sum_{i=1}^{n}{\rm expit}^{\prime}(-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}^{*})\{Y_{i}e^{-\beta_{0}A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})\}\boldsymbol{X}_{i}\right\|_{\infty}\|\boldsymbol{\gamma}^{*}-\widetilde{\boldsymbol{\gamma}}\|_{1}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right).

Thus, we have $|\Delta_{3}|=O_{\mathbb{P}}(s\log p/n)$ , and by Assumption HD5,

n^{-1}\sum_{i=1}^{n}h(\boldsymbol{D}_{i};\widehat{\beta}_{\sf\scriptscriptstyle HD},\widehat{\eta})-h(\boldsymbol{D}_{i};\widehat{\beta}_{\sf\scriptscriptstyle HD},\bar{\eta})=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right)=o_{\mathbb{P}}\left(\frac{1}{\sqrt{n}}\right),

which leads to that

n^{-1}\sum_{i=1}^{n}h(\boldsymbol{D}_{i};\widehat{\beta}_{\sf\scriptscriptstyle HD},\bar{\eta})+o_{\mathbb{P}}\left(\frac{1}{\sqrt{n}}\right)=0.

This combined with Remark 1 that $\mathbb{E}h(\boldsymbol{D};\beta_{0},\bar{\eta})=0$ under Assumption HD1, the regularity Assumption REG, and Theorem 5.21 of Van der Vaart, (2000), finally finishes proving Theorem 1.

∎

Appendix B Proof of Theorem 2

Following the general results of the DML estimator with non-linear Neyman orthogonal score presented in Section 3.3 and Theorem 3.3 of Chernozhukov et al., 2018a , we only need to verify their Assumptions 3.3 and 3.4 on our orthogonal score function $h(\boldsymbol{D};\beta,\eta)$ , specifically, Conditions A1 and A2 presented as follows.

Condition A1 (Moment condition with Neyman orthogonality).

It holds that: (a) $\mathbb{E}h(\boldsymbol{D};\beta_{0},\eta_{0})=0$ and $\mathcal{B}$ contains an interval of length $\Theta(n^{-1/2}\log n)$ centred at $\beta_{0}$ ; (b) the map $(\beta,\eta)\rightarrow\mathbb{E}h(\boldsymbol{D};\beta_{0},\eta_{0})$ is twice continuously Gateaux-differentiable; (c) $|\mathbb{E}h(\boldsymbol{D};\beta,\eta_{0})|\geq\min\{|J_{0}(\beta-\beta_{0})|,c_{0}\}$ where the parameters $\eta_{0}=\{r_{0}(\cdot),m_{0}(\cdot),\bar{\psi}\}$ , $c_{0}=\Theta(1)$ and $J_{0}=\partial_{\beta}\mathbb{E}h(\boldsymbol{D};\beta,\eta_{0})|_{\beta=\beta_{0}}=\Theta(1)$ ; (d) $h(\boldsymbol{D};\beta_{0},\eta_{0})$ obeys Neyman orthogonality, i.e. $\partial_{\eta}\mathbb{E}h(\boldsymbol{D};\beta_{0},\eta_{0})[\eta-\eta_{0}]=0$ for all $\eta\in\mathcal{E}$ where the parameter space of $\eta$ : $\mathcal{E}\subseteq\{\eta\mathrel{\mathop{\mathchar 58\relax}}\mathbb{E}|h(\boldsymbol{D};\beta_{0},\eta_{0})[\eta-\eta_{0}]|<\infty\}$ .

Condition A2 (Regularity of the score and quality of the nuisance estimators).

It holds that: (a) $\widehat{\eta}^{[\text{-}k]}$ belongs to the realization set $\mathcal{T}_{n}$ for each $k\in\{1,2,\ldots,K\}$ , with probability approaching 1 where $\mathcal{T}_{n}$ satisfies $\eta_{0}\in\mathcal{T}_{n}$ and conditions given as follows; (b) The space of $\beta$ , $\mathcal{B}$ is bounded and for each $\eta\in\mathcal{T}_{n}$ , the functional space $\mathcal{F}_{\eta}=\{h(\cdot;\beta,\eta)\mathrel{\mathop{\mathchar 58\relax}}\beta\in\mathcal{B}\}$ is measurable and its uniform covering number satisfies: there exists positive constant $R=\Theta(1)$ and $\nu=\Theta(1)$ such that

\sup_{Q}\log\mathcal{N}(\epsilon\|F_{\eta}(\cdot)\|_{Q,2},\mathcal{F}_{\eta},\|\cdot\|_{Q,2})\leq\nu\log(R/\epsilon),\quad\forall\epsilon\in(0,1],

where $F_{\eta}(\cdot)$ is a measurable envelope function for $\mathcal{F}_{\eta}$ : $F_{\eta}(\boldsymbol{D})\geq|h(\boldsymbol{D};\beta,\eta)|$ for all $\boldsymbol{D}$ and $\beta\in\mathcal{B}$ , and there exists $q>2$ such that $\|F_{\eta}(\cdot)\|_{P,q}=O(1)$ ; (c) there exists sequence $\tau_{n}$ :

	$\displaystyle\sup_{\eta=\{r,m,\psi\}\in\mathcal{T}_{n},\beta\in\mathcal{B}}\|\mathbb{E}h(\boldsymbol{D};\beta,\eta)-\mathbb{E}h(\boldsymbol{D};\beta,\{r_{0},m_{0},\psi\})\|=o(\tau_{n}),$
	$\displaystyle\sup_{\eta\in\mathcal{T}_{n},\|\beta-\beta_{0}\|\leq\tau_{n}}\mathbb{E}[h(\boldsymbol{D};\beta,\eta)-h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})]^{2}+\mathbb{E}[h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})-h(\boldsymbol{D};\beta_{0},\eta_{0})]^{2}=o(1),$
	$\displaystyle\sup_{r\in(0,1),\eta\in\mathcal{T}_{n},\|\beta-\beta_{0}\|\leq\tau_{n}}\|\partial_{r}^{2}\mathbb{E}h\{\boldsymbol{D};\beta_{0}+r(\beta-\beta_{0}),\eta_{0}+r(\eta-\{r_{0}(\cdot),m_{0}(\cdot),\psi\})\}\|=o(n^{-1/2});$

and (d) $\mathbb{E}h^{2}(\boldsymbol{D};\beta_{0},\eta_{0})=\Theta(1)$ .

We make some simplification and adaptation on the original assumptions in Chernozhukov et al., 2018a to form Conditions A1-A2, according to our own setting. The only non-trivial change made here is that in Condition A2 (c), we require

\begin{split}&\sup_{\eta=\{r,m,\psi\}\in\mathcal{T}_{n},\beta\in\mathcal{B}}|\mathbb{E}h(\boldsymbol{D};\beta,\eta)-\mathbb{E}h(\boldsymbol{D};\beta,\{r_{0},m_{0},\psi\})|=o(\tau_{n}),\\ &\sup_{r\in(0,1),\eta\in\mathcal{T}_{n},|\beta-\beta_{0}|\leq\tau_{n}}|\partial_{r}^{2}\mathbb{E}h\{\boldsymbol{D};\beta_{0}+r(\beta-\beta_{0}),\eta_{0}+r(\eta-\{r_{0}(\cdot),m_{0}(\cdot),\psi\})\}|=o(n^{-1/2});\end{split}

(A5)

instead of:

	$\displaystyle\sup_{\eta\in\mathcal{T}_{n},\beta\in\mathcal{B}}\|\mathbb{E}h(\boldsymbol{D};\beta,\eta)-\mathbb{E}h(\boldsymbol{D};\beta,\eta_{0})\|=o(\tau_{n}),$
	$\displaystyle\sup_{r\in(0,1),\eta\in\mathcal{T}_{n},\|\beta-\beta_{0}\|\leq\tau_{n}}\|\partial_{r}^{2}\mathbb{E}h\{\boldsymbol{D};\beta_{0}+r(\beta-\beta_{0}),\eta_{0}+r(\eta-\eta_{0})\}\|=o(n^{-1/2}),$

as used in Assumption 3.4 (c) of Chernozhukov et al., 2018a . The first inequality of (A5) is used by Chernozhukov et al., 2018a to derive a preliminary rate for the DML estimator: $|\widehat{\beta}_{\sf\scriptscriptstyle ML}-\beta_{0}|=o_{\mathbb{P}}(\tau_{n})$ (see their Step 1 of the proof of Lemma 6.3), and the second inequality of (A5) is used in their Step 3 the proof of Lemma 6.3 to process the second order error of

n^{-1}\sum_{k=1}^{K}\sum_{i\in\mathcal{I}_{k}}h(\boldsymbol{D}_{i};\beta_{0},\eta_{0})-h(\boldsymbol{D}_{i};\beta,\widehat{\eta}^{[\text{-}k]}),

uniformly for $\beta$ satisfying $|\beta-\beta_{0}|\leq\tau_{n}$ .

Note that our modified two assumptions are still sufficient for deriving these results. Since $\mathbb{E}h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})=0$ holds for all $\psi$ (see our Remark 1), there is actually no need to consider $\mathbb{E}h(\boldsymbol{D};\beta,\{r_{0},m_{0},\psi\})-\mathbb{E}h(\boldsymbol{D};\beta,\{r_{0},m_{0},\psi_{0}\})$ when deriving $|\widehat{\beta}_{\sf\scriptscriptstyle ML}-\beta_{0}|=o_{\mathbb{P}}(\tau_{n})$ . While for the Step 3 of Chernozhukov et al., 2018a , one can instead handle the second order error of

n^{-1}\sum_{k=1}^{K}\sum_{i\in\mathcal{I}_{k}}h(\boldsymbol{D}_{i};\beta_{0},\{r_{0},m_{0},\widehat{\psi}^{[\text{-}k]}\})-h(\boldsymbol{D}_{i};\beta,\widehat{\eta}^{[\text{-}k]}),

again using that $\mathbb{E}h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})=0$ holds for all $\psi$ , and then remove the remaining error:

n^{-1}\sum_{k=1}^{K}\sum_{i\in\mathcal{I}_{k}}h(\boldsymbol{D}_{i};\beta_{0},\{r_{0},m_{0},\psi_{0}\})-h(\boldsymbol{D}_{i};\beta_{0},\{r_{0},m_{0},\widehat{\psi}^{[\text{-}k]}\}),

through concentration based on the fact $\partial_{\psi}\mathbb{E}h(\boldsymbol{D};\beta_{0},\eta_{0})[\psi-\psi_{0}]=0$ and the second inquality of Condition A2 (c): $\sup_{\eta\in\mathcal{T}_{n},|\beta-\beta_{0}|\leq\tau_{n}}\mathbb{E}[h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})-h(\boldsymbol{D};\beta_{0},\eta_{0})]^{2}=o(1)$ .

On the other hand, this modification essentially reduces our requirement on the quality of $\widehat{\psi}^{[\text{-}k]}(\cdot)$ , as mentioned in Remark 13. As will be seen from the proof below, to fulfill the modified Condition A2 (c), we do not require $\|\widehat{\psi}^{[\text{-}k]}(\cdot)-\bar{\psi}(\cdot)\|_{P,2}=o_{\mathbb{P}}(\tau_{n})$ as one shall when following the original version of Chernozhukov et al., 2018a but only needs $\widehat{\psi}^{[\text{-}k]}(\cdot)$ to be uniformly consistent (though the former one may still be justifiable in our case since we take $\bar{\psi}(\boldsymbol{x})={\rm expit}\{-\bar{r}(\boldsymbol{x})\}$ ). We now verify Conditions A1 and A2 based on our Assumptions REG and ML1.

Proof.

Condition A1 (a) is directly given by our logistic partial model assumption 1. Condition A1 (b) is naturally satisfied as $h(\boldsymbol{D};\beta,\eta_{0}+r(\eta-\eta_{0}))$ is a twice continuously differentiable in $(\beta,r)$ . Condition A1 (c) is directly given by Assumption REG. And Condition A1 (d) holds by equation (3), combined with the model assumptions (1) and $\mathbb{E}[A|\boldsymbol{X}=\boldsymbol{x},Y=0]=m_{0}(\boldsymbol{x})$ .

By Assumption ML1 and $\bar{\psi}(\boldsymbol{x})={\rm expit}(-\bar{r}(\boldsymbol{x}))$ , there exists $\zeta_{1,n}=o(1)$ and $\zeta_{2,n}=o(n^{-1/4})$ such that $\widehat{\eta}^{[\text{-}k]}$ belongs to

	$\displaystyle\mathcal{T}_{n}=\mathrel{\mathop{\mathchar 58\relax}}\Big{\{}\eta=(r,m,\psi)\mathrel{\mathop{\mathchar 58\relax}}$	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X}}\|\psi(\boldsymbol{x})-\bar{\psi}(\boldsymbol{x})\|+\|r(\boldsymbol{x})-r_{0}(\boldsymbol{x})\|+\|m(\boldsymbol{x})-m_{0}(\boldsymbol{x})\|\leq\zeta_{1,n},$
		$\displaystyle\mbox{and}\quad\\|r(\cdot)-r_{0}(\cdot)\\|_{P,2}+\\|m(\cdot)-m_{0}(\cdot)\\|_{P,2}\leq\zeta_{2,n}\Big{\}},$

with probability approaching $1$ , for each $k\in\{1,\ldots,K\}$ . We define $\mathcal{T}_{n}$ of Condition A2 in this way such that A2 (a) is satisfied. Now we validate Condition A2 (b). By Assumption REG that $A$ and $\beta$ belong to compact sets, and $|m(\boldsymbol{x})|\leq m_{0}(\boldsymbol{x})+o(1)$ is uniformly bounded for $\eta=\{r,m,\psi\}\in\mathcal{T}_{n}$ , there exists positive $C_{1}=\Theta(1)$ such that for $\eta=\{r,m,\psi\}\in\mathcal{T}_{n}$ ,

	$\displaystyle h(\boldsymbol{D};\beta,\eta)=$	$\displaystyle\psi(\boldsymbol{X})\{Ye^{-\beta A}-(1-Y)e^{r(\boldsymbol{X})}\}\{A-m(\boldsymbol{X})\}$
	$\displaystyle\leq$	$\displaystyle\|\psi(\boldsymbol{X})Ye^{-\beta A}\{A-m(\boldsymbol{X})\}\|+\|\psi(\boldsymbol{X})e^{r(\boldsymbol{X})}(1-Y)\{A-m(\boldsymbol{X})\}\|$
	$\displaystyle\leq$	$\displaystyle C_{1}\{\psi(\boldsymbol{X})+\psi(\boldsymbol{X})e^{r(\boldsymbol{X})}\}=C_{1}\{{\rm expit}(-\psi(\boldsymbol{X}))+{\rm expit}(\psi(\boldsymbol{X}))\}\leq C_{1}+1;$
	$\displaystyle\\|\partial_{\beta}h(\boldsymbol{D};\beta,\eta)\\|_{P,2}^{2}=$	$\displaystyle\mathbb{E}[\psi(\boldsymbol{X})Ye^{-\beta A}A\{A-m(\boldsymbol{X})\}]^{2}\leq C_{1}.$

Then by Example 19.7 of Van der Vaart, (2000), Condition A2 (b) holds with $\nu=1$ and $R$ being the diameter of $\mathcal{B}$ . Note that Condition A2 (d) is again directly given by Assumption REG. It remains to verify Condition A2 (c). For each $\eta=\{r,m,\psi\}\in\mathcal{T}_{n}$ and $\beta\in\mathcal{B}$ , using the boundness of $\beta$ , $A$ , $m_{0}(\boldsymbol{x})$ $\psi(\boldsymbol{x})$ and $\psi(\boldsymbol{x})e^{r(\boldsymbol{x})}$ , there exists $C_{2}=\Theta(1)$ such that

		$\displaystyle\|\mathbb{E}h(\boldsymbol{D};\beta,\eta)-\mathbb{E}h(\boldsymbol{D};\beta,\{r_{0},m_{0},\psi\})\|$
	$\displaystyle\leq$	$\displaystyle\|\mathbb{E}\psi(\boldsymbol{X})Ye^{-\beta A}\{m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\}\|+\|\mathbb{E}\psi(\boldsymbol{X})(1-Y)e^{r(\boldsymbol{X})}\{m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\}\|$
		$\displaystyle+\|\mathbb{E}\psi(\boldsymbol{X})e^{r(\boldsymbol{X})}(1-Y)\{1-e^{r_{0}(\boldsymbol{X})-r(\boldsymbol{X})}\}\{A-m(\boldsymbol{X}_{0})\}\|$
	$\displaystyle\leq$	$\displaystyle C_{2}(\\|m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\\|_{P,2}+\\|r_{0}(\boldsymbol{X})-r(\boldsymbol{X})\\|_{P,2}+\\|r_{0}(\boldsymbol{X})-r(\boldsymbol{X})\\|_{P,2}^{2})\leq 3C_{2}\zeta_{2,n}.$

So we take $\tau_{n}=n^{-1/4}$ and by $\zeta_{2,n}=o(n^{-1/4})$ , the first inequality of Condition A2 (c) is satisfied. Again by the boundness of $\beta$ , $A$ , $m_{0}(\boldsymbol{x})$ $\psi(\boldsymbol{x})$ and $\psi(\boldsymbol{x})e^{r(\boldsymbol{x})}$ ; and $\zeta_{1,n},\tau_{n}=o(1)$ , there exists $C_{3}=\Theta(1)$ such that

		$\displaystyle\mathbb{E}[h(\boldsymbol{D};\beta,\eta)-h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})]^{2}+\mathbb{E}[h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})-h(\boldsymbol{D};\beta_{0},\eta_{0})]^{2}$
	$\displaystyle\leq$	$\displaystyle\mathbb{E}[h(\boldsymbol{D};\beta,\eta)-h(\boldsymbol{D};\beta_{0},\eta)]^{2}+\mathbb{E}[h(\boldsymbol{D};\beta_{0},\eta)-h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})]^{2}$
		$\displaystyle+\mathbb{E}[h(\boldsymbol{D};\beta_{0},\{r_{0},m_{0},\psi\})-h(\boldsymbol{D};\beta_{0},\eta_{0})]^{2}$
	$\displaystyle\leq$	$\displaystyle\mathbb{E}[\psi(\boldsymbol{X})e^{-\beta_{0}A}\{e^{(\beta_{0}-\beta)A}-1\}\{A-m(\boldsymbol{X})\}]^{2}$
		$\displaystyle+\mathbb{E}[\psi(\boldsymbol{X})e^{-\beta_{0}A}\{m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\}]^{2}+\mathbb{E}[\psi(\boldsymbol{X})e^{r(\boldsymbol{X})}\{m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\}]^{2}$
		$\displaystyle+\mathbb{E}[\psi(\boldsymbol{X})e^{r(\boldsymbol{X})}(1-Y)\{1-e^{r_{0}(\boldsymbol{X})-r(\boldsymbol{X})}\}\{A-m_{0}(\boldsymbol{X})\}]^{2}$
		$\displaystyle+\mathbb{E}[\|\bar{\psi}(\boldsymbol{X})-\psi(\boldsymbol{X})\|\{e^{-\beta_{0}A}+e^{r_{0}(\boldsymbol{X})}\}\|A-m_{0}(\boldsymbol{X})\|]^{2}$
	$\displaystyle\leq$	$\displaystyle C_{3}\sup_{a\in\mathcal{A}}\|e^{(\beta_{0}-\beta)a}-1\|+C_{3}\sup_{x\in\mathcal{X}}\{\|(1+e^{r(\boldsymbol{x})})[\bar{\psi}(\boldsymbol{x})-\psi(\boldsymbol{x})]\|+\|m(\boldsymbol{x})-m_{0}(\boldsymbol{x})\|+\|e^{r(\boldsymbol{x})-r_{0}(\boldsymbol{x})}-1\|\}$
	$\displaystyle\leq$	$\displaystyle C_{3}\sup_{a\in\mathcal{A}}\|e^{(\beta_{0}-\beta)a}-1\|+C_{3}\sup_{x\in\mathcal{X}}\{\|\bar{\psi}(\boldsymbol{x})-\psi(\boldsymbol{x})\|+\|m(\boldsymbol{x})-m_{0}(\boldsymbol{x})\|+\{{\rm expit}(r_{0}(\boldsymbol{x}))+1\}\|e^{r(\boldsymbol{x})-r_{0}(\boldsymbol{x})}-1\|\}$
	$\displaystyle\leq$	$\displaystyle C_{3}\{(e^{\tau_{n}C_{3}}-1)+2\zeta_{1,n}+2\|e^{\zeta_{1,n}}-1\|\}=o(1),$

which validates the second inequality of Condition A2 (c). At last, for each $r\in(0,1)$ , denote by $\beta^{*}=\beta_{0}+r(\beta-\beta_{0})$ , $\eta^{*}=\{r^{*},m^{*},\psi\}=\{r_{0}(\cdot),m_{0}(\cdot),\psi\}+r(\eta-\{r_{0}(\cdot),m_{0}(\cdot),\psi\})$ . Similar as above deduction, we have that there exists $C_{4}=\Theta(1)$ ,

		$\displaystyle\partial_{r}^{2}\mathbb{E}h\{\boldsymbol{D};\beta_{0}+r(\beta-\beta_{0}),\eta_{0}+r(\eta-\{r_{0}(\cdot),m_{0}(\cdot),\psi\})\}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\psi(\boldsymbol{X})Ye^{-\beta^{*}A}A(\beta-\beta_{0})\{m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\}$
		$\displaystyle+\mathbb{E}\psi(\boldsymbol{X})(1-Y)e^{r^{*}(\boldsymbol{X})}\{r_{0}(\boldsymbol{X})-r(\boldsymbol{X})\}\{m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\}$
	$\displaystyle\leq$	$\displaystyle C_{4}\|\beta-\beta_{0}\|\cdot\mathbb{E}\|m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\|+C_{4}\mathbb{E}\|r_{0}(\boldsymbol{X})-r(\boldsymbol{X})\|\cdot\mathbb{E}\|m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\|$
	$\displaystyle=$	$\displaystyle O(\\|m(\cdot)-m_{0}(\cdot)\\|_{P,2}^{2})+O\{(\beta-\beta_{0})^{2}\}+O(\\|r(\cdot)-r_{0}(\cdot)\\|_{P,2}^{2})=O(\zeta_{2,n}^{2})+o(\tau_{n}^{2})=o(n^{-1/4}).$

∎

Using the verified Conditions A1 and A2, one can follow nearly the same proof procedures as those of Theorem 3.3 and Lemma 6.3 in Chernozhukov et al., 2018a to prove our Theorem 2. The only minor difference concerning the processing of $\psi^{[\text{-}k]}$ has been presented as above. As we point out, one can handle this smoothly by first considering $n^{-1}\sum_{k=1}^{K}\sum_{i\in\mathcal{I}_{k}}h(\boldsymbol{D}_{i};\beta_{0},\{r_{0},m_{0},\widehat{\psi}^{[\text{-}k]}\})$ when deriving the initial rate and asymptotic expansion of $\widehat{\beta}_{\sf\scriptscriptstyle ML}$ as $\mathbb{E}[h(\boldsymbol{D}_{i};\beta_{0},\{r_{0},m_{0},\widehat{\psi}^{[\text{-}k]}\})|\widehat{\psi}^{[\text{-}k]}]=0$ , and finally concentrate $n^{-1}\sum_{k=1}^{K}\sum_{i\in\mathcal{I}_{k}}h(\boldsymbol{D}_{i};\beta_{0},\{r_{0},m_{0},\psi_{0}\})-h(\boldsymbol{D}_{i};\beta_{0},\{r_{0},m_{0},\widehat{\psi}^{[\text{-}k]}\})$ using that $\partial_{\psi}\mathbb{E}h(\boldsymbol{D};\beta_{0},\eta_{0})[\psi-\psi_{0}]=0$ and $\widehat{\psi}^{[\text{-}k]}(\cdot)$ is uniformly consistent.

Appendix C Justification of the FMR procedure

In this section, we derive error rates for the ML estimator $\widehat{r}^{[\text{-}k]}(\cdot)$ resulted from the FMR procedure introduced in Section 3.2. Assume that the learning algorithm $\mathscr{L}$ attains the same strong convergence and rate properties as those for $\widehat{m}^{[\text{-}k]}(\cdot)$ in Assumption ML1, i.e. for each $j\in\{1,2,\ldots,K\}$ and $k\in\{1,2,\ldots,K\}$ :

	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X},a\in\mathcal{A}}\|\widehat{M}^{[\text{-}k,\text{-}j]}(a,\boldsymbol{x})-M_{0}(a,\boldsymbol{x})\|=$	$\displaystyle o_{\mathbb{P}}(1);\quad\\|\widehat{M}^{[\text{-}k,\text{-}j]}(\cdot)-M_{0}(\cdot)\\|_{P,2}=o_{\mathbb{P}}(n^{-1/4});$
	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X}}\|\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{x})-a_{0}(\boldsymbol{x})\|=$	$\displaystyle o_{\mathbb{P}}(1);\quad\\|\widehat{a}^{[\text{-}k,\text{-}j]}(\cdot)-a_{0}(\cdot)\\|_{P,2}=o_{\mathbb{P}}(n^{-1/4});$
	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X}}\|\widehat{t}^{[\text{-}k]}(\boldsymbol{x})-t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\boldsymbol{x})\|=$	$\displaystyle o_{\mathbb{P}}(1);\quad\\|\widehat{t}^{[\text{-}k]}(\cdot)-t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\cdot)\\|_{P,2}=o_{\mathbb{P}}(n^{-1/4}),$

where $t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\boldsymbol{x})=\mathrel{\mathop{\mathchar 58\relax}}\mathbb{E}\left[{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(a,\boldsymbol{x})\}\Big{|}\boldsymbol{X}=\boldsymbol{x},\widehat{M}^{[\text{-}k,\text{-}j]}(\cdot)\right]$ . We justify as follows that

\sup_{\boldsymbol{x}\in\mathcal{X}}|\widehat{r}^{[\text{-}k]}(\boldsymbol{x})-r_{0}(\boldsymbol{x})|=o_{\mathbb{P}}(1);\quad\|\widehat{r}^{[\text{-}k]}(\cdot)-r_{0}(\cdot)\|_{P,2}=o_{\mathbb{P}}(n^{-1/4}).

First, since ${\rm logit}$ is a smooth function, it is not hard to show that

\sup_{\boldsymbol{x}\in\mathcal{X},a\in\mathcal{A}}|{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(a,\boldsymbol{x})\}-{\rm logit}\{M_{0}(a,\boldsymbol{x})\}|=o_{\mathbb{P}}(1);\quad\|{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(\cdot)\}-{\rm logit}\{M_{0}(\cdot)\}\|_{P,2}=o_{\mathbb{P}}(n^{-1/4}),

under some mild regularity conditions. Then derive the error rate of $\breve{\beta}^{[\text{-}k]}$ as follows.

		$\displaystyle\|\mathcal{I}_{\text{-}k}\|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(A_{i},\boldsymbol{X}_{i})\}\{A_{i}-\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{X}_{i})\}$
	$\displaystyle=$	$\displaystyle\|\mathcal{I}_{\text{-}k}\|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}{\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}\{A_{i}-a_{0}(\boldsymbol{X}_{i})\}$
		$\displaystyle+\|\mathcal{I}_{\text{-}k}\|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}[{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(A_{i},\boldsymbol{X}_{i})\}-{\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}]\{A_{i}-a_{0}(\boldsymbol{X}_{i})\}$
		$\displaystyle+\|\mathcal{I}_{\text{-}k}\|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}{\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}\{a_{0}(\boldsymbol{X}_{i})-\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{X}_{i})\}$
		$\displaystyle+\|\mathcal{I}_{\text{-}k}\|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}[{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(A_{i},\boldsymbol{X}_{i})\}-{\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}]\{a_{0}(\boldsymbol{X}_{i})-\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{X}_{i})\}$
	$\displaystyle=$	$\displaystyle\|\mathcal{I}_{\text{-}k}\|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}{\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}\{A_{i}-a_{0}(\boldsymbol{X}_{i})\}$
		$\displaystyle+\\|\widehat{a}^{[\text{-}k,\text{-}j]}(\cdot)-a_{0}(\cdot)\\|_{P,2}+\\|{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(\cdot)\}-{\rm logit}\{M_{0}(\cdot)\}\\|_{P,2}+O_{\mathbb{P}}(n^{-1/2})$
	$\displaystyle=$	$\displaystyle\|\mathcal{I}_{\text{-}k}\|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}{\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}\{A_{i}-a_{0}(\boldsymbol{X}_{i})\}+o_{\mathbb{P}}(n^{-1/4})+O_{\mathbb{P}}(n^{-1/2}),$

under some mild regularity conditions. Similarly, we have

\displaystyle|\mathcal{I}_{\text{-}k}|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}\{A_{i}-\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{X}_{i})\}^{2}=|\mathcal{I}_{\text{-}k}|^{-1}\sum_{j=1}^{K}\sum_{i\in\mathcal{I}_{\text{-}k,j}}\{A_{i}-a_{0}(\boldsymbol{X}_{i})\}^{2}+o_{\mathbb{P}}(n^{-1/4})+O_{\mathbb{P}}(n^{-1/2}).

And consequently, by Proposition 1 and

\breve{\beta}^{[\text{-}k]}=\frac{\sum_{i\in\mathcal{I}_{\text{-}k}}{\rm logit}\{M_{0}(A_{i},\boldsymbol{X}_{i})\}\{A_{i}-a_{0}(\boldsymbol{X}_{i})\}}{\sum_{i\in\mathcal{I}_{\text{-}k}}\{A_{i}-a_{0}(\boldsymbol{X}_{i})\}^{2}}+o_{\mathbb{P}}(n^{-1/4})=\beta_{0}+o_{\mathbb{P}}(n^{-1/4}).

Then by Assumption REG that $\beta_{0}$ and $|a_{0}(\boldsymbol{x})|$ are bounded and recall $\widehat{a}^{[\text{-}k]}(\boldsymbol{x})=K^{-1}\sum_{j=1}^{K}\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{x})$ , the estimator $\widehat{r}^{[\text{-}k]}(\cdot)$ given by equation (8) satisfies that:

		$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X}}\|\widehat{r}^{[\text{-}k]}(\boldsymbol{x})-r_{0}(\boldsymbol{x})\|$
	$\displaystyle\leq$	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X}}\|\widehat{t}^{[\text{-}k]}(\boldsymbol{x})-t_{0}(\boldsymbol{x})\|+\|\breve{\beta}^{[\text{-}k]}-\beta_{0}\|\|a_{0}(\boldsymbol{x})\|+\|\breve{\beta}^{[\text{-}k]}\|\|\widehat{a}^{[\text{-}k]}(\boldsymbol{x})-a_{0}(\boldsymbol{x})\|$
	$\displaystyle\leq$	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X},j}\|\widehat{t}^{[\text{-}k]}(\boldsymbol{x})-t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\boldsymbol{x})\|+\|t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\boldsymbol{x})-t_{0}(\boldsymbol{x})\|+o_{\mathbb{P}}(1)=o_{\mathbb{P}}(1),$

where $\sup_{\boldsymbol{x}\in\mathcal{X},j}|t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\boldsymbol{x})-t_{0}(\boldsymbol{x})|=o_{\mathbb{P}}(1)$ is a consequence of $\sup_{\boldsymbol{x}\in\mathcal{X},a\in\mathcal{A}}|{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(a,\boldsymbol{x})\}-{\rm logit}\{M_{0}(a,\boldsymbol{x})\}|=o_{\mathbb{P}}(1)$ . And

		$\displaystyle\\|\widehat{r}^{[\text{-}k]}(\cdot)-r_{0}(\cdot)\\|_{P,2}$
	$\displaystyle\leq$	$\displaystyle\\|\widehat{t}^{[\text{-}k]}(\cdot)-t_{0}(\cdot)\\|_{P,2}+\|\breve{\beta}^{[\text{-}k]}-\beta_{0}\|\\|a_{0}(\boldsymbol{x})\\|_{P,2}+\|\breve{\beta}^{[\text{-}k]}\|\\|\widehat{a}^{[\text{-}k]}(\boldsymbol{x})-a_{0}(\boldsymbol{x})\\|_{P,2}$
	$\displaystyle=$	$\displaystyle\max_{j\in\{1,\ldots,K\}}\\|t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\cdot)-t_{0}(\cdot)\\|_{P,2}+o_{\mathbb{P}}(n^{-1/4})$
	$\displaystyle=$	$\displaystyle\max_{j\in\{1,\ldots,K\}}\\|\mathbb{E}[{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(A,\boldsymbol{X})\}\|\boldsymbol{X}=\boldsymbol{x}]-\mathbb{E}[{\rm logit}\{M_{0}(A,\boldsymbol{X})\}\|\boldsymbol{X}=\boldsymbol{x}]\\|_{P,2}+o_{\mathbb{P}}(n^{-1/4})$
	$\displaystyle\leq$	$\displaystyle\max_{j\in\{1,\ldots,K\}}\\|{\rm logit}\{\widehat{M}^{[\text{-}k,\text{-}j]}(\cdot)\}-{\rm logit}\{M_{0}(\cdot)\}\\|_{P,2}+o_{\mathbb{P}}(n^{-1/4})=o_{\mathbb{P}}(n^{-1/4}).$

Thus, $\widehat{r}^{[\text{-}k]}(\boldsymbol{x})$ satisfies Assumption ML1 under our assumption for the learning tasks of $\mathscr{L}$ .

Appendix D Numerical implementation of the HD approach

We present and demonstrate the implementation procedure of our HD approach mentioned in Remark 5 that uses LASSO instead of the dantzig equation and modifies the construction procedures to make it solvable using the R-package RCAL. Let $G(u)=\int g(u)du$ , and recall that $\widetilde{\boldsymbol{\gamma}}$ is some initial estimator obtained through $\ell_{1}$ -regularized logistic regression for $Y$ versus $\{A,\boldsymbol{X}\}$ and $\widehat{\psi}(\boldsymbol{x})={\rm expit}(-\boldsymbol{x}^{\intercal}\widetilde{\boldsymbol{\gamma}})$ . Then we fit

\min_{\boldsymbol{\alpha}\in\mathbb{R}^{p}}n^{-1}\sum_{i=1}^{n}(1-Y_{i})\widehat{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}}}\left\{-A_{i}\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\alpha}+G(\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\alpha})\right\}+\lambda_{\alpha}\|\boldsymbol{\alpha}\|_{1},

(A6)

to obtain $\widehat{\boldsymbol{\alpha}}$ . It is not hard to see that the KKT (or subgradient) condition of (A6) is equivalent to the $\ell_{\infty}$ -constraint in (4). And when the link function of $g(\cdot)$ is identity (liner model) or ${\rm expit}(\cdot)$ (logistic model), (A6) can be solved using the R-package glmnet with proper specification of the sample weights, i.e. $(1-Y_{i})\widehat{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}}}$ . Then we solve

n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\beta A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}=0,

to obtain a preliminary estimator $\widetilde{\beta}$ . It can be shown that when either $r(\boldsymbol{x})$ or $m(\boldsymbol{x})$ is correctly specified, the estimator $\widetilde{\beta}$ should approach $\beta_{0}$ at the rate $O_{\mathbb{P}}\{(s\log p/n)^{1/2}\}$ , i.e. the $\ell_{2}$ errors of $\widetilde{\boldsymbol{\gamma}}$ and $\widehat{\boldsymbol{\alpha}}$ . So $\widetilde{\beta}$ provides a good enough approximation of $\beta_{0}$ that can be used for the $\ell_{1}$ -regularized (weighted) calibration regression (Tan, 2020b, ) to estimate $\boldsymbol{\gamma}$ :

{\rm min}_{\beta\in\mathbb{R},\boldsymbol{\eta}\in\mathbb{R}^{p}}n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}}}g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\{Y_{i}e^{-\widetilde{\beta}A_{i}-\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma}}+(1-Y_{i})(\widetilde{\beta}A_{i}+\boldsymbol{X}_{i}^{\intercal}\boldsymbol{\gamma})\}+\lambda_{\gamma}\|\boldsymbol{\gamma}\|_{1}.

(A7)

Again, KKT condition of (A7) corresponds to the $\ell_{\infty}$ -constraints in (5), though they are not always imposing the same moment conditions: when the nuisance model $r(\boldsymbol{x})$ is misspecified, $\widetilde{\boldsymbol{\gamma}}$ and $\widehat{\boldsymbol{\gamma}}$ typically have different limits. We use R-package RCAL to solve (A7) with the response taken as $Y_{i}$ , regressors as $\boldsymbol{X}_{i}$ , sample weight $\widehat{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widetilde{\boldsymbol{\gamma}}}g^{\prime}(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})$ and offset $\widetilde{\beta}A_{i}$ for each $i$ . Denoting the solution of (A7) as $\widehat{\boldsymbol{\gamma}}$ , we finally obtain the estimator $\widehat{\beta}_{\sf\scriptscriptstyle HD}$ by solving

n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})e^{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\widehat{\boldsymbol{\gamma}})}\{Y_{i}e^{-\beta A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}=0.

(A8)

Here the final estimating equation is asymptotically equivalent to the second row of (5) only when $r(\boldsymbol{x})$ is correctly specified ( $\widetilde{\boldsymbol{\gamma}}$ and $\widehat{\boldsymbol{\gamma}}$ have the same limiting values). When $r(\boldsymbol{x})$ is misspecified, the orthogonal score function used in (A8), denoted by $h^{\prime}(\boldsymbol{D},\beta_{0},\eta)$ , is not the same as the $h(\boldsymbol{D},\beta_{0},\eta)$ used in the main text. We shall point out that this does not hurt the Neyman orthogonality of $h^{\prime}(\boldsymbol{D},\beta_{0},\eta)$ . It is because that when $r(\boldsymbol{x})$ is misspecified but $m(\boldsymbol{x})$ is correct (by our model assumption, at least one nuisance model need to be correct), $\partial_{r}h^{\prime}(\boldsymbol{D};\beta_{0},\bar{\eta})[r-\bar{r}]$ is naturally satisfied due to the correctness of $m(\boldsymbol{x})$ and $\partial_{m}h^{\prime}(\boldsymbol{D};\beta_{0},\bar{\eta})[m-\bar{m}]$ is satisfied according to the KKT (moment) condition of (A7). When $m(\boldsymbol{x})$ is misspecified but $r(\boldsymbol{x})$ is correct, $\partial_{m}h^{\prime}(\boldsymbol{D};\beta_{0},\bar{\eta})[m-\bar{m}]$ is naturally satisfied and (A8) is asymptotically equivalent with

n^{-1}\sum_{i=1}^{n}\widehat{\psi}(\boldsymbol{X}_{i})\{Y_{i}e^{-\beta A_{i}}-(1-Y_{i})e^{\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\gamma}}}\}\{A_{i}-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}=0,

as $\widetilde{\boldsymbol{\gamma}}$ and $\widehat{\boldsymbol{\gamma}}$ approach the true $\boldsymbol{\gamma}_{0}$ and the second order errors is asymptotically negligible. So $\partial_{r}h^{\prime}(\boldsymbol{D};\beta_{0},\bar{\eta})[r-\bar{r}]$ is satisfied by (A7). Thus, our modified construction procedure does not break the theoretical guarantee of $\widehat{\beta}_{\sf\scriptscriptstyle HD}$ .

Appendix E Additional details of numerical experiments

First, we present the mean vector and covariance matrix used to generate $A$ and $\boldsymbol{X}$ in Section 5.1:

(i)

Take $\boldsymbol{\mu}_{1}=(0.4,-0.25,-0.25,0,\ldots,0)$ , $\boldsymbol{\mu}_{0}=\mathbb{0}$ , and

({\boldsymbol{\Sigma}}^{-1})_{ij}=\begin{cases}1.5&{i=j=1}\\ 1.2&{i=j\geq 2}\\ 0.2&{i=1,2\leq j\leq 5\ or\ 2\leq i\leq 5,j=1}\\ 0&\text{else}\end{cases}

(ii)

Take $\boldsymbol{\mu}_{1}=(0.4,-0.25,-0.25,0,\ldots,0)$ , $\boldsymbol{\mu}_{0}=\mathbb{0}$ ,

({\boldsymbol{\Sigma}_{1}}^{-1})_{ij}=\begin{cases}1.5&{i=j=1}\\ 1.2&{i=j\geq 2}\\ 0.2&{i=1,2\leq j\leq 5\ or\ 2\leq i\leq 5,j=1}\\ 0&\text{else}\end{cases}

and

({\boldsymbol{\Sigma}_{0}}^{-1})_{ij}=\begin{cases}1.5&{i=j=1}\\ 1.2&{i=j\geq 2}\\ 0.2&{i=1,2\leq j\leq 5\ or\ 2\leq i\leq 5,j=1}\\ 0.075&{3\leq i\leq 4,j=2\ or\ i=2,3\leq i\leq 4\ or\ i=3,j=4\ or\ i=4,j=3}\\ 0&\text{else}\end{cases}

(iii)

Take the covariance of $\boldsymbol{X}$ as

${\boldsymbol{\Sigma}}_{ij}=\begin{cases}0.5&{i=j}\\ 0.15&{i\leq 4,j\leq 4,i\neq j}\\ 0&\text{else}\end{cases}.$

Then we present the specific choice on the basis functions for data generation in Section 5.2. In specific, we take

	$\displaystyle f_{a}(\boldsymbol{x})$	$\displaystyle=\left\{\frac{1}{1+e^{x_{1}}},\frac{1}{1+e^{x_{2}}},{\rm sin}(x_{3}),{\rm cos}(x_{4}),I(x_{5}>0),I(x_{6}>0),x_{7}x_{8},x_{9}x_{10}\right\}^{\intercal};$
	$\displaystyle\boldsymbol{\zeta}_{a}$	$\displaystyle=(2,-2,1,1,0.5,-0.5,0.2,0.2)^{\intercal},$

for $a_{0}(\boldsymbol{x})$ , the conditional mean of $A$ given $\boldsymbol{X}=\boldsymbol{x}$ . And that

	$\displaystyle f_{r}(\boldsymbol{x})$	$\displaystyle=\left\{x_{1}x_{2}x_{3},x_{4}x_{5},x_{6}^{3},{\rm sin}^{2}(x_{7}),{\rm cos}(x_{8}),\frac{1}{1+x_{9}^{2}},\frac{1}{1+e^{x_{10}}},I(x_{11}>0),I(x_{12}>0)\right\}^{\intercal};$
	$\displaystyle\boldsymbol{\zeta}_{r}$	$\displaystyle=(0.1,0.1,0.1,-0.5,0.5,1,-1,0.25,-0.25)^{\intercal},$

to specify $r_{0}(\boldsymbol{x})$ .

	$\displaystyle\|\Delta_{11}\|\leq$	$\displaystyle n^{-1}\sum_{i=1}^{n}M_{1}^{2}M_{2}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\|\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\|\|\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\|$
	$\displaystyle\leq$	$\displaystyle M_{1}^{2}M_{2}\left[n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\}^{2}\cdot n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}^{2}\right]^{1/2}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{12}\|\leq$	$\displaystyle n^{-1}\sum_{i=1}^{n}M_{1}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\|\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\|\|g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\|$
	$\displaystyle\leq$	$\displaystyle M_{1}\left[n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}^{2}\cdot n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}^{2}\right]^{1/2}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{13}\|\leq$	$\displaystyle n^{-1}\sum_{i=1}^{n}M_{1}M_{2}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\}^{2}=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{14}\|\leq$	$\displaystyle n^{-1}\sum_{i=1}^{n}e^{\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\gamma}}}\|\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\gamma}}-\bar{\boldsymbol{\gamma}})\|\|g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\|=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right),~{}\mbox{similar to }\|\Delta_{12}\|.$

	$\displaystyle\|\Delta_{21}\|$	$\displaystyle=O\left(\|\widehat{\beta}_{\sf\scriptscriptstyle HD}-\beta_{0}\|\left[n^{-1}\sum_{i=1}^{n}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}^{2}\right]^{1/2}\right)=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{22}\|$	$\displaystyle=O\left(\left[n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}\bar{\boldsymbol{\gamma}}}\}\{\widehat{\psi}(\boldsymbol{X}_{i})-\bar{\psi}(\boldsymbol{X}_{i})\}^{2}\cdot n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}\bar{\boldsymbol{\gamma}}}\}\{g(\boldsymbol{X}_{i}^{\intercal}\bar{\boldsymbol{\alpha}})-g(\boldsymbol{X}_{i}^{\intercal}\widehat{\boldsymbol{\alpha}})\}^{2}\right]^{1/2}\right)$
		$\displaystyle=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{23}\|$	$\displaystyle=O\left(n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}\bar{\boldsymbol{\gamma}}}\}\{\boldsymbol{X}_{i}^{\intercal}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})\}^{2}\right)=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right).$

	$\displaystyle\|\Delta_{31}\|=$	$\displaystyle O\left(\|\widehat{\beta}_{\sf\scriptscriptstyle HD}-\beta_{0}\|\left[n^{-1}\sum_{i=1}^{n}\{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\}^{2}\right]^{1/2}\right)=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right);$
	$\displaystyle\|\Delta_{32}\|=$	$\displaystyle O\left(n^{-1}\sum_{i=1}^{n}\{1+e^{\boldsymbol{X}_{i}\bar{\boldsymbol{\gamma}}}\}\{\boldsymbol{X}_{i}^{\intercal}(\widetilde{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*})\}^{2}\right)=O_{\mathbb{P}}\left(\frac{s\log p}{n}\right).$

		$\displaystyle\|\mathbb{E}h(\boldsymbol{D};\beta,\eta)-\mathbb{E}h(\boldsymbol{D};\beta,\{r_{0},m_{0},\psi\})\|$
	$\displaystyle\leq$	$\displaystyle\|\mathbb{E}\psi(\boldsymbol{X})Ye^{-\beta A}\{m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\}\|+\|\mathbb{E}\psi(\boldsymbol{X})(1-Y)e^{r(\boldsymbol{X})}\{m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\}\|$
		$\displaystyle+\|\mathbb{E}\psi(\boldsymbol{X})e^{r(\boldsymbol{X})}(1-Y)\{1-e^{r_{0}(\boldsymbol{X})-r(\boldsymbol{X})}\}\{A-m(\boldsymbol{X}_{0})\}\|$
	$\displaystyle\leq$	$\displaystyle C_{2}(\\|m_{0}(\boldsymbol{X})-m(\boldsymbol{X})\\|_{P,2}+\\|r_{0}(\boldsymbol{X})-r(\boldsymbol{X})\\|_{P,2}+\\|r_{0}(\boldsymbol{X})-r(\boldsymbol{X})\\|_{P,2}^{2})\leq 3C_{2}\zeta_{2,n}.$

	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X},a\in\mathcal{A}}\|\widehat{M}^{[\text{-}k,\text{-}j]}(a,\boldsymbol{x})-M_{0}(a,\boldsymbol{x})\|=$	$\displaystyle o_{\mathbb{P}}(1);\quad\\|\widehat{M}^{[\text{-}k,\text{-}j]}(\cdot)-M_{0}(\cdot)\\|_{P,2}=o_{\mathbb{P}}(n^{-1/4});$
	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X}}\|\widehat{a}^{[\text{-}k,\text{-}j]}(\boldsymbol{x})-a_{0}(\boldsymbol{x})\|=$	$\displaystyle o_{\mathbb{P}}(1);\quad\\|\widehat{a}^{[\text{-}k,\text{-}j]}(\cdot)-a_{0}(\cdot)\\|_{P,2}=o_{\mathbb{P}}(n^{-1/4});$
	$\displaystyle\sup_{\boldsymbol{x}\in\mathcal{X}}\|\widehat{t}^{[\text{-}k]}(\boldsymbol{x})-t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\boldsymbol{x})\|=$	$\displaystyle o_{\mathbb{P}}(1);\quad\\|\widehat{t}^{[\text{-}k]}(\cdot)-t^{[\text{-}k,\text{-}j]}_{{\dagger}}(\cdot)\\|_{P,2}=o_{\mathbb{P}}(n^{-1/4}),$

Double/Debiased Machine Learning for Logistic Partially Linear Model

Abstract

1 Introduction

2 Neyman orthogonal score

Remark 1.

Remark 2.

3 Method

3.1 High dimensional parametric model realization

Remark 3.

Remark 4.

Remark 5.

3.2 Machine learning realization

Proposition 1.

Proof.

Remark 6.

Remark 7.

3.3 Efficiency consideration

4 Asymptotic analysis

Assumption REG (Regularity of estimating equation for β\beta).

4.1 High dimensional (parametric) setting

Assumption HD1 (Model double robustness).

Assumption HD2 (Concentration rate).

Assumption HD3 (Smooth link function).

Assumption HD4 (Risk of the ℓ1\ell_{1}-regularized estimators).

Assumption HD5 (Ultra-sparsity).

Remark 8.

Remark 9.

Remark 10.

Theorem 1.

Remark 11.

4.2 Machine learning (nonparametric) setting

Assumption ML1 (Quality of the ML nuisance estimators).

Remark 12.

Theorem 2.

Remark 13.

5 Simulation study

5.1 High dimensional (parametric) setting

5.2 Machine learning (nonparametric) setting

6 A real example: effect of EC pill on early gestation foetal

Acknowledgements

References

Appendix

Appendix A Proof of Theorem 1

Proof.

Appendix B Proof of Theorem 2

Condition A1 (Moment condition with Neyman orthogonality).

Condition A2 (Regularity of the score and quality of the nuisance estimators).

Proof.

Appendix C Justification of the FMR procedure

Appendix D Numerical implementation of the HD approach

Appendix E Additional details of numerical experiments

Double/Debiased Machine Learning for Logistic Partially
Linear Model

Assumption REG (Regularity of estimating equation for $\beta$ ).

Assumption HD4 (Risk of the $\ell_{1}$ -regularized estimators).