Selective Machine Learning of the Average Treatment Effect with an Invalid Instrumental Variable
Baoluo Sun Yifan Cui Eric Tchetgen Tchetgen
stat.ME






Victor Chernozhukov

Selective Machine Learning of the Average Treatment Effect with an Invalid Instrumental Variable
    Baoluo Sun stasb@nus.edu.sg 

       Department of Statistics and Data Science

       National University of Singapore
       Yifan Cui cuiyf@zju.edu.cn 

       Center for Data Science

       Zhejiang University
       Eric Tchetgen Tchetgen ett@wharton.upenn.edu 

       Department of Statistics and Data Science

       The Wharton School, University of Pennsylvania
    December 30, 2023
========================================================================================================================================================================================================================================================================================================================================================================================




Instrumental variable  methods have been widely used to identify causal effects in the presence of unmeasured confounding. A key identification condition known as the exclusion restriction states that the
instrument cannot have a direct effect on the outcome which is not mediated by the
exposure in view. In the health and social sciences, such an assumption is often not credible. To address this concern, we consider identification conditions of the population average treatment effect with an invalid instrumental variable which does not satisfy  the exclusion restriction, and derive the efficient influence function targeting the identifying functional under a nonparametric observed data model. We propose a novel multiply robust locally efficient estimator of the average treatment effect that is consistent in the union of  multiple parametric nuisance models, as well as a multiply debiased machine learning estimator for which the nuisance parameters are estimated using generic machine learning methods,  that effectively exploit various forms of linear or nonlinear structured sparsity in the nuisance parameter space. When one cannot be confident that any of these machine learners is consistent at sufficiently fast rates  to ensure √n-consistency for the average treatment effect, we introduce new criteria for selective machine learning which leverage the multiple robustness property in order to ensure small bias. The proposed methods are illustrated through extensive simulations  and a data analysis evaluating the causal effect of 401(k) participation on savings.



Average treatment effect, Exclusion restriction, Instrumental variable, Machine learning, Multiple robustness




§ INTRODUCTION



One of the main concerns with drawing causal inferences from observational data is the inability to categorically rule out the existence of unobserved factors  that are associated with both the exposure and outcome variables. The instrumental variable (IV) method is  widely used in the health and social sciences for identification and estimation of causal effects under potential unmeasured confounding <cit.>.   A valid IV is a pre-exposure variable that is (a) associated with treatment, (b) independent of any unmeasured confounder of the exposure-outcome relationship, and (c) has no direct causal effect on the outcome which is not fully mediated by the exposure. While the IV approach has a longstanding tradition in econometrics going back to the original works of <cit.> and <cit.> in the context of linear structural modeling, <cit.>, <cit.>, <cit.> and <cit.> formalized the approach under the potential outcomes framework <cit.> which allows one to nonparametrically define the causal estimands of interest and  clearly articulate assumptions needed to identify this effect; see recent reviews provided by <cit.>, <cit.>,  <cit.> and <cit.>.  The efficient score for the target estimands of interest in the nonparametric IV model satisfies the so called Neyman orthogonality condition <cit.>, which translates to reduced local sensitivity with respect to nuisance parameters. This allows for √n-consistent estimation of the causal estimands of interest even when the complexity of the nuisance parameter space  is no longer tractable by standard empirical process methods (e.g. Vapnik-Chervonenkis and Donsker classes) <cit.>, which represents a significant advancement in the use of machine learning methods for causal inference. 

While (b) may be ensured partly through the randomization of the IV either by design or through some natural or quasi-experiments,   the exclusion restriction (c) is not always credible  in observational studies as it requires extensive understanding of the causal mechanism by which each potential IV influences the outcome <cit.>. In randomized controlled studies with non-compliance, treatment assignment  may have a direct effect on the outcome if double-blinding is either absent or compromised, therefore rendering it invalid as an IV for the effects of treatment actually taken <cit.>.  Throughout, we shall refer to an invalid IV as a potential IV for which exclusion restriction (c) is violated. In response to this concern, there has been growing interest in the development of statistical methods
to detect and account for violation of the exclusion restriction <cit.>, primarily in a system of linear structural equation models. To the best of our knowledge, to date there has been no published work on the population average treatment effect (ATE) as a nonparametric functional targeted with an invalid IV, which prevents the use of  data-adaptive approaches such as machine learning methods for estimation. In this paper, we provide a novel, general set of sufficient conditions under which the ATE is nonparametrically identified despite the IV being invalid, without a priori restricting the nuisance parameters including the model for the conditional treatment effect given observed covariates. In the absence of covariates, identification and inference reduces to a setting studied recently by <cit.>.  Our work in this paper considerably broadens the scope of inference by allowing for potentially high dimensional covariates,  which is far more challenging than what prior literature has considered. 

For inference about the ATE, we pursue two distinct strategies for modeling the nuisance parameters: the first using standard parametric models, while the second leverages modern machine learning. In the former case we propose a multiply robust locally efficient estimator of the ATE which remains consistent under a union of multiple models, each of which restricts a separate subset of parameters indexing the observed data likelihood through low-dimensional parametric specifications. When one cannot be confident that any of these dimension-reducing models is correctly specified, we propose  flexible machine learning of nuisance parameters by selecting a learner for each nuisance parameter from an ensemble of highly adaptive candidate machine learners such as random forests, Lasso or post-Lasso and gradient boosting trees. Building upon recent work by <cit.> and <cit.>, the second main contribution of this paper is to introduce a novel framework for selective machine learning based on minimization of a certain cross-validated quadratic pseudo-risk which embodies the multiple robustness property. The proposed approach ensures that selection of a machine learning algorithm for a given nuisance function is made to minimize
bias of the ATE estimator associated with a suboptimal choice of machine learning algorithms to estimate the other nuisance functions. Our selective machine learning framework can be generally used for making inferences about a finite-dimensional functional defined on semiparametric models which admit multiply robust estimating functions; examples include multiply robust estimation in the context of longitudinal measurements with nonmonotone missingness <cit.>, randomized trials with drop-outs <cit.>, statistical interactions <cit.>, causal mediation analysis <cit.>, instrumental variable analysis <cit.> and causal inference leveraging negative controls <cit.>. 

The rest of the article is organized as follows.  In Section <ref>, we introduce the invalid IV model and provide formal identification conditions for the ATE in this setting.  We present semiparametric estimation methods in Section <ref>, and discuss the use of flexible machine learning of nuisance parameters in Section <ref>. We evaluate the finite-sample performance of these proposed methods through extensive simulation studies in Section <ref> and illustrate the approach with an application  to estimate the causal effect of 401(k) retirement programs on savings using data from the  Survey of Income and Program Participation in Section <ref>. We conclude in Section <ref> with a brief discussion.



§ PRELIMINARIES



Suppose that (O_1,...,O_n) are independent and identically distributed observations of O = (Y, A,
Z,X), where Y is an outcome variable, A is a binary treatment variable encoding the presence
(A = 1) or absence of treatment (A = 0), Z is a binary instrument and X is a set of measured baseline covariates.   To formally define the causal estimands of interest under the potential outcomes framework <cit.>, let Y(z,a) denote the potential outcome that would be observed had the instrument and exposure been set to the level z and a respectively, and let A(z) denote the potential exposure if the instrument would take value z.  We  make the fundamental causal inference assumptions of (i) no interference between units and (ii) no multiple versions of the instrument and treatment; (i) and (ii) are collectively also known as the stable-unit-treatment-value assumption (SUTVA) described in <cit.>. The potential outcomes are related to the observed data via the consistency assumptions  Y=Y(z,a) if Z=z and A=a, and A=A(z) if Z=z. We also define Y(a) to be the potential outcome had only the treatment been set to level a, which is related via the consistency assumption Y(a)=(1-Z)Y(0,a)+ZY(1,a):= Y(Z,a).

Figure <ref>(a) gives causal graph representations <cit.> of the invalid IV model considered in this paper. We assume that U contains all unmeasured common causes of A and Y, such that  conditional on (Z,X,U), the effect of A on Y is unconfounded.



    Y(z,a) A|Z=z,X,U,  for all z,a∈{0,1}.


We will also assume that Z is essentially randomized by design or through some natural experiments within strata of X <cit.>.



    Y(z,a) Z|U,X  and  Z  U |X,   for all z,a∈{0,1}.


Assumptions <ref> and <ref> may also be read (via d-separation) from the corresponding single-world intervention graph <cit.> in Figure <ref>(b). The prototypical example of an invalid IV model is a randomized study where Z is the treatment assignment while A is the treatment actually administered, which may be influenced by some latent factors U  correlated with Y. If double-blinding is either absent or compromised, then knowledge of Z may influence the post-randomization variable Y directly <cit.>.
 





 

 

 



The ATE in the overall population E{Y(1)-Y(0)}  is arguably the causal parameter of interest in many studies  for policy questions <cit.>, but
it cannot be identified under Assumptions <ref> and <ref> without further restrictions. Much of the invalid IV literature considered structural assumptions primarily in
multiple-IV settings Z=(Z_1,Z_2,...,Z_p)^ which imply the joint semiparametric partially linear model 
 
    E(Y|A,Z,X,U)   = θ^_1 Z+β A +ξ_y(X,U);
     E(A|Z,X,U)   =θ^_2 Z+ξ_a(X,U),

indexed by the parameters β∈ IR, θ_1=(θ_11,...,θ_1p)^∈ IR^p, θ_2=(θ_21,...,θ_2p)^∈ IR^p, and  the confounding effects of the measured and unmeasured confounders on the outcome and treatment are respectively encoded by the measurable and square integrable ξ_y(·) and ξ_a(·), which remain unspecified.  Under the multivariate-IV version of Assumption <ref> and correct specification of the partially linear model (<ref>), the scalar parameter β equals the population ATE. The parameter θ_1 represents direct effects of Z on Y. Identification of β (or equivalently the ATE) under (<ref>)–(<ref>) when θ_1j≠ 0 for some j∈𝒥⊆{1,2,...,p} has been an area of active research, generally by imposing additional restrictions on the nuisance parameter space of {θ_1,θ_2} <cit.>. 



 §.§ Nonparametric identification without exclusion restriction

In this paper, we consider the following generalization of (<ref>).



    E(Y|A,Z,X,U)   = θ_1(X)Z+β(X)A+ξ_y(X,U);
     E(A|Z,X,U)   =θ_2(X)Z+ξ_a(X,U),

where {β(·), θ_1(·),θ_2(·)} are unknown measurable and square integrable scalar functions of the measured covariates.

