---
abstract: 'While model selection is a well-studied topic in parametric and nonparametric regression or density estimation, model selection of possibly high dimensional nuisance parameters in semiparametric problems is far less developed. In this paper, we propose a new model selection framework for making inferences about a finite dimensional functional defined on a semiparametric model, when the latter admits a doubly robust estimating function. The class of such doubly robust functionals is quite large, including many missing data and causal inference problems. Under double robustness, the estimated functional should incur no bias if either of two nuisance parameters is evaluated at the truth while the other spans a large collection of candidate models. We introduce two model selection criteria for bias reduction of functional of interest, each based on a novel definition of pseudo-risk for the functional that embodies this double robustness property and thus may be used to select the candidate model that is nearest to fulfilling this property even when all models are wrong. Both selection criteria have a bias awareness property that selection of one nuisance parameter can be made to compensate for excessive bias due to poor learning of the other nuisance parameter. We establish an oracle property for a multi-fold cross-validation version of the new model selection criteria which states that our empirical criteria perform nearly as well as an oracle with a priori knowledge of the pseudo-risk for each candidate model. We also describe a smooth approximation to the selection criteria which allows for valid post-selection inference. Finally, we apply the approach to perform model selection of a semiparametric estimator of average treatment effect given an ensemble of candidate machine learning methods to account for confounding in a study of right heart catheterization in the intensive care unit of critically ill patients.'
author:
- |
    Yifan Cui, Eric Tchetgen Tchetgen\
    Department of Statistics, The Wharton School, University of Pennsylvania
bibliography:
- 'causal.bib'
- 'survtrees.bib'
title: 'Bias-aware model selection for machine learning of doubly robust functionals'
---

[***Keywords:*** Model Selection, Machine Learning, Doubly Robust, Influence Function, Average Treatment Effect, Cross-validation]{}

