{'fraction_non_alphanumeric': 0.05373797548708433, 'fraction_numerical': 0.0211191952732067, 'mean_word_length': 4.311762417325898, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Double/Debiased Machine Learning for Logistic Partially Linear Model

We propose double/debiased machine learning approaches to infer (at the parametric rate) 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 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 the 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).


Introduction
Consider a logistic partially linear model. Let {(Y i , A i , X i ) : i = 1, 2, . . . , n} be independent and identically distributed samples of Y ∈ {0, 1}, A ∈ R and X ∈ R p . Assume that where expit(·) = logit −1 (·), logit(a) = log{a/(1 − a)} and r 0 (·) is an unknown nuisance function of X. In an experimental or observational study with A taken as the exposure variable, Y being the binary response of interests and X representing the observed confounding variables, parameter β 0 where φ(x) is an estimation of some scalar nuisance function φ(x) affecting the asymptotic efficiency of the estimator, and α and γ are two fixed dimensional nuisance model estimators. Estimator β solved from (2) is doubly robust that it is valid when either r(x) = x ⊺ γ is correctly specified for the logistic model nonparametric component, or m(x) = g(x ⊺ α) is correctly specified for the conditional mean model m 0 (x) = E(A | Y = 0, X = x). Prior to this, the doubly robust semiparametric estimation of odds ratio was built upon p(A | X, Y = 0), the conditional density of A given X and Y = 0 (Chen, 2007;Tchetgen Tchetgen et al., 2010, e.g.), requiring a stronger model assumption 2 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 ≫ 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 β 0 at the rate O p (n −1/2 ) when either r(·) or m(·) is correctly specified. Under the ML framework, the non-linearity and unextractablity of the logit link makes it impossible to naturally estimate r 0 (·) with a learning algorithm of conditional mean 3 as trivially done under the partially linear model 1 Our 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.
2 Chen (2007) and Tchetgen Tchetgen et al. (2010) estimate the conditional distribution of A given X while Tan (2019a) only needs to specify a conditional mean model of A.
3 Without 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. (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 β 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, 2018Knaus, , 2020Yang 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,b); Ning et al. (2020); Dukes and Vansteelandt (2020) constructed ℓ 1 -regularized estimating equations with certain ℓ ∞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 ℓ 1 -regularized regression to correct for the bias of logistic LASSO estimator. They imposed a technical ultrasparsity condition on inverse of the information matrix, which has been criticized as unreasonable and unverifiable . 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 4 . 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 (·) 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 β 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.

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 . . , n and D = {Y, A, X} be a realization of D i . Motivated by equation (2) proposed and studied in Tan (2019a), we define the Neyman orthogonal score as where η = {r(·), m(·), ψ(·)} represents the whole set of nuisance functions. In analog to (2), r(·) and m(·) corresponds to the nonparametric component r 0 (·) defined in (1) and the nuisance model m 0 (x) = E(A | Y = 0, X = x), respectively. And ψ(x) is a nuisance function affecting the asymptotic variance of the estimator that may depend on r(x) and m(x) and actually corresponds to φ(x)e −r(x) with φ(x) defined by (2) 5 .