Following the tradition in the IV literature <cit.>, we focus on the canonical case of binary A and Z; the framework can be extended readily to categorical A and Z. The structural equation (<ref>) models the marginal effect of each scalar Z on the outcome and the treatment which can vary with the value of observed covariates, rather than the joint effects of all available IVs, and therefore represents a significant relaxation of the restrictions in (<ref>).  We follow  the latent IV formulation of <cit.> and formally define the no direct effect assumption or exclusion restriction as 
    E{Y(z,a)|X,U}=E{Y(z^',a)|X,U} for all z,z^',a∈{0,1},
 which explicitly incorporates the unmeasured confounder U <cit.>. Under Assumptions <ref> and <ref>, E{Y(z,a)|X,U}=E{Y(z,a)|A=a,Z=z,X,U}=E{Y|A=a,Z=z,X,U}. Therefore θ_1(X)=E{Y(1,a)|X,U}-E{Y(0,a)|X,U} encodes the population average direct effect on the outcome within levels of (X,U) for a change of the IV's value from 0 to 1 at each treatment level; exclusion restriction is violated in model (<ref>) if θ_1(X)≠ 0 for at least one value in the support of X. In addition, under Assumption <ref> and consistency, E{Y(1)-Y(0)|Z,X,U}=E{Y(1)|A=1,Z,X,U}-E{Y(0)|A=0,Z,X,U}=E{Y|A=1,Z,X,U}-E{Y|A=0,Z,X,U}=β(X) equals the conditional ATE within levels of (Z,X,U). 


A design implication of (<ref>) is that even when Z is randomized, it remains important to measure as many effect modifiers in the outcome and treatment models as possible in the hope that no residual effect modification involving U remains within strata of the measured covariates X. We show in the Appendix that (<ref>) may be  relaxed so that {β(·), θ_1(·),θ_2(·)} varies with U (albeit in restricted ways) even after controlling for X, a setting  also known as essential heterogeneity in the IV literature <cit.>.  Essential heterogeneity is more realistic  in a variety of settings. For example, the choice of medical treatment is likely influenced by idiosyncratic gains from alternative treatment  in the analysis of health-care decisions.  The direct effect of treatment assignment on the outcome may also be influenced by knowledge of such gains if double-blinding is either absent or compromised in randomized studies <cit.>. For these reasons, we focus on (<ref>) for identification and inference, although an alternative identification approach involving the nonlinear multiplicative model 

    log{p_1(X,U)/p_0(X,U)}=θ̃_2(X),

where p_z(X,U):=P(A=1|Z=z,X,U), may be used for binary treatment which rules out essential heterogeneity <cit.>. As pointed out by the reviewers, the function θ_2(X) cannot in general be variation independent of the function ξ_a(X,U) if the resulting treatment conditional mean E(A|Z,X,U) must remain in the unit interval. Nevertheless,  a variation independent parameterization of E(A|Z,X,U) is possible such that the aforementioned dependence can be encoded in a manner compatible with our identifying assumptions, by using the odds product parameterization of <cit.> detailed in Appendix B which is compatible with natural constraints of the data generating mechanism.


 Let β_0(x),  μ_0(z,x):= P(A=1|Z=z,X=x) and ε:=A-μ_0(Z,X) denote the true conditional ATE function, treatment propensity score and the treatment regression residual respectively. In what follows the residual ε serves to tease out the treatment effect via its orthogonality to direct effect component θ_1(X)Z. We show in the Appendix that under Assumptions <ref> and <ref>,

    E(ε Y|Z,X)=β_0(X)Var(A|Z,X)+ρ_0(X),

 where Var(A|Z=z,X=x):= μ_0(z,x){1-μ_0(z,x)} denotes the conditional variance of A within the subpopulation {Z=x,X=x} and ρ_0(x):= E{ε(Y-β_0(X)A)|X=x}. The conditional covariance independence restriction 
    E{ε(Y-β_0(X)A)|Z,X}=Cov{ξ_a(X,U) ,ξ_y(X,U)|Z,X}=ρ_0(X),
 holds almost surely under Assumptions <ref> and <ref>. Therefore the function ρ_0(·) may be interpreted as encoding the degree of stratum-specific unmeasured confounding. Additional regularity conditions on the  observed data law are required for identification of the ATE.
 
 
The true observed data distribution lies in the interior of the nonparametric model ℳ that satisfies Var(A|Z=1,X)-Var(A|Z=0,X)≠ 0 (heteroscedasticity), and π_0(1|X):=P(Z=1|X)∈(c,1-c) for some c ∈ (0,1/2) (positivity), almost surely.
  
 Assumption <ref> consists of observed data restrictions that are  empirically testable.  Heteroscedasticity has been widely used in prior works  as a source of identification in linear structural models without exclusion restrictions <cit.> and represents a strengthening of the traditional IV relevance assumption P(A=1|Z=1,X)-P(A=1|Z=0,X)≠ 0 in the context of binary A and Z, as it further requires P(A=1|Z=1,X)+P(A=1|Z=0,X)≠ 1 to hold almost surely. Positivity ensures that there is overlap in the distribution of baseline covariates X among Z=0 and Z=1 units so that the treatment effect within each level of X can be identified. Equation (<ref>) in conjunction with Assumptions <ref> and <ref> implies that 
  
    ATE=E{β_0(X)}=E{E(ε Y|Z=1,X)-E(ε Y|Z=0,X)/Var(A|Z=1,X)-Var(A|Z=0,X)}.

 The nonparametric representation in (<ref>) appears to be new in literature and has a form similar to the well-known Wald estimand as a ratio of differences between the two instrument groups.  Similar to the Wald estimand, estimation based on (<ref>) may be vulnerable to bias and large variance if  Var(A|Z,X) only 
weakly depends on Z. In recent work,  <cit.> proposed a measure of weak identification  relative to sample size and developed inference under a many weak invalid IVs asymptotic regime, which however requires correct specification of parametric models for all nuisance parameters. The observed data density P(O) with respect to some appropriate dominating measure factorizes as P(Y,A|Z,X)× P(Z|X) × P(X). Evaluation of (<ref>) requires knowledge of the joint density P(Y,A|Z,X). As will be shown below, identification of the ATE may be established based on some but not necessarily all of these factors. We introduce the additional notation τ_0(z,x):= E{Y-β_0(X)A|Z=z,X=x} to simplify presentation for this purpose. 






Under Assumptions <ref>–<ref>, the ATE γ:=E{Y(1)-Y(0)}=E{β_0(X)} is identified in ℳ through the following three representations, each of which involves a distinct set of nuisance parameters:


  Explicit representation (i):  
0=E{φ_1(O;π_0,μ_0)|X}-β_0(X) almost surely, where

    φ_1(O;π_0,μ_0):=(2Z-1)ε Y /π_0(Z|X)(Var(A|Z=1,X)-Var(A|Z=0,X));


  Implicit representation (ii): 
0=E{φ_2(O;π_0,β_0,τ_0)|X} almost surely, where 

    φ_2(O;π_0,β_0,τ_0):=(2Z-1)A{Y-β_0(X)A-τ_0(Z,X)}/π_0(Z|X);


  Implicit representation (iii): 0=E{φ_3(O;μ_0,β_0,ρ_0)|Z,X} almost surely, where   

    φ_3(O;μ_0,β_0,ρ_0):=ε{Y-β_0(X)A}-ρ_0(X).



 In particular, representation (i) provides a generalization of the results in <cit.> and <cit.> which both rely on a priori restrictions on the functional form of the  conditional treatment effect within strata of measured covariates,

    E{φ_1(O;π_0,μ_0)| X=x}=Cov{Z,ε Y|X=x}/Cov{Z,ε A|X=x}=β(x;η),

where η is a finite-dimensional parameter.
G-estimators developed in the context  of additive and multiplicative structural mean models <cit.> may be constructed based on unconditional forms of the equivalent restriction 

    Cov{Z,ε(Y-β(X;η)A)|X=x}=0.

No such restriction is needed in representation (i). Thus in principle we can construct the plug-in estimator γ̂=ℙ_n{φ_1(O;π̂,μ̂)}, where  ℙ_n denotes the empirical mean operator ℙ_n {G(O)} = n^-1∑_i G(O_i) and (π̂,μ̂) are nonparametric first-step estimators of (π_0,μ_0) which consists of conditional mean functions. Because we are not restricting  ℳ except for regularity conditions, nonparametric estimators of γ based on representations (i)–(iii) are in fact asymptotically equivalent with common influence function given in the following Theorem <ref>.



The efficient influence function for estimating γ in ℳ is given by 
    φ_(O;π_0,μ_0,β_0,τ_0,ρ_0)-γ,
 where

    φ_(O;π,μ,β,τ,ρ)    =(2Z-1){ε(Y-β(X)A-τ(Z,X))-ρ(X)
    }/π(Z|X){Var(A|Z=1,X)-Var(A|Z=0,X)} +β( X).

Therefore, the semiparametric efficiency bound for estimating γ in ℳ is 
    E{(φ_(O;π_0,μ_0,β_0,τ_0,ρ_0)-γ)^2}.



In most practical settings, we anticipate that X will generally be of moderate to high dimension relative
to the sample size, as analysts consider a broad collection of covariates and their functional forms in the
hope of capturing the salient features of the confounding effects. In this case, nonparametric
estimators of γ may exhibit poor finite-sample behavior due to the curse of dimensionality <cit.>. Below, we describe two distinct strategies for modeling the nuisance parameters: the first
uses standard parametric models, while the second leverages modern machine learning. 



§ MULTIPLY ROBUST ESTIMATION



Consider the working parametric models {π(z|x;η_1),μ(z,x;η_2),β(x;η_3),τ(z,x;η_4),ρ(x;η_5)} indexed by finite-dimensional parameters η=(η^_1,η^_2,η^_3,η^_4,η^_5)^. A two-step procedure to estimate the nuisance parameters is as follows:


Procedure 1.

(i) Solve the score equation 0=ℙ_n{S(A,Z,X; η_1,η_2)}  to obtain (η̂^_1,η̂^_2)^, where 
    S(A,Z,X; η_1,η_2)= ∂log P(A,Z|X; η_1,η_2 )/∂ (η^_1,η^_2)^.


(ii)  Solve 0=ℙ_n{ G(O;η̂_1,η̂_2,η_3,η_4,η_5 )} to obtain (η̂^_3,η̂^_4, η̂^_5)^, where

    G(O;η):=[[ D_3(X)2Z-1/π(Z|X;η_1){ε(η_2)(Y-β(X;η_3)A-τ(Z,X;η_4))-ρ(X;η_5)};                                 D_4(X){Y-β(X;η_3)A-τ(Z,X;η_4)};                           D_5(X){ε(η_2)(Y-β(X;η_3)A)-ρ(X;η_5)} ]],

and D_j(X) is a user-specified vector function of the same dimension as η_j for j=3,4,5.


Similar to <cit.>, in the following we propose the estimator γ̂_mr= ℙ_n{φ_eff(O;π̂,μ̂,β̂,τ̂,ρ̂)} based on the form of the efficient influence function given in Theorem <ref>, where π̂=π(·;η̂_1 ), μ̂=μ(·;η̂_2 ), β̂=β(·;η̂_3 ), τ̂=τ(·;η̂_4 ) and ρ̂=ρ(·;η̂_5 ). Let η^∗ denote the probability limit of η̂. Because the two-step estimator η̂ may be viewed as solving the joint moment equation 0=ℙ_n{G̃(O;η)} where G̃(O;γ)={S^(A,Z,X;η_1,η_2), G^(O;η)}^ <cit.>, the following result holds by invoking the n^-1/2 asymptotic expansion for η̂-η^∗, allowing for model misspecification <cit.>.




Under standard regularity conditions for method of moments estimation  <cit.>, η^∗ is the unique solution to E{G̃(O;η)}=0. Furthermore, γ̂_mr is a consistent and asymptotically normal (CAN) estimator of γ^∗_mr=E{φ_(O;η^∗)},

    √(n)(γ̂_mr-γ^∗_mr)N(0,Σ),

where Σ=E[{φ_mr(O;η^∗)-γ^∗_mr}^2] and

    φ_mr(O;η^∗)=φ_(O;η^∗) -E{∂φ_(O;η)/∂η^|_η=η^∗}×[E{∂G̃(O;η)/∂η|_η=η^∗}]^-1G̃(O;η^∗).



Let π^∗=π(·;η^∗_1 ), μ^∗=μ(·;η^∗_2 ), β^∗=β(·;η^∗_3), τ^∗=τ(·;η^∗_4 ) and ρ^∗=ρ(·;η^∗_5 ) denote the probability limits under the (possibly misspecified) working models. Based on the three distinct sets of nuisance parameters characterized in Theorem <ref>, the efficient influence function has the multiple robustness property that γ=E{φ_eff(O;π^∗, μ^∗,β^∗,τ^∗,ρ^∗)} if at least one of the following holds: (i) (π^∗, μ^∗)=(π_0,μ_0); (ii) (π^∗,β^∗,τ^∗)=(π_0,β_0,τ_0) and (iii) (μ^∗,β^∗, ρ^∗)=(μ_0,β_0,ρ_0). This suggests that γ̂_mr is a CAN estimator of γ  under one, but not necessarily more than one, of the following three different sets of model
assumptions: 


	
  ℳ_1: models for (π_0,μ_0) are correctly specified;
           
  ℳ_2: models for (π_0,β_0,τ_0) are correctly specified;
           
  ℳ_3: models for (μ_0,β_0,ρ_0) are correctly specified.




γ̂_mr is a CAN estimator of γ in the union model ℳ_union=∪_k=1^3 ℳ_k. Furthermore, γ̂_mr attains the semiparametric efficiency bound in ℳ at the intersection submodel ∩_k=1^3 ℳ_k where all the working models are correctly specified.

Following a theorem due to <cit.>, γ̂_mr can be shown to also attain the semiparametric efficiency bound for  ℳ_union at the intersection submodel ∩_k=1^3 ℳ_k. Because the nuisance parameters in each of ℳ_1, ℳ_2 and ℳ_3 are variation independent of each other, multiply robust estimation gives the analyst three genuine opportunities to obtain valid inferences about γ, even under partial misspecification of the observed data models. 



 §.§ Comparison with existing estimators


Under the particular specification {β(x;η_3),τ(z,x;η_4),ρ(x;η_5)}={0,0,0},  γ̂_mr reduces to the semiparametric plug-in estimator γ̂_1=ℙ_n{φ_1(O;π̂,μ̂)} which is CAN only in ℳ_1. An appealing feature of γ̂_1 is that the nuisance parameters can all be estimated in step (i) of procedure 1 without involving outcome data, and therefore mitigates potential for “data-dredging” exercises <cit.>. However, γ̂_1 is neither multiply robust nor locally efficient.  Furthermore, γ̂_mr is expected to be more efficient than γ̂_1, since the latter fails to incorporate information from (A,Z,X) which may be predictive of the outcome values. Such efficiency considerations are analogous to related results on covariate adjustment in completely randomized experiments with full compliance <cit.>.

Based on moment condition (<ref>),  <cit.> proposed the covariate-adjusted “Mendelian Randomization G-Estimation under No
Interaction with Unmeasured Selection” (MR GENIUS) estimator η̂^g_3 which solves

    0= ℙ_n{D_3(X)(2Z-1)ε(η̂_2)(Y-β(X;η_3)A)/π(Z|X;η̂_1)}.

 It is straightforward to verify that γ̂_g=ℙ_n{β(X;η̂_3^g)} is a CAN estimator of γ under the model assumption
  
	
  ℳ_1^': models for (π_0,μ_0, β_0)  are correct.

 Interestingly, under Assumptions <ref>–<ref> and no unmeasured confounding given (Z,X), i.e., if either U A|Z,X or U  Y|A,Z,X, the G-estimator η̂_3 of <cit.>  solves

    0= ℙ_n{D_3(X)ε(η̂_2)(Y-β(X;η_3)A)},

and γ̂=ℙ_n{β(X;η̂_3)} is a CAN estimator of γ.

 Similar to standard G-estimation, a more efficient MR GENIUS estimator η̃^g_3 may be obtained as the joint solution to

    0= ℙ_n
    {[ D_3(X)(2Z-1)ε(η̂_2)(Y-β(X;η_3)A-τ(Z,X;η_4))/π(Z|X;η̂_1);                          D_4(X)(Y-β(X;η_3)A-τ(Z,X;η_4));                                                         ]},

 where information about the association between (Z,X) and Y is incorporated via an additional working model for τ_0(·). <cit.> considered semiparametric estimation based on moment restrictions similar to (<ref>) but with nonparametric plug-ins for the nuisance parameters (η_1,η_2). The resulting estimator γ̃_g=ℙ_n{β(X;η̃^g_3)} is doubly robust in the union model ℳ^'_1 ∪ℳ_2, which is in turn a submodel of ℳ_union.






§ FLEXIBLE ESTIMATION OF NUISANCE PARAMETERS


With high-dimensional X, various flexible and data-adaptive statistical or machine learning  methods may be adopted to estimate the nuisance parameters η_0=(π_0,μ_0,β_0,τ_0,ρ_0),  including random forests, Lasso, neural nets, boosting or their ensembles. Recent work by <cit.> show that √n-consistent estimation of  γ is possible even when the complexity of the nuisance parameters is not tractable by standard empirical process theory (e.g. Vapnik-Chervonenkis and Donsker classes). Let (I_k)_k=1^K be a K-fold random partition of the observation indices {1,2,...,n}. For each k, let π̂(k), μ̂(k), β̂(k), τ̂(k) and ρ̂(k)  be  learners of the nuisance parameters that are constructed using all observations not in I_k based on the following two-step procedure.


Procedure 2.

(i) Obtain π̂(k) and μ̂(k) by machine learning of the conditional mean functions (π_0,μ_0). 

(ii) Given π̂(k) and μ̂(k), obtain β̂(k), τ̂(k) and ρ̂(k)  sequentially by machine learning based on the following conditional mean relationships: β(x; π,μ)=E{φ_1(O;π,μ)| X=x}, τ(z,x; β)=E{Y-β(X)A| Z=z, X=x} and ρ(x;μ,β)=E{ε(μ)(Y-β(X)A)| X=x}.


The cross-fitted debiased machine learning (DML) estimator of γ is

    γ̂_dml=1/n∑^K_k=1∑_i∈ I_kφ_eff(O_i;π̂(k), μ̂(k), β̂(k), τ̂(k) , ρ̂(k)).

 By definition, the efficient influence function φ_eff(O;η) satisfies the Neyman orthogonality condition <cit.>, as all first order influence functions admit second order bias <cit.>. Under general regularity conditions established by <cit.>, γ̂_dml is CAN if all the nuisance parameters are estimated with mean-squared error rates diminishing faster than n^-1/4. Such rates are achievable  for many highly data-adaptive machine learning  methods, including LASSO <cit.>, gradient boosting trees <cit.>, random forests <cit.> or ensembles of these methods.  We note that in low-dimensional settings, Lemma <ref> shows that γ̂_mr is CAN even when some of the nuisance models is misspecified by invoking the usual n^-1/2 asymptotic expansion <cit.>, which is not applicable when nuisance parameters are estimated via machine learning  methods. Therefore while methods such as DML and CV-TMLE  <cit.> with machine learning remain consistent when various strict subsets of nuisance parameter learners (e.g. π̂ and μ̂) are consistent due to the multiple robustness property of the efficient influence function φ_eff(O;η), they generally require consistent estimation of all nuisance parameters in order to obtain valid confidence intervals.





 §.§ Selective machine learning of multiply robust functionals


The performance of DML estimators is intimately related to the choice of the nuisance parameter learners, even when the latter includes flexible machine learning or other nonparametric data adaptive methods. For this reason, generally one would like to learn adaptively from data and avoid choosing models ex ante. The task of model selection of parametric nuisance models was recently considered by <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.> in specific semiparametric doubly robust estimation settings.  A related strand of work is CV-TMLE which can provide notable improvements by incorporating an ensemble of semiparametric or nonparametric methods. Nonetheless, the above methods primarily focused on optimal estimation of nuisance parameters, but not bias reduction of the functional ultimately of interest. This latter task is considerably more challenging since the risk of a nonparametric functional does not typically admit an unbiased estimator and therefore may not be minimized without excessive error. <cit.> proposed a novel model selection criteria for bias reduction in estimating  nonparametric functionals of interest, based on minimization of a cross-validated empirical quadratic pseudo-risk  in the context of doubly robust estimating functions. In this paper we propose to extend their work to the multiply robust setting.


Consider the collection of candidate parametric or nonparametric learners 
    ℒ={η̂(α)=η̂(α_1,α_2,α_3,α_4,α_5)=(π̂_α_1, μ̂_α_2,β̂_α_3,τ̂_α_4,ρ̂_α_5): 1≤α_j≤ r_j  for j=1,...,5},
 with probability limits {η(α)=η(α_1,α_2,α_3,α_4,α_5): 1≤α_j≤ r_j  for j=1,...,5}. Suppose one of the candidate learners η̂(α̌)∈ℒ is consistent so that η(α̌)=(π_0,μ_0,β_0,τ_0,ρ_0).  The proposed procedure relies crucially on the following  two sets of mean zero implications due to the multiply robust property of the efficient influence function, that for all 1≤α_j≤ r_j, 1≤α^'_j≤ r_j, j∈{1,...,5},


    0   =E{φ̅_eff(α̌_1,α̌_2,α̌_3,α̌_4,α̌_5)-φ̅_eff(α̌_1,α̌_2,α_3,α_4,α_5)};
    
    0   =E{φ̅_eff(α̌_1,α̌_2,α̌_3,α̌_4,α̌_5)-φ̅_eff(α̌_1,α_2,α̌_3,α̌_4,α_5)};
    
    0   =E{φ̅_eff(α̌_1,α̌_2,α̌_3,α̌_4,α̌_5)-φ̅_eff(α_1,α̌_2,α̌_3,α_4,α̌_5)},

and

    0   =E{φ̅_eff(α̌_1,α̌_2,α^'_3,α^'_4,α^'_5)-φ̅_eff(α̌_1,α̌_2,α_3,α_4,α_5)};
    
    0   =E{φ̅_eff(α̌_1,α^'_2,α̌_3,α̌_4,α^'_5)-φ̅_eff(α̌_1,α_2,α̌_3,α̌_4,α_5)};
    
    0   =E{φ̅_eff(α^'_1,α̌_2,α̌_3,α^'_4,α̌_5)-φ̅_eff(α_1,α̌_2,α̌_3,α_4,α̌_5)},

where φ̅_eff(α):= φ_eff(O;η(α)). To ease presentation, we introduce the sets 𝒞={1,2,3,4,5}, 𝒞_1={3,4,5}, 𝒞_2={2,5} and 𝒞_3={1,4}  which index the nuisance learner components. Let α_+k denote the  counters for the nuisance learner components indexed by the elements in 𝒞_k, e.g. α_+1=(α_3,α_4,α_5). Similarly, let α_-k denote the counters for the nuisance learner components  indexed by the elements in 𝒞 but not in 𝒞_k, e.g. α_-1=(α_1,α_2). Then (<ref>) and (<ref>) may be restated more concisely as 

    0=E{φ̅_eff(α̌_-k,α̌_+k)-φ̅_eff(α̌_-k,α_+k)},

and

    0=E{φ̅_eff(α̌_-k,α_+k)-φ̅_eff(α̌_-k,α^'_+k)},

respectively, for all α_+k, α^'_+k∈𝒜_k:={α_+k:1≤α_j≤ r_j, j∈𝒞_k}, k∈{1,2,3}. These mean zero conditions suggest perturbing the learners indexed by α_+k and using some measure of the resulting
spread as a basis for selecting between the  learners indexed by α_-k. Towards this end we introduce two different norms to define the spread or pseudo-risk. The first type is
given by the overall maximum squared bias (i.e., change in the estimated functional) induced by perturbing one distinct set of learners at
a time while holding the remaining ones fixed. For an arbitrary learner η̂(α^∗)∈ℒ,  we define the  minimax pseudo-risk ℛ^(1)(α^∗)=max_k∈{1,2,3}Λ^(1)_k(α), where 

    Λ^(1)_k(α^∗)=max_α_+k∈𝒜_k [E{φ̅_eff(α^∗_-k,α^∗_+k)-φ̅_eff(α^∗_-k,α_+k)}]^2,  for k=1,2,3.

The second type is given
by the sum of three maximum squared bias terms, each capturing the bias induced by
perturbing a distinct set of learners. We define the  mixed minimax pseudo-risk ℛ^(2)(α^∗)=∑_k=1^3 Λ^(2)_k(α^∗), where
   
    Λ^(2)_k(α^∗)=max_α_+k,α^'_+k∈𝒜_k [E{φ̅_eff(α^∗_-k,α_+k)-φ̅_eff(α^∗_-k,α_+k^')}]^2,  for k=1,2,3.

For instance, suppose we have 2 candidate learners for each of the 5 nuisance parameters, i.e., ℒ={η̂(α)=η̂(α_1,α_2,α_3,α_4,α_5): 1≤α_j≤ 2  for j=1,...,5}. Then the  minimax pseudo-risk for the learner η̂(1,1,1,1,1)∈ℒ  is 
 
    ℛ^(1)(1,1,1,1,1)=   max(max_α_3=1,2;α_4=1,2;α_5=1,2 [E{φ̅_eff(1,1,1,1,1)-φ̅_eff(1,1,α_3,α_4,α_5)}]^2,
       max_α_2=1,2; α_5=1,2[E{φ̅_eff(1,1,1,1,1)-φ̅_eff(1,α_2,1,1,α_5)}]^2,
       max_α_1=1,2; α_4=1,2[E{φ̅_eff(1,1,1,1,1)-φ̅_eff(α_1,1,1,α_4,1)}]^2),

and its mixed minimax pseudo-risk is 
  
    ℛ^(2)(1,1,1,1,1)=   max_α_3=1,2; α_4=1,2; α_5=1,2
    α^'_3=1,2; α^'_4=1,2; α^'_5=1,2 [E{φ̅_eff(1,1,α_3,α_4,α_5)-φ̅_eff(1,1,α^'_3,α^'_4,α^'_5)}]^2 
       +max_α_2=1,2; α_5=1,2 
    α^'_2=1,2; α^'_5=1,2[E{φ̅_eff(1,α_2,1,1,α_5)-φ̅_eff(1,α^'_2,1,1,α^'_5)}]^2
       +max_α_1=1,2; α_4=1,2 
    α^'_1=1,2; α^'_4=1,2  [E{φ̅_eff(α_1,1,1,α_4,1)-φ̅_eff(α^'_1,1,1,α^'_4,1)}]^2.

 The pseudo-risks for the remaining 2^5-1 learners in ℒ are evaluated similarly. The population version of minimax learners are defined as  {_αℛ^(1)(α)} and {_αℛ^(2)(α)}  respectively. 





 

 §.§ Multi-fold cross-validated selection


 We repeatedly split the data into a training set and a validation set S times to avoid overfitting in selecting the minimax learners. For the s-th split where s∈{1,2,...,S}, let {I^m_s}_m=0,1 be a random bipartition of the observation indices {1,2,...,n}. We use the training sample {1≤ i≤ n: i∈ I^s_0} to construct the estimators {η̂(α;s): 1≤α_j≤ r_j  for j=1,...,5} based on procedure 2. For each fixed learner η̂(α^∗)∈ℒ, the validation sample is used to evaluate
 
    Λ^(1)_k(α^∗)   =max_α_+k∈𝒜_k1/S∑_s=1^S[ℙ^1_s{ϕ_s(α^∗;α^∗_-k,α_+k)}]^2,    k=1,2,3,