Introduction {#sec:intro}
============

Model selection is a well-studied topic in statistics, econometrics and machine learning. In fact, methods for model selection and corresponding theory abound in these disciplines, although primarily in settings of parametric and nonparametric regression and density estimation ([@akaike1974new; @BIC1978; @vuong1989likelihood; @zhang1993model; @wand1994kernel; @fan1995jrssb; @ruppert1995; @HALL1996165; @lasso; @yang2000; @birge2001gaussian; @fan2001variable; @rao; @wegkamp2003; @ruppert2003semiparametric; @efron2004least; @zhao2006lasso; @BIRGE2006497; @candes2007dantzig; @celisse2014; @belloni2014] and many others). Model selection methods are far less developed in settings where one aims to make inferences about a finite dimensional, pathwise differentiable functional defined on a semiparametric model. Model selection for the purpose of estimating such a functional may involve selection of an infinite dimensional parameter, say a nonparametric regression for the purpose of more accurate estimation of functional in view, which can be considerably more challenging than selecting a regression model strictly for the purpose of prediction. This is because whereas the latter admits a risk, e.g., mean squared error loss, that can be estimated unbiasedly and therefore can be minimized with small error, the risk of a semiparametric functional will typically not admit an unbiased estimator and therefore may not be minimized without excessive error. This is an important gap in both model selection and semiparametric theory which this paper aims to address.

Specifically, we propose a novel approach for model selection of a functional defined on a semiparametric model, in settings where inferences about the targeted functional involves infinite dimensional nuisance parameters, and the functional of scientific interest admits a doubly robust estimating function. Doubly robust inference has received considerable interest in the past few years across multiple disciplines including Statistics, Epidemiology and Econometrics [@robins1994; @Rotnitzky1998; @Scharfstein1999; @robins2000; @robins2001comment; @unified; @Lunceford2004; @bang2005; @tan2006; @cao2009; @ett2010dr; @funk2011; @rotnitzky2012dr; @han2013; @FARRELL20151; @2015biasreduce; @vermeulen2016adaptive; @Chernozhukov2018; @rotnitzky2019mix; @tan2019; @fulcher2017robust]. An estimator is said to be doubly robust if it remains consistent if one of two nuisance parameters needed for estimation is consistent, even if both are not necessarily consistent. The class of functionals that admit doubly robust estimators is quite rich, and includes estimation of pathwise differentiable functionals in missing data problems under missing at random assumptions, and also in more complex settings where missingness process might be not at random. Several problems in causal inference also admit doubly robust estimating equations, the most prominent of which is the average treatment effect under assumptions that include positivity, consistency and no unmeasured confounding [@Scharfstein1999; @robins2000]. All of these functionals are members of a large class of doubly robust functionals recently studied by [@robins2008HOIF] in a unified theory of first and higher order influence functions.

The literature on double robustness combined with machine learning methods is rapidly expanding [@van2010collaborative; @van2011targeted; @belloni2014res; @FARRELL20151; @belloni2017Econometrica; @robins2017; @Chernozhukov2018; @athey2018jrssb; @van2018targeted; @dukes2018high; @rotnitzky2019mix; @tan2019model]. A well-documented advantage of using doubly robust influence functions is that flexible machine learning or other nonparametric data adaptive methods may generally be used to estimate high dimensional nuisance parameters such that valid inferences may be obtained about the functional of interest provided that estimated nuisance parameters have mean squared error of order smaller than $n^{-1/4}$, which can be considerably slower than converge rates attained by parametric models [@robins2017; @Chernozhukov2018]. As in practice, one cannot be certain that any model is either correctly specified or estimated with small bias, model selection remains important even in the advent of doubly robust estimation including for methods that leverage machine learning. Clearly, the performance of doubly robust semiparametric estimators is intimately related to the performance of estimators of its nuisance parameters, a task towards which model selection is paramount.

This paper aims at the selection of an optimal estimator for the functional $\psi(\theta)$ in a class of doubly robust functionals where $\theta$ is a parameter (possibly infinite dimensional) indexing the observed data law within a semiparametric/nonparametric model. Given a large collection of doubly robust estimators $\Psi_{K}=\{\widehat \psi_k: k=1,\ldots,K\}$ of size $K$ (which may grow with sample size) indexed by candidate estimators of nuisance parameters, we ultimately wish to identify an estimator that minimizes the risk associated with a measurable loss function. A natural choice would be to try to select the estimator that minimizes the mean squared error $E(\widehat \psi_{k,\tk}-\psi)^2$. However, it is clear that this cannot be done empirically in a straightforward fashion, as an unbiased estimator of the mean squared error (even up to a constant shift) is generally not available, so that model selection becomes challenging.

In this paper, we propose two novel model selectors each based on minimization of a certain cross-validated quadratic pseudo-risk for a large class of doubly robust functionals. The proposed pseudo-risk embodies the idea of double robustness: The first kind of pseudo-risk is given by the overall maximum squared bias (i.e., change in the estimated functional) at a given candidate estimator, induced by perturbing one nuisance parameter at the time over candidate models holding the other one fixed; The second proposed pseudo-risk is given by the sum of two maximum squared bias quantities, each capturing the bias induced by perturbing a single nuisance parameter only. As we establish both procedures are guaranteed to recover a consistent estimator for the functional whenever consistent estimators of nuisance parameters are available, with corresponding pseudo-risk converging to zero. However, even when all models are wrong, as in many practical settings where parametric models are used, and therefore all candidate estimators are inconsistent with pseudo-risk bounded away from zero, a minimizer of pseudo-risk nevertheless corresponds to a choice of models that is least sensitive to perturbations, i.e., misspecification of either nuisance parameter. Both selection criteria have a bias awareness property that selection of one nuisance parameter is aware and therefore may be made to compensate for excessive bias due to poor learning of the other nuisance parameter. We find such awareness may be key to bias reduction in context of machine learning of doubly robust functionals.

Our cross-validation scheme is akin to that of [@vanderlann2003cross] and [@vaart2006cv; @vaart2006oracle], who formally established that such a scheme can perform nearly as well as an oracle with access to underlying data generating mechanism, in selecting an optimal estimator in settings such as nonparametric regression or density estimation. In contrast, we aim to perform model selection for a pathwise differentiable functional of such nonparametric regression or density function, and therefore to minimize a risk function for the functional; a different task which generally proves to be more challenging. For each split of the observed sample, a training sample is used to estimate each candidate model of the nuisance parameters. The validation subsample is then used to construct corresponding candidate estimators of functional $\psi$, and subsequently, to estimate the pseudo-risk of each candidate estimator conditional on the training sample. The optimal model is selected by minimizing multi-fold cross-validated pseudo-risk over the set of candidate nuisance models. To our knowledge, this is the first model selection result for doubly robust functionals which aims directly at bias reduction of the functional. Significant amounts of work have been devoted to improving performance of doubly robust estimators [@bang2005; @tan2006; @cao2009; @tan2010; @rotnitzky2012dr; @FARRELL20151; @2015biasreduce; @vermeulen2016adaptive; @van2018targeted; @dukes2018high; @smucler2019unifying; @rotnitzky2019mix; @tan2019; @bradic2019sparsity] from a variety of perspectives, however, none have considered model selection for the underlying functional, over a generic collection of candidate nuisance parameter models that may include classical parametric, semiparametric and nonparametric estimators, as well as modern highly data adaptive machine learning estimators. The task of model selection of parametric nuisance models for specific semiparametric doubly robust problems was recently considered by [@han2013; @chan2013; @HAN2014101; @han2014jasa; @chan2014; @duan2017; @chen2017; @liu2019], although, their goal differs from ours as they aim to select parametric nuisance models that best approximate each nuisance model, which may generally conflict with selecting the nuisance models that minimize a well-defined pseudo-risk of the targeted functional, especially when as often the case in practice, all candidate models are wrong.

A related targeted maximum likelihood learning approach for model selection in functional estimation, known as cross-validated targeted maximum likelihood estimation [@zheng2010asymptotic; @van2011targeted; @van2018targeted] can provide notable improvements on the above methods by allowing the use of an ensemble of semiparametric or nonparametric methods, including modern machine learning for flexible estimation of nuisance parameters, still the ensemble learning is targeted at optimal estimation of nuisance parameters, not bias reduction of the functional ultimately of interest. Another state of the art approach recently proposed incorporates modern machine learning in functional estimation via double debiased machine learning (DDML) [@Chernozhukov2018]; however the approach uses a single machine learning algorithm for estimating each nuisance parameter, and does not leverage model selection targeted at the functional of interest. In comparison, as we will show, our approach ensures that selection of one nuisance model is made to minimize bias due to possible misspecification of the other, such bias awareness for the functional of interest endows the proposed model selection procedure with additional robustness.

The proposed approach is generic, in the sense that it allows the space of candidate models/learners to be quite large ($K_1\times K_2$ of order $c^{n^\gamma}$ for any constants $c>0$ and $\gamma<1$), and arbitrary in the sense of including parametric, semiparametric, nonparametric, as well as modern machine learning highly adaptive estimators. Importantly, our results are completely agnostic as to whether the collection of models includes a correct model for nuisance parameters, in the sense that our procedure will select the nuisance models that optimize our doubly robust criteria. Another aspect in which our approach is generic is that it does not depend on a particular choice of doubly robust estimator of a given functional. In the sense, the approach may be used with say doubly robust targeted maximum likelihood learning to construct an ensemble of doubly robust targeted maximum likelihood estimators, each of which based on different estimators of nuisance parameters. As discussed in Section \[sec:dr\], a very general class of doubly robust functionals are considered here. The purpose of considering a broad class is to demonstrate the flexibility of our method for various functionals that are generally of interest. Several functionals, such as the expected conditional covariance, marginal mean of an outcome subject to missingness as well as the closely related marginal mean of a counterfactual outcome are within our class. As a running example to develop the proposed methodology, throughout we consider the average treatment effect under unconfoundedness as the target of inference.

In settings where all candidate estimators of the functional are regular and asymptotically linear, although not necessarily consistent, we propose a smooth approximation of the proposed criteria which allows for valid post-selection inference. In case the selected model fails to be consistent for the functional of interest, because all candidate models fail to consistently estimate nuisance parameters, valid inference can nevertheless be obtained for the approximate functional that minimizes a population version of the proposed doubly-robust-inspired pseudo-risk function, whenever such approximate functional is well-defined. Confidence intervals can then be constructed either using an estimate of asymptotic variance of smooth selected estimator of the functional based on a standard application of the delta method, or via the nonparametric bootstrap.

The paper is organized as follows: in Section \[sec:dr\], we introduce the general class of doubly robust functionals and give specific examples of interest within this class. In Section \[sec:prelim\], we introduce the problem setting, and demonstrate the main challenge of model selection. Section \[sec:selection\] is devoted to developing our proposed selection criteria. Utilizing these criteria, we propose a general cross-validation scheme to construct empirical minimax functional selectors in Section \[sec:cf\]. In Section \[sec:theory\], we use powerful exponential inequalities for tails of extreme of second order U-statistics to establish a risk bound for the cross-validated minimax and mixed-minimax criteria. The risk bound firmly establishes that our empirical criteria select a pair of nuisance models which performs nearly as well as the pair of models selected by an oracle with access to the law that generated the data. In Section \[sec:simulations\], we present simulation studies to evaluate the performance of the proposed approach in a range of settings. In Section \[sec:softmax\], we describe a smooth approximation to the cross-validated pseudo-risk minimizer which allows for post-selection inferences. In Section \[sec:real\], we illustrate the proposed methods by studying the effectiveness of right heart catheterization in the insentive care unit (ICU) of critically ill patients. Details of proofs are given in the appendices.

A class of doubly robust functionals\[sec:dr\]
==============================================

Suppose we observe $n$ i.i.d. samples $\mathcal O \equiv \{O_i,i=1,\cdots,n\}$ from a law $F_0$, belonging to a model $\mathcal M=\{F_\theta: \theta \in \Theta\}$, where $\Theta$ may be infinite dimensional. We are interested in inference about a functional $\psi(\theta)=\psi(F_\theta)$ on $\cal M$ for a large class of functionals known to admit a doubly robust first order influence function as defined in [@robins2016tr].

Suppose that $\theta = \theta_1 \times \theta_2$, where $\times$ denotes Cartesian product, $\theta_1 \in \Theta_1$ governs the marginal law of $X$ which is a $d$-dimensional subset of variables in $O$, and $\theta_2  \in \Theta_2$ governs the conditional distribution of $O|X$. \[as1\]

An influence function is a fundamental object of statistical theory that allows one to characterize a wide range of estimators and their efficiency. The influence function of a regular and asymptotically linear estimator $\widehat \psi$ of $\psi (\theta)$, $\theta \in \cal M$, is a random variable $IF(\theta)\equiv IF(O;\theta)$ which captures the first order asymptotic behavior of $\widehat \psi$, such that ${n}^{1/2}\{\widehat \psi-\psi(\theta)\}=n^{-1/2} \sum_{i=1}^n IF(O_i;\theta) + o_p(1)$. The set of influence functions of all regular and asymptotically linear estimators of a given functional $\psi(\theta)$ on $\cal M$ is contained in the Hilbert subspace of mean zero random variables $U\equiv u(O;\theta)$ that solve the following equation, $$d(\psi(\theta_t)/dt|_{t=0} =E\{US\},$$ for all regular parametric submodels of $\cal M$, $F_{\theta_t}$, $t \in (-\epsilon,\epsilon)$ with $F(\theta_0)=F_0$, and $S$ the score function of $f(O; \theta_t)$ at $t = 0$ [@newey1990semiparametric; @bickel1993efficient; @van2000asymptotic; @tsiatis2007semiparametric]. Once one has identified the influence function of a given estimator, one knows its asymptotic distribution, and can construct corresponding confidence intervals for the target parameter. We now describe a large class of doubly robust influence functions.

The parameter $\theta_2$ contains components $b:\mathbbm{R}^d\rightarrow \mathbbm{R}$ and $p:\mathbbm{R}^d\rightarrow \mathbbm{R}$, such that the functional $\psi(\theta)$ of interest has a first order influence function $IF(\theta)= H(b,p)-\psi(\theta)$, where $$\begin{aligned}
\label{eq:H}
H(b,p) \equiv b(X)p(X)h_1(O) + b(X)h_2(O) + p(X)h_3(O) + h_4(O),\end{aligned}$$ and $h_1,\ldots,h_4$ are measurable functions of $O$. \[as2\]

$\Theta_{2b} \times \Theta_{2p} \subseteq \Theta_2$, where $\Theta_{2b}$ and $\Theta_{2p}$ are the parameter spaces for the functions $b$ and $p$. Furthermore, the sets $\Theta_{2b}$ and $\Theta_{2p}$ are dense in $L_2(F_0)$ at each $\theta_1\in \Theta_1$, where $L_2(F_0)$ is the Hilbert space of all functions with finite variance. \[as3\]

[@robins2016tr] point out that Assumptions \[as1\]-\[as3\] imply the following double robustness property, $$\begin{aligned}
\label{eq:dr0}
E_\theta[H(b^*,p^*)]-E_\theta[H(b,p)]=E[(b(X)- b^*(X))(p(X)- p^*(X))h_1(O)],\end{aligned}$$ for all $(b^*,p^*)\in \Theta_{2b} \times \Theta_{2p}$. In which case $E[H(b^*,p^*)]=\psi$ if either $b^*=b$ or $p^*=p$. Examples of functionals within this class include:

(Expected product of conditional expectations) Suppose we observe $O=(A,Y,X)$, where $A$ and $Y$ are univariate random variables. Let $\psi(\theta)= E_\theta[p(X)b(X)] $, where $b(X)=E_\theta[Y|X]$ and $p(X)=E_\theta[A|X]$ are a priori unrestricted. In this nonparametric model, the first order influence function of $\psi$ is given by $$IF(\theta)= p(X)b(X)+p(X)\{Y-b(X)\} +b(X)\{A-p(X)\}-\psi(\theta),$$ so $h_1(O)=-1, h_2(O)=A, h_3(O)=Y, h_4(O)=0$.

(Expected conditional covariance) Suppose $O=(A,Y,X)$, where $A$ and $Y$ are univariate random variables. Let $\psi(\theta)= E_\theta [Cov_\theta (Y,A|X)]= E_\theta[AY]- E_\theta[p(X)b(X)]$, where $b(X)=E_\theta[Y|X]$ and $p(X)=E_\theta[A|X]$. In this model, the first order influence function is $$IF(\theta)= AY- \big[ p(X)b(X)-\psi(\theta)+p(X)\{Y-b(X)\} +b(X)\{A-p(X)\} \big]-\psi(\theta),$$ so $h_1(O)=1, h_2(O)=-A, h_3(O)=-Y, h_4(O)=AY$.

As pointed out by [@robins2008HOIF; @robins2016tr; @robins2017], inference about expected conditional covariance is key to obtaining valid inferences about $\beta$ in the widely used semiparametric regression model $E(Y|A,X)=\beta A+b(X)$, where $b(X)$ is unrestricted [@robins2008HOIF].

(Missing at random) Suppose $O=(AY,A,X)$, where $A$ is the binary missing indicator, and $X$ is a $d$-dimensional vector of fully observed continuous covariates. We assume $Y$ is missing at random, i.e., $A{\rotatebox[origin=c]{90}{$\models$}}Y|X$. Let $b(X) = E(Y|A=1,X)$ be the outcome model and $\Pr(A = 1|X)>0$. The parameter of interest $\psi(\theta)$ is the marginal mean of $Y$. In this model, the first order influence function is $$IF(\theta)= Ap(X)\{Y-b(X)\}+b(X)-\psi(\theta),$$ where $p(X) = 1/\Pr(A = 1|X)$. So $h_1(O)=-A, h_2(O)=1, h_3(O)=AY, h_4(O)=0$.

(Missing not at random) We consider the setting in the last example allowing for missing not at random. We assume that $\Pr(A=1|X,Y)=[1+\exp\{-[\gamma(X)+\alpha Y]\} ]^{-1}$, where $\gamma(X)$ is an unknown function and $\alpha$ is a known constant. The marginal mean of $Y$ is again of interest and given by $\psi(\theta)=E_\theta (AY[1+\exp\{-[\gamma(X)+\alpha Y]\} ])$. [@robins2001comment] derived the first order influence function of $\psi$, $$IF(\theta)=A[1+\exp\{-\alpha Y\}p(X)][Y-b(X)]+b(X)-\psi(\theta),$$ where $b(X)=E[ Y\exp\{ -\alpha Y\} |A=1,X]/E[ \exp\{-\alpha Y\}|A=1,X]$ and $p(X)=\exp\{-\gamma(X)\}$. So $h_1(O)= -A \exp\{-\alpha Y\}, h_2(O)=1-A, h_3(O)=AY\exp\{-\alpha Y\}$, $h_4(O)=AY. $

(Average treatment effect) Suppose we observe $O=(A,Y,X)$, where $A$ is a binary treatment taking values in $\{0,1\}$, $Y$ is a univariate response, and $X$ is a collection of covariates. We wish to make inferences about the average treatment effect $E\left\{ Y_{1}-Y_{0}\right\}$, where $Y_{1}$ and $Y_{0}$ are potential outcomes.

Three important assumptions are sufficient for identification of the average treatment effect from the observed data. First, we make the consistency assumption that $Y = Y_A$ almost surely. This assumption essentially states that one observes $Y_a$ only if the treatment $a$ is equal to a subject’s actual treatment assignment $A$. The next assumption is known as ignorability [@10.2307/2335942], which requires that there are no unmeasured confounders for the effects of $A$ on $Y$, i.e., for both $a\in \{0,1\}$, $Y_a {\rotatebox[origin=c]{90}{$\models$}}A|X$. Finally, we assume that $\pi(a|X=x)=\Pr(A = a|X=x)>0$ for $a\in \{0,1\}$ if $f(x)>0$. This positivity assumption basically states that any subject with an observed value of $x$ has a positive probability of receiving both values of the treatment.

Under these three identifying conditions, functional $\psi_0(\theta) = E[E(Y |A = 1, X) - E(Y |A = 0, X)]$ is the average effect of treatment on the outcome. The first order influence function of this functional is $$\begin{aligned}
IF(\theta)=\frac{\left( -1\right) ^{1-A}}{\pi\left( A|X\right) }Y
-\left\{ \frac{\left( -1\right) ^{1-A}}{\pi\left( A|X\right) }E(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}E(Y|a,X)\right\}-\psi_0(\theta).
\label{eq:dr}\end{aligned}$$ In fact, this model has four nuisance parameters. Note that we can rewrite the influence function as $$\frac{A}{\pi(1|X)}\{Y-E(Y|1,X)\}+E(Y|1,X) - \left[\frac{1-A}{\pi(0|X)}\{Y-E(Y|0,X)\}+E(Y|0,X) \right]-\psi_0(\theta).$$ Then $IF(\theta)$ can be viewed as a difference of two influence functions of similar form as the missing at random (MAR), where $p^{(1)}(X)=1/\pi(1|X),$ $b^{(1)}(X)= E(Y|1,X)$, $p^{(2)}(X)=1/\pi(0|X),$ $b^{(2)}(X)= E(Y|0,X)$, $h^{(1)}_1(O)=-A, h^{(1)}_2(O)=1, h^{(1)}_3(O)=AY, h^{(1)}_4(O)=0$, $h^{(2)}_1(O)=-(1-A), h^{(2)}_2(O)=1, h^{(2)}_3(O)=(1-A)Y, h^{(2)}_4(O)=0$.

\[remark:rm2\] [@rotnitzky2019mix] study a more general class of doubly robust influence functions, which admit the following “mixed bias property”: For each $\theta$, there exist functions $c(X;\theta)$ and $d(X;\theta)$ such that for any $\theta'$, $$\begin{aligned}
\psi(\theta') - \psi(\theta) + E_\theta (IF(\theta')) = E_\theta[ s_{ab}(O) \{ c(X,\theta')-c(X,\theta)  \}\{ d(X,\theta') - d(X,\theta) \}  ],
\label{eq:mix}\end{aligned}$$ where $s_{ab}$ is a known function not depending on $\theta$. Note that the selection procedure proposed in Section \[sec:select\] extends to this richer class of doubly robust influence functions which includes both the classes of [@riesz2018] and of doubly robust functionals described above [@robins2016tr]. In fact, all that is required by the proposed approach is the influence function has mean zero when either nuisance parameter is evaluated at the truth. As will be discussed later, the approach can readily be extended to multiply robust influence functions in the sense of [@tchetgentchetgen2012; @wang2018bounded; @Caleb2019; @Shi2019MultiplyRC; @sun2019multiple].

The practical implication of double robustness is that the asymptotic bias of an estimator obtained by solving $\PP_n \widehat {IF}(\widehat \psi)=\PP_n IF(\widehat p, \widehat b,\widehat \psi) = 0$ is guaranteed to be zero provided either but not necessarily both $\widehat p$ is consistent for $p$ or $\widehat b$ is consistent for $b$. Despite this local robustness property, in practice one may be unable to ensure that either model is consistent, and even when using nonparametric models, that the resulting bias is small. For this reason, model selection over a class of candidate estimators may be essential to optimize performance in practical settings.

Challenges of model selection for doubly robust inference\[sec:prelim\]
=======================================================================

Hereinafter, in order to ground ideas, we focus the presentation to the average treatment effect functional of Example 2.5. It is well known that $$\begin{aligned}
\psi_0  &=&E\left[ E\left( Y|A=1,X\right) -E\left( Y|A=0,X\right) \right]  \\
&=&E\left( \frac{\left( -1\right) ^{1-A}}{\pi\left( A|X\right) }Y\right)  \\
&=&E\left(
\begin{array}{c}
\frac{\left( -1\right) ^{1-A}}{\pi\left( A|X\right) }Y
-\left\{ \frac{\left( -1\right) ^{1-A}}{\pi\left( A|X\right) }E(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}E(Y|a,X)\right\}
\end{array}\right).\end{aligned}$$

The first representation is known as outcome regression as it depends on the regression of $Y$ on $\left( A,X\right) $; the second is inverse probability weighting with weights depending on the propensity score [@10.2307/2335942]; the third representation is known as doubly robust as it relies on outcome regression or propensity score model to be correct but not necessarily both. In fact, the doubly robust representation based on the efficient influence function of $\psi _{0}$ is given by Equation , which will be used as our estimating equation for $\psi_0$ for the proposed model selection.

In order to describe the inherent challenges of performing model selection for $\psi_0$, consider the sample splitting scheme whereby a random half of the sample is used to construct $\widehat \pi_k(A|X)$ and $\widehat E_\tk(Y|A,X)$, while the other half is used to obtain the doubly robust estimator $\widehat \psi_{k,\tk}$. Consider the goal of selecting a pair of models $(k, \tk)$ that minimizes the mean squared error $E[(\widehat \psi_{k, \tk}-\psi_0)^2|\text{Training sample}]$ = bias$^2(\widehat \psi_{k,\tk})$ + variance$(\widehat \psi_{k,\tk})$, where bias$^2(\widehat \psi_{k,\tk})$ is given by Equation  As we expect the variance term to be of order $1/n$ conditional on training sample, we may focus primarily on minimizing the squared bias. As no unbiased estimator of $\text{bias}^2(\widehat \psi_{k,\tk})$ exists, minimizing $\text{bias}^2(\widehat \psi_{k,\tk})$ will generally not be possible without incurring excessive bias. Hereafter, for a given split of sample we shall refer to $\arg \min_{k,\tk} \text{bias}^2(\widehat \psi_{k,\tk})$ as “squared bias minimizer”, which depends on the true data generating law (through $\pi(A|X)$ and $E(Y|A,X)$), and therefore may not be accurately estimated even in large samples. In the next section, we propose alternative criteria for selecting an estimator with a certain optimality condition that is nearly attainable empirically.

Recall that consistent estimators of the propensity score and outcome regression are not necessarily contained as candidates for model selection, so the minimal squared bias may not necessarily converge to zero asymptotically; nevertheless, it will do so when at least one nuisance parameter is estimated consistently. Furthermore, as we formally establish in Sections \[sec:theory1\] and \[sec:theory2\] and illustrate in our simulations, when a library of flexible machine learning estimators is used to estimate nuisance parameters, the approach proposed in the next section behaves nearly as well as an oracle that selects the estimator with smallest average squared bias, which vanishes at least as fast as any given choice of machine learners. This is quite remarkable as the proposed approach avoids directly estimating the squared bias.

Model selection via a minimax cross-validation\[sec:select\]
============================================================

Minimax criteria for model selection \[sec:selection\]
------------------------------------------------------

In this section, we consider alternative selection criteria which avoid estimating and directly minimizing Equation . Suppose that we have candidate models $\pi_k =\left( A|X\right), k \in \mathcal K_1 \equiv \{1,\cdots,K_1\}$ for the propensity score and $E_\tk\left( Y|A,X\right), \tk \in \mathcal K_2 \equiv \{1,\cdots,K_2\}$ for the outcome model, respectively. We begin by describing the population version of our minimax criteria, i.e., we focus on $\pi_{k}\left( A|X\right)$ and $E_{\tk}(Y|a,X)$, the asymptotic limits of $\widehat \pi_{k}\left( A|X\right)$ and $\widehat E_{\tk}(Y|a,X)$. We will introduce the cross-validated estimator in Section \[sec:cf\]. For each pair of candidate models $(k_1, \tk_1)$, we have $$\begin{aligned}
\psi _{k_1,\tk_1} = E\left(\frac{\left( -1\right) ^{1-A}}{\pi_{k_1}\left( A|X\right) }Y
-\left\{ \frac{\left( -1\right) ^{1-A}}{\pi_{k_1}\left( A|X\right) }E_{\tk_1}(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}E_{\tk_1}(Y|a,X)\right\}\right).
\label{eq:if2}\end{aligned}$$ The working models for the propensity score and outcome model could be parametric, semiparametric or nonparametric. A simple parametric case may entail positing that $\pi_k \left( A|X\right)$ and $E_\tk\left( Y|A,X\right)$ are chosen to be the following regression models $$\begin{aligned}
\text{logit}\Pr \left( A=1|X\right) &=&\alpha _{k,0}+\alpha _{k,1}^{T}h_{k}(X),\label{model1}\\
E\left( Y|A,X\right) &=&\beta _{\tk,0}+ \beta _{\tk,1}^{T}h_{\tk}(X)+\beta _{\tk,2}^{T}g_{\tk}(X)A+\beta _{\tk,3}A, \label{model2}\end{aligned}$$ for dictionary $\{h_k, h_\tk, g_\tk: k \in \mathcal K_1;\tk \in \mathcal K_2\}$. Subsequently, [$$\begin{aligned}
\psi_{k,\tk}=E \bigg( \frac{\left( -1\right) ^{1-A}}{\pi \left( A|X;\alpha _{k}\right) }\Big\{
Y -\beta _{\tk,0}-\beta _{\tk,1}^{T} h_{\tk}(X)
- \beta _{\tk,2}^{T} g_{\tk}(X)  A -\beta _{\tk,3} A
\Big\}
+ \beta _{\tk,2}^{T} g_{\tk}(X)+ \beta _{\tk,3} \bigg).\end{aligned}$$]{}

Recall that the doubly robust estimator  which depends on both unknown functions has zero bias if either one contains the truth. Motivated by this observation, we define the following perturbation of a fixed index pair $(k_1,\tk_1)$, $$\begin{aligned}
\per(k,\tk; k_1,\tk_1) \equiv (\psi_{k,\tk}- \psi_{k_1,\tk_1})^2. \label{definition:per}\end{aligned}$$

The perturbations defined above have the following forms.

\[lemma:1\] $$\begin{aligned}
\per( k_1,\tk; k_1,\tk_1) = E\left[\sum_{a}\left( -1\right) ^{1-a}  (\frac{\pi\left(
a|X\right) }{\pi_{k_1}\left( a|X\right) }-1) (E_{\tk_1}(Y|a,X)-E_{\tk}(Y|a,X))\right]^2,\end{aligned}$$

$$\begin{aligned}
\per( k,\tk_1; k_1,\tk_1) = E\left[\sum_{a}\left( -1\right) ^{1-a}  (\frac{\pi\left(
a|X\right) }{\pi_k\left( a|X\right) }- \frac{\pi\left(
a|X\right) }{\pi_{k_1}\left( a|X\right) } ) (E(Y|a,X)-E_{\tk_1}(Y|a,X))\right]^2.\end{aligned}$$

Subsequently, for each fixed pair $(k_1,\tk_1)$, we only consider perturbations over pairs $(k,\tk)$ with either $k=k_1$ or $\tk=\tk_1$, and evaluate the perturbation of $ \psi_{k,\tk}$ at $\psi_{k_1,\tk_1}$ as $$\text{per}(k,\tk; k_1,\tk_1)=
    \begin{cases}
      \text{per}(k_1,\tk; k_1,\tk_1) & \text{if}\ k=k_1, \\
      \text{per}(k,\tk_1; k_1,\tk_1) & \text{if}\ \tk=\tk_1,\\
      0 & \text{otherwise}.
    \end{cases}
\label{bias}$$ We may define the pseudo-risk $$\begin{aligned}
B^{(1)}_{k_1,\tk_1}=\max_{\substack{k=k_1,\tk \in \mathcal K_2;\\ \tk=\tk_1, k\in \mathcal K_1}} \text{per}(k,\tk; k_1,\tk_1),\end{aligned}$$ which measures the maximum change of underlying functional at a candidate selected model $(k_1,\tk_1)$ induced by perturbing one of the nuisance parameters at the time, and holding the other fixed. We call this a pseudo-risk because unlike a standard definition of risk (e.g., mean squared error) which is typically defined in terms of the data generating mechanism and a given candidate model/estimator, the proposed definition is in terms of all candidate models/estimators. We also consider the following pseudo-risk, $$\begin{aligned}
B^{(2)}_{k_1,\tk_1}=\max_{\tk_0 \in \mathcal K_2} \max_{\tk \in \mathcal K_2} \text{per}(k_1,\tk; k_1,\tk_0) + \max_{k_0\in \mathcal K_1} \max_{k\in \mathcal K_1} \text{per}(k,\tk_1; k_0,\tk_1).\end{aligned}$$ Evaluating the above perturbation for each pair $(k_1,\tk_1)$ gives $K_1 \times K_2$ pseudo-risk values $B_{k_1,\tk_1}, k_1\in \mathcal K_1;\tk_2\in \mathcal K_2$. Finally, we define $$\begin{aligned}
\arg\min_{(k_1,\tk_1)} B^{(1)}_{k_1,\tk_1},~~ \arg\min_{(k_1,\tk_1)} B^{(2)}_{k_1,\tk_1},\end{aligned}$$ as population version of selected models, respectively. we refer to $B^{(1)}_{k_1,\tk_1}$ as population minimax pseudo-risk and $B^{(2)}_{k_1,\tk_1}$ as population mixed-minimax pseudo-risk.

One may also consider the following alternative criterion: $$\begin{aligned}
\begin{cases}
 \dot{k} &=\arg\min_{k_1}  \max_{\tilde k} \text{per}( k_1,\tk; k_1,\dot{\tilde k}),\\
 \dot{\tilde k} &= \arg\min_{\tilde k_1}  \max_{k} \text{per}( k,\tk_1; \dot{k},\tk_1).\\
 \end{cases}\end{aligned}$$ However, there may not exist such pair $(\dot{k},\dot{\tilde k})$. The proposed minimax criteria solve the optimization jointly and avoid this difficulty.

There may be different ways to define pseudo-risk using different norms, e.g., the first kind $B^{(1)}_{k_1,\tk_1}$ can also be defined as $$\max_{\tk \in \mathcal K_2} \text{per}(k_1,\tk; k_1,\tk_1) + \max_{ k\in \mathcal K_1} \text{per}(k,\tk_1; k_1,\tk_1);$$ the second kind $B^{(2)}_{k_1,\tk_1}$ can also be defined as $$\sum_{\tk_0 \in \mathcal K_2} \max_{\tk \in \mathcal K_2} \text{per}(k_1,\tk; k_1,\tk_0) + \sum_{k_0\in \mathcal K_1} \max_{k\in \mathcal K_1} \text{per}(k,\tk_1; k_0,\tk_1).$$ The proposed two criteria are representatives of two novel ideas: The first type is given by the overall maximum squared bias (i.e., change in the estimated functional) at a given candidate estimator, induced by perturbing one nuisance parameter at the time over candidate models holding the other one fixed; The second type is given by the sum of two maximum squared bias terms, each capturing the bias induced by perturbing a single nuisance parameter only.

The second mixed-minimax criterion has a doubly robust property, i.e., $\psi_{\arg\min_{(k_1,\tk_1)} B^{(2)}_{k_1,\tk_1}}$ has zero bias if either nuisance model is consistently estimated by at least one candidate learner.

Multi-fold cross-validated estimator\[sec:cf\]
----------------------------------------------

Following [@vaart2006oracle], we avoid overfitting in implementing an empirical minimax selector by considering a multi-fold cross-validation scheme which repeatedly randomly splits the data $\mathcal O$ into two subsamples: a training set $\mathcal O^{0s}$ and a test set $\mathcal O^{1s}$, where $s$ refers to $s$-th split. The splits may be either deterministic or random without loss of generality. In the following, we consider random splits, whereby we let $T^s = (T_1^s,\cdots,T_n^s)\in \{0,1\}^n$ denote a random vector independent of $O_1,\cdots, O_n$. If $T_i^s=0$, $O_i$ belongs to the $s$-th training sample $\mathcal O^{0s}$; otherwise it belongs to the $s$-th test sample $\mathcal O^{1s}$, $s=1, \ldots,S$. For each $s$ and $(k,\tk)$, our construction uses the training samples to construct estimators $\widehat \pi_k^s$ and $\widehat E_\tk^s(Y|A,X)$. The validation sample is then used to estimate the perturbation defined in Equation , $$\begin{aligned}
\widehat {\text{per}}(k,\tk; k_1,\tk_1) =& \frac{1}{S}\sum_{s=1}^S \left[  (\widehat \psi_{k,\tk}^s  - \widehat \psi_{k_1,\tk_1}^s )
\right]^{2},
$$ where $$\widehat \psi_{k,\tk}^s = \PP^1_s
\left\{ \frac{\left( -1\right) ^{1-A}}{\widehat \pi^s_k\left( A|X\right) }\left\{
\begin{array}{c}
Y -  \widehat E^s_{\tk}(Y| X,A)
\end{array}\right\} +\sum_a (-1)^{1-a}\widehat E^s_\tk(Y|a,X) \right\},$$ $$\widehat \psi_{k_1,\tk_1}^s = \PP^1_s
\left\{ \frac{\left( -1\right) ^{1-A}}{\widehat \pi^s_{k_1}\left( A|X\right) }\left\{
\begin{array}{c}
Y -  \widehat E^s_{\tk_1}(Y| X,A)
\end{array}\right\} +\sum_a (-1)^{1-a}\widehat E^s_{\tk_1}(Y|a,X) \right\},$$ $$\PP^j_s= \frac{1}{\#\{1\leq i\leq n:T_i^s=j\}}\sum_{i:T_i^s=j} \delta_{X_i}, \quad j=0,1,$$ and $\delta_X$ is the Dirac measure. We then select the empirical minimizers of $$\widehat B^{(1)}_{k_1,\tk_1} = \max_{\substack{k=k_1,\tk \in \mathcal K_2;\\ \tk=\tk_1, k\in \mathcal K_1}} \widehat {\text{per}}(k,\tk; k_1,\tk_1)$$ and $$\widehat B^{(2)}_{k_1,\tk_1} = \max_{\tk_0 \in \mathcal K_2} \max_{\tk \in \mathcal K_2} \widehat {\text{per}}(k_1,\tk; k_1,\tk_0) + \max_{k_0\in \mathcal K_1} \max_{k\in \mathcal K_1} \widehat {\text{per}}(k,\tk_1; k_0,\tk_1)$$ among all candidate pairs $(k_1,\tk_1)$ as our final models, i.e., $$\begin{aligned}
 (k^\dagger,\tk^\dagger) = \arg \min_{(k_1,\tk_1)}\widehat B^{(1)}_{k_1,\tk_1}, \label{eq:dagger} \\
  (k^\diamond,\tk^\diamond) = \arg \min_{(k_1,\tk_1)}\widehat B^{(2)}_{k_1,\tk_1}.
 \end{aligned}$$ The final estimators are $$\widehat \psi_{k^\dagger,\tk^\dagger}=\frac{1}{S}\sum_{s=1}^S \widehat \psi^s_{k^\dagger,\tk^\dagger},~~\widehat \psi_{k^\diamond,\tk^\diamond}=\frac{1}{S}\sum_{s=1}^S \widehat \psi^s_{k^\diamond,\tk^\diamond}.$$ We provide a high-level Algorithm 1 for the proposed selection procedure. We also define the cross-validated oracle selectors $$\begin{aligned}
 (k^\star,\tk^\star) = \arg \min_{(k_1,\tk_1)} B^{(1)}_{k_1,\tk_1}, \label{eq:defstar1}\\
  (k^\circ,\tk^\circ) = \arg \min_{(k_1,\tk_1)} B^{(2)}_{k_1,\tk_1}, \label{eq:defstar2}
 \end{aligned}$$ where in a slight abuse of notation, we define the cross-validated pseudo-risk $$B^{(1)}_{k_1,\tk_1} = \max_{\substack{k=k_1,\tk \in \mathcal K_2;\\ \tk=\tk_1, k\in \mathcal K_1}} \frac{1}{S}\sum_{s=1}^S \left[  (\psi_{k,\tk}^s  - \psi_{k_1,\tk_1}^s )
\right]^{2},$$ $$B^{(2)}_{k_1,\tk_1} = \max_{\tk_0 \in \mathcal K_2} \max_{\tk \in \mathcal K_2} \frac{1}{S}\sum_{s=1}^S \left[  (\psi_{k_1,\tk}^s  - \psi_{k_1,\tk_0}^s )
\right]^{2} + \max_{k_0\in \mathcal K_1} \max_{k\in \mathcal K_1} \frac{1}{S}\sum_{s=1}^S \left[  (\psi_{k,\tk_1}^s  - \psi_{k_0,\tk_1}^s )
\right]^{2},$$ $$\psi_{k,\tk}^s = \PP^1
\left\{ \frac{\left( -1\right) ^{1-A}}{\widehat \pi^s_k\left( A|X\right) }\left\{
\begin{array}{c}
Y -  \widehat E^s_{\tk}(Y| X,A)
\end{array}\right\} +\sum_a (-1)^{1-a}\widehat E^s_\tk(Y|a,X) \right\},$$ $$\psi_{k_1,\tk_1}^s = \PP^1
\left\{ \frac{\left( -1\right) ^{1-A}}{\widehat \pi^s_{k_1}\left( A|X\right) }\left\{
\begin{array}{c}
Y -  \widehat E^s_{\tk_1}(Y| X,A)
\end{array}\right\} +\sum_a (-1)^{1-a}\widehat E^s_{\tk_1}(Y|a,X) \right\},$$ and $\PP^1$ denotes the true measure of $\PP_s^1$.

For each pair $(k_1,\tk_1)$, average the perturbations over the splits and obtain $$\widehat {\text{per}}(k,\tk; k_1,\tk_1) = \frac{1}{S}\sum_{s=1}^S \left[  (\widehat \psi_{k,\tk}^s  - \widehat \psi_{k_1,\tk_1}^s )
\right]^{2},$$ where $k=k_1,\tk \in \mathcal K_2$ or $\tk=\tk_1, k\in \mathcal K_1$ Calculate $$\widehat B^{(1)}_{k_1,\tk_1}=\max_{\substack{k=k_1,\tk \in \mathcal K_2;\\ \tk=\tk_1, k\in \mathcal K_1}} \wper(k,\tk,k_1,\tk_1),$$ $$\widehat B^{(2)}_{k_1,\tk_1} = \max_{\tk_0 \in \mathcal K_2} \max_{\tk \in \mathcal K_2} \widehat {\text{per}}(k_1,\tk; k_1,\tk_0) + \max_{k_0\in \mathcal K_1} \max_{k\in \mathcal K_1} \widehat {\text{per}}(k,\tk_1; k_0,\tk_1)$$ for each pair $(k_1,\tk_1)$ Pick $(k^\dagger,\tk^\dagger)=\arg\min_{(k,\tk)} \widehat B^{(1)}_{k,\tk}$, $(k^\diamond,\tk^\diamond)=\arg\min_{(k,\tk)} \widehat B^{(2)}_{k,\tk}$ as our selected models, and obtain the estimations of the parameter over the splits $$\widehat \psi_{k^\dagger,\tk^\dagger}=\frac{1}{S}\sum_{s=1}^S \widehat \psi^s_{k^\dagger,\tk^\dagger},~~\widehat \psi_{k^\diamond,\tk^\diamond}=\frac{1}{S}\sum_{s=1}^S \widehat \psi^s_{k^\diamond,\tk^\diamond};$$\
**Return** $(k^\dagger,\tk^\dagger)$, $(k^\diamond,\tk^\diamond)$ and $\widehat \psi_{k^\dagger,\tk^\dagger}$, $\widehat \psi_{k^\diamond,\tk^\diamond}$

Theoretical results\[sec:theory\]
=================================

Optimality of oracle selectors\[sec:theory1\]
---------------------------------------------

In this section, we establish certain optimality properties of the minimax and mixed-minimax oracle pseudo-risk selectors defined by Equations  and respectively. As we will later show by establishing excess risk bounds relating empirical selectors to their oracle counterparts, these optimality results imply near optimal behavior of the corresponding empirical (cross-validated) selector. In this vein, focusing on a functional of the doubly robust class with property , consider the collections of learners for $p$ and $b$ obtained from an independent sample of size $n$: $$\begin{aligned}
\mathcal{C}_{p} = \left\{ \widehat{p}_{1},\ldots,\widehat{p}_{K_{1}}\right\};~~
\mathcal{C}_{b} = \left\{ \widehat{b}_{1},\ldots,\widehat{b}_{K_{2}}\right\}.\end{aligned}$$For the purpose of inference, in the following, our analysis is conditional on $\mathcal{C}_{p}$ and $\mathcal{C}_{b}$. Suppose further that these learners satisfy the following assumptions.

\[asm:1\] Given any $\epsilon>0$, there exist constants $C_p, C_b>1$ and sufficiently large $n_0$ such that for $n>n_0$, $$\begin{aligned}
\frac{1}{C_{p}}\nu _{j}\leq \left\vert \widehat{p}_{j}(x)-p\left( x\right)
\right\vert \leq C_{p}\nu _{j},~~ j\in \mathcal K_1,\\
\frac{1}{C_{b}}\omega _{j}\leq \left\vert \widehat{b}_{j}(x)-b\left(
x\right) \right\vert \leq C_{b}\omega _{j},~~ j\in \mathcal K_2,\end{aligned}$$for any $x$ with probability $1-\epsilon$, where $\nu _{j}$ and $\omega_{j}$ depend on $n$.

In the following we write $a_n \lesssim b_n$ when there exists a constant $C>0$ such that $a_n \leq Cb_n$ for sufficiently large $n$. Without loss of generality, suppose that$$\nu _{\min }=\min_{j}\left\{ \nu _{j}:j \in \mathcal K_1\right\} \lesssim \omega
_{\min }=\min_{j}\left\{ \omega _{j}:j \in \mathcal K_2\right\}.$$ Let $$\begin{aligned}
\nu _{\max } =\max_{j}\left\{ \nu _{j}:j \in \mathcal K_1 \right\},~~ \omega _{\max } =\max_{j}\left\{ \omega _{j}:j \in \mathcal K_2\right\}.\end{aligned}$$

\[asm:2\] We assume that $\lim_{n\rightarrow \infty }v_{\max } < \infty ; ~~
\lim_{n\rightarrow \infty }\omega _{\max } < \infty .$

\[asm:3\] Suppose $$(p(X)- \widehat p_i(X))(b(X)- \widehat b_j(X))E[h_1(O)|X],$$ is continuous with respect to $X$ for $i\in \mathcal K_1; j\in \mathcal K_2$. Furthermore, suppose that the support of $X$ is closed and bounded.

Assumption 5.1 essentially states that the bias of $\widehat p_j$ and $\widehat b_j$ is eventually exactly of order $v_j$ and $w_j$ with large probability. Note that $\widehat p_j$ and $\widehat b_j$ may not necessarily be consistent i.e., $v_j$ and $w_j$ may converge to a positive constant. Assumption \[asm:2\] guarantees the bias of each learner does not diverge. Note also that Assumption \[asm:3\] need only hold for $i$ and $j$ such that $v_i=v_{\text{min}}$ and $w_j=w_{\text{max}}$ for Lemma \[lemma:rate\] given below to hold for the minimax bias selector. Let $\psi_{k^\star,\tk^\star} =\PP^1 \{H(\widehat{p}_{k^\star},\widehat{b}_{\tk^\star})\}$ and $\psi_{k^\circ,\tk^\circ} =\PP^1 \{H(\widehat{p}_{k^\circ},\widehat{b}_{\tk^\circ})\}$, where $H(\cdot,\cdot)$ is defined in Equation , $(k^\star,\tk^\star)$ and $(k^\circ,\tk^\circ)$ are defined in Equations  and . We have the following lemma.

\[lemma:rate\] Under Assumptions \[asm:1\]-\[asm:3\], we have that the bias of the minimax oracle selector is of the order of:$$\left\vert \psi _{k^\star,\tk^\star}-\psi_0 \right\vert =O_{P}\left( \frac{\nu _{\max }}{\omega _{\max }}\omega _{\min }^{2}\right),$$while the bias of the mixed-minimax oracle selector is of the order of: $$\left\vert \psi _{k^\circ,\tk^\circ}-\psi_0
\right\vert =O_{P}\left( \nu _{\min }\omega _{\min }\right).$$

The above lemma implies that in the event $\nu _{\max }/\omega _{\max }~ {\rightarrow }$ as $n \rightarrow \infty$ in probability, as would be the case, say in a setting where at least one machine learner of both $p$ and $b$ fails to be consistent, the bias of the oracle minimax selector $\left\vert \psi _{k^\star,\tk^{\star}}-\psi_0
\right\vert $ is of order of the maximum (comparing learners of $b$ to those of $p)$ of the minimum (across learners for $b$ and $p,$ respectively) squared bias, that is the maximin squared bias of learners of $b$ and $p$ which under our assumptions is equal to $\omega _{\min }^{2}$. In this case, the minimax selector provides a guaranty for adaptation only up to the least favorable optimal learner across nuisance parameters, and therefore may fail to fully leverage the fact that a fast learner of one nuisance parameter may compensate for a slower learner of another. In contrast, the mixed-minimax selector can leverage such a gap to improve estimation rate, so that in the above scenario, its rate of estimation would be $\nu _{\min }\omega _{\min
}\leq \omega _{\min }^{2}$ with equality only if $\nu _{\min }=\omega _{\min
},$ that is if one can learn $b$ and $p$ equally well.

Excess risk bound of the proposed minimax selector\[sec:theory2\]
-----------------------------------------------------------------

In this section, we focus on our first estimator, however, analogous results hold for the second estimator. We first introduce some notation used to study the excess risk bound of $\widehat \psi_{k^\dagger,\tk^\dagger}$. Define $U^s_{(k,\tk)}(k_1,\tk_1)$ as $$\begin{aligned}
&\frac{\left( -1\right) ^{1-A}}{\widehat \pi_k^s \left(
A|X\right) }Y -\left\{ \frac{\left( -1\right) ^{1-A}}{\widehat \pi_k^s \left( A|X\right) }\widehat E_\tk^s(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}\widehat E_\tk^s(Y|a,X)\right\} \nonumber \\
-& \frac{\left( -1\right) ^{1-A}}{\widehat \pi_{k_1}^s \left(
A|X\right) }Y -\left\{ \frac{\left( -1\right) ^{1-A}}{\widehat \pi_{k_1}^s \left( A|X\right) }\widehat E_{\tk_1}^s(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}\widehat E_{\tk_1}^s(Y|a,X)\right\}.
\label{eq:U}\end{aligned}$$ Based on our minimax selection criterion, Equation  and are equivalently expressed as $$(k^\dagger,\tk^\dagger) = \arg\min_{(k_1,\tk_1)} \max_{\substack{k=k_1,\tk \in \mathcal K_2;\\ \tk=\tk_1, k\in \mathcal K_1}} \frac{1}{S}\sum_{s=1}^{S} [\PP^1_s\{U^s_{(k,\tk)}(k_1,\tk_1) \}]^2,$$ $$(k^\star,\tk^\star) = \arg\min_{(k_1,\tk_1)} \max_{\substack{k=k_1,\tk \in \mathcal K_2;\\ \tk=\tk_1, k\in \mathcal K_1}} \frac{1}{S} \sum_{s=1}^{S}  [\PP^1 \{U^s_{(k,\tk)}(k_1,\tk_1) \} ]^2.$$

Next, we derive a risk bound for empirically selected model $(k^\dagger,\tk^\dagger)$ which states that its risk is not much bigger than the risk provided by the oracle selected model $(k^\star,\tk^\star)$. For this purpose, it is convenient to make the following boundedness assumption.

\[asm:positivity\] (1) $\pi(a|X)\geq M_1$ and $\widehat \pi_{k}(a|X)\geq M_1$ almost surely for $a=0,1$, $k \in \{1,\ldots,K_1\}$, and some $0<M_1<1$. (2) $|Y| \leq M_2$ and $|\widehat E_{\tk}(Y|A,X)| \leq M_2$ almost surely for $\tk \in \{1,\ldots,K_2\}$, and some $M_2>0$.

Suppose Assumptions \[asm:positivity\] holds, then we have that $$\begin{aligned}
& \frac{1}{S}\sum_{s=1}^S \PP^0[\PP^1 \{U^s_{(k^\ddagger,\tk^\ddagger)}(k^\dagger,\tk^\dagger)\}]^2 \\
\leq & \frac{1+2\delta}{S}\sum_{s=1}^S \PP^0[\PP^1 \{U^s_{(k^\sstar,\tk^\sstar)}(k^\star,\tk^\star)\}]^2 + O\left(\frac{(1+\delta)\log(1+K_1^2 \times K_2^2)}{n^{1/p}} \left(\frac{1+\delta}{\delta}\right)^{(2-p)/p} \right),\end{aligned}$$ for any $\delta>0$, and $1\leq p\leq 2$, where $$(k^\ddagger,\tk^\ddagger)=\arg\max_{\substack{k=k^\dagger,\tk \in \mathcal K_2;\\ \tk=\tk^\dagger, k\in \mathcal K_1}} \frac{1}{S} \sum_{s=1}^S [\PP^1_s\{U^s_{(k,\tk)}(k^\dagger,\tk^\dagger)\}]^2,$$ $$(k^\sstar,\tk^\sstar)=\arg\max_{\substack{k=k^\star,\tk \in \mathcal K_2;\\ \tk=\tk^\star, k\in \mathcal K_1}} \frac{1}{S}\sum_{s=1}^S [\PP^1 \{U^s_{(k,\tk)}(k^\star,\tk^\star)\}]^2,$$ and $\PP^0$ denotes the true measure of $\PP_s^0$. \[thm:2\]

The proof of this result is based on a finite sample inequality for $$\begin{aligned}
 \frac{1}{S}\sum_{s=1}^S \PP^0[\PP^1 \{U^s_{(k^\ddagger,\tk^\ddagger)}(k^\dagger,\tk^\dagger)\}]^2 - \frac{1}{S}\sum_{s=1}^S \PP^0[\PP^1 \{U^s_{(k^\sstar,\tk^\sstar)}(k^\star,\tk^\star)\}]^2,\end{aligned}$$ the excess pseudo-risk of our selected model compared to the oracle’s, which requires an extension of Lemma 8.1 of [@vaart2006oracle] to handle second order U-statistics. We obtained such an extension by making use of a powerful exponential inequality for the tail probability of the maximum of a large number of second order U-statistics derived by [@10.1007/978-1-4612-1358-1_2]. Note that Theorem \[thm:2\] generalizes to the doubly robust functionals in the class of [@rotnitzky2019mix], with Equation  replaced by $IF^s_{k,\tk}(\widehat \psi^s_{k_1,\tk_1})-IF^s_{k_1,\tk_1}(\widehat \psi^s_{k_1,\tk_1})$ in the definition of $U^s_{(k,\tk)}(k_1,\tk_1)$, where $IF^s_{k,\tk}(\widehat \psi^s_{k_1,\tk_1})$ is an influence function of $\psi$ evaluated at nuisance parameters $(k,\tk)$ and $\widehat \psi^s_{k_1,\tk_1}$ solves $\PP_1^s IF^s_{k_1,\tk_1}(\widehat \psi^s_{k_1,\tk_1})=0$ (see Algorithm 2 in the Appendix for details).

The bound given in Theorem \[thm:2\] is valid for any $\delta>0$, such that the error incurred by empirical risk is of order $n^{-1}$ for any fixed $\delta$ if $p=1$, suggesting in this case that our cross-validated selector performs nearly as well as a oracle selector with access to the true pseudo-risk. Theorem \[thm:2\] is most of interest in nonparametric/machine learning setting where pseudo-risk can be of order substantially larger than $O(n^{-1})$ in which can the error made in selecting optimal learner is negligible relative to risk. By allowing $\delta_n \rightarrow 0$, the choice $p=2$ gives an error of order $n^{-1/2}$, which may be substantially larger. If we write the bound in the form of [$$\begin{aligned}
\frac{\sum_{s=1}^S \PP^0[\PP^1 \{U^s_{(k^\ddagger,\tk^\ddagger)}(k^\dagger,\tk^\dagger)\}]^2}{\sum_{s=1}^S \PP^0[\PP^1 \{U^s_{(k^\sstar,\tk^\sstar)}(k^\star,\tk^\star)\}]^2} \leq (1+2\delta) + O\left(\frac{n^{-1/p}S(1+\delta)^{2/p}\log(1+K_1^2 \times K_2^2)}{{\sum_{s=1}^S \PP^0[\PP^1 \{U^s_{(k^\sstar,\tk^\sstar)}(k^\star,\tk^\star)\}]^2 \delta^{(2-p)/p}}}  \right),\end{aligned}$$ ]{}the risk ratio converges to 1 with the remainder of order $n^{-1/2}$. Furthermore, the derived excess risk bound holds for as many as $c^{n^\gamma}$ models ($K_1\times K_2$) for any $\gamma<1$ and $c>0$.

Suppose that $\widehat \pi_k(X)\rightarrow \pi(X)$ and $\widehat E_\tk(Y|A,X) \rightarrow E(Y|A,X)$ in probability for some pair $(k,\tk)$, $[\PP^1 \{U^s_{(k^\sstar,\tk^\sstar)}(k^\star,\tk^\star)\}]^2$ converges to zero as $n \rightarrow \infty$ by Lemma \[lemma:1\], otherwise $\lim_{n\rightarrow \infty}[\PP^1 \{U^s_{(k^\sstar,\tk^\sstar)}(k^\star,\tk^\star)\}]^2>0$ and we select model/estimator that is nearest to satisfying the double robustness property.

An theorem analogous to Theorem \[thm:2\] holds for the mixed-minimax selector. Although details are omitted, the proof is essentially the same.

Simulation studies\[sec:simulations\]
=====================================

In Section \[sec:simu1\], we first consider two different settings in context of (possibly misspecified) parametric models as illustrative examples of the proposed approach. Next, in Section \[sec:simu3\], we evaluate the proposed approach in the context of primary interest, where modern machine learning methods are used to estimate nuisance parameters; and compare to double debiased machine learning (DDML) [@Chernozhukov2018] for each candidate machine learner of nuisance parameters. For each setting, we use three-fold cross-validation, i.e., $S=3$ with $T_i$’s from Bernoulli(0.5).

Illustrative examples with parametric candidate models {#sec:simu1}
------------------------------------------------------

Consider the following five-variate functional forms [@zhao2017selective] as potential candidate parametric models, $$\begin{aligned}
f_1(x) &= \Big((x_1-0.5)^2, (x_2-0.5)^2, \ldots, (x_{5}-0.5)^2\Big)^T,\\
f_2(x) &= \Big((x_1-0.5)^3, (x_2-0.5)^3, \ldots, (x_{5}-0.5)^3\Big)^T,\\
f_3(x) &= \Big(x_1,x_2, \ldots, x_{5}\Big)^T,\\
f_4(x) &= \Big(\frac{1}{1+\exp(-20(x_1-0.5))}, \ldots, \frac{1}{1+\exp(-20(x_5-0.5))}\Big)^T.\end{aligned}$$

In the first simulation setting, the true model is not included as candidate model and all working models are misspecified, whereas in the second setting, the true model is included as a candidate model. Each simulation was repeated 500 times. For each setting, covariates $X_i$’s were independently generated from a uniform distribution, and the noise of outcome was normal with mean 0 and standard deviation $0.1$. In the first scenario, the data were generated from $$\begin{aligned}
\text{logit}\Pr \left( A=1|X\right) &=&(1,-1,1,-1,1)f_1(X),\\
E\left( Y|A,X\right) &=&1 + \mathbbm{1}^Tf_1(X)+ \mathbbm{1}^Tf_1(X) A+ A,\end{aligned}$$ where $ \mathbbm{1}=(1,1,1,1,1)^T$. In the second scenario, the data were generated from $$\begin{aligned}
\text{logit}\Pr \left( A=1|X\right) & =& (1,-1,1,-1,1)f_2(X),\\
E\left( Y|A,X\right) &= & 1 + \mathbbm{1}^Tf_2(X)+ \mathbbm{1}^Tf_2(X) A+ A.\end{aligned}$$ For both scenarios, we used $\{f_2,f_3,f_4\}$ as candidate models of $g$ and $h$ specified in Equations  and .

The squared bias of $\widehat \psi$ for both scenarios is shown in Figures \[fig:1\]-\[fig:3\], respectively. (bias$^2 \times 10^{-4}$) is displayed. “Minimizer” refers to “squared bias minimizer” given by Equation ; “Oracle1” and “Oracle2” refer to the oracle minimax and mixed-minimax selectors evaluated by Lemma \[lemma:1\], respectively; “Proposed1” and “Proposed2” refer to the proposed minimax and mixed-minimax selectors, respectively; “Separate” refers to the more conventional practical approach which performs model selection separately for each nuisance parameter via AIC [@akaike1974new]; “Truth” refers to using the underlying true candidate models for estimation.

In the first scenario, because the true model is not a candidate model, there is a notable gap between the squared bias of estimators obtained from working models and those estimated directly from true models. In addition, because the “squared bias minimizer” minimizes the squared bias, it naturally has smaller Monte Carlo squared bias than the proposed criteria. Note that we do not expect “squared bias minimizer” and the proposed approach to perform similarly even asymptotically because they minimize different objective functions, and recall that the former may not be attainable in this specific setting. The proposed method has smaller bias than selecting models separately. However, this may not always be the case because the proposed method does not necessarily minimize the squared bias directly. Additional simulations in the Appendix illustrate this point. In the second scenario, the gap observed in the first scenario vanishes asymptotically, and both proposed methods perform nearly as well as both oracle and “squared bias minimizer” selectors in large samples as illustrated in Figure \[fig:3\]. Mixed-minimax selector appears to perform somewhat better than minimax selector in small to moderate samples.

![Squared bias of Scenario 1[]{data-label="fig:1"}](biasplot1.pdf "fig:"){width="4in" height="2.5in"}\

![Squared bias of Scenario 2[]{data-label="fig:3"}](biasplot2.pdf "fig:"){width="4in" height="2.5in"}\

Model selection with machine learners {#sec:simu3}
-------------------------------------

Finally, we report simulation results for the setting of primarily interest, where various machine learners are used to form candidate estimators of nuisance parameters. We considered the following machine learning methods. For the propensity score model: 1. Logistic regression with $l_1$ regularization [@lasso; @friedman2010regularization]; 2. Classification random forests [@Breiman2001; @randomForest]; 3. Gradient boosting trees [@friedman2001greedy; @gbm2019]. For the outcome model: 1. Lasso [@lasso; @friedman2010regularization]; 2. Regression random forests [@Breiman2001; @randomForest]; 3. Gradient boosting trees. Data were generated from $$\begin{aligned}
\text{logit}\Pr \left( A=1|X\right) &=&(1,-1,1,-1,1)^Tf_4(X),\\
E\left( Y|A,X\right) &=&1 + \mathbbm{1}^Tf_4(X)+ \mathbbm{1}^Tf_4(X) A+ A.\end{aligned}$$ The outcome error term was normal with mean 0 and standard deviation $0.25$.

Implementing candidate estimators required selecting corresponding tuning parameters: The penalty $\lambda_n$ for Lasso was chosen using 10-fold cross-validation over the pre-specified grid $[10^{10},\ldots,10^{-2}]$; For gradient boosting trees [@gbm2019], all parameters were tuned using a 4-fold cross-validation over the following grid: `ntrees`=\[100,300\], `depth`=\[1,2,3,4\], `shrinkag`=\[0.001,0.01,0.1\]; We used the default values of minimum node size (1 for classification, 5 for regression), and number of trees (500) for random forest [@randomForest], while the number of variables randomly sampled at each split, i.e., `mtry`, was tuned by `tuneRF` function [@randomForest].

We compared the proposed estimators with three DDML estimators using Lasso, random forests, and gradient boosting trees to estimate nuisance parameters respectively [@Chernozhukov2018]. Each DDML was estimated by cross-fitting [@Chernozhukov2018], i.e., 1) using training data (random half of sample) to estimate nuisance parameters and validation data to obtain $\widehat \psi_1$; 2) swaping the role of training and validation dataset to obtain $\widehat \psi_2$; 3) computing the estimator as the average $\widehat \psi_{\text{CF}} = (\widehat \psi_1 + \widehat \psi_2)/2$. The squared bias of $\widehat \psi$ of different methods are shown in Figure \[fig:ml\]. (bias$^2 \times 10^{-4}$) is displayed. “LASSO” refers to using logistic regression with $l_1$ regularization for propensity score model, and standard Lasso for outcome model; “RF” refers to using classification forests for propensity score model, and regression forests for outcome model; “GBT” refers to using gradient boosting classification tree for propensity score model, and gradient boosting regression tree for outcome model. Both proposed estimators have smallest bias across sample sizes, and there is a notable gap between the proposed estimators and those estimated by DDML. It is not surprising that Lasso has the largest bias because the working models are not correctly specified. This confirms our earlier claim that combined with flexible nonparametric/machine learning methods, our proposed approach can in finite samples yield smaller squared bias than any given machine learning estimator, without directly estimating the corresponding squared bias.

![Squared bias of the proposed estimator and different machine learners[]{data-label="fig:ml"}](bias1.pdf "fig:"){width="4in" height="2.5in"}\

A smooth-max approach to post-selection approximate inference {#sec:softmax}
=============================================================

In this section, we propose a novel smooth-max approach as smooth approximation to proposed minimax and mixed-minimax model selection criteria. Such smooth approximation provides a genuine opportunity to perform valid post-selection inference, appropriately accounting for uncertainty in both selecting and estimating nuisance parameters. It is well known that the following smooth-max function $$\begin{aligned}
\Gamma(\tau) = \frac{1}{\tau} \log \sum_{i=1}^{m} \exp(\tau z_i),\end{aligned}$$ converges to $ \max(z_1,\ldots,z_m)$ as $\tau \rightarrow \infty$, where $z_1,\ldots,z_m$ are positive real numbers. Similarly, we define $$\begin{aligned}
\Gamma_{k_1,\tk_1}(\tau) = \frac{1}{\tau} \log \sum_{\substack{k=k_1,\tk \in \mathcal K_2;\\ \tk=\tk_1, k\in \mathcal K_1}} \exp\{\frac{\tau}{S}\sum_{s=1}^{S}(\widehat \psi^s_{k,\tk} -\widehat \psi^s_{k_1,\tk_1})^2 \},\end{aligned}$$ for minimax selector and $$\begin{aligned}
\Gamma_{k_1,\tk_1}(\tau) =& \frac{1}{\tau} \log \sum_{\tk,\tk_0 \in \mathcal K_2} \exp\{\frac{\tau}{S}\sum_{s=1}^{S}(\widehat \psi^s_{k_1,\tk} -\widehat \psi^s_{k_1,\tk_0})^2 \} \\+& \frac{1}{\tau} \log \sum_{ k,k_0 \in \mathcal K_1} \exp\{\frac{\tau}{S}\sum_{s=1}^{S}(\widehat \psi^s_{k,\tk_1} -\widehat \psi^s_{k_0,\tk_1})^2 \},\end{aligned}$$ for mixed-minimax selector. Note that $$\begin{aligned}
 \Gamma_{k_1,\tk_1}(\tau) \rightarrow  \max_{\substack{k=k_1,\tk \in \mathcal K_2;\\ \tk=\tk_1, k\in \mathcal K_1}} \frac{1}{S}\sum_{s=1}^S (\widehat \psi^s_{k,\tk} -\widehat \psi^s_{k_1,\tk_1})^2,
 \label{eq:limit1}
 \end{aligned}$$ and $$\begin{aligned}
 \Gamma_{k_1,\tk_1}(\tau) \rightarrow  \max_{\tk,\tk_0 \in \mathcal K_2} \frac{1}{S}\sum_{s=1}^S (\widehat \psi^s_{k_1,\tk} -\widehat \psi^s_{k_1,\tk_0})^2 + \max_{ k,k_0 \in \mathcal K_1} \frac{1}{S}\sum_{s=1}^S (\widehat \psi^s_{k,\tk_1} -\widehat \psi^s_{k_0,\tk_1})^2,
 \label{eq:limit2}
 \end{aligned}$$ as $\tau \rightarrow \infty$.

Recall that our goal is to select the model $(k_1,\tk_1)$ minimizing the right hand side of Equations  or , which is then used to estimate $\psi_0$. A smooth-max approximation to this selection process is given by $$\begin{aligned}
\widehat \psi(\tau) = \sum_{(k_1,\tk_1)} p_{k_1,\tk_1}(\tau) \widehat \psi_{k_1,\tk_1},\end{aligned}$$ where $$\begin{aligned}
p_{k_1,\tk_1}(\tau)
=  \frac{\exp\{ \tau [\sum_{(k,\tk)} \Gamma_{k,\tk}(\tau)-\Gamma_{k_1,\tk_1}(\tau)] \}}{ \sum_{(k',\tk')} \exp\{ \tau [\sum_{(k,\tk)} \Gamma_{k,\tk}(\tau)-\Gamma_{k',\tk'}(\tau)] \} },\end{aligned}$$ and $$\begin{aligned}
\widehat \psi_{k_1,\tk_1}=\frac{1}{S}\sum_{s=1}^S \widehat \psi^s_{k_1,\tk_1}.\end{aligned}$$ Because $\widehat \psi(\tau)$ is a smooth transformation of $\widehat \psi_{k_1,\tk_1}, k_1=1,\ldots,K_1, \tk_1=1,\ldots,K_2$, justification of this smooth approximation is that $p_{k_1, \tk_1}(\tau) \rightarrow 1$ as $\tau \rightarrow \infty$ if $(k_1, \tk_1)=\arg\min_{k,\tk}\Gamma_{k,\tk}(\tau)$, otherwise $p_{k_1, \tk_1}(\tau) \rightarrow 0$. Given a fixed $\tau$, inference can be carried out with the nonparametric bootstrap or sandwich asymptotic variance estimators using the delta method. In following section, we use this smooth approximation for post-selection inference in a data application.

The above claim certainly holds when $K_1$ and $K_2$ are bounded, and may be established even when $K_1$ and $K_2$ grow with $n$ by high-dimensional central limit theorem [@chernozhukov2017] although we do not formally prove this here. Assuming $K_1$ and $K_2$ are bounded is not a real limitation, particularly in machine learning settings, as most practical applications similar to that in the next section will likely be limited to a small to moderate number of machine learners. Furthermore, unlike in previous section, our proposed approach for post-selection inference technically requires that each $\widehat \psi_{k,\tk}$ admits a first order influence function, which is a condition that as previously mentioned may still hold even when candidate estimators include nonparametric regression or flexible machine learning methods provided that these estimators are consistent at rates faster than $n^{-1/4}$ [@robins2017; @Chernozhukov2018].

To conclude this section, we briefly discuss selection of $\tau$. The choice of $\tau$ essentially determines how well the smooth-max function approximates the minimax estimator as captured by the following inequality, $$\begin{aligned}
\max\{z_1,\ldots,z_m\}\leq \frac{1}{\tau} \log \sum_{i=1}^m \exp(\tau z_i) \leq \frac{1}{\tau} \log m + \max\{z_1,\ldots,z_m\},\end{aligned}$$ which holds for any positive real numbers $z_1,\ldots,z_m$. Thus, the approximation error of $\Gamma_{k,\tk}(\tau)$ is controlled by $\epsilon = \log m/\tau$.

Data Analysis\[sec:real\]
=========================

In this section, similarly to [@tan2006; @2015biasreduce; @tan2019model; @tan2019regularized], we reanalyze data from the Study to Understand Prognoses and Preferences for Outcomes and Risks of Treatments (SUPPORT) to evaluate the effectiveness of right heart catheterization (RHC) in the intensive care unit of critically ill patients. At the time of the study by [@5c6af36c0fb64cfcbb482d75c2bc7ff1], RHC was thought to lead to better patient outcomes by many physicians. [@5c6af36c0fb64cfcbb482d75c2bc7ff1] found that RHC leads to lower survival compared to not performing RHC.

We consider the effect of RHC on 30-day survival. Data are available on 5735 individuals, 2184 treated and 3551 controls. In total, 3817 patients survived and 1918 died within 30 days. To estimate the additive treatment effect $\psi_0 = E\{Y_1-Y_0\}$, 72 covariates were used to adjust for potential confounding [@Hirano2001]. We posited $K_1=3$ candidate models/estimators for the propensity score model including all 72 covariates: 1. Logistic regression with $l_1$ regularization; 2. Classification random forests; 3. Gradient boosting trees. We posited $K_2=6$ candidate estimators for the outcome model $E_\tk\left( Y|A,X\right)$ with 72 covariates: 1. Linear regression; 2. Logistic regression; 3. Lasso; 4. Logistic regression with $l_1$ regularization; 5. Regression random forests; 6. Classification random forests. The proposed selection procedure was implemented with three-fold cross-validation. Tuning parameters were selected as in Section \[sec:simu3\].

The proposed minimax selection selected gradient boosting trees estimator of propensity score model and logistic regression with $l_1$ regularization for the outcome model. The proposed mixed-minimax selection selected classification random forest estimator of propensity score model and regression random forest for the outcome model. The estimated causal effect of RHC was $-0.0548$ and $-0.0476$ for the minimax and mixed-minimax criteria, respectively, while the point estimate obtained by smooth-max approach was $-0.0528$ and $-0.0483$ which, as expected, is close to the minimax point estimates, respectively. Results were somewhat smaller than other improved estimators considered by [@2015biasreduce], who did not perform model selection. Specifically, the TMLE with default super learner [@van2011targeted] gave $\widehat \psi_{\text{TMLE-SL}} = - 0.0586$; the bias reduced doubly robust estimation with linear and logit link gave $\widehat \psi_{\text{BR},\text{lin}}=-0.0612$ and $\widehat \psi_{\text{BR},\text{logit}}=-0.0610$, respectively; the calibrated likelihood estimator [@tan2010] gave $\widehat \psi_{\text{TAN}}=-0.0622$. Our estimates suggest that previous estimates may still be subject to a small amount of misspecification bias. To obtain valid confidence intervals, we applied the proposed smooth-max approach with error tolerance $\epsilon =0.002$. Smooth-max based 95% confidence intervals by nonparametric bootstrap with 200 replications were estimated as $(-0.1041,-0.0301)$ and $(-0.1083,-0.0277)$ for the minimax and mixed-minimax criteria, respectively, which are slightly wider than other improved estimators considered in [@2015biasreduce], e.g., the targeted maximum likelihood estimation with default super learner gave $(-0.0877, -0.0295)$; the bias reduced doubly robust estimator with linear and logit link gave $(-0.0889,-0.0335)$ and $(-0.0879,-0.0340)$, respectively; the calibrated likelihood estimator gave $(-0.0924,-0.0319)$. This is not surprising because we consider a richer class of models and formally account for such selection step, potentially resulting in smaller bias and more accurate confidence intervals.

Discussion\[sec:discussion\]
============================

We have proposed a general model selection approach to estimate a functional $\psi(\theta)$ in a general class of doubly robust functionals which admit an estimating equation that is unbiased if at least one of two nuisance parameters is correctly specified. The proposed method works by selecting the candidate model based on minimax or mixed-minimax criterion of pseudo-risk defined in terms of the doubly robust estimating equation. A straightforward, cross-validation scheme was proposed to estimate the pseudo-risk. While, throughout the paper, we have described and evaluated the proposed selection procedure primarily in estimating average treatment effect as a running example, in the appendix, all results are extended to the more general class of mixed-bias functionals of [@rotnitzky2019mix]. Extensive simulation studies and a real data example on the effectiveness of right heart catheterization in the intensive care unit of critically ill patients were also presented to illustrate the proposed approach.

As mentioned in Remark \[remark:rm2\], our selection criteria extend to multiply robust influence functions in the sense of [@tchetgentchetgen2012; @wang2018bounded; @Caleb2019; @Shi2019MultiplyRC; @sun2019multiple], where three or more nuisance parameters are needed to evaluate the influence function, however, the influence function remains unbiased if all but one of the nuisance parameters are evaluated at the truth. Briefly, in such setting, the minimax criterion entails the maximum squared change of the functional over all perturbations of one nuisance parameter holding the others fixed. The mixed-minimax selector likewise generalizes. We expect our theoretical results to extend to this setting, an application of which is in progress [@sun2019multiple].

The proposed methods may be improved or extended in multiple ways. The choice of the criterion could be more flexible and one may use a different norm rather than $L_\infty$ norm, e.g., $L_2$ or $L_1$ norm. Another potential extension of our method is in the direction of statistical inference. It would be both interesting and important to derive the exact asymptotic distribution of the proposed estimators, as originally described in Section \[sec:select\], instead of relying on a smooth approximation. Finally, in principle one could develop a stacked generalization of our proposed approach by forming linear combinations of various candidate estimators of nuisance parameters [@wolpert1992stacked; @breiman1996bagging]. An optimal estimator could then be obtained by minimizing the pseudo-risk for the functional of interest with respect to the weights. Clearly, the current minimax approach explores only a small set of possible values for such weights, i.e., all values that have unit mass at one candidate model and zero elsewhere, and therefore may be sup-optimal relative to the stacked generalization. Because candidate learners may yield estimates that are highly correlated, one may need a form of regularization to ensure good performance, such as restricting set of weights to a finite support (in the spirit of [@van2007super]), or alternatively penalizing the pseudo-risk. We are currently investigating theoretical properties of such stacked minimax functional learning which we plan to report elsewhere.

Acknowledgment {#acknowledgment .unnumbered}
==============

We thank James Robins, Andrea Rotnitzky, and Weijie Su for helpful discussions and suggestions. We thank Karel Vermeulen and Stijn Vansteelandt for providing the dataset. This research is supported in part by U.S. National Institutes of Health grants.

Proofs
======

In this section, we present proofs of the theoretical results.

[**Proof of Lemma \[lemma:1\].**]{}

We have that

$$\begin{aligned}
\left[ \psi_{k_1,\tk} -  \psi_{k_1,\tk_1}\right]^2 &= E \Bigg[\frac{\left( -1\right) ^{1-A}}{\pi_{k_1}\left(
A|X\right) }Y -\left\{ \frac{\left( -1\right) ^{1-A}}{\pi_{k_1} \left( A|X\right) }
E_\tk(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}E_\tk(Y|a,X)\right\}  \\
&-  \frac{\left( -1\right) ^{1-A}}{\pi_{k_1}\left(
A|X\right) }Y + \left\{ \frac{\left( -1\right) ^{1-A}}{\pi_{k_1}\left( A|X\right) }
E_{\tk_1}(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}E_{\tk_1}(Y|a,X)\right\} \Bigg]^2\\
&= E\Bigg[\frac{\left( -1\right) ^{1-A}}{\pi_{k_1}\left(
A|X\right) } \left\{ E_{\tk_1}(Y|A,X)-E_{\tk}(Y|A,X) \right\} \\
& + \sum_{a}\left( -1\right) ^{1-a}\left\{ E_\tk(Y|a,X)-E_{\tk_1}(Y|a,X)\right\} \Bigg]^2\\
& = E\left[\sum_{a}\left( -1\right) ^{1-a}  \left\{\frac{\pi\left(
a|X\right) }{\pi_{k_1}\left( a|X\right) }-1\right\} \left\{E_{\tk_1}(Y|a,X)-E_{\tk}(Y|a,X)\right\}\right]^2.\end{aligned}$$

$$\begin{aligned}
\left[\psi_{k,\tk_1} - \psi_{k_1,\tk_1}\right]^2 &= E \Bigg[\frac{\left( -1\right) ^{1-A}}{\pi_{k}\left(
A|X\right) }Y -\left\{ \frac{\left( -1\right) ^{1-A}}{\pi_{k} \left( A|X\right) }
E_{\tk_1}(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}E_{\tk_1}(Y|a,X)\right\}  \\
&-  \frac{\left( -1\right) ^{1-A}}{\pi_{k_1}\left(
A|X\right) }Y + \left\{ \frac{\left( -1\right) ^{1-A}}{\pi_{k_1}\left( A|X\right) }
E_{\tk_1}(Y|A,X)-\sum_{a}\left( -1\right) ^{1-a}E_{\tk_1}(Y|a,X)\right\} \Bigg]^2\\
&= E\Bigg[\left\{\frac{\left( -1\right) ^{1-A}}{\pi_{k}\left(
A|X\right)}-\frac{\left( -1\right)^{1-A}}{\pi_{k_1}\left(
A|X\right)}\right\} E(Y|A,X) \\
&- \left\{\frac{\left( -1\right) ^{1-A}}{\pi_{k}\left(
A|X\right)}-\frac{\left( -1\right)^{1-A}}{\pi_{k_1}\left(
A|X\right)}\right\} E_{\tk_1}(Y|A,X)\Bigg]^2\\
&= E\left[\sum_{a}\left( -1\right) ^{1-a}  \left\{\frac{\pi\left(
a|X\right) }{\pi_k\left( a|X\right) }- \frac{\pi\left(
a|X\right) }{\pi_{k_1}\left( a|X\right) } \right\} \left\{E(Y|a,X)-E_{\tk_1}(Y|a,X)\right\}\right]^2.\quad \Box\end{aligned}$$

[**Proof of Lemma \[lemma:rate\]**]{}

Without loss of generality, we focus on $S=1$ and the proof of cross-validated oracle selector follows similarly. Let $$k^{\max }=\arg \max_{k\in \mathcal K_1} \text{per}\left( k,\tk^\star,k^\star,\tk^\star\right) ^{1/2},$$ $$\tk^{\max}=\arg \max_{\tk\in \mathcal K_2} \text{per}\left( k^\star,\tk,k^\star,\tk^\star\right) ^{1/2}.$$ By Assumptions \[asm:1\]-\[asm:3\] and the mean value theorem, there exist constants $C_{0},C_{1},$ and a value $x_{0}$ in the support of $X$ such that$$\begin{aligned}
&&\max_{\substack{ k=k^\star,\tk\in \mathcal K_2 \\ \tk = \tk^\star,k\in \mathcal K_1}}\text{per}\left(
k,\tk,k^\star,\tk^\star\right) ^{1/2}\\
&=&\max_{\tk=\tk^\star,k\in \mathcal K_1}\text{per}\left(
k,\tk,k^\star,\tk^\star \right) ^{1/2} \\
&=&\left\vert \left\{ \int \left( \widehat{p}_{k^\star}(x)-\widehat{p}_{k^{\max }}\left( x\right) \right) \left( \widehat{b}_{\tk^\star}(x)-b(x)\right) E[h_1(O)|X=x] \right\} dF\left( x\right) \right\vert  \\
&=&\left\vert \left\{ \int \left( \widehat{p}_{k^\star}(x)-\widehat{p}_{k^{\max }}\left( x\right) \right) \left( \widehat{b}_{1}(x)-b(x)\right)E[h_1(O)|X=x] \right\} dF\left( x\right) \right\vert  \\
&=&\left\vert \left( \widehat{p}_{k^\star}(x_{0})-\widehat{p}_{k^{\max
}}\left( x_0 \right) \right) \left( \widehat{b}_{1}(x_{0})-b(x_{0})\right) E[h_1(O)|X=x_0]
\right\vert  \\
&=&C_{0}\nu _{\max }\omega _{\min } \\
&\geq &
\max_{k=k^\star,\tk\in \mathcal K_2}\text{per}\left(
k,\tk,k^\star,\tk^\star\right) ^{1/2}\\
 &=&\left\vert \left\{ \int
\left( \widehat{p}_{k^\star}(x)-p\left( x\right) \right) \left( \widehat{b}_{1}(x)-\widehat{b}_{\tk^{\max }}(x)\right) E[h_1(O)|X=x] \right\} dF\left( x\right)
\right\vert  \\
&=&C_{1}\nu _{k^\star}\omega _{\max }.\end{aligned}$$Therefore $$\nu _{k^\star}\leq \frac{C_{0}}{C_{1}}\nu _{\max }\frac{\omega _{\min }}{\omega _{\max }}.$$Implied by Equation  and the mean value theorem, there exists a positive constant $C_{2}$ and a value $x^{\ast }$ in the support of $X$ such that $$\begin{aligned}
\left\vert \psi _{k^\star,\tk^\star}-\psi_0 \right\vert  &=&\left\vert \int
\left( \widehat{p}_{k^\star}(x)-p\left( x\right) \right) \left( \widehat{b}_{1}(x)-b(x)\right) E[h_1(O)|X=x] dF\left( x\right) \right\vert  \\
&=&C_{2}\left\vert \left( \widehat{p}_{k^\star}(x^{\ast })-p\left( x^{\ast
}\right) \right) \left( \widehat{b}_{1}(x^{\ast })-b(x^{\ast })\right) E[h_1(O)|X=x^*]
\right\vert  \\
&\lesssim &\nu_{k^\star}\omega _{\min } \\
&\lesssim &\frac{\nu _{\max }}{\omega _{\max }}\omega _{\min }^{2}.\end{aligned}$$For $\psi _{k^\circ,\tk^\circ}$, it is straightforward to verify that, by Equation  and the mean value theorem, there exists a positive constant $\overline{C}_{2}$ and a value $\overline{x}$ in the support of $X$ such that $$\begin{aligned}
\left\vert \psi _{k^\circ,\tk^\circ}-\psi_0
\right\vert  &=&\overline{C}_{2}\left\vert \left( \widehat{p}_{k^\circ}(\overline{x})-p\left( \overline{x}\right) \right) \left( \widehat{b}_{\tk^\circ}(\overline{x})-b(\overline{x})\right) \right\vert  \\
&=&O_{P}\left( \nu _{\min }\omega _{\min }\right).\end{aligned}$$

[**Proof of Theorem \[thm:2\].**]{}

By the definition of our estimator, $$\max_{\substack{k=k^\dagger,\tk \in \mathcal K_2;\\ \tk=\tk^\dagger, k\in \mathcal K_1}}  \sum_{s=1}^{S} [\PP^1_s\{U^s_{(k,\tk)}(k^\dagger,\tk^\dagger)\}]^2 \leq  \max_{\substack{k=k^\star,\tk \in \mathcal K_2;\\ \tk=\tk^\star, k\in \mathcal K_1}} \sum_{s=1}^{S} [\PP^1_s(U^s_{(k,\tk)}(k^\star,\tk^\star))]^2.$$ So we have that, $$\sum_{s=1}^{S} [\PP^1_s\{U^s_{(k^\ddagger,\tk^\ddagger)}(k^\dagger,\tk^\dagger)\}]^2 \leq \sum_{s=1}^{S}  [\PP^1_s(U^s_{(k^\vartriangle,\tk^\vartriangle)}(k^\star,\tk^\star))]^2,$$ where $$(k^\ddagger,\tk^\ddagger)=\arg\max_{\substack{k=k^\dagger,\tk \in \mathcal K_2;\\ \tk=\tk^\dagger, k\in \mathcal K_1}} \sum_{s=1}^{S} [\PP^1_s(U^s_{(k,\tk)}(k^\dagger,\tk^\dagger))]^2,$$ and $$(k^\vartriangle,\tk^\vartriangle)=\arg\max_{\substack{k=k^\star,\tk \in \mathcal K_2;\\ \tk=\tk^\star, k\in \mathcal K_1}} \sum_{s=1}^{S} [\PP^1_s(U^s_{(k,\tk)}(k^\star,\tk^\star))]^2.$$

We denote $n^s_j = \#\{1\leq i\leq n:T_i^s=j\},$ $j=0,1$. For simplicity, in a slight abuse of notation in the following, we use $(k,\tk)$ instead of $(k^\ddagger,\tk^\ddagger), (k^\vartriangle,\tk^\vartriangle)$ for the sub-indices of $U$, i.e., $U^s_{(k,\tk)}(k^\dagger,\tk^\dagger)$ denotes $U^s_{(k^\ddagger,\tk^\ddagger)}(k^\dagger,\tk^\dagger)$; $U^s_{(k,\tk)}(k^\star,\tk^\star)$ denotes $U^s_{(k^\vartriangle,\tk^\vartriangle)}(k^\star,\tk^\star)$. By simple algebra, we have $$\begin{aligned}
& [\PP^1_s (U^s_{(k,\tk)}(k^\dagger,\tk^\dagger))]^2 \\ =& \frac{1}{{n_1^s}^2} \sum_{i,j} [U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j ]\\
= &\frac{1}{{n_1^s}^2} \sum_{i,j} \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big] \Big)\\
+ & \frac{1}{{n_1^s}^2} \sum_{i,j}  \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big]   U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j\\
+ & \frac{1}{{n_1^s}^2} \sum_{i,j}  \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big]\end{aligned}$$ $$\begin{aligned}
= &   \frac{1}{{n_1^s}^2} \sum_{i,j}  \Big(\PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big] \Big) \\
+ & \frac{2}{{n_1^s}^2} \sum_{i,j}  \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big]\\
+ & \frac{1}{{n_1^s}^2} \sum_{i,j} \Big ( U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big] \Big),\end{aligned}$$ where $U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i$ denotes the estimating equation evaluated at $i$-th observation. Thus, $[\PP^1_s (U^s_{(k,\tk)}(k^\dagger,\tk^\dagger))]^2$ further equals to $$\begin{aligned}
& [\PP^1_s (U^s_{(k,\tk)}(k^\dagger,\tk^\dagger))]^2\\
= &  [\PP^1  \big(  U^s_{k,\tk} ({k^\dagger,\tk^\dagger}) \big)]^2 \\
+ & \frac{2}{{n_1^s}^2} \sum_{i,j}  \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1 \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big]\\
+ & \frac{1}{{n_1^s}^2} \sum_{i,j} \Big ( U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big] \Big).\end{aligned}$$ A similar decomposition holds for $[\PP_s^1(U^s_{(k,\tk)}({k^\star,\tk^\star}))]^2$. By definition of our estimator, for any $\delta>0$, we have that $$\begin{aligned}
& \sum_{s=1}^S [\PP^1 (U^s_{(k,\tk)}({k^\dagger,\tk^\dagger}))]^2 \\
\leq &(1+2\delta)\sum_{s=1}^S[\PP^1 (U^s_{(k,\tk)}({k^\star,\tk^\star}))]^2\\
 + & \frac{1}{\sqrt{\nn}} \Big\{ (1+\delta)\sqrt{\nn} \sum_{s=1}^S\Big[  \big(\PP_s^1(U^s_{(k,\tk)}({k^\star,\tk^\star})) \big)^2 -\big(  \PP^1 (U^s_{(k,\tk)}({k^\star,\tk^\star})) \big)^2 \Big] \\ &-  \delta \sqrt{\nn}\sum_{s=1}^S \Big[ \PP^1 (U^s_{(k,\tk)}({k^\star,\tk^\star}))  \Big]^2   \Big\}\\