Method
In this section, we separately present our specific construction procedures for HD parametric and ML nonparametric realization of the debiased/DML estimator for β 0 , based on the Neyman orthogonal score derived in Section 2.
First, we obtain γ as some initial estimator for γ that converges to some limiting parameter γ * equaling to the true model parameter γ 0 when r(x) is correctly specified: r(x) = x ⊺ γ 0 . And let ψ(x) be some estimator of the nuisance function ψ(x) depending on γ with its limiting function beingψ(x), whose choice will be discussed in Section 3.3 later. According to (3), we obtain α through the dantzig moment equation: where λ α is a tuning parameter controlling the regularization bias. Finally, we obtain the nuisance estimator γ and the targeted HD estimator β HD simultaneously from: Let the limits of { α, γ} be {ᾱ,γ}, andη = {r(·),m(·),ψ(·)} wherer(x) = x ⊺γ ,m(x) = g(x ⊺ᾱ ) andψ(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 4. Similar the HD partially linear (or log-linear) setting studied in Dukes and Vansteelandt (2020), estimating equation for the nuisance parameter γ in our framework involves the unknown β. Unlike their construction procedure that plug-in β 0 as every β ∈ R to estimate γ and invert the resulted score-test p-values for interval estimation of β 0 , we solve for β and γ jointly from (5) with the moment equation for β being doubly robust, as demonstrated in Remark 1. Compared with their method, ours is more friendly in computation and implementation, additionally provides n −1/2consistent point estimator of β 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 ℓ ∞ -norm constraints in (4) and (5), to obtain the estimators α and γ, which has equivalent theoretical properties as the dantzig equations. 6 Due to the non-convexity 7 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.

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: I 1 , I 2 , . . . , I K of equal size, to assist removing the first order (over-fitting) bias through concentration. Then the cross-fitted estimating equation for β is constructed as where . . , n} \ I k and independent of the samples in I k . Now we present the specific construction procedures for r [-k] (·) and m [-k] . . , n} with some response R i and covariates C i , and outputs an estimation of E[R i | C i , i ∈ I]. We outline as follows our approach utilizing L to estimate the nuisance functions. Corresponding to the definition of m 0 (·), it can be estimated by: Compared to the partially linear setting in Chernozhukov et al. (2018a), estimation of r 0 (·) with L is more sophisticated since it is defined through a non-linear form: P(Y = 1 | A, X) = expit{β 0 A + r 0 (X)}. One could modify some ML approaches, e.g. neural network 8 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 L to estimate r 0 (·). Our method is motivated by a simple proposition: Proof. For any β ∈ R, we have Further randomly each split I -k into K folds I -k,1 , . . . I -k,K , of equal size and denote by I -k,-j = I -k \ I -k,j . Motivated by Proposition 1, we first estimate the "full" model M 0 (A, X) with I -k,-j as: Then we fit the (cross-fitted) least square regression to obtain: as an estimator approaching β 0 at certain rate typically larger than n −1/2 . Then r 0 (·) could be identified through r 0 ( Then the estimator of r 0 (·) is given by: Alternatively, one can estimate r 0 (·) through inspired by the moment condition that is sufficient to identify r 0 (·): We refer the estimation step forβ [-k] and r [-k] (·) introduced above as "refitting", and the whole procedure as FMR since we "refit" the least square problem (7) and ML models L to estimate r 0 (·) with the initially estimated full model logit{M 0 (A i , X i )} as a psuedo-outcome. Finally, we solve (6) based on η [-k] to obtain the DML estimator β ML .
Remark 6. We further use cross-fitting in FMR to avoid over-fitting of the models M [-k,-j] (·) and a [-k,-j] (·) when they are used to obtain the estimatorsβ [-k] Remark 7. The FMR implicitly assumes that L should perform similarly well on different learning objects with the covariates set as either X or (A, X ⊺ ) ⊺ . 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 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).

Efficiency consideration
The nuisance function ψ(·) in our framework is included and chosen in consideration of estimation efficiency. Tan (2019a) proposed and studied two options for φ(·) used and defined in (2), with the corresponding function ψ(·) taken as: It was shown that when both nuisance models are correctly specified, the estimator solved with the weight ψ opt (·) achieves the minimum asymptotic variance. However, computation of ψ opt (·) involves numerical integration with respect to X given Y = 0, making it sometimes inconvenient to implement. So Tan (2019a) proposed a simplified but reasonable choice ψ simp (x) defined as above that is obtained by evaluating ψ opt (x) at β 0 = 0. In the following theoretical and numerical studies,