where  ϕ_s(α;α^'):= φ_eff(O;η̂(α;s))-φ_eff(O;η̂(α^';s)) and ℙ^m_s:= 1/#{1≤ i≤ n: i∈ I^m_s}∑_i∈ I^m_sδ_O_i for m=0,1, with δ_O denoting the Dirac measure. The empirical terms {Λ^(2)_k(α^∗)}_k=1,2,3 may be evaluated similarly. We select the minimizers of the empirical pseudo-risks ℛ^(1)(α)=max_k∈{1,2,3}Λ^(1)_k(α) and ℛ^(2)(α)=∑_k=1^3 Λ^(2)_k(α) as our nuisance parameter learners. Let α̂^(ℓ)=_αℛ^(ℓ)(α) for ℓ=1,2 respectively. The two proposed selective machine learning  (SML) estimators of γ are given by 

    γ̂^(ℓ)_sml=1/S∑_s=1^Sℙ^1_s{φ_eff(η̂(α̂^(ℓ);s))},
 for ℓ=1,2. We provide a high-level Algorithm <ref> for the proposed selective machine learning  procedure in Appendix C.
 


 §.§ Excess risk bound of the proposed selectors


 We derive risk bounds for the empirically selected minimax learners α̂^(1), α̂^(2) and show that their risks are not much bigger than the risks provided by the respective oracle selected learners α^(1)=_α{max_k∈{1,2,3}Λ̇^(1)_k(α)} and α^(2)=_α{∑_k=1^3Λ̇^(2)_k(α)}, where
 
    Λ̇^(1)_k(α^∗)   =max_α_+k∈𝒜_k1/S∑_s=1^S[ℙ^1{ϕ_s(α^∗;α^∗_-k,α_+k)}]^2;
    Λ̇^(2)_k(α^∗)   =max_α_+k,α^'_+k∈𝒜_k1/S∑_s=1^S[ℙ^1{ϕ_s(α^∗_-k,α_+k;α^∗_-k,α^'_+k)}]^2,