-& \frac{1}{\sqrt{\nn}} \Big\{ (1+\delta)\sqrt{\nn} \sum_{s=1}^S \Big[  \big(\PP_s^1(U^s_{(k,\tk)}({k^\dagger,\tk^\dagger})) \big)^2 -\big(  \PP^1(U^s_{(k,\tk)}({k^\dagger,\tk^\dagger})) \big)^2 \Big] \\ &+ \delta \sqrt{\nn} \sum_{s=1}^S \Big[ \PP^1 (U^s_{(k,\tk)}({k^\dagger,\tk^\dagger}))  \Big]^2   \Big\}.\\\end{aligned}$$ Combined with the decomposition of $[\PP_s^1(U^s_{(k,\tk)}({k^\dagger,\tk^\dagger}))]^2$ and $[\PP_s^1(U^s_{(k,\tk)}({k^\star,\tk^\star}))]^2$, we further have that $$\begin{aligned}
&\sum_{s=1}^S [\PP^1(U^s_{(k,\tk)}({k^\dagger,\tk^\dagger}))]^2 \\
\leq &(1+2\delta)\sum_{s=1}^S [\PP^1(U^s_{(k,\tk)}({k^\star,\tk^\star}))]^2\\
 + & \frac{1}{\sqrt{n_1^s}} \Big\{ (1+\delta)\sqrt{n^s_1} \sum_{s=1}^S \Big[ \frac{2}{n_1^s} \sum_{i}  \Big (  U^s_{k,\tk} ({k^\star,\tk^\star})_i - \PP^1 \big[  U^s_{k,\tk} ({k^\star,\tk^\star})_i \big] \Big) \PP^1  \big[  U^s_{k,\tk} (\widehat \psi_{k^\star,\tk^\star}) \big] \\
 + & \frac{1}{{n_1^s}^2} \sum_{i, j} \Big ( U^s_{k,\tk} ({k^\star,\tk^\star})_i - \PP^1  \big[  U^s_{k,\tk} ({k^\star,\tk^\star})_i \big] \Big) \Big (  U^s_{k,\tk} ({k^\star,\tk^\star})_j- \PP^1  \big[  U^s_{k,\tk} ({k^\star,\tk^\star})_j \big] \Big)\Big]  \\ -&  \delta \sqrt{n_1^s}\sum_{s=1}^S \Big[ \PP^1(U^s_{(k,\tk)}({k^\star,\tk^\star}))  \Big]^2  \Big\}\\
  - & \frac{1}{\sqrt{n_1^s}} \Big\{ (1+\delta)\sqrt{n^s_1}\sum_{s=1}^S \Big[  \frac{2}{n_1^s} \sum_{i}  \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger}) \big] \\
 + & \frac{1}{{n_1^s}^2} \sum_{i, j} \Big ( U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \Big (  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j - \PP^1  \big[  U^s_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big] \Big)\Big]   \\+&  \delta \sqrt{n_1^s} \sum_{s=1}^S \Big[ \PP^1(U^s_{(k,\tk)}({k^\dagger,\tk^\dagger}))  \Big]^2  \Big\}.
 \end{aligned}$$