Asymptotic analysis
, Ω(α n ) and Θ(α n ) represent the sequences growing at a smaller, equal/smaller, larger, equal/larger and equal rate of α n , respectively. And let o P , O P , ω P , Ω P and Θ P be the corresponding rates with probability approaching 1 as n → ∞. Let X ⊆ R p be the domain of X. First, we introduce the regularity condition for β 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 β HD and β ML in Sections 4.1 and 4.2 respectively.
Assumption HD1 (Model double robustness). At least one of the following conditions hold: Assumption HD2 (Concentration rate). It holds that To accommodate the notations of both HD and ML, we usem(·) andr(·) to represent the limiting models defined as g(x ⊺ᾱ ) and x ⊺γ under the HD setting and just the true models m0(·) and r0(·) under ML.
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 0 by Remark 3. Then their maximum norms can be controlled by O P {(log p/n) 1/2 } as assumed in HD2, when the covariates 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 ℓ 1 -regularized estimation introduced in Candes et al. (2007) (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, 202 Note that (4) and (5) involves the estimators γ or α obtained beforehand. This will require some additional effort on removing the "plug-in" errors of γ or α when deriving the risk rates for α and γ, 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∈{1,...,n} |g(X ⊺ i α)| = O P (1) imposed in HD4 is not a standard assumption but is mild and very likely to hold: sup x∈X |g(x ⊺ᾱ )| = O(1) according to Assumption REG so we only need g(x ⊺ α) − g(x ⊺ᾱ ) to be O 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 β HD in Theorem 2 and its proof in Appendix A.
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 . 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 P(Y = 1 | A, X) = expit{β 0 A + X ⊺ γ 0 } is correctly specified. As is argued by Xia et al. (2020), assuming the information matrix of the logistic model has an ultra-sparse 10 inverse is crucial to ensure the desirable parametric rate of the debiased LASSO estimator for β 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 E . This assumption has two advantages. First, it accommodates nonlinear link function g(·) 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: and A | X, Y = 0 follows a gaussian linear model with the coefficient α 0 determined by Σ −1 . Therefore, our sparsity assumptions on α 0 and γ 0 actually assumes the data generation parameters Σ −1 and µ 1 − µ 0 to be sparse, which are more explainable and verifiable in practice.

Machine learning (nonparametric) setting
Define that f (·) Q,q =: f (U ) Q,q =: { |f (u)| q dQ(u)} 1/q for any real number q > 0, function f (·), random variables U and probability measure Q. Let P denote the probability measure of the observed D. We assume that K = Θ(1) and introduce the following assumption.
Assumption ML1 (Quality of the ML nuisance estimators). For each k ∈ {1, 2, . . . , K}, 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 (·) and m 0 (·) are uniformly consistent and their mean squared errors (MSE) achieve the rate o 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 r [-k] (·) and m [-k] (·) to be o P (n −1/2 ). In Appendix C, we provide justification for our proposed FMR procedure to derive that the resulted r [-k] (·) satisfies Assumption ML1 as long as the learning algorithm 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 β ML in Theorem 2 with its proof found in Appendix B.
achieves the same strong convergence and rate properties as r [-k] (·) under Assumption ML1. While generally speaking, uniform consistency of ψ [-k] (·) is sufficient for the desirable conclusion in Theorem 2 so our framework accommodates more flexible choices on ψ(x), for example, ψ opt (x) as introduced in Section 3.3. We demonstrate this point during the proof in Appendix B.

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.
(iii) First generate X ∼ N (0, Σ) with Σ ∈ R p×p given in Appendix E. Then we generate A given X from a gaussian linear model with unit variance and conditional mean We realize configurations (i)-(iii) with the sample size n = 1000, 1500 or 2000 separately and the dimension of X fixed as p = 200. Under all these settings, we specify the nuisance models as: r(X) = X ⊺ γ and m(X) = X ⊺ α. Then both nuisance models are correctly specified under (i), only m(X) is correctly specified under (ii), and only r(X) is correct under (iii). Note that we could not extract the explicit form of m 0 (X) under (iii) since A is generated conditional on X without fixed Y = 0. While we still expect the linear model m(X) = X ⊺ α is misspecified under (iii) as there is non-linear terms included in E(A|X). Implementing details of our HD approach are presented in Appendix D. Specifically, all the tuning parameters in ℓ 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 β 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 β 0 (= 0.5) when n = 1000, and at most 12% and 4% of the β 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 ±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.

Machine learning (nonparametric) setting
To study our proposed method under the ML setting, we let Σ ∈ R p×p with Σ ii = 1; Σ ij = 0.2 for i = j, and generate N (0, Σ) random vectors and truncate their entries by (−2, 2) to obtain X. We then generate A from the gaussian model given X with unit variance and conditional mean a 0 (X) = ζ ⊺ a f a (X) where f a (x) is a non-linear basis function of x including various types of effects (interaction, indicator and trigonometric function, etc.), defined in Appendix E and ζ a represents its loading coefficients also given in Appendix E. Finally, we set β 0 = 1, r 0 (X) = ζ ⊺ r f r (X) (see Appendix E), and generate Y with P(Y = 1 | A, X) = expit{β 0 A + r 0 (X)}. We fix p = 20 and set n = 1000 or n = 2000 separately as two settings.
To estimate the nuisance function r 0 (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 L , we consider four ML methods and a hybrid method of the ML models introduced as follows.
(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 β 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. 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 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 (1) denote the indicator for the status of early gestation foetal death in each individual record and Y (2) indicate giving new birth (pregnant and did not incur foetal death). Assume that where r (1) 0 (·) and r (2) 0 (·) are two unknown functions. We are interested in inferring the two parameters β (1) 0 and β (2) 0 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 β (1) 0 , 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 (1) and Y (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 (1) and Y (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 β (1) 0 (abortion) has n (1) = 5, 824 samples. While we only take a subset with n (2) = 10, 000 samples for β (2) 0 so that our algorithms will not require excessive computation times.
For our HD approach, we let X be the p = 175 dimensional bases joining Z, all the interaction terms of Z and the three-dimensional natural splines of all the continuous variable in Z. And we specify the nuisance functions as m(x) = expit(x ⊺ α (ℓ) ), and r(x) = x ⊺ γ (ℓ) for j = 1, 2. We take X = 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 11 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 β (1) 0 and β (2) 0 , respectively. For β (1) 0 , log-OR of early gestation foetal death to EC pill, point estimations of all methods are positive and around −0.175 ± 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: β (1) 0 = 0 at level 0.05. This is because NN outputs slightly worse prediction model for A | 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 β (2) 0 as well. For β (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 12 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 Z and perform our analysis from the individual level granting us to have larger sample sizes. 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.
As sup i∈{1,...,n} |X ⊺ i ( γ −γ)| = O P (1) by Assumption HD4, there exists M 1 = O(1) such that with probability approaching 1, And by Assumptions REG and HD4, there exists M 2 = Θ(1) that sup i∈{1,...,n} |A i − g(X ⊺ iᾱ )| + |A i − g(X ⊺ i α)| ≤ M 2 . Consequently, by Assumptions HD4 and boundness of ψ(·), we have By Assumptions HD2 and HD4, Thus, we have |∆ 1 | = O P (s log p/n). For ∆ 2 , we have Again using Assumptions REG and HD4, there exists M 3 = Θ(1) such that And by Assumption HD3 and the mean value theorem, for each i, there exits t i lying between X ⊺ i α and X ⊺ iᾱ such that These combined with (A1), (A2) and Assumptions REG and HD4 lead to that And by Assumptions HD2 and HD4, So we also have |∆ 2 | = O P (s log p/n). For ∆ 3 , we have Again using the mean value theorem and the fact that Then by (A2), (A3), and Assumptions REG and HD4, we have And by Assumptions HD2 and HD4, Thus, we have |∆ 3 | = O P (s log p/n), and by Assumption HD5, which leads to that This combined with Remark 1 that Eh(D; β 0 ,η) = 0 under Assumption HD1, the regularity Assumption REG, and Theorem 5.21 of Van der Vaart (2000), finally finishes proving Theorem 1.

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(D; β, η), specifically, Conditions A1 and A2 presented as follows.
Condition A1 (Moment condition with Neyman orthogonality). It holds that: (a) Eh(D; β 0 , η 0 ) = 0 and B contains an interval of length Θ(n −1/2 log n) centred at Condition A2 (Regularity of the score and quality of the nuisance estimators). It holds that: (a) η [-k] belongs to the realization set T n for each k ∈ {1, 2, . . . , K}, with probability approaching 1 where T n satisfies η 0 ∈ T n and conditions given as follows; (b) The space of β, B is bounded and for each η ∈ T n , the functional space F η = {h(·; β, η) : β ∈ B} is measurable and its uniform covering number satisfies: there exists positive constant R = Θ(1) and ν = Θ(1) such that where F η (·) is a measurable envelope function for F η : F η (D) ≥ |h(D; β, η)| for all D and β ∈ B, and there exists q > 2 such that F η (·) P,q = O(1); (c) there exists sequence τ n :  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: | β ML −β 0 | = o P (τ 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 uniformly for β satisfying |β − β 0 | ≤ τ n . Note that our modified two assumptions are still sufficient for deriving these results. Since Eh(D; β 0 , {r 0 , m 0 , ψ}) = 0 holds for all ψ (see our Remark 1), there is actually no need to consider Eh(D; β, {r 0 , m 0 , ψ}) − Eh(D; β, {r 0 , m 0 , ψ 0 }) when deriving | β ML − β 0 | = o P (τ n ). While for the Step 3 of Chernozhukov et al. (2018a), one can instead handle the second order error of again using that Eh(D; β 0 , {r 0 , m 0 , ψ}) = 0 holds for all ψ, and then remove the remaining error: On the other hand, this modification essentially reduces our requirement on the quality of ψ [-k] (·), as mentioned in Remark 13. As will be seen from the proof below, to fulfill the modified Condition A2 (c), we do not require ψ [-k] (·) −ψ(·) P,2 = o P (τ n ) as one shall when following the original version of Chernozhukov et al. (2018a) but only needs ψ [-k] (·) to be uniformly consistent (though the former one may still be justifiable in our case since we takeψ(x) = expit{−r(x)}). We now verify Conditions A1 and A2 based on our Assumptions REG and ML1.

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) = g(u)du, and recall that γ is some initial estimator obtained through ℓ 1 -regularized logistic regression for Y versus {A, X} and ψ(x) = expit(−x ⊺ γ). Then we fit to obtain α. It is not hard to see that the KKT (or subgradient) condition of (A6) is equivalent to the ℓ ∞ -constraint in (4). And when the link function of g(·) is identity (liner model) or expit(·) (logistic model), (A6) can be solved using the R-package glmnet with proper specification of the sample weights, i.e. (1 − Y i ) ψ(X i )e X ⊺ i γ . Then we solve to obtain a preliminary estimator β. It can be shown that when either r(x) or m(x) is correctly specified, the estimator β should approach β 0 at the rate O P {(s log p/n) 1/2 }, i.e. the ℓ 2 errors of γ and α. So β provides a good enough approximation of β 0 that can be used for the ℓ 1 -regularized (weighted) calibration regression (Tan, 2020b) to estimate γ: Again, KKT condition of (A7) corresponds to the ℓ ∞ -constraints in (5), though they are not always imposing the same moment conditions: when the nuisance model r(x) is misspecified, γ and γ typically have different limits. We use R-package RCAL to solve (A7) with the response taken as Y i , regressors as X i , sample weight ψ(X i )e X ⊺ i γ g ′ (X ⊺ i α) and offset βA i for each i. Denoting the solution of (A7) as γ, we finally obtain the estimator β HD by solving Here the final estimating equation is asymptotically equivalent to the second row of (5) only when r(x) is correctly specified ( γ and γ have the same limiting values). When r(x) is misspecified, the orthogonal score function used in (A8), denoted by h ′ (D, β 0 , η), is not the same as the h(D, β 0 , η) used in the main text. We shall point out that this does not hurt the Neyman orthogonality of h ′ (D, β 0 , η). It is because that when r(x) is misspecified but m(x) is correct (by our model assumption, at least one nuisance model need to be correct), ∂ r h ′ (D; β 0 ,η)[r −r] is naturally satisfied due to the correctness of m(x) and ∂ m h ′ (D; β 0 ,η)[m−m] is satisfied according to the KKT (moment) condition of (A7). When m(x) is misspecified but r(x) is correct, ∂ m h ′ (D; β 0 ,η)[m −m] is naturally satisfied and (A8) is asymptotically equivalent with as γ and γ approach the true γ 0 and the second order errors is asymptotically negligible. So ∂ r h ′ (D; β 0 ,η)[r −r] is satisfied by (A7). Thus, our modified construction procedure does not break the theoretical guarantee of β HD .