and ℙ^1 denotes the true measure of ℙ^1_s. 




Suppose the nuisance parameter learners satisfy the boundedness conditions (i) P(c≤π̂_α_1(1|X)≤ 1-c)=1 and P(|Var(A|Z=1,X;α_2)-Var(A|Z=0,X;α_2)|> 0)=1 for 1≤α_1≤ r_1, 1≤α_2≤ r_2 and some c>0, where Var(A|Z,X;α_2):=μ̂_α_2(Z,X){1-μ̂_α_2(Z,X)}; (ii) P(|β̂_α_3(X)|≤ M)=1,  P(|α̂_α_4(Z,X)|≤ M)=1 and P(|ρ̂_α_5(X)|≤ M)=1 for 1≤α_3≤ r_3, 1≤α_4≤ r_4, 1≤α_5≤ r_5 and some M>0. Then we have that  

    ℙ^0{ℛ^(1)(α̂^(1))}≤ 
    (1+2ϵ)ℙ^0{ℛ̅^(1)(α^(1))}
       +(1+ϵ/n^1/q)(1+ϵ/ϵ)^(2-q)/qClog{1+(r_3r_4r_5)^2(r_2r_5)^2(r_1r_4)^2},

for any ϵ>0, 1≤ q≤ 2, and some constant C, where ℙ^0 denotes the expectation with respect to training data, 
    ℛ^(1)(α̂^(1))=  max_k∈{1,2,3}1/S∑_s=1^S [ℙ^1{ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2;
    ℛ̅^(1)(α^(1))= max_k∈{1,2,3}1/S∑_s=1^S [ℙ^1{ϕ_s(α^(1); α^(1)_-k,α̅_+k)}]^2,

and for k=1,2,3,

    α̃_+k   =_α_+k∈𝒜_k1/S∑_s=1^S[ℙ^1_s{ϕ_s(α̂^(1); α̂^(1)_-k,α_+k)}]^2;
      α̅_+k=_α_+k∈𝒜_k1/S∑_s=1^S[ℙ^1{ϕ_s(α^(1); α^(1)_-k,α_+k)}]^2.

 Analogous results hold for the mixed minimax selected learner,
 
    ℙ^0{ℛ^(2)(α̂^(2))}≤ 
    (1+2ϵ)ℙ^0{ℛ̅^(2)(α^(2))}
       +(1+ϵ/n^1/q)(1+ϵ/ϵ)^(2-q)/q∑_k=1^3 C_klog[1+(r_3r_4r_5)^{1+I(k=1)}(r_2r_5)^{1+I(k=2)}(r_1r_4)^{1+I(k=3)}] ,

for any ϵ>0, 1≤ q≤ 2, and some constants C_1, C_2 and C_3, where I(·) is the indicator function,

    ℛ^(2)(α̂^(2))= 1/S∑_k=1^3  ∑_s=1^S [ℙ^1{ϕ_s(α̂^(2)_-k,α̃_+k; α̂^(2)_-k,α̃^'_+k)}]^2;
    ℛ̅^(2)(α^(2))= 1/S∑_k=1^3  ∑_s=1^S [ℙ^1{ϕ_s(α^(2)_-k,α̅_+k; α^(2)_-k,α̅^'_+k)}]^2,

and for k=1,2,3,

    (α̃_+k,α̃^'_+k)   =_α_+k,α^'_+k∈𝒜_k1/S∑_s=1^S[ℙ^1_s{ϕ_s(α̂^(2)_-k,α_+k ; α̂^(2)_-k,α^'_+k)}]^2;
    
    (α̅_+k,α̅^'_+k)   =_α_+k,α^'_+k∈𝒜_k1/S∑_s=1^S[ℙ^1{ϕ_s(α^(2)_-k,α_+k ; α^(2)_-k,α^'_+k)}]^2.


The bound given in Theorem <ref> extends the risk bound established  in  <cit.> to the multiply robust setting, and shows that  the error incurred by the empirical risk is of order n^-1 for any fixed ϵ if q=1. Therefore, although the proposed selection procedure is able to incorporate both parametric and nonparametric candidate learners, it is  of most interest in machine learning  settings  where the pseudo-risk can be of order substantially larger than O(n^-1), so that the error made in selecting the cross-validated minimax learner is negligible relative to its risk and the proposed selector performs nearly as well as a oracle
selector with access to the true pseudo-risk. It is also of interest in such settings to compare  the proposed SML estimators with machine learning  estimators using ensemble methods such as super learner <cit.> which selects through  cross-validation the optimal combination from a library of candidate learners to estimate each nuisance parameter separately; we investigate their empirical performances via a simulation study in the next section.


The proposed approach is completely agnostic as to whether the collection ℒ includes a consistent  learner of all the nuisance parameters. Indeed if none of them are consistent there is no estimator of γ that can still be consistent, and the proposed  approach is mostly geared towards identifying the learner that minimizes the minimax pseudo-risks for a given data set. Standard machine learning methods such as DML do not have this built-in data-adaptive feature. To illustrate the implications of this selection procedure, we note that the bias of a DML estimator of γ evaluated with the nuisance parameter learner η̂(α)=(π̂_α_1, μ̂_α_2,β̂_α_3,τ̂_α_4,ρ̂_α_5) chosen ex ante is typically of the order  
    O_p{n^-1/2+||π̂_α_1-π_0||_2||ρ̂_α_5-ρ_0||_2+||μ̂_α_2-μ_0||_2(||β̂_α_3-β_0||_2+||τ̂_α_4-τ_0||_2) },

 which  depends crucially  on products of the learners' estimation errors. Because  Λ̇^(2)_1(α) captures the maximum squared bias in the estimated functional induced by
perturbing only the learners indexed by (α_3,α_4,α_5), its minimizer corresponds to learners indexed by (α_1,α_2) with smallest bias. <cit.> provided formal proof of a related result. The  mixed minimax pseudo-risk  represents a natural extension of this idea as the sum of three maximum squared bias terms, each capturing the bias induced by perturbing a distinct set of learners. Due to the dependence across cross-validation samples, formal machine learning  post selection inference is challenging and the subject of ongoing research.




§ SIMULATION STUDIES


In this section, we investigate the finite-sample properties of the proposed estimators  under a variety of settings. Baseline covariates X=(X_1,...,X_5)^T are generated from independent standard uniform distributions. We consider the functional form X^⋆_k=[1+exp{-20(X_k-.5)}]^-1 for k=1,...,5. The unmeasured confounder U is generated from a truncated normal distribution in the interval [-.5,.5] with mean 0 and variance .25+.5X^⋆_1+.15X^⋆_2-.1X^⋆_3-.1X^⋆_4+.1X^⋆_5. Conditional on (U,X), the invalid instrument Z, treatment A and outcome Y are generated from the models

    P(Z=1|U,X)   ={1+exp(-0.8-X^⋆_1+.2X^⋆_2+.2X^⋆_3+.2X^⋆_4-.1X^⋆_5)}^-1;
    
    P(A=1|Z,U,X)   ={1+exp(2-1.5Z-.6X^⋆_1+.2X^⋆_2+.2X^⋆_3+.1X^⋆_4-.1X^⋆_5)}^-1+κ_1 U;
    
    E(Y|A,Z,U,X)   =-2+(2X^⋆_1+.5X^⋆_2+.5X^⋆_3) A+2X^⋆_1+.5X^⋆_2+.2X^⋆_3+.1X^⋆_4+.1X^⋆_5
       =-2Z+κ_2 U,

where (κ_1,κ_2)=(.1,1) and the outcome error term followed standard normal distribution. We are interested in estimating γ=E(2X^⋆_1+.5X^⋆_2+.5X^⋆_3)=1.5 based on the generated data for (Y,A,Z,X). 

                      


 §.§ Semiparametric estimators


We implement the five semiparametric estimators γ̂, γ̂_1, γ̂_g, γ̃_g and γ̂_mr  using  the R package  <cit.>, and evaluate their performances in situations where some models may be misspecified. A particular working model is misspecified when the quadratic functional form  X^⋆⋆_k=(X_k-0.5)^2 is used in place of X^⋆_k, k=1,...,5. Specifically, we report results from the following four scenarios:
 
	
  𝒮_0: All models are correctly specified;
	
  𝒮_1: models for (π_0,μ_0)  are correct, but models for (δ_0,τ_0,ρ_0) are misspecified;
	
  𝒮_2: models for (π_0,β_0,τ_0) are correct, but models for (μ_0,ρ_0) are misspecified;
	
  𝒮_3:   models for (μ_0,β_0,ρ_0)  are correct, but models for (π_0,τ_0) are misspecified.

Table <ref> summarizes the results  based on 1000 repeated simulations with sample size n=2000 or 4000. Standard errors are obtained using the empirical sandwich estimator for generalized method of moments <cit.>.  The g-estimator γ̂ which does not account for unmeasured confounding shows notable bias relative to its standard error, with coverage below nominal level in all scenarios. In agreement with theory, γ̂_1 has negligible bias and coverage proportions close to nominal levels in scenarios {𝒮_j}_j=0,1, γ̂_g only in 𝒮_0, γ̃_1  in  {𝒮_j}_j=0,2, and γ̂_mr in  {𝒮_j}_j=0,1,2,3, confirming its multiple robustness property.  The estimators  γ̂_mr and γ̂_1 perform similarly to each other in terms of absolute bias, variance and coverage in 𝒮_1, but γ̂_mr yields smaller variance than γ̂_1 in 𝒮_0 where all models are correct. 
 
 
 
 




 §.§ Machine learning estimators


We implement the DML estimators γ̂_LASSO, γ̂_RF, γ̂_GBM and γ̂_SL with covariates X and K=2, whereby the nuisance parameters were estimated  with (i) LASSO <cit.>, (ii) classification or regression random forests <cit.>, (iii) gradient boosting machines <cit.> or the ensemble method super learner  based on a library
consisting of (i), (ii) and (iii), using the R packages  <cit.>,  <cit.>,  <cit.> or  <cit.> respectively. In addition, we implement the proposed SML estimators with covariates X by minimizing the empirical quadratic pseudo-risks over the candidate learners {η̂(α): 1≤α_j≤ 3  for j=1,...,5} with covariates X and S=2, whereby each learner component is based on (i), (ii) or (iii). Table <ref> summarizes the results  based on 1000 repeated simulations  with sample size n=2000 or 4000.   Because the functional form of the regressors is misspecified, γ̂_LASSO has noticeable bias, although it has the smallest Monte Carlo standard error. The bias of γ̂_RF decreases with increasing sample size, and becomes negligible at n=4000. The DML estimator γ̂_GBM has considerably large bias
due to outliers when n =2000, but its bias decreases when n=4000. Remarkably, without access to the true underlying data generating mechanism or the performance of individual DML estimators, the proposed SML estimators  nearly attain the minimum absolute bias at n=4000.  The mixed minimax SML estimator γ̂^(2)_sml tends to be more efficient than γ̂^(1)_sml, in agreement with previous simulation results for doubly robust functionals <cit.>.
 
 




§ APPLICATION


The causal relationship between 401(k) retirement programs and savings has been a subject of considerable interest in economics <cit.>. The main concern with  causal inference based on observational data is that program participation is not randomly assigned, but rather are self-selected by individuals. Potential unmeasured confounders U such as individual preferences may affect both program participation and savings. Thus, estimation of the effects of tax-deferred retirement programs may be biased even after controlling for observed covariates <cit.>. <cit.> proposed 401(k) eligibility as an instrument for program participation. If individuals made employment decisions based on income and within jobs classified by income categories, whether or not a firm offers a 401(k) plan can essentially be viewed as randomized conditional on income and other measured covariates since eligibility is determined by employers. However, 401(k) eligibility may also interact  with some unobserved heterogeneity at the firm's level in  influencing savings other than through 401(k) participation <cit.>.

  In this section, we illustrate the proposed methods by reanalyzing the data from the 1991 Survey of Income and Program Participation (n=9,915) used in <cit.>. The treatment variable A is a binary indicator of participation in a 401(k) plan and Z is a binary indicator of 401(k) eligibility. In this dataset, 37% are eligible for 401(k) programs and 26% participated. The outcomes of interest are net financial assets and net non-401(k) financial assets  in 1991 (US dollars). Following <cit.>, <cit.> and <cit.>, the vector of measured covariates X includes an intercept, family size, indicators for marital status, two-earner status, defined benefit pension status, IRA participation status, homeownership status, four categories of number of years of education, five categories of age and seven income categories. We  consider as benchmarks the two models

    E(Y|A,Z,X,U)   =β(X) A+ξ_y(X);
    
    E(Y|A,Z,X,U)   =β(X)A+ξ_y(X,U),

which are special cases of the outcome structural equation in Assumption <ref>. The former model holds when the effect of A on Y is unconfounded conditional on X, and yields the observed data model E(Y|A,X)=β(X) A+τ(X). If we specify the parametric models β(X;η_3)=η_3 and τ(X;η_4)=η^_4 X, the parameters (η_3,η_4)  indexing the conditional mean model E(Y|A,X;η_3,η_4) may be estimated via  ordinary least squares. On the other hand, the latter model in (<ref>) holds in the presence of unmeasured confounding if Z is a valid instrument that satisfies exclusion restriction. The observed data model E(Y|Z,X)=β(X)μ(Z,X)+τ(X) may be estimated using two-stage instrumental variable estimation under an additional model for the propensity score, which we specify as μ(Z,X;η_2)={1+exp(-η^_2 (Z,X^)^)}^-1. We denote the resulting semiparametric ATE estimators as γ̂_ols and γ̂_tsiv respectively. For comparison, we implement the proposed semiparametric estimators under the same parametric models for δ(X) and μ(Z,X), as well as the additional models π(1|X; η_1)={1+exp(-η^_1 X)}^-1, τ(Z,X;η_4)=η^_4 (Z,X^)^ and ρ(X;η_5)=η^_5 X. The results are summarized in Table <ref>.
 
 
	


 §.§ Effect of 401(k) participation on  net financial assets

The point estimate of γ̂_tsiv is noticeably smaller than that of γ̂_ols, which  is consistent with the results in  <cit.> and suggests that unmeasured confounding generates an upward-biased estimate of the effect of 401(k) participation on savings. The proposed estimators yield significant and uniformly positive point estimates that are close to that of γ̂_tsiv, which provide further evidence that 401(k) participation increases net financial assets even when exclusion restriction may be implausible. This similarity between γ̂_tsiv and the proposed estimators is due to the near zero point estimate for the coefficient corresponding to the main effect of Z in η_4, which encodes the direct effect of  401(k) eligibility on net financial assets. Compared to the semiparametric estimators, the proposed SML estimators  implemented with covariates X and whereby each learner component is based on LASSO, random forests or gradient boosting machines yield similar positive effect estimates, but with noticeably larger nominal standard errors. This difference in efficiency between semiparametric and nonparametric data-adaptive estimation agrees with the simulation results in Section <ref>.



 §.§ Effect of 401(k) participation on  net non-401(k)  financial assets
 
 Consistent with the findings in  <cit.>, the  point estimate of γ̂_tsiv is negative and noticeably lower than that of  γ̂_ols. On the other hand, the proposed estimators yield uniformly positive point estimates, which suggests that 401(k) participation does not crowd out  non-401(k) savings, although the estimates are not statistically significant. This difference may be partially explained by the noticeably negative point estimate for the direct effect of 401(k) eligibility on  non-401(k) savings, indicating asset substitution in 401(k) eligible firms which obscured the effect of 401(k) participation.   This may happen, for example, if 401(k) eligible employees have access to improved financial education
and advice on 401(k) plans at the firm level which encourage  asset substitution. 


	





§ DISCUSSION

 

 There are several improvements and extensions for future work. The finite sample performance of the proposed semiparametric estimators can be improved in terms of efficiency <cit.> and bias <cit.>. Efficiency can also potentially be improved by incorporating a priori knowledge such as degree of exclusion restriction violation or using IVs that are known to be valid in conjunction with the invalid ones.  Lastly, multiple invalid weak IVs can be incorporated by adopting the generalized method of moments approach <cit.>.



The authors would like to thank the Action Editor and two anonymous referees for many
constructive comments which greatly improved the paper. Baoluo Sun's work is supported by the National University of Singapore Start-Up Grant R-155-000-203-133. Eric Tchetgen Tchetgen's work is funded by NIH grants R01AI27271, R01CA222147, R01AG065276 and R01GM139926. 








§ APPENDIX A. PROOFS



 








 §.§ Proof of Theorem <ref>


 <cit.> considered the following generalization of Assumption <ref>  in which both the treatment effects and treatment choice may depend on U:
 

 
 Assumption 3^'

    E(Y|A,Z,X,U)   = θ_1(X,U)Z+β(X,U)A+ξ_y(X,U);
     E(A|Z,X,U)   =θ_2(X,U)Z+ξ_a(X,U),

where {β(·), θ_1(·),θ_2(·)} are unknown measurable and square integrable functions of both measured and unmeasured confounders.



This situation is also known as essential heterogeneity in the econometrics literature <cit.>. Following the proof of Lemma 3.1 in <cit.>, the covariate-specific equality 

    E(ε Y|Z=z,X=x)   =β_0(x)Var(A|Z=z,X=x)+ρ_0(x)
       +∑_j=1,2m_0j(z,x)ψ_0j(x)+∑_k=0,1,2m_1k(z,x)ψ_1k(x),

holds under Assumptions <ref> and 3^', where β_0(x):=E{β(X,U)|X=x} and

    ρ_0(x)   :=Cov{ξ_y(U,X),ξ_a(U,X)|X=x}
    ψ_01(x)   :=Cov{β(U,X),ξ_a(U,X)|X=x};
    ψ_02(x)   :=Cov{θ_1(U,X),ξ_a(U,X)|X=x};
    ψ_10(x)   :=Cov{ξ_y(U,X),θ_2(U,X)|X=x};
    ψ_11(x)   :=Cov{β(U,X),θ_2(U,X)|X=x};
    ψ_12(x)   :=Cov{θ_1(U,X),θ_2(U,X)|X=x}.

It is straightforward to verify that equation (<ref>) holds if for all j=1,2 and k=0,1,2

    ψ_0j(X)=0  and ψ_1k(X)=0,

almost surely. The orthogonality condition (<ref>) does not rule out non-linear forms of essential heterogeneity  <cit.>. In particular, if β_0(X,U)=β_0(X) and θ_j(X,U)=θ_j(X) for j=1,2 almost surely (i.e., Assumption <ref> holds), then (<ref>) holds. Assumption <ref> ensures that (<ref>) is well-defined while Assumption <ref> imbues the identifying functional therein with causal interpretation as the ATE. The rest of the proof below follows from (<ref>).


  §.§.§ Proof of explicit representation (i)


    E{φ_1(O;π_0,μ_0)|X}=E(ε Y |Z=1,X)-E(ε Y |Z=0,X)/Var(A|Z=1,X)-Var(A|Z=0,X)=β_0(X).




  §.§.§ Proof of implicit representation (ii)


    E{φ_2(O;π_0,β_0,τ_0)|X}   =E{(2Z-1)ρ_0(X)/π_0(Z|X)| X}=0.




  §.§.§ Proof of implicit representation (iii)


    E{φ_3(O;μ_0,β_0,ρ_0)|Z,X}   =ρ_0(X)-ρ_0(X)=0.


  

 §.§ Proof of Theorem <ref>


We follow closely the semiparametric efficiency theory of <cit.> and <cit.>. Consider a parametric submodel for the law of the observed data, 
    f_t(o)=f_t(y|a,z,x)μ_t(z,x)^a{1-μ_t(z,x)}^1-aπ_t(x)^z{1-π_t(x)}^1-zf_t(x),
 where μ_t(z,x):= P_t(A=1|Z=z,X=x) and π_t(x):= P_t(Z=1|X=x). The score function S_t(o) is given by S_t(y|a,z,x) +S_t(a|z,x)+ S_t(z|x) + S_t(x), where S_t(y|a,z,x)=∂log f_t(y|a,z,x)/∂ t, S_t(a|z,x)=a-μ_t(z,x)/μ_t(z,x){1-μ_t(z,x)}∂μ_t(z,x)/∂ t,  S_t(z|x)=z-π_t(x)/π_t(x){1-π_t(x)}∂π_t(x)/∂ t and S_t(x)=log f_t(x)/∂ t. A representation of the tangent space is therefore given by 
    𝒯={S_t(y|a,z,x)+{a-μ_t(z,x)}ϱ_1,t(z,x)+{z-π_t(x)}ϱ_2,t(x)+ϱ_3,t(x) },
 where {ϱ_j,t(·)}_j=1,2,3 are arbitrary square-integrable functions. Pathwise differentiability follows if we can find a random element G_t(O)∈𝒯 such that it satisfies E_t(G_t)=0 and 
∂γ_t/∂ t=E_t{G_t(O)S_t(O)}, where E_t{h(O)}=∫ h(o) dF_t.  We make use of the following equalities in the proof, that for any arbitrary square-integrable functions ϱ_1(A,Z,X), ϱ_2(Z,X) and ϱ_3(X),

    E_t{ϱ_1(A,Z,X)S_t(Y|A,Z,X)}  = 0; 
          E_t{ϱ_2(Z,X)S_t(A|Z,X)}  =0;    E_t{ϱ_2(X,Z)ε_t }= 0; 
          E_t{ϱ_3(X)S_t(Z|X)}  =0;     E_t [ϱ_3(X){Z-π_t(Z|X)}]= 0;
          E_t{φ_1(O;π_t,μ_t)| X}=β_t(X),

where (<ref>) follows from the proof of Theorem <ref>. We start with representation (i) of Theorem <ref>. To ease notation, let σ^2_t(z,x):=μ_t(z,x){1-μ_t(z,x)} denote the conditional variance. Differentiating the right hand side of γ_t=E_t{φ_1(O;π_t,μ_t)} under the integral with respect to t yields

    ∇ _tγ _t   =   ∇ _tE_t{2Z-1 /
    π_t(Z|X)ε_t Y /{σ _t^2(1,X)-σ _t^2(0,X)}}
       =   E_t{2Z-1 /
    π_t(Z|X)ε_t YS_t( O) /{σ _t^2(1,X)-σ _t^2(0,X)}}
          -E_t{2Z-1 /
    π_t(Z|X)E_t( AS_t(A|Z,X)|Z,X) Y
    /{σ _t^2(1,X)-σ _t^2(0,X)}}
          -E_t{2Z-1 /
    π_t(Z|X)ε_t  Y∇ _t{σ _t^2(1,X)-σ _t^2(0,X)}/{σ _t^2(1,X)-σ _t^2(0,X)} ^2}
          -E_t{2Z-1 /
    π_t(Z|X)ε_t Y S_t(Z|X)  /{σ _t^2(1,X)-σ _t^2(0,X)}}
       :=   ℒ_1-ℒ_2-ℒ_3-ℒ_4.


We consider the terms ℒ_1 to ℒ_4 separately:


    ℒ_1    =    E_t{2Z-1 /
    π_t(Z|X)ε_t YS_t( O) /{σ _t^2(1,X)-σ _t^2(0,X)}
     };
    ℒ_2    =    E{2Z-1 /
    π_t(Z|X)E_t( AS_t(A|Z,X)|Z,X) Y
    /{σ _t^2(1,X)-σ _t^2(0,X)}}
       =    E{2Z-1 /
    π_t(Z|X)E_t( AS_t(A|Z,X)|Z,X) E_t(Y|Z,X)
    /{σ _t^2(1,X)-σ _t^2(0,X)}}
       =    E{2Z-1 /
    π_t(Z|X) AS_t(A|Z,X) E_t(Y|Z,X)
    /{σ _t^2(1,X)-σ _t^2(0,X)}}
       =    E{2Z-1 /
    π_t(Z|X)ε_t S_t(A|Z,X) E_t(Y|Z,X)
    /{σ _t^2(1,X)-σ _t^2(0,X)}}  by (<ref>)
       =    E{2Z-1 /
    π_t(Z|X)ε_t E_t(Y|Z,X) S_t(O)
    /{σ _t^2(1,X)-σ _t^2(0,X)}};   by (<ref>; <ref>)
    ℒ_3    =   E_t{2Z-1 /
    π_t(Z|X)ε_t  Y∇ _t{σ
    _t^2(1,X)-σ _t^2(0,X)}/{σ_t
    ^2(1,X)-σ_t^2(0,X)} ^2}
       =   E_t{2Z-1 /
    π_t(Z|X)ε_t  Y[1-μ_t(1,X)]E_t[AS_t(A|Z=1,X)|Z=1,X ]  /{σ
    ^2_t(1,X)-σ ^2_t(0,X)} ^2}
          -E_t{2Z-1 /
    π_t(Z|X)ε_t  Yμ_t(1,X)E_t[AS_t(A|Z=1,X)|Z=1,X ] /{σ
    ^2_t(1,X)-σ ^2_t(0,X)} ^2}
          +E_t{2Z-1 /
    π_t(Z|X)ε_t  Y[1-μ_t(0,X)]E_t[AS_t(A|Z=0,X)|Z=0,X ] /{σ
    ^2_t(1,X)-σ ^2_t(0,X)} ^2}
          -E_t{2Z-1 /
    π_t(Z|X)ε_t Yμ_t(0,X)E_t[AS_t(A|Z=0,X)|Z=0,X ] /{σ_t
    ^2(1,X)-σ_t^2(0,X)} ^2}
       =   E_t{β_t( X) /{σ_t
    ^2(1,X)-σ_t ^2(0,X)}Z/π_t( Z|X) 
     ε_t σ ^2_t(Z,X) S_t(O)  /μ_t(Z,X)}
          -E_t{β_t( X) /{σ
    ^2_t(1,X)-σ ^2_t(0,X)}Z/π_t( Z|X) 
     ε_t  σ ^2_t(Z,X) S_t(O)/( 1-μ_t(Z,X)) }
          -E_t{β_t( X) /{σ_t
    ^2(1,X)-σ_t ^2(0,X)}1-Z/π_t( Z|X) 
     ε_t σ ^2_t(Z,X) S_t(O)/μ_t (Z,X)}
          +E_t{β_t( X) /{σ_t
    ^2(1,X)-σ_t ^2(0,X)}1-Z/π_t( Z|X) 
     ε_t  σ ^2_t(Z,X) S_t(O)/[ 1-μ_t (Z,X)] }  by (<ref>, <ref>, <ref>)
       =   E_t{2Z-1 /
    π_t(Z|X)β_t( X) /{σ_t
    ^2(1,X)-σ_t ^2(0,X)}
     ε_t  σ ^2_t(Z,X) S_t(O)/μ_t (Z,X)}
          -E_t{2Z-1 /
    π_t(Z|X)β_t( X) /{σ_t
    ^2(1,X)-σ_t ^2(0,X)}
     ε_t σ ^2_t(Z,X) S_t(O)/[ 1-μ_t (Z,X)] };
    ℒ_4   =   E_t{2Z-1 /π_t(Z|X)ε_tYS_t(Z|X)   /{σ
    ^2(1,X)-σ ^2(0,X)}}
       =   E_t{{[ E_t[ 2Z-1 /π_t(Z|X)ε_t /{σ _t^2(1,X)-σ
    _t^2(0,X)}Y | Z,X];   -E_t[ 2Z-1 /π_t(Z|X)ε_t/{σ _t^2(1,X)-σ
    _t^2(0,X)}Y | X] ]
    } S_t(Z|X)}  by (<ref>)
       =   E_t{{[  E_t[ 2Z-1/π_t(Z|X)ε_t/{σ _t^2(1,X)-σ
    _t^2(0,X)}Y | Z,X]; -E_t[ 2Z-1 /π_t(Z|X)ε_t /{σ _t^2(1,X)-σ
    _t^2(0,X)}Y | X] ]
    } S_t(O)}  by (<ref>, <ref>, <ref>)
       =   E_t{{2Z-1 /π_t(Z|X)[σ _t^2(Z,X)β_t( X) /{σ _t^2(1,X)-σ _t^2(0,X)}
    +ρ_t(X)/{σ _t^2(1,X)-σ
    _t^2(0,X)}]-β_t(X)}S_t(O)},

where the last equality holds from the proof of Theorem  <ref> as well as identity (<ref>). Combining the terms ℒ_1 to ℒ_4 yields

    ∇ _tγ_t   =   E{[2Z-1 /π_t(Z|X)ε_t [ Y-β_t( X) μ_t(Z,X) -τ_t(Z,X) ] /{σ
    ^2_t(1,X)-σ ^2_t(0,X)}
          -2Z-1 /π_t(Z|X)σ _t^2(Z,X)/{σ _t^2(1,X)-σ _t^2(0,X)}β_t( X)
    +β_t( X)  
          -(2Z-1)ρ_t(X)/π_t(Z|X){σ _t^2(1,X)-σ
    _t^2(0,X)}
          - 2Z-1/π_t(Z|X)ε_tβ_t( X)[1-μ_t(Z,X)] /{σ
    ^2_t(1,X)-σ ^2_t(0,X)}
          + 2Z-1 /π_t(Z|X)ε_t β_t( X)μ_t(Z,X)/{σ
    ^2_t(1,X)-σ ^2_t(0,X)}]S_t(O)}
       =   E{[2Z-1 /π_t(Z|X)ε_t[ Y-β_t( X) A -τ_t(Z,X) ] /{σ
    ^2_t(1,X)-σ ^2_t(0,X)}
          + 2Z-1 /π_t(Z|X)ε_tβ_t( X)A/{σ
    ^2_t(1,X)-σ_t^2(0,X)}
          -2Z-1 /π_t(Z|X)σ _t^2(Z,X)/{σ _t^2(1,X)-σ _t^2(0,X)}β_t( X)
    +β_t( X)  
          -(2Z-1) ρ_t(X)/π_t(Z|X){σ _t^2(1,X)-σ
    _t^2(0,X)}
          - 2Z-1/π_t(Z|X)ε_t β_t( X)[1-μ_t(Z,X)] /{σ
    ^2_t(1,X)-σ ^2_t(0,X)}]S_t(O)}
       =   E{[2Z-1 /π_t(Z|X)ε_t [ Y-β_t( X) A -τ_t(Z,X)] /{σ
    ^2_t(1,X)-σ ^2_t(0,X)}
          + 2Z-1 /π_t(Z|X)σ _t^2(Z,X)/{σ _t^2(1,X)-σ _t^2(0,X)}β_t( X)
          -2Z-1 /π_t(Z|X)σ _t^2(Z,X)/{σ _t^2(1,X)-σ _t^2(0,X)}β_t( X)
    +β_t( X)  
          -(2Z-1) ρ_t(X)/π_t(Z|X){σ _t^2(1,X)-σ
    _t^2(0,X)}]S_t(O)}
       =   E{[2Z-1 /π_t(Z|X)ε_t [ Y-β_t( X) A -τ_t(Z,X) ] /{σ
    ^2_t(1,X)-σ ^2_t(0,X)}
          -(2Z-1)ρ_t(X)/π_t(Z|X){σ _t^2(1,X)-σ
    _t^2(0,X)} +β_t( X)   ]S_t(O)}
       =   E{[2Z-1/π_t(Z|X)ε_t [ Y-β_t( X) A -τ_t(Z,X) ] /{σ
     _t^2(1,X)-σ  _t^2(0,X)}
          -(2Z-1)ρ_t(X)/π_t(Z|X){σ _t^2(1,X)-σ
    _t^2(0,X)}+β_t( X) -γ_t  ]S_t(O)}
       :=    E[{φ_eff(O;π_t,μ_t,β_t,τ_t,ρ_t)-γ_t}S_t(O)].

We can readily verify that φ_eff(O;π_t,μ_t,β_t,τ_t,ρ_t)-γ_t ∈𝒯. It follows that  
    φ_eff(O;π_0,μ_0,β_0,τ_0,ρ_0)-γ
 is the efficient influence function for estimating γ in ℳ by Theorem 3.1 of  <cit.>.
   

    
    
  

 §.§ Proof of Lemma <ref>


It suffices to show that γ=E{φ_eff(O;π^∗, μ^∗,β^∗,τ^∗,ρ^∗)}=E{φ_eff(O;η^∗)} if at least one of the following holds: (i) Suppose (π^∗, μ^∗)=(π_0,μ_0). Then

    E{φ_eff(O;η^∗)}   =   E[
    2Z-1/π_0(Z|X)ε{Y-β_0(X)A}/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)
          -2Z-1/π_0(Z|X)ε{β^∗(X)-β_0(X)}A/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)
          -2Z-1/π_0(Z|X)ετ^∗(Z,X)/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)
          -2Z-1/π_0(Z|X)ρ^∗(X)/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)+E{β^∗(X)}]
       =   E[
    2Z-1/π_0(Z|X)ρ_0(X)/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)
          -2Z-1/π_0(Z|X)σ^2(Z,X;μ_0){β^∗(X)-β_0(X)}/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)+β^∗(X)]
       =   E[
    2Z-1/π_0(Z|X)ρ_0(X)/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)
          -{β^∗(X)-β_0(X)}+β^∗(X)]
       =   γ.


(ii) Next, suppose (π^∗,β^∗,τ^∗)=(π_0,β_0,τ_0). Then

    E{φ_eff(O;η^∗)}   =   E[
    2Z-1/π_0(Z|X)ε{ Y-β_0(X)A}/σ ^2(1,X;μ^∗)-σ ^2(0,X;μ^∗)
          -2Z-1/π_0(Z|X){μ^∗(Z,X)-μ_0( Z,X)}{ Y-β_0(X)A-τ_0(Z,X)}/σ ^2(1,X;μ^∗)-σ ^2(0,X;μ^∗)
          -2Z-1/π_0(Z|X)ρ^∗(X)/σ ^2(1,X;μ^∗)-σ ^2(0,X;μ^∗)+β_0(X)]
       =   E[
    2Z-1/π_0(Z|X)ρ_0(X)/σ ^2(1,X;μ^∗)-σ ^2(0,X;μ^∗)+β_0(X)]
       =   γ.


(iii) Finally, suppose (μ^∗,β^∗, ρ^∗)=(μ_0,β_0,ρ_0). Then

    E{φ_eff(O;η^∗)}   =   E[
    2Z-1/π^∗(Z|X)ε{ Y-β_0(X)A}-ρ_0(X)/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)
          -2Z-1/π^∗(Z|X)ετ^∗(Z,X)/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)+β_0(X)]
       =   E[
    2Z-1/π^∗(Z|X)ρ_0(X)-ρ_0(X)/σ ^2(1,X;μ_0)-σ ^2(0,X;μ_0)+β_0(X)]
       =   γ.


The last claim in Lemma  <ref> follows by noting that under the intersection submodel ∩_k=1^3ℳ_k,  E{∂φ_eff(O;η)/∂η^|_η=η^∗}=0 by  Neyman orthogonality <cit.>.



 §.§ Proof of Theorem <ref>



To ease notation, let n_m^s:= #{1≤ i≤ n: i∈ I^s_m} for m=0,1. By the definition of the proposed model selector,

    max_k∈{1,2,3}, α^'_+k∈𝒜_k1/S∑_s=1^S[ℙ^1_s{ϕ_s(α̂^(1); α̂^(1)_-k,α^'_+k)}]^2≤max_k∈{1,2,3}, α^'_+k∈𝒜_k1/S∑_s=1^S[ℙ^1_s{ϕ_s(α^(1); α^(1)_-k,α^'_+k)}]^2.

It follows that

    ℛ^(1)(α̂^(1)) ≤max_k∈{1,2,3}1/S∑_s=1^S[ℙ^1_s{ϕ_s(α^(1); α^(1)_-k,α̇_+k)}]^2,

where α̇_+k=_α^'_+k∈𝒜_k1/S∑_s=1^S[ℙ^1_s{ϕ_s(α^(1); α^(1)_-k,α^'_+k)}]^2. By simple algebra, we have that for k∈{1,2,3},

    [^1_s {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2=    1/n_1^s^2∑_i,j[^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i }^1 {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j}] 
       + 2/n_1^s^2∑_i,j [  ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i - ^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i}] ^1 {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j }
       + 1/n_1^s^2∑_i,j [ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i - ^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i}]×
       + [ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j- ^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j}],

where ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i denotes the estimating equation evaluated at i-th observation. Thus,

    [^1_s {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2=  [^1 {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2
       + 2/n_1^s^2∑_i,j [ ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i - ^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i}] ^1 {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j }
        + 1/n_1^s^2∑_i,j [ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i - ^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i}]×
       +[ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j- ^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j}].

The same decomposition holds for [^1_s {ϕ_s(α^(1); α^(1)_-k,α̇_+k)}]^2. By definition of our estimator, for any ϵ>0, we have that

    ℛ^(1)(α̂^(1)) ≤  (1+2ϵ)max_k1/S∑_s=1^S  [^1 {ϕ_s(α^(1); α^(1)_-k,α̇_+k)}]^2
    
    +   { (1+ϵ)max_k1/S∑_s=1^S(   [^1_s {ϕ_s(α^(1); α^(1)_-k,α̇_+k)}]^2-  [^1 {ϕ_s(α^(1); α^(1)_-k,α̇_+k)}]^2 ) 
       -ϵmax_k1/S∑_s=1^S [^1{ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2   }
    
    -   { (1+ϵ)max_k1/S∑_s=1^S (   [^1_s {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2-  [^1 {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2 ) 
        + ϵmax_k1/S∑_s=1^S [^1{ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2    }.

Combined with the decomposition of [^1_s {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2 and [^1_s {ϕ_s(α^(1); α^(1)_-k,α̇_+k)}]^2, we further have that

    ℛ^(1)(α̂^(1)) ≤ (1+2ϵ)max_k1/S∑_s=1^S  [^1 {ϕ_s(α^(1); α^(1)_-k,α̇_+k)}]^2
    
    +    { (1+ϵ)  max_k1/S∑_s=1^S [ 2/n_1^s∑_i  ( ϕ_s(α^(1); α^(1)_-k,α̇_+k)_i- ^1 {ϕ_s(α^(1); α^(1)_-k,α̇_+k)_i})^1 {ϕ_s(α^(1); α^(1)_-k,α̇_+k)}
    
     +    1/(n_1^s)^2∑_i,j ( ϕ_s(α^(1); α^(1)_-k,α̇_+k)_i - ^1  {ϕ_s(α^(1); α^(1)_-k,α̇_+k)_i })×
        ( ϕ_s(α^(1); α^(1)_-k,α̇_+k)_j - ^1  {ϕ_s(α^(1); α^(1)_-k,α̇_+k)_j })] -ϵmax_k1/S∑_s=1^S [^1 {ϕ_s(α^(1); α^(1)_-k,α̇_+k)}]^2  }
    
     -    { (1+ϵ) max_k1/S∑_s=1^S [  2/n_1^s∑_i  ( ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i- ^1 {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i})^1 {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}
    
      +    1/(n_1^s)^2∑_i,j (ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i - ^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_i }) ×
        (ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j - ^1  {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)_j })] + ϵmax_k1/S∑_s=1^S  [^1 {ϕ_s(α̂^(1); α̂^(1)_-k,α̃_+k)}]^2   }.

Note that the only assumption on {I_s^m}_m=0,1 is its stochastic independence of the observations, we omit sup-index s hereinafter.
Because the maximum of sum is at most the sum of maxima, we deal with the first order and second order terms separately.
By Lemma 2.2 in <cit.>, we further have the following bounds for the first order term,

    ^0 [max_k∈{1,2,3}, α, α^'_+k∈𝒜_k{2(1+ϵ)√()/∑_i ( ϕ(α; α_-k,α^'_+k)_i - ^1   {ϕ(α; α_-k,α^'_+k)_i })  
       ×^1   {ϕ(α; α_-k,α^'_+k)}- ϵ√()[^1   {ϕ(α; α_-k,α^'_+k)} ]^2 }] 
    ≤   ^0  (  16(1+ ϵ)/n_1^1/q-1/2log{1+(r_3r_4r_5)^2(r_2r_5)^2(r_1r_4)^2 }×
       max_k∈{1,2,3}, α, α^'_+k∈𝒜_k [||  ϕ(α; α_-k,α^'_+k)^1   {ϕ(α; α_-k,α^'_+k)} ||_∞/^1-1/q
    
    +   (3^1  [ϕ(α; α_-k,α^'_+k) ^1   {ϕ(α; α_-k,α^'_+k)} ]^2 2^(1-q)(1+δ)^(2-q)/ϵ^2-q [ ^1   {ϕ(α; α_-k,α^'_+k)}]^2-2q)^1/q ]  ):= (I),

and

    ^0[max_k∈{1,2,3}, α, α^'_+k∈𝒜_k{-[ 2(1+ϵ)√()/∑_i ( ϕ(α; α_-k,α^'_+k)_i - ^1   {ϕ(α; α_-k,α^'_+k)_i}) ×
       ^1   {ϕ(α; α_-k,α^'_+k)}+ ϵ√()[^1   {ϕ(α; α_-k,α^'_+k)} ]^2 }]  ≤ (I),

where the maximum is taken over sets that fixes α_-k as α and varies α^'_+k, for k∈{1,2,3}. Thus, ℙ^0{ℛ^(1)(α̂^(1))} is further bounded by

    ℙ^0{ℛ^(1)(α̂^(1))}
       ≤ (1+2ϵ)^0(max_k[^1{ϕ(α^(1); α^(1)_-k,α̇_+k)}]^2) +  2/√()(I)
    
     +    ℙ^0[max_k{1+ϵ/n_1^2∑_i,j ( ϕ(α^(1); α^(1)_-k,α̇_+k)_i - ^1  {ϕ(α^(1); α^(1)_-k,α̇_+k)_i})×
        (ϕ(α^(1); α^(1)_-k,α̇_+k)_j - ^1  {ϕ(α^(1); α^(1)_-k,α̇_+k)_j}) 
    
      -   1+ϵ/n_1^2∑_i,j (ϕ(α̂^(1); α̂^(1)_-k,α̃_+k)_i- ^1  {ϕ(α̂^(1); α̂^(1)_-k,α̃_+k)_i})×
        (ϕ(α̂^(1); α̂^(1)_-k,α̃_+k)_j - ^1  {ϕ(α̂^(1); α̂^(1)_-k,α̃_+k)_j})}].

Following Lemmas 4 and 5 of <cit.>, the U-statistics are bounded and we have the following excess risk bound,	

    ℙ^0{ℛ^(1)(α̂^(1))}
       ≤ (1+2ϵ)^0(max_k[^1{ϕ(α^(1); α^(1)_-k,α̇_+k)}]^2) 
    
    +     (1+ϵ)C {( 2M/n_1^2log(1+M (r_3r_4r_5)^2(r_2r_5)^2(r_1r_4)^2 /2 ) )^1/2
       + 2M/n_1log(1+M (r_3r_4r_5)^2(r_2r_5)^2(r_1r_4)^2 /2 ) 
       +  4M^3/2/n_1^3/2log^3/2(1+M (r_3r_4r_5)^2(r_2r_5)^2(r_1r_4)^2 /2+D_0 )+
       4 M^2/n_1^2log^2(1+M (r_3r_4r_5)^2(r_2r_5)^2(r_1r_4)^2 /2+D_1 ) }
    
    +    ^0  (  16(1+ ϵ)/n_1^1/q-1/2log{1+(r_3r_4r_5)^2(r_2r_5)^2(r_1r_4)^2 }×
       max_k∈{1,2,3}, α, α^'_+k∈𝒜_k [||  ϕ(α; α_-k,α^'_+k)^1   {ϕ(α; α_-k,α^'_+k)} ||_∞/^1-1/q
    
    +   (3^1  [ϕ(α; α_-k,α^'_+k) ^1   {ϕ(α; α_-k,α^'_+k)} ]^2 2^(1-q)(1+δ)^(2-q)/ϵ^2-q [ ^1   {ϕ(α; α_-k,α^'_+k)}]^2-2q)^1/q ]  ),
	
where C, M, D_0, and D_1 are some universal constants. Finally, recall that for the term 
    (1+2ϵ)^0(max_k[^1{ϕ(α^(1); α^(1)_-k,α̇_+k)}]^2),
 α̇_+k is  chosen corresponding to α^(1) under measure _s^1. It is further bounded by  (1+2ϵ)ℙ^0{ℛ̅^(1)(α^(1))}, where α̅_+k is chosen corresponding to α^(1) under true measure ^1. This completes the proof for the risk bound of the minimax estimator α̂^(1).  Although details are omitted, the proof for the risk bound of the mixed minimax estimator α̂^(2) is essentially the same. 	
  
      0.2in




§ APPENDIX B.

For binary treatment, the baseline risk ξ_a(X,U)=p_0(X,U) and the risk difference θ_2(X)=p_1(X,U)-p_0(X,U) are variation dependent, where p_z(X,U):=P(A=1|Z=z,X,U). <cit.> showed that the log odds product

    ω_1(X,U):=log[p_0(X,U)p_1(X,U)/{1-p_0(X,U)}{1-p_1(X,U)}],
 is variation independent of ω_2(X):=arctanh{θ_2(X)}, which yields the appropriate linear structural model  
P(A=1|Z,X,U;ω_1,ω_2)=θ_2(X;ω_2)Z+ξ_a(X,U; ω_1,ω_2) for binary treatment in (<ref>), where

    θ_2(X;ω_2)   =tanh{ω_2(X)};
    ξ_a(X,U; ω_1,ω_2)   =e^ω_1(2-θ_2)+θ_2-√({e^ω_1(θ_2-2)-θ_2}^2+4e^ω_1(1-θ_2)(1-e^ω_1))/2(e^ω_1-1).






§ APPENDIX C.





mycommfont

KwInputInput                
KwOutputOutput              