Note that the only assumption on $s$ 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 \[lemma:0\], and note that given Assumption \[asm:positivity\], the estimating equation is bounded, we further have the following bound for the first order term, $$\begin{aligned}
& \PP^0 \max_{k,\tk,k_1,\tk_1} \Big\{ \frac{2(1+\delta)\sqrt \nn}{\nn} \sum_{i}  \Big (  U_{k,\tk} ({k_1,\tk_1})_i - \PP^1  \big[  U_{k,\tk} ({k_1,\tk_1})_i \big] \Big) \PP^1  \big(  U_{k,\tk} ({k_1,\tk_1}) \big) \\
- & \delta \sqrt \nn \Big[ \PP^1  \big(  U_{k,\tk} ({k_1,\tk_1}) \big)  \Big]^2 \Big\} \\
\leq &  \PP^0 \frac{16(1+\delta)}{\nn^{1/p-1/2}} \log(1+K^2_1K^2_2)\max_{k,\tk,k_1,\tk_1}\bigg [\frac{||  U_{k,\tk} ({k_1,\tk_1})  \PP^1  \big[  U_{k,\tk} ({k_1,\tk_1}) \big] ||_\infty}{\nn^{1-1/p}}\\
 +& \bigg(\frac{3 \PP^1 \Big(U_{k,\tk} ({k_1,\tk_1})\PP^1  \big[  U_{k,\tk} ({k_1,\tk_1}) \big] \Big)^2 2^{(1-p)}(1+\delta)^{(2-p)}}{\delta^{2-p} \Big[ \PP^1  \big(  U_{k,\tk} ({k_1,\tk_1}) \big) \Big]^{4-2p}} \bigg)^{1/p} \bigg ] = (I),\end{aligned}$$ and $$\begin{aligned}
& \PP^0 \max_{k,\tk,k_1,\tk_1} \Big\{-\Big[ \frac{2(1+\delta)\sqrt \nn}{\nn} \sum_{i}  \Big (  U_{k,\tk} ({k_1,\tk_1})_i - \PP^1  \big[  U_{k,\tk} ({k_1,\tk_1})_i \big] \Big) \PP^1  \big(  U_{k,\tk} ({k_1,\tk_1}) \big)\\
+& \delta \sqrt \nn \Big[ \PP^1  \big(  U_{k,\tk} ({k_1,\tk_1}) \big) \Big]^2 \Big] \Big\}
  \leq (I),\end{aligned}$$ where $\max$ is taken on the set $\{ k_1\in \{1,\ldots,K_1\},\tk_1 \in \{1,\ldots,K_2\}, k=k_1~\text{or}~ \tk = \tk_1\}$. Thus, $\PP^0 [\PP^1 (U_{(k,\tk)}(k^\dagger,\tk^\dagger))]^2$ is further bounded by $$\begin{aligned}
&\PP^0 [\PP^1 (U_{(k,\tk)}({k^\dagger,\tk^\dagger}))]^2  \\
\leq &(1+2\delta)\PP^0 [\PP^1 (U_{(k,\tk)}({k^\star,\tk^\star}))]^2 +  \frac{2}{\sqrt \nn}(I)\\  + &  \frac{1+\delta}{\nn^{2}} \PP^0 \sum_{i, j} \Big ( U_{k,\tk} ({k^\star,\tk^\star})_i - \PP^1   \big[  U_{k,\tk} ({k^\star,\tk^\star})_i \big] \Big) \Big (  U_{k,\tk} ({k^\star,\tk^\star})_j - \PP^1   \big[  U_{k,\tk} ({k^\star,\tk^\star})_j \big] \Big)  \\
  -& \frac{1+\delta}{\nn^{2}} \PP^0  \sum_{i, j} \Big ( U_{k,\tk} ({k^\dagger,\tk^\dagger})_i - \PP^1   \big[  U_{k,\tk} ({k^\dagger,\tk^\dagger})_i \big] \Big) \Big(  U_{k,\tk} ({k^\dagger,\tk^\dagger})_j - \PP^1   \big[  U_{k,\tk} ({k^\dagger,\tk^\dagger})_j \big] \Big).\end{aligned}$$

By Lemma \[lemma:2\] and \[lemma:3\], the U-statistics are bounded and we have the bound of the risk of our selector $({k^\dagger,\tk^\dagger})$, $$\begin{aligned}
&\PP^0[\PP^1(U_{(k^\ddagger,\tk^\ddagger)}({k^\dagger,\tk^\dagger}))]^2 \\
\leq &(1+2\delta)\PP^0[\PP^1(U_{(k^\vartriangle,\tk^\vartriangle)}({k^\star,\tk^\star}))]^2 \\
+ & (1+\delta)C \Bigg\{ \left( \frac{2M}{n_1^2} \log(1+\frac{MK^2_1 K^2_2}{2} ) \right)^{1/2} + \frac{2M}{n_1} \log(1+\frac{MK^2_1 K^2_2}{2} ) \\
+ & \frac{4M^{3/2}}{n_1^{3/2}} \log^{3/2}(1+\frac{MK^2_1K^2_2}{2}+D_0)+ \frac{4 M^2}{n_1^2} \log^{2}(1+\frac{MK^2_1K^2_2}{2}+D_1) \Bigg\}\\
+ & \PP^0 \frac{16(1+\delta)}{\nn^{1/p}} \log(1+K^2_1K^2_2)\max_{k,\tk,k_1,\tk_1}\bigg [\frac{||   U_{k,\tk} ({k_1,\tk_1})\PP^1   \big[  U_{k,\tk} ({k_1,\tk_1}) \big] ||_\infty}{\nn^{1-1/p}} \\
+& \bigg(\frac{3\PP^1  \Big(U^2_{k,\tk} ({k_1,\tk_1}) \Big) 2^{(1-p)}(1+\delta)^{(2-p)}}{\delta^{2-p} \Big[ \PP^1   \big(  U_{k,\tk} ({k_1,\tk_1}) \big) \Big]^{2-2p}} \bigg)^{1/p} \bigg ],\\
   \end{aligned}$$ where $C$, $M$, $D_0$, and $D_1$ are some universal constants.

Finally, recall that for the term $(1+2\delta)\PP^0[\PP^1 (U_{(k^\vartriangle,\tk^\vartriangle)}({k^\star,\tk^\star}))]^2$, $(k^\vartriangle,\tk^\vartriangle)$ is chosen corresponding to $(k^\star,\tk^\star)$ under measure $\PP_s^1$. It is further bounded by $(1+2\delta)\PP^0[\PP^1 (U_{(k^\sstar,\tk^\sstar)}({k^\star,\tk^\star}))]^2$, where $(k^\sstar,\tk^\sstar)$ is chosen corresponding to $(k^\star,\tk^\star)$ under true measure $\PP^1$. $\Box$

[(Lemma 2.2 in [@vaart2006oracle])]{}\[lemma:0\] Assume that $E f\geq 0$ for every $f \in \cal F$. Then for any $1 \leq p \leq 2$ and $\delta>0$, we have that $$\begin{aligned}
E\max_{f \in \mathcal F}( \GG_n - \delta \sqrt{n} E) f \leq \frac{8}{n^{1/p-1/2}} \log(1+\#{\cal F}) \max_{f\in \mathcal F} [ \frac{M(f)}{n^{1-1/p}}+(\frac{v(f)}{(\delta E f)^{2-p}})^{1/p} ],\end{aligned}$$ where $\GG_n$ is the empirical process of the $n$ i.i.d. observations, and $(M(f),v(f))$ is any pair of Bernstein numbers of measurable function $f$ such that $$M(f)^2 E\Big(\exp\{ |f|/M(f)\} -1-|f|/M(f)\Big)\leq 1/2v(f).$$ Furthermore, if $f$ is uniformly bounded, then $(||f||_\infty,1.5E f^2)$ is a pair of Bernstein number.

$$\begin{aligned}
& \Pr\Bigg\{ \frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2})\geq x \Bigg\}\\
\leq & \frac{M}{2}\exp\Bigg\{ -\frac{1}{M} \min\bigg( \frac{x^2 n_1^2}{Eh^2}, \frac{xn_1}{||h||_{L_2\rightarrow L_2}} ,\frac{x^{2/3} n_1}{[( ||E_{O_1}h^2||_\infty + ||E_{O_2}h^2||_\infty)]^{1/3}},\frac{x^{1/2}n_1}{||h||_\infty^{1/2}} \bigg) \Bigg\},\end{aligned}$$

and $$\begin{aligned}
& \Pr \Bigg\{ -\frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2})\geq x \Bigg\}\\
\leq & \frac{M}{2}\exp\Bigg\{ -\frac{1}{M} \min\bigg( \frac{x^2 n_1^2}{Eh^2}, \frac{xn_1}{||h||_{L_2\rightarrow L_2}} ,\frac{x^{2/3} n_1}{[( ||E_{O_1}h^2||_\infty + ||E_{O_2}h^2||_\infty)]^{1/3}},\frac{x^{1/2}n_1}{||h||_\infty^{1/2}} \bigg) \Bigg\}.\end{aligned}$$ where $M$ is some universal constant, and $||h||_{L_2\rightarrow L_2}$ is defined as $$||h||_{L_2\rightarrow L_2}=\sup\{ E [h(O_1,O_2)a(O_1)c(O_2)]: E(a^2(O_1))\leq 1, E(c^2(O_2))\leq 1 \}.$$ \[lemma:2\]

[**Proof of Lemma \[lemma:2\].**]{}

The inequality follows directly from the Corollary 3.4 in [@10.1007/978-1-4612-1358-1_2].

$$\begin{aligned}
& E\max_{h\in \cal H}\Bigg\{  \frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2}) \Bigg\}\\
\leq & \left( 2M \log(1+\frac{M \# \cal H}{2} ) \max_{h} \frac{Eh^2}{ n_1^2} \right)^{1/2} \\
+ & 2M \log(1+\frac{M\# \cal H}{2} ) \max_{h} \frac{||h||_{L_2\rightarrow L_2}}{n_1 } \\
+ & 4 M^{3/2} \log^{3/2}(1+\frac{M\# \cal H}{2}+D_0) \max_{h}  \frac{||E_{O_1}h^2||_\infty^{1/2}}{n_1^{3/2}} \\
+ & 4 M^2 \log^{2}(1+\frac{M\# \cal H}{2}+D_1) \max_{h}\frac{ ||h||_\infty}{n_1^2},\end{aligned}$$

and the same result holds for $E\max_{h\in \cal H} \Bigg\{-  \frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2}) \Bigg\}$. \[lemma:3\]

[**Proof of Lemma \[lemma:3\].**]{}

Denote $$\begin{aligned}
\omega_{n_1}(x)=\min\bigg( \frac{x^2 n_1^2}{Eh^2}, \frac{xn_1}{||h||_{L_2\rightarrow L_2}} ,\frac{x^{2/3} n_1}{[ 2||E_{O_1}h^2||_\infty]^{1/3}},\frac{x^{1/2}n_1}{||h||_\infty^{1/2}} \bigg),\end{aligned}$$ and the following four sets:

$\Omega_{1,n_1}(h)=\{x: \omega_{n_1}(x)= \frac{ x^2 n_1^2}{ Eh^2} \}$

$\Omega_{2,n_1}(h)=\{x: \omega_{n_1}(x)= \frac{xn_1}{||h||_{L_2\rightarrow L_2}} \}$

$\Omega_{3,n_1}(h)=\{x: \omega_{n_1}(x)= \frac{x^{2/3}n_1}{[ 2||E_{O_1}h^2||_\infty]^{1/3}} \}$

$\Omega_{4,n_1}(h)=\{x: \omega_{n_1}(x)= \frac{ x^{1/2}n_1}{||h||_\infty^{1/2}} \}$

Then $$\begin{aligned}
&\Pr\Bigg\{ \bigg( \frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2}) \bigg) I_{\Omega_{1,n_1}}\geq x \Bigg\}
\leq  \frac{1}{2}M\exp\Bigg\{- \frac{1}{M} \frac{ x^2 n_1^2}{ Eh^2} \Bigg\} ,\\
&\Pr\Bigg\{ \bigg( \frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2})\bigg) I_{\Omega_{2,n_1}}\geq x \Bigg\}\leq \frac{1}{2}M\exp\Bigg\{- \frac{1}{M} \frac{xn_1}{||h||_{L_2\rightarrow L_2}} \Bigg\},\\
& \Pr\Bigg\{ \bigg( \frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2}) \bigg) I_{\Omega_{3,n_1}}\geq x \Bigg\}\leq   \frac{1}{2}M\exp\Bigg\{- \frac{1}{M} \frac{x^{2/3}n_1}{[ 2||E_{O_1}h^2||_\infty]^{1/3}} \Bigg\},\\
&\Pr\Bigg\{ \bigg( \frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2})\bigg) I_{\Omega_{4,n_1}}\geq x \Bigg\} \leq   \frac{1}{2}M\exp\Bigg\{- \frac{1}{M} \frac{x^{1/2}n_1}{||h||_\infty^{1/2}}\Bigg\}.\end{aligned}$$

Then by the above inequalities and Lemma 8.1 in [@vaart2006oracle], $$\begin{aligned}
& E\max_{h}\Bigg\{  \frac{1}{n_1^2} \sum_{1\leq i_1, i_2 \leq n_1} h(O_{i_1},O_{i_2}) \Bigg\}\\
\leq & \left( 2M \log(1+\frac{M\# \cal H}{2} ) \max_{h} \frac{Eh^2}{ n_1^2} \right)^{1/2} \\
+ & 2M \log(1+\frac{M\# \cal H}{2} ) \max_{h} \frac{||h||_{L_2\rightarrow L_2}}{n_1 } \\
+ & 4 M^{3/2} \log^{3/2}(1+\frac{M\# \cal H}{2}+D_0) \max_{h}  \frac{||E_{O_1}h^2||_\infty^{1/2}}{n_1^{3/2}} \\
+ & 4 M^2 \log^{2}(1+\frac{M\# \cal H}{2}+D_1) \max_{h}\frac{ ||h||_\infty}{n_1^2},\end{aligned}$$ where $D_0$ and $D_1$ are some universal constants. $\Box$\

For each pair $(k_1,\tk_1)$, average the perturbations over the splits and obtain $$\widehat {\text{per}}(k,\tk; k_1,\tk_1) = \frac{1}{S}\sum_{s=1}^S \left[  \PP_s^1( IF^s_{k,\tk}(\widehat \psi^s_{k_1,\tk_1}) - IF^s_{k_1,\tk_1}(\widehat \psi^s_{k_1,\tk_1}) )
\right]^{2},$$ where $k=k_1,\tk \in \{1\cdots, K_2\}$ or $\tk=\tk_1, k\in \{1\cdots, K_1\}$ Calculate $$\widehat B^{(1)}_{k_1,\tk_1}=\max_{\substack{k=k_1,\tk \in \{1\cdots, K_2\};\\ \tk=\tk_1, k\in \{1\cdots, K_1\}}} \wper(k,\tk,k_1,\tk_1),$$ $$\widehat B^{(2)}_{k_1,\tk_1} =  \max_{\tk_0 \in \mathcal K_2} \max_{\tk \in \mathcal K_2} \widehat {\text{per}}(k_1,\tk; k_1,\tk_0) + \max_{k_0\in \mathcal K_1} \max_{k\in \mathcal K_1} \widehat {\text{per}}(k,\tk_1; k_0,\tk_1)$$ for each pair $(k_1,\tk_1)$ Pick $(k^\dagger,\tk^\dagger)=\arg\min_{(k,\tk)} \widehat B^{(1)}_{k,\tk}$, $(k^\diamond,\tk^\diamond)=\arg\min_{(k,\tk)} \widehat B^{(2)}_{k,\tk}$ as our selected models, and obtain the estimations of the parameter over the splits $$\widehat \psi_{k^\dagger,\tk^\dagger}=\frac{1}{S}\sum_{s=1}^S \widehat \psi^s_{k^\dagger,\tk^\dagger},~~\widehat \psi_{k^\diamond,\tk^\diamond}=\frac{1}{S}\sum_{s=1}^S \widehat \psi^s_{k^\diamond,\tk^\diamond};$$\
**Return** $(k^\dagger,\tk^\dagger)$, $(k^\diamond,\tk^\diamond)$ and $\widehat \psi_{k^\dagger,\tk^\dagger}$, $\widehat \psi_{k^\diamond,\tk^\diamond}$

[**Additional simulations.**]{} In this section, we present more simulations in which all models are misspecified. In both scenarios, the data were generated from $$\begin{aligned}
\text{logit}\Pr \left( A=1|X\right) &=&(1,-1,1,-1,1)f_1(X),\\
E\left( Y|A,X\right) &=&1 + \mathbbm{1}^Tf_1(X)+ \mathbbm{1}^Tf_1(X) A+ A.\end{aligned}$$ For the first scenario, we used $\{f_2,f_3,f_4, f_5\}$ as candidate models of $g$ and $h$ specified in Equations  and , where $$f_5(x) = \Big(\text{cos}(\pi x_1), \ldots, \text{cos}(\pi x_5)\Big)^T.$$ For the second scenario, we used $\{f_2,f_3,f_4, f_6\}$ as candidate models of $g$ and $h$ specified in Equations  and , where $$f_6(x) = \Big(\text{cos}(\pi x_1/2), \ldots, \text{cos}(\pi x_5/2)\Big)^T.$$

The squared bias of $\widehat \psi$ for each scenario is shown in Figures \[fig:11\]-\[fig:33\], respectively. We see that when all models are misspecified, the performance of each method depends on the class of candidate models. While in Scenario 1, Figure \[fig:11\] is similar to Figure \[fig:1\], and the proposed methods have smaller bias; In Scenario 2, conducting model selection separately performs similar to and sometimes better than the proposed methods.

![Squared bias of Scenario 1[]{data-label="fig:11"}](biasplot3.pdf "fig:"){width="4in" height="2.5in"}\

![Squared bias of Scenario 2[]{data-label="fig:33"}](biasplot4.pdf "fig:"){width="4in" height="2.5in"}\
