Robust Q-learning
Ashkan Ertefaie James R. McKay David Oslin Robert L. Strawderman
stat.ME stat.ML


#1 1




1
  
  
  Robust Q-learning
    Ashkan Ertefaie^1, James R. McKay^2, David Oslin^3,
    Robert L. Strawderman^1  

 

    ^1Department of Biostatistics and Computational Biology, University of Rochester

    ^2Center on the Continuum of Care in the Addictions, Department of Psychiatry, University of Pennsylvania

    ^3Philadelphia Veterans Administration Medical Center, and Treatment Research Center and Center for Studies of Addictions, Department of Psychiatry, University of Pennsylvania 

    December 30, 2023
====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

 


0
  

    Robust Q-learning


   



Q-learning is a regression-based approach that is widely used to formalize the development of an optimal dynamic treatment strategy. Finite dimensional working models are typically used to estimate certain nuisance parameters, and misspecification of these working models can result in residual confounding and/or efficiency loss. We propose a robust Q-learning approach which allows estimating such nuisance parameters using data-adaptive techniques. We study the asymptotic behavior of our estimators and provide simulation studies that highlight the need for and usefulness of
the proposed method in practice. We use the data from the “Extending Treatment Effectiveness of Naltrexone” multi-stage randomized 
trial to illustrate our proposed methods.  



Keywords: Cross-fitting, Data-adaptive techniques, Dynamic treatment strategies, Residual confounding  


1.5 




§ INTRODUCTION

A dynamic treatment strategy is a sequence of decision rules that maps individual characteristics to a treatment option at each decision point (i.e., a specific point in time in which a treatment is to be considered or altered). An optimal dynamic treatment strategy seeks to make these decisions to maximize a particular expected health outcome <cit.>. 
This is similar to clinical decision making whereby care providers tailor the type/dose of treatment over the course of clinical care based on ongoing information regarding patient progress in treatment. 

The main goal of precision medicine (i.e., developing an effective dynamic treatment strategy) is to use patient characteristics to inform a personalized treatment plan as a sequence of decision rules that leads to the best possible health outcome for each patient  <cit.>. Q-learning is a reinforcement learning algorithm that is widely used to estimate an optimal dynamic treatment strategy using data from multi-stage randomized clinical  trials or observational studies <cit.>.  Starting with the final study stage, Q-learning finds the treatment option that optimizes the desired expected outcome. Fixing the optimally-chosen treatment at the final stage,  Q-learning moves backward to the immediately preceding stage and searches for a treatment option assuming that future treatments will be optimized. The process continues until the first stage is reached. This  backward induction procedure is designed to avoid treatment options that appear to be optimal in the short term but may lead to a less desirable long-term outcome <cit.>.  Similar to other model-based approaches, model misspecification can seriously affect the result of Q-learning and the problem exacerbates substantially as the number of stages increases. Specifically, it can lead to residual confounding and suboptimal  dynamic treatment strategies <cit.>. An alternative approach is A-learning; this backward induction strategy is also model-based, 
hence subject to the possibility of misspecification, but imposes somewhat less restrictive regression model assumptions through modeling
only the contrasts between treatments and the propensity of treatment assignment given the observed patient history <cit.>. However, this extra robustness comes at the price substantially reduced efficiency; indeed, Q-learning may lead to a parameter estimate that  can be up to 170% more 
efficient than A-learning <cit.>. Policy learning methods are another class of methods for estimating an optimal dynamic treatment regime that circumvent the need for the conditional outcome models by directly optimizing the expected outcome among a class of rules <cit.>. Despite this appealing feature, policy learning methods are inefficient and  fail to provide reasonable inference for the parameter estimates that define non-smooth decision rules  (e.g., indicator or max operators) due
to slow rates of convergence <cit.>.


 Doubly robust estimators that are based on modeling both the treatment and outcome processes have also been proposed for policy learning and for structural nested models <cit.>. In general terms, consistency of the doubly robust estimators is guaranteed as long as either the treatment assignment mechanism or the postulated conditional mean outcome models are correctly specified, and semiparametric efficiency follows when both models are correctly specified <cit.>.  However, 
while doubly robust estimators give two routes for consistent estimation, the performance of these estimators depends critically on the modeling choice for the indicated treatment and mean outcome parameters. In practice, finite-dimensional models are used that are too restrictive and likely to be misspecified <cit.>. <cit.> showed that doubly robust estimators can have poor performance when both models are misspecified. To mitigate this problem, bias reduction techniques have been proposed <cit.>.  An alternative is to use flexible learning-based methods that may reduce the chance of inconsistency <cit.>.

In this paper, we consider the problem of Q-learning for the setting of a two-stage dynamic binary treatment choice regime. 
A typical approach involves postulating linear models for both the first and second stage Q-functions that have
stage-specific main effects involving pretreatment variables and interactions of these variables with the stage-specific treatment choice 
(e.g., ); however, only the interaction terms in these models directly influence the corresponding treatment decision functions 
that are optimized as part of the Q-learning process. Due to the nature of backward induction, the first stage model is likely to be a complicated 
function of the relevant covariates; misspecification of the main effect models in either stage can induce non-ignorable residual confounding.

To increase robustness, we therefore consider the indicated main effects as unknown nuisance parameters and we adapt an approach originally
proposed in <cit.> for partial linear models that allows us to eliminate these parameters from the Q-functions. In particular, these hard-to-estimate
parameters are replaced with models for the treatment assignment probability and mean outcome given pretreatment covariates that can be more easily estimated using, for example, nonparametric regression methods or related statistical learning methods (e.g., random forests). The resulting transformation process leads to consistent and asymptotically normal estimators of the so-called first- and second-stage blip functions that are robust to misspecification of the main effect nuisance parameters.



§ NOTATION AND FORMULATION


Consider a two-stage study where binary treatment decisions are made at each time point. 
Let 𝒪 = (X_1,A_1,X_2,A_2,Y) follow some probability distribution P_0 and 
suppose we observe N independent, identically distributed trajectories of 𝒪.
The vector X_1 ∈𝒳_1 ⊂^p_1 consists of all available baseline covariates measured 
before treatment at the first decision point A_1∈{0,1} and the vector X_2 ∈𝒳_2 ⊂^p_2 
consists of all available intermediate covariates measured before treatment at the second decision point A_2∈{0,1}. 
For notational convenience we define  
S^0_i=(X_1i^⊤,A_1i,X_2i^⊤)^⊤∈𝒮⊂^p_1+p_2+1 and 
 W^0_i=X_1i∈𝒳_1 ⊂^p_1.



For later use, we also define variables S_i and W_i, i=1,…,n; respectively,
each represents some finite dimensional function of the variables in S^0_i and  W^0_i.

Hence, knowledge of S^0_i and W^0_i respectively implies knowledge of S_i and W_i;
however, the reverse may not hold.

The observed outcome Y ∈ (measured after A_2) is assumed continuous, with a larger value of Y indicating  a better clinical outcome. 




§ A ROBUST FORMULATION OF Q-LEARNING


The outcome Y is assumed to satisfy the model

    Y_i = η_2(S^0_i) + A_2i·Δ_2(S^0_i) + ϵ_2i, 
       ( ϵ_2i | S^0_i, A_2i ) = 0,  i=1,…,N,

where the deterministic, unknown real-valued functions η_2(·) and Δ_2(·) are defined on 𝒮.
Because the treatment variable A_2 is binary,  the additive error model  (<ref>) 
places no parametric constraints on the conditional mean function.

In Q-learning,  backward induction is used to characterize the optimal dynamic treatment regime. 
Define the second stage Q-function 

    Q_2(s^0, a_2) = ( Y  | S^0 = s^0, A_2 = a_2 ) =  η_2(s^0) + a_2 Δ_2(s^0);

this measures “quality” when treatment a_2 is assigned 
to a patient with characteristics s^0 at the second stage <cit.>.
Similarly, define the first stage Q-function as

    Q_1(w^0, a_1) = {max_a_2 Q_2(S^0, a_2)  | W^0 = w^0, A_1 = a_1 }

where, for a_2 ∈{0,1}, 

    max_a_2 Q_2(s^0, a_2) = η_2(s^0) + Δ_2(s^0) I{Δ_2(s^0) > 0 }.

Analogously to Q_2(s^0, a_2), Q_1(w^0, a_1) measures “quality” when treatment a_1 is assigned to a patient with 
characteristics w^0 at baseline, assuming the optimal treatment choice is also made in the second stage.  
Because A_1 is binary, Q_1(w^0, a_1) can without loss of generality be written as η_1(w^0) +  a_1 Δ_1(w^0),
where the real-valued functions η_1(·) and Δ_1(·) are defined on 𝒳_1. 
It then follows that

    max_a_1 Q_1(w^0, a_1) = η_1(w^0) + Δ_1(w^0) I{Δ_1(w^0) > 0 }.

Taken together, 


the optimal dynamic treatment regime is given by d^opt(s^0)={d^opt_1(w^0),d^opt_2(s^0)}, where 
d^opt_1(w^0)=I{Δ_1(w^0)>0} and d^opt_2(s^0)=I{Δ_2(s^0)>0}. We note that, under standard causal assumptions
and formulated appropriately, Δ_1(w^0) and Δ_2(s^0) are
commonly referred to as the first and second stage blip functions.

A widely used convention in the literature on dynamic  treatment regimes is to respectively model 
Q_2(s^0, a_2; b_2) =  b^⊤_20s_20 + a_2 b^⊤_21s_21
and
Q_1(w^0, a_1; b_1) =  b^⊤_10w_10 + a_1 b^⊤_11w_11,
where s_2k, k=0,1 are sets of variables derived from s^0 (i.e., a realization of S^0) and
w_1k, k=0,1 are sets of variables derived from w^0 (i.e., a realization of W^0).
These model formulations impose restrictive assumptions on both η_j(·) and Δ_j(·), j = 1,2.
However, the decision functions of interest only depend on the Δ_j(·), j = 1,2, or
per the indicated linear models, on the b_j1s. Because misspecification of the models for
the nuisance parameters η_j(·), j =1,2 can induce residual confounding and affect the 
causal interpretation of the interaction terms, we propose a novel modification of the Q-learning
approach that eliminates the need to directly model η_j(·), j=1,2. In particular, 
we adapt techniques originally introduced by  <cit.> for root-n consistent
inference in semiparametric regression models, specifically partially linear models, to the
 Q-learning problem.
 
 

 §.§ Regression model for Stage 2 decision function

 
Define μ_2Y(s^0) = ( Y | S^0 = s^0 ) and μ_2A(s^0) = ( A_2 | S^0 = s^0 );
for now, we will proceed as if these two functions are known. Under 
(<ref>), or equivalently (<ref>),
μ_2Y(s^0) = η_2(s^0) + μ_2A(s^0)  Δ_2(s^0), implying that

    Y_i -  μ_2Y( S^0_i) = { A_2i - μ_2A(S^0_i) }·Δ_2(S^0_i) + ϵ_2i, 
       ( ϵ_2i | S^0_i, A_2i ) = 0,  i=1,…,N.

Define the centered second stage Q-function 

    Q_2c(s^0, a_2) = { Y -  μ_2Y(S^0) | S^0 = s^0, A_2 = a_2 } = 
    { a_2 - μ_2A(s^0) }Δ_2(s^0);

observe that Q_2c(s^0, a_2) = Q_2(s^0, a_2) - μ_2Y(s^0). 
Although Q_2c(s^0, a_2) differs from Q_2(s^0, a_2), the action
a_2 that maximizes Q_2c(s^0, a_2) is the same as that which maximizes
Q_2(s^0, a_2) because μ_2Y(s^0) does not depend on the value of a_2.

With η_2(s^0) eliminated from Q_2c(s^0, a_2; β_2), we 
propose to model (<ref>) via 

    Q_2c( s^0, a_2; β_2 ) = { a_2 - μ_2A(s^0) }s^⊤β_2,

where s denotes the realization of S
(i.e., some finite set of variables derived from S^0).
Because 

    { Y -  μ_2Y(S^0) - Q_2c(S^0, A_2; β_2) | S^0 = s^0, A_2 = a_2 } = 
    { a_2 - μ_2A(s^0) }{Δ_2(s^0) - s^⊤β_2 },

it is easily shown that

    β^*_2 
        =    β_2min [ { Y -  μ_2Y(S^0) - Q_2c(S^0, A_2; β_2) }^2 ] 
        =    β_2min [ { A_2 - μ_2A(S^0) }^2  {Δ_2(S^0) - S^⊤β_2  }^2 ].

The second expression shows that S^⊤β^*_2 is the best (weighted) linear predictor of  Δ_2(S^0).
For data (Y_i,S^0_i),i = 1,…,N, the above developments further show that one can estimate β^*_2 by 
β̃_2N, calculated as the minimizer of

    ∑_i=1^N { Y_i -  μ_2Y(S^0_i) - Q_2c(S^0_i, A_2i; β_2) }^2.

Finally, maximizing Q_2c( s^0, a_2; β^*_2 )
for a_2 ∈{0,1} gives d^opt_2(s; β^*_2) = I{s^⊤β^*_2 >0}
as the optimal model-based treatment decision in Stage 2, with 
a corresponding estimated decision rule of d̃^opt_2(s) = I{s^⊤β̃_2N >0}
(i.e., assuming μ_2Y(·) and μ_2A(·) are known).

The calculations above evidently rely on the availability of μ_2Y(·) and  μ_2A(·). The assumption that
μ_2Y(·) is known is particularly unrealistic; hence, in Section <ref>, we establish the properties of the 
corresponding least squares estimator when these functions are estimated from the available
data using suitable consistent nonparametric estimators, such as those derived from random forests <cit.>
or Super Learner <cit.>.



 §.§ Regression model for Stage 1 decision function



By construction, the relevant first stage Q-function depends on the model for  (<ref>).
Using the partially linear model of Section <ref>, the model-based analog of the second stage Q-function (<ref>)
is given by

    Q_2(s^0, a_2; β^*_2) = μ_2Y(s^0) + Q_2c(s^0, a_2; β^*_2),

where the second term is defined in (<ref>). In view of (<ref>), we therefore re-define the first stage
Q-function of interest as

    Q_1(w^0, a_1) = {max_a_2 Q_2(S^0, a_2; β^*_2)  | W^0 = w^0, A_1 = a_1 }.

Without loss of generality, and arguing similarly to the previous section, 
the fact that A_1 is binary means 

    Q_1(w^0, a_1) = η_1(w^0) + a_1 Δ_1(w^0)

is a saturated nonparametric model for (<ref>); here, a_1 ∈{0,1} and
the real-valued functions η_1(·) and Δ_1(·) are defined on 𝒳_1
and are not necessarily assumed to be the same functions that were initially used to define Q_1(w^0, a_1) a
the beginning of Section <ref>.

Define

    Ỹ^† = max_a_2 Q_2(S^0, a_2; ^*_2);

then, it is not difficult to show that

    Ỹ^†    =    μ_2Y(S^0) + S^⊤^*_2 
    { I(S^⊤^*_2  > 0 ) - μ_2A(S^0) }.

Under (<ref>), 
( Ỹ^† |  W=w^0, A_1 = a_1 ) = η_1(w^0) + a_1 Δ_1(w^0);
similarly,

    μ_1Y(w^0) =  (Ỹ^† |  W=w^0 )
    = η_1(w^0) + μ_1A(w^0)  Δ_1(w^0),

where μ_1A(w^0) = ( A_1 | W=w^0).  Similarly to (<ref>), 
 we can write

    Q_1c(w^0, a_1) = {Ỹ^† -  μ_1Y(W^0) | W^0 = w^0, A_1 = a_1 } = 
    { a_1 - μ_1A(w^0) }Δ_1(w^0)

and, considering  (<ref>), can model (<ref>) via

    Q_1c(w^0, a_1; β_1 ) = { a_1 - μ_1A(w^0) }w^⊤β_1,

where w is defined analogously to s.
Together, these results imply that

    {Ỹ^† -  μ_1Y(W^0) - Q_1c(W^0, A_1; β_1) | W^0 = w^0, A_1 = a_1 } = 
    { a_1 - μ_1A(w^0) }{Δ_1(w^0) - w^⊤β_1 }.

Similarly to the second stage problem, it now follows that

    β^*_1 
        =    β_1min [ {Ỹ^† -  μ_1Y(W^0) - Q_1c(W^0, A_1; β_1) }^2 ] 
        =    β_1min [ { A_1 - μ_1A(W^0) }^2  {Δ_1(W^0) - W^⊤β_1  }^2 ],

the latter implying that W^⊤β^*_1 is the best (weighted) linear predictor of  Δ_1(W^0). 
For data (Y_i,S^0_i),i = 1,…,N, 
and assuming that  β^*_2, μ_jY(·) and μ_jA(·), j = 1,2 are all known,
the above developments further imply that one can estimate β^*_1 using

    β̃_1N = β_1min∑_i=1^N {Ỹ^†_i -  μ_1Y(W^0_i) - Q_1c(W^0_i, A_1i; β_1) }^2.

Parallel to the second stage problem, 
the optimal model-based treatment decision in Stage 1, assuming
the optimal model-based treatment is also given in Stage 2, would be
d^opt_1(w; β^*_1) = I{w^⊤β^*_1 >0},
and may be estimated by d̃^opt_1(w) = 
I{w^⊤β̃_1N >0}.

Of course, none of β^*_2, μ_jY(·) and possibly μ_jA(·), j = 1,2 are known in practice; in Section <ref>, we 
establish the properties of the corresponding least squares estimator when these quantities are all
estimated.

Remark: As an alternative to (<ref>), one can substitute 

    Ỹ = Y + S^⊤β^*_2{ I(S^⊤β^*_2 > 0 ) - A_2 }

for (<ref>) when calculating β̃_1N;
this follows directly from the equalities
( Ỹ | W^0 = w^0, A_1 = a_1 ) = ( Ỹ^†  | W^0 = w^0, A_1 = a_1 )
and 
( Ỹ  | W^0 = w^0 ) = (Ỹ^† | W^0 = w^0 ). 





§ ROBUST Q-LEARNING: ESTIMATION IN PRACTICE AND CORRESPONDING THEORY



The two-stage procedure described in the previous section leads to a class of decision rules indexed by finite-dimensional parameter vectors, that is,
d^opt_1(W;β_1) =I( W^⊤β_1>0) and d^opt_1(S;β_2)=I( S^⊤β_2>0) 
<cit.>. 

Although not explicit in prior developments, the variable sets S and W are each assumed to contain a column of ones, 
so that the main effects of treatment at each stage can be included as part of the decision rule.

The proposed approach eliminates the nuisance parameters η_j(·), j=1,2 from the first and 
the second stage decision rules at the expense of introducing the four additional unknown functions  μ_jY(·) and μ_jA(·), j = 1,2. The advantage
of the proposed  approach is that the indicated functions depend on observables and can be easily estimated using any nonparametric regression or statistical learning method having sufficiently good prediction performance. Importantly, in the case of a sequentially randomized clinical trial, the functions μ_jA(·), j = 1,2 are known and correct models are easily formulated.



 §.§ Estimation in practice


The  developments in the next two subsections assume that the original sample, with N elements independently and identically distributed
as P_0, has been randomly split into two disjoint and independent samples, say D_I_n and D_I^c_n, with 
n = O(N) (e.g., n = N/2) and where I_n and its complement I^c_n form a partition of the index set {1,…,N}.
The induced nuisance parameters μ̂_2Y(·), μ̂_2A(·), μ̂_1Y(·), and μ̂_1A(·)
are to be estimated as described earlier using the data in D_I^c_n;
the finite dimensional parameters of interest are then estimated using the data D_I_n,
treating μ̂_2Y(·), μ̂_2A(·), μ̂_1Y(·), and μ̂_1A(·)
as known functions. 

As described above, the use of such sample-splitting is a particularly simple form of cross-fitting and can be generalized
easily <cit.>; our use of sample splitting as described above will be sufficient to establish the main ideas for
both estimation and asymptotics without unnecessarily complicating notation. 
Generalization to cross-fitting is straightforward and will be discussed at the end of Section <ref>.

Let μ̂_2Y(·), μ̂_1Y(·), μ̂_2A(·), and μ̂_1A(·) denote suitable estimates of 
μ_jY(·) and μ_jA(·), j = 1,2 derived from the data in D_I^c_n.
Backward induction, implemented as described earlier with obvious modifications, 
 can be used estimate the optimal dynamic treatment regime. In particular, for the second stage, we compute

    β̂_2n= 
    β_2min ∑_i ∈I_n[ Y_i-μ̂_2Y(S^0_i) - {A_2i-μ̂_2A(S^0_i)}·S_i^⊤β_2    ]^2.

To estimate the first stage parameters, we first calculate the estimated first stage pseudo-outcome

    Ỹ̂_i
    =Y_i + S_i^⊤β̂_2n{I(S_i^⊤β̂_2n  >0) -A_2i};

and then compute

    β̂_1n= β_1min ∑_i ∈I_n[ Ỹ̂_̂̃̂ĩ̂-μ̂_1Y(W^0_i) - {A_1i-μ̂_1A(W^0_i)}·W_i^⊤β_1    ]^2.

The notation in (<ref>) and (<ref>) emphasizes the fact that the nuisance parameters
μ_2Y(·), μ_2A(·), μ_1Y(·), and μ_1A(·) are estimated using
the outcome and full set of either second and first stage covariates, whereas the linear specifications used for modeling
the centered Q-functions might not use all available covariate information.

As defined, the pseudo-outcomes are non-smooth functions of the data, hence so is β̂_1n; this can cause non-regularity problems for 
β̂_1n <cit.>. In particular, when Pr(| S^⊤β^*_2 | =0)> 0, 
i.e., there exists a strata of the covariates S used to model the Q-function that occurs with positive probability and 
for which treatment is neither beneficial nor harmful, the estimators of first stage regression coefficients become non-regular 
due to the non-differentiability of the indicator function in  the definition of the pseudo-outcome. 

The proposed Q-learning models essentially utilize the propensity score regression approach of <cit.> to eliminate the problem 
of mismodeling hard-to-estimate infinite-dimensional parameters (i.e., η_j(·), j=1,2) on the estimators of the β^*_js. The resulting 
estimator of β^*_j is consistent and asymptotically normal under suitable conditions on μ̂_jA(.) and μ̂_jY(.), j=1,2.  
In particular, the estimate of β^*_j is robust to misspecification of μ_jY(·) provided that μ_jA(·) is consistently 
estimated where β^*_j represents the parameters of the best linear approximation of the unknown Δ_j(·).

In practice, we recommend using ensemble learning methods such as Super Learner <cit.> for estimating 
both μ_jA(·) and μ_jY(·). Asymptotically, Super Learner performs as well as the best convex combination of the 
base learners in the chosen library, in the sense of minimizing the difference in risk compared to the corresponding oracle estimator. 
Moreover, the size of the library can grow at a polynomial rate compared with the sample size without affecting its oracle performance 
 <cit.>.  For these reasons, it is recommended that the library consist of a 
 large and diverse set of regression modeling procedures (i.e., nonparametric, semiparametric, parametric).
Importantly, these theoretical results only imply that Super Learner can match the performance of the (unknown) best possible convex 
combination of choices in the specified library. Thus, consistency is not guaranteed unless the corresponding oracle estimator is consistent 
and converges sufficiently fast. However, with the use of a sufficiently flexible library, Super Learner clearly improves one's ability to 
construct a consistent estimator because  it eliminates the need to select and subsequently rely on a single method of estimation.
 
 
 


 §.§ Theoretical results


To further simplify notation, let Δ_2i = Δ_2(S^0_i), Δ_1i = Δ_1(W^0_i), 
μ_2Ai = μ_2A(S^0_i), μ̂_2Ai = μ̂_2A(S^0_i), μ_1Ai = μ_1A(W^0_i), 
and μ̂_1Ai = μ̂_1A(W^0_i).  In addition, with
x^⊗ 2 = xx^⊤ for any vector x,
define the matrices

    _2n = 1/n∑_i ∈I_n   (A_2i- μ_2Ai)^2  S^⊗ 2_i  
      _2n = 1/n∑_i ∈I_n   (A_2i- μ̂_2Ai)^2  S^⊗ 2_i


    _1n = 1/n∑_i ∈I_n   (A_1i- μ̂_1Ai)^2  W^⊗ 2_i  
      _1n = 1/n∑_i ∈I_n   (A_1i- μ̂_1Ai)^2  W^⊗ 2_i.


Let x_q denote the usual Q- norm of a vector x  for q=1,2,∞. Also, for Z ∼ P for 
some probability measure P, suppose f(·) is any real-valued, P-measurable function; then, we  
define the L^2(P) norm of f(·) as f(Z) _P,2 = {∫ f(ω)^2 d P(ω) }^1/2. For a 
real-valued function h(s^0; D_I^c_n) defined for s^0 ∈𝒮 whose calculation 
may depend on the data contained in D_I^c_n, we can also define the random norm
h(S^0; D_I^c_n) _P_0,2 as the square-root of

    h(S^0; D_I^c_n)  ^2_P_0,2
    = {  h(S^0) ^2_ℙ_n,2| D_I^c_n},

where ℙ_n denotes the empirical measure on D_I_n.


Our results are established under the following assumptions.

 
(i)   The support of W^0 and the conditional treatment effect Δ_1(W^0) are uniformly bounded; 
(ii)  the support of S^0 and the conditional treatment effect Δ_2(S^0) are uniformly bounded; and,
the supports of S and W are uniformly bounded.


 
  (i) μ̂_1A(W^0; D_I^c_n) - μ_1A(W^0) ^2_P_0,2 = o_p(n^-1/2); 
  (ii) μ̂_2A(S^0; D_I^c_n) - μ_2A(S^0) ^2_P_0,2 = o_p(n^-1/2).


  
  (i) μ̂_1Y(W^0; D_I^c_n) - μ_1Y(W^0) ^2_P_0,2 = o_p(1); 
  (ii) μ̂_2Y(S^0; D_I^c_n) - μ_2Y(S^0) ^2_P_0,2 = o_p(1).


   
(i) μ̂_1Y(W^0; D_I^c_n) - μ_1Y(W^0) _P_0,2μ̂_1A(W^0; D_I^c_n) - μ_1A(W^0) _P_0,2 = o_p(n^-1/2);

(ii) μ̂_2Y(S^0; D_I^c_n) - μ_2Y(S^0) _P_0,2μ̂_2A(S^0; D_I^c_n) - μ_2A(S^0) _P_0,2 = o_p(n^-1/2)


 

There exists 1 ≤ n_0 < ∞ such that _jn and _jn, j = 1,2 are positive definite for n ≥ n_0.


 
 
P( | S_1^⊤β^*_2 | = 0 ) = 0.


Assumption <ref> requires no discussion. Assumptions <ref> - <ref> 
impose reasonable conditions on the estimators of the nuisance parameters estimated using cross fitting that are satisified by many 
machine learning algorithms; see <cit.> for further discussion.

Assumption <ref> imposes reasonable conditions on the existence and uniqueness of the least squares estimators
(<ref>) and (<ref>). Assumption <ref> combined with independent, identically distributed sampling ensures
that the limiting matrices

    _2  =  {( A_2 | S^0 ) S^⊗ 2}  _1  =  {( A_1 | W^0 ) W^⊗ 2}

both exist and are positive definite.

Finally, Assumption <ref> is imposed to avoid non-regular asymptotic behavior
in the first stage least squares estimator (<ref>).  Inferences for the parameters that define the
estimated optimal dynamic treatment regime 
d̂^opt_2(s) = I{s^⊤β̂_2n >0}
and
d̂^opt_1(w) = I{w^⊤β̂_1n >0}
can now be derived using the results in the following theorem.


Suppose that Assumptions <ref> - <ref> hold. 


  (a) Let β̂_2n be given by (<ref>). Then,
√(n) ( β̂_2n  - β^*_2 ) 
  d→ N( 0,   ^-1_2 _2 ^-1_2 )
where the matrices _2  =  {( A_2 | S^0 ) S^⊗ 2} 
and _2 =  ( J_2^⊗ 2 ) for

    J_2 = 
    {A_2-  μ_2A(S^0)}S[
    Y - μ_2Y(S^0)  - { A_2-  μ_2A(S^0) }S^⊤β^*_2].


  (b) Let β̂_1n be given by (<ref>). Then,
√(n) ( β̂_1n  - β^*_1 ) 
  d→ N( 0,   ^-1_1 _1 ^-1_1 )
where the matrices _1  =  {( A_1 | W^0 ) W^⊗ 2}
and _1 = { ( J_1 + K^-1_2 J_2) ^⊗ 2} for

    J_1 = {A_1-  μ_1A(W^0)}W[
    Ỹ- μ_1A(W^0) - { A_1 - μ_1A(W^0) }W^⊤β^*_1]

and

    K = [ { A_1-  μ_1A(W^0) }{ I(S^⊤β^*_2  > 0 ) - A_2}WS^⊤].



The following corollary to Theorem <ref> shows
that Assumption <ref> is not required to establish the results
in part (b) in certain settings, in contrast to the standard form of Q-learning. 


Suppose Assumptions <ref> - <ref> hold. In addition, suppose 
( A_1i - μ_1Ai|W_i, S_i,  I_n ) = 0, i ∈I_n.
Then, Theorem <ref>, part (b) remains true.



The set of variables S used for modeling the second stage decision rule are those thought to be potential effect modifiers
for the second stage treatment assigment A_2; hence, a sufficient condition for the Corollary to hold is that A_1 is 
independent of S, conditionally on the set of pre-treatment covariates W included in the first stage model. Note
that this does not preclude the possibility that A_1 can affect variables in S^0 that are not part of S.
Due to the way in which A_1i - μ_1Ai enters the estimating equation of <cit.>, 
it is unclear whether their approach avoids non-regularity under the same conditions as Corollary <ref> even
when the propensity model μ_1A(W^0) has been correctly specified. 



 §.§ Generalization to cross-fitting


Sample splitting, as used in the previous two sections, does not make use of the full sample of N observations
to estimate the finite-dimensional regression parameters, and this can negatively impact efficiency. We now
describe an alternative approach, cross-fitting, that uses the full sample to estimate the desired target parameters.

Suppose that N = n K for some integer n and some integer K ≥ 2. Using an extension of previous
notation, we first randomly split the original sample into 
disjoint (hence independent) samples (D_I_n,k)_k=1^K such that the size of each sample is n=N/K and 
I_n,k, k = 1,…,K partition the indices {1,…,N}.  
Analogously to before, define I_n,k^c as the set of sample indices that are not included in I_n,k; that is, 
I_n,k^c = {1,2,⋯,N}∖I_n,k, k=1,…,K. Then, for each k=1,2,⋯,K, estimate the 
nuisance parameters  μ_2Y(·), μ_2A(·), μ_1Y(·), and μ_1A(·) using the data in D_I_n,k^c;
we respectively denote these estimators 
μ̂_2Y(·;D_I_n,k^c), μ̂_2A(·;D_I_n,k^c), μ̂_1Y(·;D_I_n,k^c), 
and μ̂_1A(·;D_I_n,k^c), k=1,… K. 
Finally, we define

    β̂_2n= 
    β_2min ∑_k=1^K ∑_i ∈I_n,k[ Y_i-μ̂_2Y(S^0_i;D_I_n,k^c) - {A_2i-μ̂_2A(S^0_i;D_I_n,k^c)}·S_i^⊤β_2    ]^2

and

    β̂_1n= β_1min ∑_k=1^K ∑_i ∈I_n,k[ Ỹ̂_̂̃̂ĩ̂-μ̂_1Y(W^0_i;D_I_n,k^c) - {A_1i-μ̂_1A(W^0_i;D_I_n,k^c)}·W_i^⊤β_1    ]^2.

This form of cross-fitting essentially corresponds to “DML2” as described in <cit.>. Like sample splitting,  cross-fitting helps
to guarantee that some of the remainder terms in the asymptotic linearity expansion converge to zero at an appropriately fast rate. However, in contrast
to sample splitting, cross-fitting as described above is also capable of asymptotically achieving the same efficiency as in the case 
where estimators of the regression parameters are computed using all N observations (i.e., with μ_jA(·) and μ_jY(·) j=1,2 being known).




§ SIMULATION STUDIES


 We examined the performance of our proposed Q-learning method under different simulation scenarios 
 with various functional complexities and degrees of non-regularity (i.e., violation of Assumption <ref>). 
 
 
The main simulation in the regular setting uses the following data generation mechanism. 
Let X_1=(X_11,X_12,X_13,X_14,X_15)^⊤ be a 5-dimensional vector of baseline covariates independently
generated and uniformly distributed on [-0.5,0.5]. 

Let  X_2 = (X_21,X_22,X_23,X_24,X_25)^⊤,  X_2l = X_1l + U_l,  l=1,2,3; 
X_24 = 0.35 X_15 + U_4; and, X_25 = U_5, where U_l, l =1,…,5 are independent and uniformly
distributed on [-0.5,0.5]. 


It is assumed that only nonresponders to the first stage treatment will receive the second stage treatment. This nonresponse indicator 
R equals 1 if X_24 is less than its median value and is 0 otherwise. 

Finally, the first and second stage treatments A_j are generated from a Bernoulli distribution with success probability 
μ_jA(·) = [ 1+exp{-λ_jA(·) } ]^-1, where  λ_jA(·)
depends on either S^0 = (X^⊤_1,A_1, X^⊤_2)^⊤ (j=2) or 
W^0 = X_1 (j=1); see Section <ref>.



 §.§ Performance: regular setting
 


In this case, we consider performance for models that satisfy Assumption <ref>.

To implement our proposed method, we used the  package  <cit.> to estimate 
μ_1Y(·), μ_2Y(·), μ_1A(·), and μ_2A(·).  The library used for  included 
generalized linear models (i.e., ), generalized additive models (i.e., ; ), multivariate adaptive 
regression splines (i.e., ; ),
random forests (i.e., ; ),  and  
support vector machines (i.e.,  from the  package ; );
estimation was implemented with all tuning parameters set to their respective default values.  

This simulation study uses four different functional forms for the treatment assignment model
μ_jA(·) = [ 1+exp{-λ_jA(·) } ]^-1, j=1,2:
 

  * Randomized: λ_jA(·) =0

  * Linear: λ_jA(·) =2X_j1+2X_j2+X_j3+0.1X_j4+0.1X_j5

  * Quadratic: λ_jA(·) = 1.4 {(X_j1-0.5)^2+(X_j2-0.5)^2+0.6(X_j3-0.5)^2+0.5(X_j4-0.5)^2+0.5(X_j5-0.5)^2+X_j1+X_j2+0.6X_j3+0.5X_j4+0.5X_j5-2}.

  * InterQuad: λ_jA(·) = 1.4 {(X_j1-0.5)^2+(X_j2-0.5)^2+0.6(X_j3-0.5)^2+0.5(X_j4-0.5)^2+0.5(X_j5-0.5)^2+X_j1+X_j2+0.6X_j3+0.5X_j4+0.5X_j5+X_j1X_j2-2}.

The Randomized model corresponds to a SMART-like trial where simple randomization is used at baseline and then
simple re-randomization occurs among the set of non-responders at the first stage.  The randomization model, part of the trial design,
is therefore known and the inclusion of an appropriate  model in the Super Learner library should ensure that 
μ_jA(·) can be consistently estimated at the usual parametric rate. 

The other three settings are meant to correspond to increasingly complex observational data settings,  where the “assignment” mechanism by 
which patients follow a particular treatment regimen is covariate-dependent, not randomized, and is not considered to be known by design. 

Hence, the analyst cannot knowingly select a correctly specified parametric model a priori. 

In the case of the Linear model, the inclusion of a  model in the Super Learner library again 
ensures that μ_jA(·) can be consistently estimated at the usual parametric rate. For the
other two models, the inclusion of methods such as  and  will help to mitigate, 
but not necessarily eliminate, the possibility of inconsistent estimation. 
These observations highlight the importance of using flexible methods when modeling μ_jA(·), j = 1,2, 
particularly in observational data settings.

The outcome models are given by :


  * Linear^R: Y= _1^⊤_1+_2^⊤_2+A_1_1^⊤_1+ A_2 R _2^⊤_2+ϵ 
where _1=_2=(1,0.1,0.1,0.1,0.1)^⊤, _1=(0,0,0,0,0)^⊤ and _2=(1,1,0,0,0)^⊤; 

  * FGS^R: Y= f(_1)+f(_2)+A_1_1^⊤_1+A_2 R g(_2)+ϵ where _1=(0,0,0,0,0)^⊤ and
for x = (x_1,x_2,x_3,x_4,x_5)^⊤, we set  g(x ) = 2 sin(π x_1 x_2) + 2
  (x_2-0.5)^2 and

    f(x) = -1.5 + sin(π x_1 x_2) + 2
      (x_3-0.5)^2 +  x_4 + 1.5 x_1/|x_2|+|x_3| +2 x_1 (x_2+x_3).


The noise variable ϵ is generated from N(0,σ=0.5).  

In connecting the above Linear^R outcome model specification with earlier notation, we have  
S^0 = (X_1^⊤ ,A_1,X_2^⊤)^⊤, η_2(S^0) = _1^⊤_1+_2^⊤_2 + A_1_1^⊤_1 
and Δ_2(S^0) = R _2^⊤_2, where R  is a function of X_24 only; we further have W^0 = X_1. 

We respectively use S = R (1,X_21, X_22,X_23)^⊤ and W = (1,X_11,X_12)^⊤ for modeling the relevant Q-functions.

In this case, the target of estimation ^*_2 = (0,θ_21,θ_22,θ_23)^⊤ and
it can additionally be shown that (<ref>) coincides with (<ref>).

However, for the FGS^R outcome model, η_2(S^0) = f(_1)+f(_2)+A_1_1^⊤_1 and 
Δ_2(S^0)  = R g(_2); here, (<ref>) does not coincide with (<ref>) since the 
linear parametric specification used in the latter is not
equal to Δ_2(S).  In this case S^⊤^*_2 still exists as the best linear 
projection of Δ_2(S^0) on to the linear space spanned by S; however, 
its value for this simulation study must be determined  numerically (e.g., through simulation). 


In general, it is not similarly straightforward to characterize the functions η_1(W^0) and Δ_1(W^0), or the value of ^*_1 
in the first stage models, without appealing to numerical methods. However, in the current simulation setting, the value of ^*_1 can 
be determined exactly for both the Linear^R and FGS^R outcome model specifications. Specifically, 
neither model involves an interaction between A_1 and A_2; more generally, there is no correlation between A_1 and the
second stage variables A_2 and _2. As a result, the linear term  A_1_1^⊤_1 that appears in both the 
Linear^R and FGS^R outcome model specifications accurately describes the interaction between 
treatment A_1 and _1 in the true first stage Q-function (i.e., Δ_1(W^0) = W^⊤_1 = _1^⊤_1 = 0).
It follows that ^*_1 = 0 and hence that expression (<ref>) also coincides with (<ref>).
 





In our main simulation study, there are 8 possible model combinations represented by the outcome and treatment assignment models,
and within each setting we compare the performance of the proposed method for estimating _j^*,j=1,2
to the standard form of Q-learning (Q_N,N) and also to the weighted least squares (dWOLS_N,N) estimator 
proposed by <cit.>. The subscripts on these latter two estimators denote the fact that 
standard errors would normally be calculated using the N-out-of-N bootstrap (i.e., in the regular setting). 
In the case of dWOLS_N,N, linear models are used for the relevant Q-function model specification and logistic 
regression models are used for estimating the treatment assignment probabilities.

The estimation of _j^*,j=1,2  is not subject to residual confounding
bias for any of the proposed methods under the Linear^R outcome model specification.

However, there is a possibility of such bias under the FGS^R in the case of Q_N,N and
dWOLS_N,N. To be more specific, residual confounding bias under the FGS^R outcome model
is expected for Q_N,N regardless of the treatment assignment model.


For dWOLS_N,N, the Randomized and Linear first and second stage treatment assignment models are 
correctly specified and easily modeled. Hence, under the FGS^R outcome model specification, a significant potential for bias 
arises only under the  Quadratic or InterQuad treatment assignment rules.  

For the proposed method, residual confounding bias when estimating _j^* is not anticipated provided 
that μ_jA(·), j=1,2 are sufficiently well-estimated.

We generate 500 datasets of size 2000 to examine the performance of our proposed method
and use cross-fitting as described in Section <ref> with K=2 to estimate the desired target parameters.
Tables <ref> and <ref> show the empirical absolute bias and standard deviations of the 
second and first stage parameter estimates (i.e., standard errors). The values of β^*_2,1≈ 0 and β^*_2,2≈ -2 
are determined by simulation.  As expected, standard Q-learning performs poorly except under the Randomized  
treatment assignment model.  The proposed method and dWOLS_N,N also perform similarly well under the 
Randomized and Linear treatment assignment models for estimating the first and second stage parameters. 
However, under the  FGS^R outcome model, the proposed method exhibits similar biases and substantially smaller
standard errors.

For the Quadratic and InterQuad treatment assignment mechanism, both of which are mis-modeled 
in the case of dWOLS_N,N, the corresponding estimators show substantial bias in some of the parameters, 
whereas those for the proposed method remain comparatively low.
For example, under the  InterQuad treatment assignment model and FGS^R outcome model,  the proposed method 
respectively results in estimators for β^*_21 and β^*_22 with absolute biases of 0.070 and 0.019; in contrast, 
those for the dWOLS_N,N estimators are 0.758 and 0.455, respectively.  We again see a substantial reduction in standard 
errors; in this same example, the standard errors under the proposed method are 0.612 and 0.271, whereas for dWOLS_N,N these
are respectively 0.916 and 0.402, the degree of reduction exceeding 30%. Overall, the proposed method is observed to be
more robust,  typically producing less biased estimators with smaller standard errors compared with the other two approaches.  

The performance of our proposed method was also assessed using smaller sample sizes. Tables S1-S6 in the supplementary material respectively show the 
results for _2^* and _1^* with sample sizes of 1000, 500, and 250. Overall, the proposed method continues to outperform both Q_N,N and dWOLS_N,N, particularly when the underlying treatment assignment and the outcome models are both nonlinear (i.e., settings in which bias
can be expected for both Q_N,N and dWOLS_N,N.). However, the performance of the proposed method is also affected by sample size. 
For example, under the Linear treatment assignment and FGS^R  outcome model with a sample size of N=250, the proposed method 
shows unacceptably high bias when estimating β^*_2 when compared to dWOLS_N,N; see Table S5. We conjecture that this 
occurs because the information available for estimating the second stage propensity model is limited to the set
non-responders (i.e., 50% of the sample) that are re-randomized.  The value functions for the estimated rules in all cases were
calculated for all sample sizes and show that the proposed method, followed by dWOLS, typically results in value functions that are closest to optimal; 
see Section 9.4 of the supplementary materials.

The supplementary material includes simulation results in which  is replaced by alternative data adaptive techniques. Specifically, in 
Tables S7 and S8 in the supplementary material, the columns RF-RF and GAM-GAM represent modeling approaches in which 
and   are used for both the marginalized outcome (i.e., μ_1Y(·) and μ_2Y(·)) and treatment assignment models  
(i.e., μ_1A(·) and μ_2A(·)). The column RF-GAM instead uses  for the outcome model and 
for the treatment assignment model.  Comparing these results with those summarized in Tables <ref> and <ref> 
shows that the use of  improves performance.







 §.§ Performance: non-regular setting
 


The treatment assignment models considered here are respectively Randomized, Linear and InterQuad, 
defined as in Section <ref>.

Additionally, define 
X̃_2=(X̃_21,X̃_22,X_23,X_24,X_25)^⊤ where 
X̃_21 is generated from a Bernoulli distribution with 
success probability [ 1+exp{-(2X_11+2X_12-1)} ]^-1, 
X̃_22 is generated from a Bernoulli distribution with 
success probability [ 1+exp{-(2X_12+X_21-1)} ]^-1, 
X_23 = U_1,
X_24 = 0.35 X_15 + U_2,
and X_25 = U_3,
where U_l, l =1,…,3 are independent and uniformly
distributed on [-0.5,0.5]. 

We consider the following outcome models:


  * Linear^NR,ϖ: Y= _1^⊤_1+ _2^⊤_2+A_1_1^⊤_1+  A_2 · (θ_2  R  X̃_21)+ϵ where _1=_2=(1,0.1,0.1,0.1,0.1)^⊤, _1=(0,0,0,0,0)^⊤ and θ_2 = 2 ϖ.

  * Non-linear^NR,ϖ: Y= f(_1)+A_1_1^⊤_1+  A_2 · (θ_2 R  X̃_21)+ϵ where _1=(0,0,0,0,0)^⊤,
θ_2 = 2 ϖ and, for x = (x_1,x_2,x_3,x_4,x_5)^⊤, we set

    f(x) = -1.5 + sin(π x_1 x_2) + 2
      (x_3-0.5)^2 +  x_4 + 1.5 x_1/|x_2|+|x_3| +2 x_1 (x_2+x_3).

   
The noise variable ϵ is generated from N(0,σ=0.5) and the constant ϖ∈{0,1} specifies the degree of non-regularity,
as will be discussed further below. 

In the above models, S^0 = (X^⊤_1,A_1,X̃^⊤_2)^⊤, W^0 = X_1,
Δ_2(S^0) = θ_2 R  X̃_21, η_2(S_0) is determined by the remaining model 
terms, and R  is a function of X_24 only.
The second and first stage Q-functions are respectively modeled as linear functions of
S= R (1, X̃_21, X̃_22, X_23)^⊤ and W = (1, X_11, X_12)^⊤. 
For both models, it is not difficult to show that β^*_2 = (0,2 ϖ,0,0)^⊤ and that
^*_1 = 0.


In both scenarios, for each subject i, the first-stage pseudo outcome is defined as in 
(<ref>) and estimated by substituting in β̂_2n for
β^*_2.   The construction of the pseudo-outcome, specifically the projection 
S^⊤^*_2, violates  Assumption <ref>. 
 
In particular,  ϖ=0 corresponds to no second-stage effect modifier, implying that 
P( |S^⊤^*_2| = 0) = 1 because ^*_2 = 0.
Setting ϖ=1 instead implies that there is no second-stage treatment effect
when R X̃_21 = 0, and a reasonably strong
effect when R X̃_21 = 1; in this case,
0 < P( |S^⊤^*_2| = 0) < 1. 
However, the conditions of Corollary <ref> hold in each case because
S does not include A_1, resulting in regular asymptotic behavior for the proposed method. 

Because these simulations focus on coverage rather than bias and standard error, we simulate 1000 datasets of size N=2000.
In the non-regular setting considered here, neither Q_N,N nor  dWLOS_N,N can necessarily be expected to perform well;
hence, we compare our proposed method to a modified version of standard Q-learning and doubly robust weighted least squares
in which the first stage confidence intervals are respectively constructed using a m-out-of-N bootstrap technique
as developed in <cit.> (i.e., Q_m,N^κ) and <cit.> (i.e.,
dWLOS_m,N^κ). 

In both of these approaches, the tuning parameter κ∈ [0,1) determines
the bootstrap sample size m; here, κ=0.05.   

Table <ref> summarizes the results; for comparison, results obtained using the N-out-of-N bootstrap in the first stage
are provided in Table S9 in the supplementary material. 

In these tables, empirical coverages that are significantly over or under the nominal level 0.95 are indicated with a dagger,
with significance being assessed using a binomial test.


The performance of both Q_N,N and Q_m,N^κ=0.05 relies heavily on the correct 
specification of the outcome model. In those cases where both methods are observed to 
exhibit reasonable performance, Tables <ref> and S9 respectively show that 
Q_m,N^κ=0.05 typically over-covers whereas 
Q_N,N either under-covers or has close to nominal coverage; 
in contrast, when Q_N,N is observed to under-cover to a very significant 
extent, so does Q_m,N^κ=0.05.


In general, both dWLOS_m,N^κ=0.05 and the proposed method lead to significant improvements in performance. 

Indeed, the proposed method produces valid confidence intervals with coverages close to the nominal level throughout
Tables  <ref> and S9. This can be readily explained by the fact that each setting satisfies the assumptions of Corollary <ref> despite 
violating Assumption <ref>.

Similarly, we see that dWLOS_m,N^κ=0.05 performs reasonably well regardless of the outcome model for both the Randomized
and Linear treatment assignment models, since in these two cases the latter can be consistently estimated at a parametric rate.
However, compared to the proposed method, the coverages tend to be slightly conservative, with longer confidence intervals. In these same
cases, dWLOS_N,N also performs reasonably, though does have a tendency to under-cover.
Under the InterQuad treatment assignment model, the performance of both dWLOS_m,N^κ=0.05 and
dWLOS_N,N declines due to misspecification of the treatment assignment model,
and in the Non-linear^NR,ϖ setting, also the outcome model.

For example,  dWOLS_m,N^κ=0.05 exhibits coverage rates as low as 87%. We conjecture that the combination of 
non-regularity, model  misspecification and residual confounding are the main reasons for the poor performance of Q_N,N,
Q_m,N^κ=0.05 and, where observed to be poor, both dWLOS_m,N^κ=0.05 and dWLOS_N,N.

In comparing the two approaches to bootstrapping for both standard Q-learning and dWOLS, our results 
further suggest that tuning m differently (i.e., increasing m) may result in better agreement with the nominal 
coverage level in cases where the relevant models are appropriately specified.



Finally, we conducted a related simulation study in which both Assumption <ref> and the conditions of Corollary <ref>
are violated. Unlike the simulation settings above, this example considers a case where the first and second stage treatments
interact with each other. This modified study is described in Section 9.3 of the supplementary document, where we compare
the proposed approach with dWLOS_m,N^κ=0.05; the results are summarized in Table S10. Overall, the methods
perform as expected. In particular,  the proposed method demonstrates either nominal or modest undercoverage for the first stage  regression
parameters and dWLOS_m,N^κ=0.05 demonstrates conservative coverage except in cases where 
the required conditions for consistency are violated.








§ APPLICATION


We use the data from the Extending Treatment Effectiveness of Naltrexone (ExTENd) clinical trial to illustrate our method. Naltrexone  (NTX) is an opioid receptor antagonist used in the prevention of relapse to alcoholism. Even though NTX has been shown to be efficacious in those that adhere to treatment, its use by clinicians has been limited, at least in some cases, because  adherence rates are often negatively impacted by the fact that NTX diminishes the pleasurable effects of alcohol use.



In the ExTENd  study (Figure <ref>), at the first decision stage, patients were randomized to one of two definitions of non-response while receiving NTX: 
(1) Stringent: a patient is a non-responder if (s)he has two or more heavy drinking days in the first 8 weeks (A_1=0); (2) Lenient: a patient is a non-responder 
if (s)he has five or more heavy drinking days in the first 8 weeks (A_1=1). At the second decision stage, the treatment assignment mechanism depends on response status.

Specifically, define A_2 = 1 if the current treatment (NTX) is augmented, and zero otherwise; in addition, we let
R̅ denote the indicator of response to treatment. Then, among responders (R̅=1), patients are randomized 
(with equal probability) to augment NTX with telephone disease management (NTX+TDM; A_2 = 1) or to maintain NTX alone (A_2 = 0). For 
non-responders (R̅=0), patients are instead randomized (with equal probability) to augment NTX with combined behavioral intervention 
(NTX+CBI; A_2 = 1) or to CBI alone (A_2 = 0). In the latter case, maintenance on NTX alone is replaced with an alternative treatment
due to non-response.


The primary outcome is the proportion of abstinence days over 24 weeks. The list of baseline and time varying variables that are used in our analyses are given in Table <ref>. There are multiple measurements of time-varying variables during the first stage. We denote the average of these variables as mcs_1, pacs_1, pdhd_1, and apc_1. 

Standardized differences in means for each covariate (i.e.,  differences in means divided by the corresponding pooled standard deviation) were used to check the covariate balance across the treatment groups. Figure <ref> indicates that there is a good balance of baseline covariates across the levels of A_1 (circle). However, we see some imbalance across the levels of second stage treatment options. This is more evident in the non-responder group (triangle). Absolute standardized differences exceeding 0.1 or 0.2 are respectively referred to as mild and substantial imbalance, and can potentially induce bias in the evaluation of effect modifiers 
if not taken into account <cit.>. In this figure, sdApc_1 and sdPdhd_1 respectively represent the standard deviation of the 
indicated variables during the first stage.






We analyzed the data using the proposed method, dWLOS_m,N^κ and  Q-learning (Q_m,N) approaches; 
the results are summarized in Table <ref>. The latter two methods use the N-out-of-N bootstrap
for calculating standard errors in the second stage model and the m-out-of-N bootstrap
for calculating standard errors in the first stage model.

Referring to earlier notation, the first stage covariate vector W^0 consists of the predictors
gender, race, alcyr_0, intox_0, and ocds_0, and the second stage covariate vector S^0 consists
of all the predictors listed in Table <ref>, along with response status.

First stage regression models are fit using W = W^0. The description of the second stage model
predictor S is more involved.  

Specifically, let Z = (Z_1, …, Z_6)^⊤ contain the variables  gender,  A_1, intox_0, ocds_0, pacs_1, and mcs_1;
then, we define S = (R̅, 1-R̅, R̅ Z_1, R̅ Z_2, (1-R̅) Z^⊤)^⊤. 

As specified, the second stage model allows the set of possible effect modifiers to differ between responders and non-responders,
with some overlap in the case of gender and A_1.

We used  to estimate μ_1Y(𝐰^0) and μ_2Y(𝐬^0),  employing the same library as 
we did in the simulation study and respectively using W^0 and S^0 for general confounding control.

In view of the fact that the randomization mechanism is known, and mostly successful in view of the overall degree of balance 
observed in Figure <ref>, the treatment propensities μ_1A(𝐰^0) and μ_2A(𝐬^0)
are estimated using logistic regression models. Specifically, the former is estimated as a function of gender, and the 
latter is estimated using gender, response status, and the interaction between gender and response status.

The parameters of the Q-functions used by dWOLS are assumed to follow linear models (i.e., including the main effects). 
Similarly, for standard Q-learning, linear working models respectively replace η_2(S^0) and η_1(W^0).  


As shown in Table <ref>, the signs of all predictor effects are the same for all methods, though magnitudes and 
confidence intervals differ. None of the effect modifiers in the second stage are deemed statistically significant among 
responders using any of the 3 methods. The proposed Q-learning method suggests that both ocds_0 and mcs_1 
are significant effect modifiers of A_2 among non-responders; specifically, individuals with higher ocsd_0 and mcs_1 
would benefit from CBI. Similarly, dWOLS identifies mcs_1 as a significant effect modifier among non-responders, whereas 
none of the effect modifiers are identified as significant using standard Q-learning.

For the first stage model, the proposed Q-learning method shows that the years of drinking (i.e., alcyr_0) and gender 
significantly modify the effect of A_1. In particular,  female individuals and those with more years of drinking would benefit from 
being treated under a stringent definition of non-response. This makes sense because, for example, individuals with more years 
of drinking at baseline likely have a higher craving for alcohol and require more immediate attention and rescue treatments 
(i.e., A_2 for non-responders). 
 
In contrast, neither standard Q-learning nor dWLOS detects any effect modifiers.  With the exception of
the interaction between A_1 and gender, the first stage point estimates are rather similar across the 3 methods,
highlighting the fact that the differences in significance stem from the tighter confidence intervals obtained
using the proposed methods (i.e., compared to those produced using the m-out-of-N bootstrap).   
 
Indeed, for both stages, the dWLOS and standard Q-learning methods yield point estimates that are mostly similar 
to each other.

These similarities are expected for two reasons. First, under successful randomization, we would not generally expect misspecification 
of the functional form of the main effects in the Q-function to bias the estimate of interaction terms (i.e., β^*_2). Second, the 
linear models being used in the Q-models are identical in both cases; the only difference is that parameter estimation is carried out using 
weighted versus unweighted least squares.   

Figure <ref> demonstrates the presence of random confounding among non-responders 
in the second stage for gender, ethnicity, A_1, pacs_0 and pacs_1, whereas there is good balance 
among the first stage predictors. Comparing dWLOS and standard Q-learning, we see that the largest differences 
in point estimates occur in the second stage model among non-responders for A_1 and pacs_1.


Comparing the proposed method to both dWOLS and standard Q-learning, we observe somewhat greater 
disparity in point estimates. These differences occur primarily among non-responders in the second stage,
and include interactions between A_2 and each of gender,  A_1 and pacs_1; as noted above, 
gender,  A_1 and pacs_1 all demonstrate substantial imbalance among non-responders in the second 
stage. The largest difference among regression coefficients in the first stage model occurs for gender, consistent 
with the disparities observed in the second stage model as well as propagation of those differences through 
the backward induction process.  We conjecture that modeling the 
true main effects (i.e., μ_jY(·), j=1,2) using Super Learner may help to reduce small sample biases 
when compared to the more restrictive linear models used by both dWOLS and standard Q-learning.
 
 






§ DISCUSSION


Much of the current work on Q-learning continues to involve parametric working models despite the fact that
finite-dimensional models are generally too restrictive to permit consistent estimation of nuisance parameters. 
We proposed a robust Q-learning approach where the working models need not all be linear and, specifically, 
where the main effects that do not influence the optimal decision rules are estimated using data-adaptive approaches.
Our simulation studies highlight the value of our proposed approach compared with existing Q-learning methods. The 
proposed method also performed relatively well in simulations when key regularity assumption 
(i.e., Assumption <ref>) is violated; however, we cannot expect this in all scenarios, as the underlying theory 
and simulation results show  otherwise.

An important advantage of the proposed method is that it does not suffer from the curse of dimensionality, 
as it produces root-n consistent estimators even when μ̂_jA(·) and μ̂_jY(·) (j=1,2)
are estimated at  rates slower than root-n. This important property facilitates the use of nonparametric methods
like Super Learner for estimating these unknown functions, substantially reducing the chance of model misspecification. 
A second important feature of the proposed approach is that consistent estimation of the treatment models leads 
to consistent estimation of the blip function parameters, whether or not these models or those for μ_jY(·), j=1,2 
are correctly specified. However, the proposed estimators are not doubly robust, in that we require that the 
μ_jA(·)s are consistently estimated at a sufficiently fast rate. True double robustness 
under (<ref>) for ^*_2 requires that one either consistently estimates μ_2A(S^0) or the treatment-free conditional 
mean model E[ Y - A_2 S^⊤^*_2 | S^0] = η_2(S^0); similarly, for ^*_1,
one must either consistently estimate μ_1A(W^0) or E[Ỹ - A_1 W^⊤^*_1 | W^0] = η_1(W^0).
Because the expectation operator is linear, correct specification of both treatment-free models essentially relies on
both μ_jY(·) and μ_jA(·), j=1,2 being correctly specified. This limitation on the practicality of finding
a truly doubly robust estimator applies  to the proposed approach as well as that taken in <cit.>. 
Further research on doubly robust estimation in this class of problems is merited.

Although data-adaptive estimation methods reduce the risk of inconsistency, there is still a chance that one or more nuisance parameters will be
estimated inconsistently. Further research is needed to study the behavior of the proposed methods under inconsistent estimation of a nuisance parameter. 
In particular, <cit.> showed that when nuisance parameters are estimated using data-adaptive approaches, inconsistently estimating 
one nuisance parameter may lead to an irregular estimator having a convergence rate slower than root-n. These authors proposed a targeted minimum 
loss-based approach to resolve the issue <cit.>. Generalization of the method of <cit.> to a multi-stage  decision making 
process would be an interesting topic for future research. Studying the asymptotic behavior of an appropriate version of the bootstrap in our proposed Q-learning 
method is also of interest as it can potentially resolve the non-regularity issues in settings where both Assumption <ref> fails and  Corollary <ref>  fail <cit.>. Finally, in practice, there are  often many candidate variables to be considered when constructing a decision rule. The inclusion of spurious variables in these analyses can substantially reduce the quality of the estimated decision rules. Although one can adapt the proposed methods to obtain regularized estimators of the target parameters, valid post-selection inference remains a challenge and merits further research <cit.>. 


Biometrika












§ SUPPLEMENTARY MATERIAL TO ROBUST Q-LEARNING


Let Z ∼ P for some probability measure P and suppose f(·) is any real-valued, P-measurable function; then, we  
define the L^2(P) norm of f(·) as f(Z) _P,2 = {∫ f(ω)^2 d P(ω) }^1/2. 
In addition, let x_q denote the usual q- norm of a vector x  for q=1,2,∞.
The following general lemmas will be helpful  in our proofs.


Let A_n and B_n be sequences of random vectors, n ≥ 1.  Let ϵ > 0 be arbitrary and, for any vector norm, 
suppose that lim_n →∞ P(  A_n  > ϵ | B_n)  = 0. Then, 
lim_n →∞ P(  A_n  > ϵ)  = 0. By Chebyshev's inequality, a sufficient condition for proving that lim_n →∞ P(  A_n  > ϵ)  = 0 
is that lim_n →∞ E(  A_n ^q |  B_n)  = 0 from some q ≥ 1. 


The above lemma essentially repeats Lemma 6.1 in <cit.>

and will not be proved here.  The following lemma is a direct consequence of a well-known result and also has an easy proof; see, for example, 
<cit.>.


Let M_n and M̂_n be two sequences of square matrices and let · be any proper matrix norm. 
Suppose  there exists n_0 < ∞ such that  (i) M^-1_n and M_n exist for n ≥ n_0, with 
 0 < C_1 ≤M^-1_n ≤ C_2 < ∞; 
and, (ii) M̂_n -  M_n ≤ (2 M^-1_n )^-1.
Then, 

    M̂^-1_n -  M^-1_n ≤ 2  C^2_2 M̂_n -  M_n .




We will also have need of the following lemma.

Let B_1,…,B_N be independent, identically distributed vectors from P_0, where B_i ∈ℬ⊂^d. 
Let I_n be a randomly chosen subset of the integers 1,…,N of length n = O(N) and let its complement I^c_n have N-n = O(N) elements. 
Let F_I_n and F_I^c_n be the corresponding disjoint subsets of B_1,…,B_N. Let γ_j: ℬ→, j=1,2 and let γ̂_j(·; F_I^c_n) be an estimator of γ_j(·) derived from the data F_I^c_n. 
Finally, define 

    G_n,N  =  1/n∑_i ∈I_n{γ̂_1(B_i; F_I^c_n) - γ_1(B_i) }{γ̂_2(B_i; F_I^c_n) - γ_2(B_i) }h( B_i )

where h(B_i ) is any finite dimensional vector- or matrix-valued function of B_i such that P(h(B_i ) _∞≤ C) = 1
for some C < ∞. Then,

    G_n,N_∞≤ C γ̂_1(B; F_I^c_n) - γ_1(B) _ℙ_n,2γ̂_2(B; F_I^c_n) - γ_2(B) _ℙ_n,2,

where ℙ_n is the empirical measure on F_I_n. Moreover, for j=1,2
 define 

    γ̂_j(B; F_I^c_n) - γ_j(B) ^2_P_0,2 
    = {γ̂_j(B; F_I^c_n) - γ_j(B) ^2_ℙ_n,2| F_I^c_n}

and suppose (<ref>) 
is o_p(N^-a_j), where a_j ≥ 0. Then,
G_n,N_∞ = o_p(N^-(a_1+a_2)/2).



Let r_ij(F_I^c_n) = γ̂_j(B_i; F_I^c_n) - γ_j(B_i) for i=1,…,n and
j=1,2. Under the assumption that
P(h(B_i ) _∞≤ C) = 1, 
the triangle and Cauchy-Schwarz equalities imply

    G_n,N_∞   ≤   1/n∑_i ∈I_n| r_i1(F_I^c_n) r_i2(F_I^c_n)  |
     h( B_i ) _∞
       ≤    C
    [ 1/n∑_i ∈I_n{ r_i1(F_I^c_n) }^2 ]^1/2[ 1/n∑_i ∈I_n{ r_i2(F_I^c_n) }^2 ]^1/2,

the representation (<ref>) now following immediately from the definition of the L^2(Q) norm given earlier 
using Q = ℙ_n. 

To establish that G_n,N_∞ = o_p(N^-(a_1+a_2)), 
we first use Markov's inequality: for any ϵ > 0,

    P ( G_n,N_∞ > ϵ | F_I^c_n) ≤ϵ^-1( G_n,N_∞  | F_I^c_n).

Using (<ref>) and the Cauchy-Schwarz inequality again, it follows that

    P ( G_n,N_∞ > ϵ | F_I^c_n)    ≤   C/ϵ{∏_j=1^2 γ̂_j(B; F_I^c_n) - γ_j(B) _ℙ_n,2| F_I^c_n}
       ≤   C/ϵ∏_j=1^2
    {γ̂_j(B; F_I^c_n) - γ_j(B) ^2_ℙ_n,2| F_I^c_n}^1/2
        =    C/ϵ∏_j=1^2 γ̂_j(B; F_I^c_n) - γ_j(B) _P_0,2

the last result following directly from  (<ref>).

Under the stated assumptions, the right-hand side is now seen to be
o_p(N^-(a_1+a_2)/2), as desired. 


The statement and proof of Lemma <ref> employs a simple form of sample splitting in which the unknown function γ_j(·) is
estimated by γ̂_j(·; F_I^c_n) from a sample F_I^c_n that is independent of the
B_is (i.e., F_I_n) appearing in the calculation of  (<ref>). Importantly, Lemma
<ref> does not preclude the possibility that γ_1(·) = γ_2(·) and
γ̂_1(·) = γ̂_2(·); in this case, (<ref>) reduces to

    L_n,m  =  1/n∑_i ∈I_n{γ̂(B_i; F_I^c_n) - γ(B_i) }^2 h( B_i )

and G_n,N_∞ = o_p(N^-a) for some a ≥ 0 provided that
 
    γ̂(B; F_I^c_n) - γ(B) ^2_P_0,2 = o_p(N^-a).

A related lemma now follows.


Let (R_1,B_1),…, (R_n,B_N) be independent, identically distributed vectors from P_0, where B_i ∈ℬ⊂^d
and E_i ∈. Suppose ( R_i | B_i ) = 0 and  ( R_i | B_i ) = θ^2_i ∈ (0,C_1] for i = 1 … N and
a constant C_1 < ∞.
Let I_n be a randomly chosen subset of the integers 1,…,N of length n = O(N) and let its complement 
I^c_n have N-n elements.
Let F_I_n and F_I^c_n be the corresponding disjoint subsets of (R_1,B_1),…,(R_N,B_N). 
Let γ: ℬ→ and let γ̂(·; F_I^c_n) be an estimator of γ(·) 
derived from the data F_I^c_n.  Finally, define 

    L_n,N  =  1/n∑_i ∈I_n R_i 
    {γ̂(B_i; F_I^c_n) - γ(B_i) }h( B_i )

where h(B_i ) is any finite dimensional vector-valued function of B_i 
such that P(h(B_i ) _2≤ C_2) = 1 for C_2 < ∞. 
Suppose that

    γ̂(B; F_I^c_n) - γ(B) ^2_P_0,2 
    = {γ̂(B; F_I^c_n) - γ(B) ^2_ℙ_n,2| F_I^c_n}

is o_p(N^-a) where a ≥ 0. Then,
L_n,N_∞ = o_p(N^-(1+a)/2).



The proof  relies on a variant of Chebyshev's inequality. Let

    L_n,N,j = 1/n∑_i ∈I_n R_i 
    {γ̂(B_i; F_I^c_n) - γ(B_i) }h_j( B_i )

be the j^th element of L_n,N.  Let ℬ_n = { I_n,  ( B_k, k ∈ I_n) }. Then,
it is easy to show that

    ( L_n,N,j| F_I^c_n) 
     = {( L_n,N,j| F_I^c_n,  ℬ_n ) 
    | F_I^c_n} = 0;

this follows from calculating the inner expectation on the right-hand-side and using
the assumption that ( R_i | B_i ) = 0 for every i. 
Using a similar conditioning argument, 

    ( L_n,N,j| F_I^c_n)     =    {( L_n,N,j| F_I^c_n,  ℬ_n ) | F_I^c_n}.

Straightforward calculations now show

    ( L_n,N,j| F_I^c_n,  ℬ_n )
        =    1/n^2∑_i ∈I_n{γ̂(B_i; F_I^c_n) - γ(B_i) } ^2
    h^2_j( B_i ) (R_i | B_i ),

implying that

    ( L_n,N,j| F_I^c_n) 
        =    [ 1/n^2∑_i ∈I_n{γ̂(B_i; F_I^c_n) - γ(B_i) }^2 
    h^2_j( B_i ) θ^2_i | F_I^c_n] 
       ≤   C_1 C_2/n[ 1/n∑_i ∈I_n{γ̂(B_i; F_I^c_n) - γ(B_i) }^2 | F_I^c_n] 
        =    C_1 C_2/nγ̂(B; F_I^c_n) - γ(B) ^2_P_0,2
        =     o_p(N^-(1+a)),

the last step following from the assumptions on (<ref>) made in the statement of
the lemma and the fact that n = O(N).  Using a vector form of Chebyshev's inequality,
it can then be shown that L_n,N_2 = o_p(N^-(1+a)/2); since
L_n,N_∞≤L_n,N_2, the stated result follows.





§ PROOF OF THEOREM 1


To review our main assumptions, we assume that we observe N independently identically distributed trajectories of 
(X_1,A_1,X_2,A_2,Y) ∼ P_0. The vector X_1 ∈𝒳_1 ⊂^p_1 
consists of baseline covariates measured before treatment at the first decision point A_1∈{0,1} and the vector 
X_2 ∈𝒳_2 ⊂^p_2 consists of intermediate covariates measured before treatment at the second decision point A_2∈{0,1}. 
For notational convenience we define  S^0_i=(X_1i^⊤,A_1i,X_2i^⊤)^⊤∈𝒮⊂^p_1+p_2+1 and 
 W^0_i=X_1i∈𝒳_1 ⊂^p_1. 
 
We will also have need to define the variables S_i and W_i, i=1,…,N; respectively,
each represents some finite dimensional function of the variables in S^0_i and  W^0_i.
We note that knowledge of S^0_i and W^0_i respectively implies knowledge of S_i and W_i;
however, the reverse may not hold.

The observed outcome Y ∈ (measured after A_2) is assumed continuous, with a larger value of Y indicating  a better clinical outcome. 


The developments below assume that the original sample, with elements independently and identically distributed
as P_0, has been split into two independent samples, say
D_I_n and D_I^c_n, being respectively of sizes n = O(N) and N-n = O(N). 
The nuisance parameters
μ̂_2Y(·), μ̂_2A(·), μ̂_1Y(·), and μ̂_1A(·)
are estimated using the data in D_I^c_n;
the finite dimensional parameters of interest are then estimated using the data D_I_n,
treating μ̂_2Y(·), μ̂_2A(·), μ̂_1Y(·), and μ̂_1A(·)
as if they were known functions. 

As developed here, our use of sample-splitting is a simple form of cross-fitting and can be generalized easily 
to make better use of the full sample

<cit.>; the simpler form used here suffices to establish

the main ideas of the proofs. Lemmas <ref> and <ref> play
an important role in several of the proofs; since n = O(N), statements
of the form o_p(N^-a) and o_p(n^-1) are equivalent, we 
use the latter to emphasize that the technical arguments rely on 
sample splitting, where a sample D_I_n of
size n is used to estimate the finite dimensional parameters
of interest.


To simplify notation, where needed all calculations implicitly condition
on the set of selected indices I_n. Using notation from the main paper, let 
Δ_2i = Δ_2(S^0_i), Δ_1i = Δ_1(W^0_i), 
μ_2Ai = μ_2A(S^0_i), μ̂_2Ai = μ̂_2A(S^0_i), μ_1Ai = μ_1A(W^0_i), 
and μ̂_1Ai = μ̂_1A(W^0_i).  In addition, as in the main paper, we define the matrices

    _2n = 1/n∑_i ∈I_n   (A_2i- μ_2Ai)^2  S^⊗ 2_i  
      _2n = 1/n∑_i ∈I_n   (A_2i- μ̂_2Ai)^2  S^⊗ 2_i


    _1n = 1/n∑_i ∈I_n   (A_1i- μ̂_1Ai)^2  W^⊗ 2_i  
      _1n = 1/n∑_i ∈I_n   (A_1i- μ̂_1Ai)^2  W^⊗ 2_i

where x^⊗ 2 = xx^⊤ for any vector x.


We make the following assumptions.

 
(i)   The support of W^0 and the conditional treatment effect Δ_1(W^0) are uniformly bounded; 
(ii)  the support of S^0 and the conditional treatment effect Δ_2(S^0) are uniformly bounded; and,
the supports of S and W are uniformly bounded.


 
  (i) μ̂_1A(W^0; D_I^c_n) - μ_1A(W^0) ^2_P_0,2 = o_p(n^-1/2); 
  (ii) μ̂_2A(S^0; D_I^c_n) - μ_2A(S^0) ^2_P_0,2 = o_p(n^-1/2).


  
  (i) μ̂_1Y(W^0; D_I^c_n) - μ_1Y(W^0) ^2_P_0,2 = o_p(1); 
  (ii) μ̂_2Y(S^0; D_I^c_n) - μ_2Y(S^0) ^2_P_0,2 = o_p(1).


   
(i) μ̂_1Y(W^0; D_I^c_n) - μ_1Y(W^0) _P_0,2μ̂_1A(W^0; D_I^c_n) - μ_1A(W^0) _P_0,2 = o_p(n^-1/2);

(ii) μ̂_2Y(S^0; D_I^c_n) - μ_2Y(S^0) _P_0,2μ̂_2A(S^0; D_I^c_n) - μ_2A(S^0) _P_0,2 = o_p(n^-1/2)


 

There exists 1 ≤ n_0 < ∞ such that _jn and _jn, j = 1,2 are positive definite for n ≥ n_0.


 
 
P( | S_1^⊤β^*_2 | = 0 ) = 0.


We prove this Theorem with help from the following lemma. 
 
Suppose Assumptions <ref>, <ref>,  and
<ref> hold. Let d_2 = dim(S) and d_1 = dim(W).
Then β̃^*_2n-β^*_2n_∞=o_p(n^-1/2) and β̃^*_1n-β^*_1n_∞=o_p(n^-1/2), where 
  
    β̃^*_2n   = _β_2∈^d_2∑_i ∈I_n{ A_2i- μ̂_2A(S^0_i)}^2 {Δ_2(S^0_i) -   S^⊤_iβ_2}^2, 
    β^*_2n   = _β_2∈^d_2∑_i ∈I_n{ A_2i- μ_2A(S^0_i)}^2 {Δ_2(S^0_i) -   S^⊤_iβ_2}^2, 
    β̃^*_1n   = _β_1∈^d_1∑_i ∈I_n{ A_1i- μ̂_1A(W_i)}^2 {Δ_1(W^0_i) -   W^⊤_iβ_1}^2, 
    β^*_1n   = _β_1∈^d_1∑_i ∈I_n{ A_1i- μ_1A(W_i)}^2 {Δ_1(W^0_i) -   W^⊤_iβ_1}^2.

 





Below, we will prove that β̃^*_2n-β^*_2n_∞=o_p(n^-1/2); the result that
β̃^*_1n-β^*_1n_∞=o_p(n^-1/2) follows from essentially identical arguments. 

Using the definitions of β̃^*_2n and β^*_2n and assuming n is large enough so that
Assumption <ref> holds, straightforward algebra shows 

    β̃^*_2n-β^*_2n =    
       (_2n^-1- _2n^-1) {1/n∑_i ∈I_n   (A_2i-  μ_2Ai)^2  S_i Δ_2i}
        +_2n^-1{1/n∑_i ∈I_n   ( μ_2Ai- μ̂_2Ai)^2  S_i Δ_2i}.

Taking norms and using the triangle inequality, it can be shown that

    √(n)β̃^*_2n-β^*_2n_∞≤√(n)_2n^-1- _2n^-1_∞  (A_n + B_n) + √(n)_2n^-1_∞  B_n

where

    A_n     = 1/n∑_i ∈I_n   (A_2i-  μ_2Ai)^2  S_i Δ_2i_∞
    
    B_n     =   1/n∑_i ∈I_n   ( μ̂_2Ai- μ_2Ai)^2 S_i Δ_2i_∞.



Suppose that _2n - _2n_∞ = o_p( n^-1/2).  Then, by Lemma <ref>
and Assumption <ref>, we have for n sufficiently large that

    _2n^-1 - _2n^-1_∞≤ K_1 _2n - _2n_∞

for any constant K_1 such that 2  _2n^-1^2_∞≤ K_1.
It can be seen that 

    A_n p→[  {A_2 - μ_2A(S^0)}^2 SΔ_2(S^0) ] _∞
    = [   SΔ_2(S^0) ( A_2 | S^0) ] _∞ < ∞.

Since cross-fitting is used to estimate 
μ̂_2A(·), we also see that B_n is an example of (<ref>);
hence, using Lemma <ref> and Assumption <ref>, 
it follows that B_n = o_p( n^-1/2). It follows from these results and (<ref>) 
that (<ref>) is o_p(1). 

In order to prove that _2n - _2n_∞ = o_p(n^-1/2), we begin by writing
_2n - _2n_∞≤ C_n + 2 D_n, where

    C_n     =  1/n∑_i ∈I_n   ( μ_2Ai-  μ̂_2Ai)^2  S^⊗ 2_i_∞
    
    D_n     =  1/n∑_i ∈I_n (A_2i-   μ_2Ai)(  μ̂_2Ai - μ_2Ai)  S^⊗ 2_i_∞.

Again, because  cross-fitting is used to estimate 
μ̂_2A(·),  we can see that C_n is also an example of (<ref>)
and it follows by previously stated arguments that C_n = o_p( n^-1/2). 
In order to establish the behavior of D_n, we first note that
D_n = max_j=1,…,d_2∑_k=1^d_2| H_njk|
where d_2 is finite and

    H_njk =  1/n∑_i ∈I_n (A_2i-   μ_2Ai)(  μ̂_2Ai - μ_2Ai)  S_ik S_ij.

It suffices to establish the behavior of H_njk. First, using the definition
of μ_2Ai = μ_2A(S^0_i) = ( A_2i | S^0_i) and the fact that
μ̂_2Ai = μ̂_2A(S^0_i) where μ̂_2A(·) is estimated
from data D_I^c_n that is independent of S^0_i ∈D_I_n for each i,
it is easy to see that ( H_njk | S^0_1, …,  S^0_n, D_I^c_n) = 0 and hence that
( H_njk | D_I^c_n) = 0. Using these same properties, it is also easily
shown that

    ( H_njk | D_I^c_n) = 1/n^2∑_i ∈I_n{ 
     (  μ̂_2Ai - μ_2Ai)^2  
     ( A_2i | S^0_i)
     (S_ik S_ij)^2 | D_I^c_n}.

Under Assumption <ref>, we can find a constant K_2 < ∞ such that

    ( H_njk | D_I^c_n)    ≤   K_2/n[ 1/n∑_i ∈I_n{μ̂_2A(S^0_i) - μ_2A(S^0_i) }^2  | D_I^c_n] 
        =    K_2/n{μ̂_2A(S^0; D_I^c_n) - μ_2A(S^0) ^2_ℙ_n,2| D_I^c_n}
        =    K_2/nμ̂_2A(S^0; D_I^c_n) - μ_2A(S^0) ^2_P_0,2,

By Chebyshev's inequality, for all ϵ > 0 we then have 

    P( n^1/2| H_njk|  > ϵ| D_I^c_n) ≤K_2/ϵ^2μ̂_2A(S^0; D_I^c_n) - μ_2A(S^0) ^2_P_0,2,

where the right-hand side is o_p(n^-1/2) by Assumption <ref>.
Lemma <ref> now implies that
H_njk = o_p(n^-1/2) and hence that D_n = o_p(n^-1/2).  Therefore, 
_2n - _2n_∞≤ C_n + 2 D_n = o_p(n^-1/2), proving the desired result.
A similar argument shows √(n)β̃^*_1n-β^*_1n_∞ = o_p(1).


As in the main paper, we define

    β^*_2   =   _β_2∈^d_2[ 
    ( A_2 | S^0 ) {Δ_2(S^0) -   S^⊤β_2}^2 ], 
    _2     =    {( A_2 | S^0 ) S^⊗ 2}
    β^*_1   =    _β_1∈^d_1[ 
    ( A_1 | W^0 ) {Δ_1(W^0) -   W^⊤β_1}^2 ], 
    _1     =    {( A_1 | W^0 ) W^⊗ 2}.

As n →∞, it is easy to see that β^*_2n converges to 

    β^*_2 = _2^-1{( A_2 | S^0 ) SΔ_2(S^0) };

similarly, β^*_1n converges to

    β^*_1 = _1^-1{( A_1 | W^0 ) WΔ_1(W^0) }.


With these preliminaries in place, we can now prove the main result.



We desire to show that β̂_2n is an asymptotically linear estimator of β_2^* with the
claimed influence function, where β^*_2 is defined in (<ref>) and

    β̂_2n   = _β_2∈^d_2∑_i ∈I_n[ Y_i-μ̂_2Y(S^0_i) - {A_2i-μ̂_2A(S^0_i)}·S_i^⊤β_2    ]^2.

In view of Lemma <ref>, we can proceed by establishing asymptotically linear representations for both
√(n)( β̂_2n - β̃^*_2n) and
√(n)(β^*_2n -  β^*_2); combined, these will lead to
that for √(n)( β̂_2n -  β^*_2).

Recalling notation introduced earlier, it is easy to show that

    √(n) ( β̂_2n    - β̃^*_2n ) = 
       _2n^-1[
    1/√(n)∑_i ∈I_n   
    (A_2i-  μ̂_2Ai)  S_i 
    { Y_i-μ̂_2Yi -(A_2i-  μ̂_2Ai)  Δ_2i}]


By adding and subtracting terms and using the model assumptions on Y_i, one can write

    Y_i-μ̂_2Yi
       -(A_2i-  μ̂_2Ai)  Δ_2i = ϵ_2i + (μ_2Yi - μ̂_2Yi) + (μ̂_2Ai - μ_2Ai) Δ_2i.

The decomposition (<ref>) implies that the term in the square brackets on the right-hand side of (<ref>)
can be decomposed into six terms:

    1/√(n)∑_i ∈I_n   
    (A_2i-   μ̂_2Ai)      S_i
    { Y_i-μ̂_2Yi -(A_2i-   μ_2Ai)  Δ_2i}
        =  1/√(n)∑_i ∈I_nϵ_2i (A_2i-  μ_2Ai)  S_i  
        -  1/√(n)∑_i ∈I_n   (A_2i-  μ_2Ai)   (μ̂_2Yi - μ_2Yi) S_i 
        +  1/√(n)∑_i ∈I_n   (A_2i-  μ_2Ai)   (μ̂_2Ai - μ_2Ai) Δ_2iS_i 
        -  1/√(n)∑_i ∈I_nϵ_2i (μ̂_2Ai -  μ_2Ai)    S_i
        +  1/√(n)∑_i ∈I_n   (μ̂_2Ai -  μ_2Ai)  (μ̂_2Yi - μ_2Yi)  S_i  
        -  1/√(n)∑_i ∈I_n   (μ̂_2Ai -  μ̂_2Ai)^2  Δ_2iS_i .

Under the assumptions of this theorem, the central limit theorem establishes the asymptotic normality 
of (<ref>), which is O_p(1). The terms (<ref>)-(<ref>) are each seen
to be examples to which Lemma <ref> applies; under Assumptions
<ref>-<ref>, it follows that each term is o_p(1). 
The terms (<ref>) and  (<ref>) are both seen
to be examples to which Lemma <ref> applies; again, under Assumptions
<ref>-<ref>, each term is o_p(1). Because
_2n^-1 - ^-1_2n_∞ = o_p(1) it follows that
√(n) ( β̂_2n  - β̃^*_2n ) can be written

    √(n) ( β̂_2n  - β̃^*_2n ) = 
     1/√(n)∑_i ∈I_n^-1_2nϵ_2i (A_2i-  μ_2Ai)  S_i  + o_p(1).


Turning to √(n)(β^*_2n -  β^*_2), we can write

    √(n)(β^*_2n -  β^*_2)       =  
    √(n)[ { n^-1∑_i ∈I_n^-1_2n  (A_2i-  μ_2Ai)^2 Δ_2iS_i }  - β^*_2] 
        = 
    1/√(n)^-1_2n∑_i ∈I_n   (A_2i-  μ_2Ai)^2  S_i ( Δ_2i  - S^⊤_i β^*_2).

Hence, using (<ref>) and (<ref>) and collecting terms,

    √(n) ( β̂_2n  -  β^*_2 )  = 
    1/√(n)^-1_2n∑_i ∈I_n  
     (A_2i-  μ_2Ai)  S_i  H_2i
    + o_p(1)

where H_2i =  ϵ_2i + 
 (A_2i-  μ_2Ai)  ( Δ_2i  - S^⊤_i β^*_2).
Using the fact that 
Y_i - μ_2Yi = ϵ_2i + (A_2i-  μ_2Ai)  Δ_2i,
we may write

    H_2i    =     Y_i - μ_2Yi - (A_2i-  μ_2Ai)  Δ_2i + 
     (A_2i-  μ_2Ai)  ( Δ_2i  - S^⊤_i β^*_2) 
        =    Y_i - μ_2Yi - (A_2i-  μ_2Ai)  S^⊤_i β^*_2.

Consequently,

    √(n) ( β̂_2n  - β^*_2 ) 
      =    1/√(n)∑_i ∈I_n_2in +o_p(1)

where

    _2in  =
     ^-1_2n (A_2i-  μ_2Ai)  S_i  
     {
     Y_i - μ_2Yi - (A_2i-  μ_2Ai)  S^⊤_i β^*_2}.

Since
_2n^-1 - ^-1_2_∞ = o_p(1),
it now follows that

    √(n) ( β̂_2n  - β^*_2 ) 
      =    1/√(n)∑_i ∈I_n_2i +o_p(1),

where

    _2i  =
     ^-1_2 
     (A_2i-  μ_2Ai)  S_i  
     {
     Y_i - μ_2Yi - (A_2i-  μ_2Ai)  S^⊤_i β^*_2}

has mean zero and variance ^-1_2 Q_2 ^-1_2 where
_2 =  ( J_2^⊗ 2 ) and

    J_2 = 
    {A_2-  μ_2A(S^0)}S[
    Y - μ_2Y(S^0)  - { A_2-  μ_2A(S^0) }S^⊤β^*_2].





As in part (a), we need to show that β̂_1n is an asymptotically linear estimator of β_1^* with a 
certain influence function, where β^*_1 is defined in (<ref>) and

    β̂_1n   = _β_1∈^d_1∑_i ∈I_n[  - μ̂_1Y(W^0_i) - {A_1i-μ̂_1A(W^0_i)}·W_i^⊤β_1    ]^2,

where  is calculated as

    = Y_i + S_i^⊤β̂_2n{ I(S_i^⊤β̂_2n  > 0 ) - A_2i}.



Proceeding similarly to the proof of part (a), we will establish asymptotically linear representations for both
√(n)( β̂_1n - β̃^*_1n) and
√(n)(β^*_1n -  β^*_1); combining these will
provide the claimed influence function for √(n)( β̂_1n -  β^*_1).

We begin with √(n)( β̂_1n - β̃^*_1n).
Define ϵ_1i =  - ( |W^0_i, A_1i),
where  is given by

    = Y_i + S_i^⊤β^*_2{ I(S_i^⊤β^*_2  > 0 ) - A_2i};

by construction, 
ϵ_1i  =  - μ_1Ai - (A_1i - μ_1Ai) Δ_1i
and
(ϵ_1i|W^0_i, A_1i)=0. 

In addition, let δ̂_i =  - .  Similarly to the proof
in part (a), we can decompose √(n)( β̂_1n - β̃^*_1n) into the
sum of several terms:




    √(n)( β̂_1n -    β̃^*_1n) =  
       _1n^-1{1/√(n)∑_i ∈I_n (A_1i -
      μ_1Ai) (  ϵ_1i  + δ̂_i   )    W_i  }
        - _1n^-1{1/√(n)∑_i ∈I_n (A_1i - μ_1Ai)  (μ̂_1Yi- μ_1Yi) W_i }
         + _1n^-1{1/√(n)∑_i ∈I_n (A_1i -  μ_1Ai) (μ̂_1Ai- μ_1Ai) Δ_1iW_i  }
        - _1n^-1{1/√(n)∑_i ∈I_n ( μ̂_1Ai- μ_1Ai)   ϵ_1iW_i }
        - _1n^-1{1/√(n)∑_i ∈I_n ( μ̂_1Ai- μ_1Ai)  δ̂_i     W_i }
        -   _1n^-1{1/√(n)∑_i ∈I_n ( μ̂_1Ai- μ_1Ai)  (μ̂_1Yi- μ_1Yi)    W_i  }
        +  _1n^-1{1/√(n)∑_i ∈I_n (μ̂_1Ai- μ_1Ai)^2  Δ_1iW_i
          }


Assuming that μ̂_1A(·) and μ̂_1Y(·) are estimated similarly to μ̂_2A(·) and μ̂_2Y(·)
(i.e., meaning, sample splitting has been used) and in view of the fact that _1n is a consistent estimator of
_1, Lemma  <ref> implies that the terms (<ref>), (<ref>), and (<ref>) are all o_p(1)
under Assumptions <ref> – <ref>; similarly,
Lemma  <ref> implies that  the terms (<ref>) and (<ref>)  are also o_p(1)
under these same assumptions. To complete this part of the proof, we must 
therefore establish the asymptotic behavior of (<ref>) and (<ref>), both of
which depend on the asymptotic behavior of δ̂_i =  - . 
The terms (<ref>) and (<ref>) isolate the potential for non-regular behavior; however,
as we will see, Assumption <ref> is only needed for controlling such behavior  in (<ref>).

To determine the asymptotic behavior of (<ref>), let D_ i = I(S_i^⊤β^*_2  > 0 ) - A_2i and

    R̂_ni = I( S_i^⊤β̂_2n  > 0) - I(S_i^⊤β^*_2  > 0)

for i ∈I_n.

Algebra now shows

    1/√(n)∑_i ∈I_n  (A_1i - μ_1Ai)  W_i      (  ϵ_1i  + δ̂_i   )   = 
       1/√(n)∑_i ∈I_n  (A_1i - μ_1Ai)  W_i  {ϵ_1i  + S_i^⊤  (β̂_2n-β^*_2)  D_i  }
         + 1/√(n)∑_i ∈I_n  (A_1i - μ_1Ai)  W_i  S_i^⊤β^*_2R̂_ni
         + 1/√(n)∑_i ∈I_n  (A_1i - μ_1Ai)  W_i  S_i^⊤  ( β̂_2n-β^*_2)  R̂_ni.

Although (<ref>) and (<ref>) can be easily combined, treating these two terms separately turns out to be advantageous.
We first consider (<ref>). Note that R̂_ni∈{-1,0,1} and, importantly, that | R̂_ni | ≤ R_ni, where
R_ni = I{ 0 ≤ | S_i^⊤β^*_2 | ≤ | S_i^⊤  ( β̂_2n-β^*_2) | }.
It follows that

    | S_i^⊤β^*_2 | R_ni≤ | S_i^⊤  ( β̂_2n-β^*_2) | R_ni,

an inequality that is trivially true when R_ni = 0 and true for R_ni = 1 in view of its definition.
Consequently, considering the j^th element of the vector W_i, we have

    | ∑_i ∈I_n  (A_1i - μ_1Ai)  W_ijS_i^⊤β^*_2R̂_ni|
        ≤   ∑_i ∈I_n|  (A_1i - μ_1Ai) W_ij|   |  S_i^⊤β^*_2|   |  R̂_ni| 
       ≤   ∑_i ∈I_n|  (A_1i - μ_1Ai)  W_ij|   |  S_i^⊤β^*_2|   R_ni
       ≤   ∑_i ∈I_n|  A_1i - μ_1Ai| |  W_ij|   |  S_i^⊤  ( β̂_2n-β^*_2)   |   R_ni
       ≤    4 C
     ∑_i ∈I_n|  S_i^⊤  ( β̂_2n-β^*_2)   |   R_ni,

the last step following from the fact that A_1i is binary, μ̂_1A(·) ∈ [0,1], and |W_ij| is bounded, say, by 
a finite constant C. 

Considering (<ref>), a similar calculation shows that

    |
     ∑_i ∈I_n  (A_1i - μ_1Ai)  W_i  S_i^⊤  ( β̂_2n-β^*_2)  R̂_ni|
        ≤    
     4 C
     ∑_i ∈I_n|  S_i^⊤  ( β̂_2n-β^*_2)   |   R_ni.

Therefore, 

    ≤    
    8 C 1/√(n)∑_i ∈I_n|  S_i^⊤  ( β̂_2n-β^*_2)   |   R_ni
       ≤    8 C √(n)[ 1/n∑_i ∈I_n{S_i^⊤  ( β̂_2n-β^*_2)   }^2 ]^1/2[ 1/n∑_i ∈I_n   R_ni^2 ]^1/2
        =     8 C λ_1n√(n)β̂_2n-β^*_2_2  R̅_n^1/2,

where ·_2 is the usual Euclidean vector norm and λ_1n > 0 is the square root of the maximum eigenvalue of 
n^-1∑_i ∈I_nS^⊗ 2_i.
Because √(n)β̂_2n-β^*_2_2 = O_p(1) and λ_1n converges to a finite constant as
n →∞ under our assumptions, it follows that (<ref>) + (<ref>) is o_p(1) if R̅_n = o_p(1). However,
Markov's inequality implies that

    P( R̅_n > α|S^0_1,…, S^0_n) ≤ (n α)^-1∑_i ∈I_n(R_ni|S^0_1,…, S^0_n
    )

for any α > 0, where
(R_ni|S^0_1,…, S^0_n)  =  
P{ R_ni = 1 |S^0_1,…, S^0_n }.
Letting

    I^*_i = I( | S_i^⊤β^*_2 | = 0 ),

an easy conditioning argument shows 
(R_ni|S^0_1,…, S^0_n ) = 
I^*_i + (1-I^*_i) P_ni where

    P_ni =
    P{  | S_i^⊤  ( β̂_2n-β^*_2) | ≥ k_i |  (S^0_1,…, S^0_n),   k_i > 0}

for k_i = | S_i^⊤β^*_2 |.
Letting n →∞, the fact that β̂_2np→β^*_2
as n →∞ implies P_ni→ 0 for each i; hence,

    lim_n →∞
    P( R̅_n > α|S^0_1,…, S^0_n) ≤1/n αlim_n →∞∑_i ∈I_n  I^*_i.

However, under our assumptions, 

    lim_n →∞1/n∑_i ∈I_n  I^*_i p→ P( | S_1^⊤β^*_2 | = 0 )

and it follows from Assumption <ref> that  R̅_n p→ 0.

Because _1n^-1 - ^-1_1n_∞ = o_p(1), we have now shown that

    = _1n^-1[ 1/√(n)∑_i ∈I_n  (A_1i - μ_1Ai)  W_i  {ϵ_1i  + S_i^⊤  (β̂_2n-β^*_2)  D_i  } + o_p(1) ]

under Assumptions <ref> – <ref>, where we recall the notation D_ i = I(S_i^⊤β^*_2  > 0 ) - A_2i.

To establish (<ref>), observe that we may similarly decompose it as above, leading to

    1/√(n)∑_i ∈I_n  (  μ̂_1Ai    - μ_1Ai)  W_i    δ̂_i      = 
       1/√(n)∑_i ∈I_n  ( μ̂_1Ai - μ_1Ai) W_i  ϵ_1i
         + 1/√(n)∑_i ∈I_n  ( μ̂_1Ai - μ_1Ai)  W_i  S_i^⊤β^*_2R̂_ni
         + 1/√(n)∑_i ∈I_n ( μ̂_1Ai - μ_1Ai)  W_i  S_i^⊤   ( β̂_2n-β^*_2)  R̂_ni
         + 1/√(n)∑_i ∈I_n  ( μ̂_1Ai - μ_1Ai) W_i  S_i^⊤  ( β̂_2n-β^*_2)  D_i.

The term (<ref>) can be handled using Lemma  <ref>. The remaining terms can be handled similarly to (<ref>) and (<ref>);
however, the required decomposition of terms differs some and, importantly, can make use of Assumption <ref>.
In particular, establishing the behavior of (<ref>)-(<ref>) can be done under 
Assumptions <ref> – <ref>, without additionally imposing Assumption <ref>, 
showing that any effect of non-regularity is limited to the behavior of (<ref>) and (<ref>), or equivalently,

    1/√(n)∑_i ∈I_n  (A_1i - μ_1Ai)  W_i  S_i^⊤β̂_2nR̂_ni.





The above proof establishes an asymptotic linear representation for
√(n)( β̂_1n - β̃^*_1n).

Turning to √(n)(β^*_1n -  β^*_1), we can write

    √(n)(β^*_1n -  β^*_1)   
        = 
    1/√(n)^-1_1n∑_i ∈I_n   (A_1i-  μ_1Ai)^2  W_i ( Δ_1i  - W^⊤_i β^*_1).

Hence, using (<ref>) and (<ref>) and collecting terms, it follows that

    √(n) ( β̂_1n  -  β^*_1 ) = 
    1/√(n)^-1_1n∑_i ∈I_n  
     (A_1i-  μ_1Ai)  W_i  H_1i
    + o_p(1)

where 
H_1i =  ϵ_1i + 
 (A_1i-  μ_1Ai)  ( Δ_1i  - W^⊤_i β^*_1)
 + S^⊤_i ( β̂_2n - β^*_2 ) D_i.
 Because 
 
    ϵ_1i  =  - μ_1Ai - (A_1i - μ_1Ai) Δ_1i

we have

    H_1i =   - μ_1Ai - (A_1i - μ_1Ai) W^⊤_i β^*_1
     + S^⊤_i ( β̂_2n - β^*_2 ) 
     { I(S_i^⊤β^*_2  > 0 ) - A_2i}.

 Using the fact that _1n^-1 - ^-1_1_∞ = o_p(1), it follows that

    √(n) ( β̂_1n  -  β^*_1 ) = 
    1/√(n)^-1_1∑_i ∈I_n  
     (A_1i-  μ_1Ai)  W_i  H_1i
    + o_p(1).

Defining

    K_n = 1/n∑_i ∈I_n  
     (A_1i-  μ_1Ai) { I(S_i^⊤β^*_2  > 0 ) - A_2i}W_i  S^⊤_i

and letting K denote its limit in probability, the results from part (a), in particular (<ref>), 
now imply that

    √(n) ( β̂_1n  -  β^*_1 ) = 
    1/√(n)∑_i ∈I_n_1i + o_p(1).

for

    _1i = ^-1_1[  (A_1i-  μ_1Ai)  W_i  
     { - μ_1Ai - (A_1i - μ_1Ai) W^⊤_i β^*_1}
     + K _2i].

This representation result implies 
√(n) ( β̂_1n  - β^*_1 ) 
  d→ N( 0,   ^-1_1 _1 ^-1_1 )
where we define the matrix _1 = { ( J_1 + K^-1_2 J_2) ^⊗ 2}, 
J_2 is given in (<ref>), 

    J_1 = {A_1-  μ_1A(W^0)}W[
    Ỹ- μ_1A(W^0) - { A_1 - μ_1A(W^0) }W^⊤β^*_1]

and

    K = [ { A_1-  μ_1A(W^0) }{ I(S^⊤β^*_2  > 0 ) - A_2}WS^⊤].


 
 


As established in the proof of Theorem 1, the regularity Assumption <ref> is imposed
only to control the potentially non-regular behavior of the terms (<ref>) and (<ref>). The origin
of this non-regular behavior is the dependence of each term on

    R̂_ni = I( S_i^⊤β̂_2n  > 0) - I(S_i^⊤β^*_2  > 0), i ∈I_n.

In view of the proof of Theorem 1, establishing that each of  (<ref>) and (<ref>) is o_p(1) is sufficient
to prove the corollary as stated.

To simply the proof of these results, let β̂^(-i)_2n be the least squares estimator based
on the subset of subjects that excludes subject i, and define

    R̂̃̂_ni = I( S_i^⊤β̂^(-i)_2n  > 0) - I(S_i^⊤β^*_2  > 0), i ∈I_n

and also

    R̃_ni = I{ 0 ≤ | S_i^⊤β^*_2 | ≤ | S_i^⊤  ( β̂^(-i)_2n-β^*_2) | };

similarly to before, | R̂̃̂_ni | ≤R̃_ni.

Then, considering (<ref>) with  R̂_ni replaced by R̂̃̂_ni, we may write

    E[ ]
       =  1/√(n)∑_i ∈I_n  E(  W_i  S_i^⊤β^*_2R̂̃̂_ni  E[ A_1i - μ_1Ai|W_i , S_i,   I_n, β̂^(-i)_2n] |I_n ).

In view of the definition of β̂^(-i)_2n, we have

    E[ A_1i - μ_1Ai|W_i , S_i,   I_n, β̂^(-i)_2n]
    = E[ A_1i - μ_1Ai|W_i , S_i,   I_n ] = 0,

the last equality following by assumption. Therefore, it follows that 
E[ |I_n ]= 0.
Arguing similarly and using the conditional variance formula,

    var[ |I_n ] = 1/n∑_i ∈I_n 
     E (  W^⊗ 2_i  ( S_i^⊤β^*_2)^2 [ A_1i - μ_1Ai]^2 R̂̃̂^2_ni|I_n ).

Let V_kj denote the (k,j) element of the matrix on the right-hand side of this last expression. Then, 
under Assumptions <ref> – <ref>, it can be shown that 

    V_kj≤C/n∑_i ∈I_n 
     E ( ( S_i^⊤β^*_2 )^2  R̃_ni|I_n )

for some finite constant C > 0. Similarly to the proof of Theorem 1, the fact that

    | S_i^⊤β^*_2 | R̃_ni≤ | S_i^⊤  ( β̂^(-i)_2n-β^*_2) | R̃_ni

now implies the inequality

    V_kj≤C/n∑_i ∈I_n 
     E (  [ S_i^⊤  ( β̂^(-i)_2n-β^*_2) ]^2   R̃_ni|I_n ),

from which it follows that

    V_kj   ≤   C/n∑_i ∈I_n   E (    (β̂^(-i)_2n-β^*_2)^⊤S^⊗ 2_i   ( β̂^(-i)_2n-β^*_2)    R̃_ni|I_n ),
       ≤    
     C λ^2_1n 
      E (    β̂^(-i)_2n-β^*_2^2_2 |I_n ),

where  λ^2_1n > 0 is the maximum eigenvalue of n^-1∑_i ∈I_nS^⊗ 2_i.
Because the right-hand side of (<ref>) goes to zero as n →∞ for all (k,j), it follows that (<ref>) is o_p(1). A similar argument
establishes that (<ref>) is o_p(1). 

 
 
 
 
 


§ ADDITIONAL SIMULATION STUDIES



 §.§ Various sample sizes
 

In this section, we complement the main simulation study by examining the performance of our proposed method under various sample sizes. 
We generate 500 datasets with sample sizes N of 1000, 500, and 250 using the same generative model as in Section 5.1 in the paper,
and except in Section 2.2, use the same methods to estimate μ_1Y(·), μ_2Y(·), μ_1A(·), and μ_2A(·)
as in the main simulation study.

Tables <ref>-<ref> in this document respectively summarize the results for sample sizes N of 1000, 500, and
250 sample sizes. Overall, our method outperforms both Q_N,N and dWOLS_N,N for all sample sizes, particularly when the underlying 
treatment assignment and the outcome models are nonlinear. More specifically,  when the postulated parametric models for the nuisance 
parameters are correctly specified, the bias of the dWOLS_N,N estimators are comparable to the estimators obtained by our proposed method,
but the latter usually has a substantially smaller standard error. However, when the parametric models for the outcome and the treatment assignment models 
are both misspecified, the dWOLS_N,N estimators exhibit large biases. 

When N=250, under the linear treatment assignment model and FGS^R outcome model, the proposed method shows larger bias in estimating 
β_2 than the dWOLS_N,N (Table <ref>). We conjecture that this occurs because there are only a relatively small number of 
units that are rerandomized at stage 2 (i.e., 50%), and that this subsequently affects the performance of Super Learner. 











































 §.§ Performance under alternative nonparametric estimation methods


In this section, we assess the performance of the proposed method when the nuisance parameters are instead estimated using either random forests or generalized additive models, recalling that both were included as part of the library used by Super Learner. As in the main simulation study, we generated 500 datasets of size N=2000 using the same generative model as in Section <ref>. In Table <ref> and <ref>, columns RF-RF and GAM-GAM respectively represent cases where a random forest and a generalized additive model are used for both the treatment assignment (i.e., μ_1A(·) and μ_2A(·)) and the outcome models (i.e., μ_1Y(·) and μ_2Y(·)). The column RF-GAM represents a case where a random forest is used for the outcome model and a generalized additive model is used for the treatment assignment model.  Overall, as the complexity of the outcome or treatment assignment model increases, modeling the nuisance parameters using a random forest results in estimators that are less biased compared with those obtained by the generalized additive model. However, comparing these results with those presented in Tables 1 and 2, we observe that Super Learner typically performs better than these other data adaptive methods.















 §.§ Additional simulation results for the non-regular setting


In Section 5.2 of the main paper, we showed that the proposed method can provide valid inferences in non-regular settings 
when the assumptions of Corollary 1 are satisfied. Table <ref> supplements the results of Table 3 for the main simulation study 
by providing additional comparisons with Q_N,N and dWOLS_N,N. 








































Below, we also consider a simulation that corresponds to a non-regular setting when the assumptions of Corollary 1 and
Assumption 6 are both violated. We consider the same general set-up as in Section 5.2, but modify the outcome models as described below:


  * Linear^NR,ϖ: Y= _1^⊤_1+ _2^⊤_2 +A_1_1^⊤_1+ A_2 R (θ_21 A_1+ θ_22X̆_21 )+ϵ where _1=(1,0.1,0.1,0.1,0.1)^⊤, _2=(1,0.1,0.1,0.1)^⊤, _1=(0,0,0,0,0)^⊤, and _2 = 2 (ϖ,ϖ)^⊤;  

  * Non-linear^NR,ϖ: Y= f(_1)+A_1_1^⊤_1+   A_2 R (θ_21 A_1+ θ_22X̆_21)+ϵ where _1=(0,0,0,0,0)^⊤,
_2 = 2 (ϖ,ϖ)^⊤,  and
for x = (x_1,x_2,x_3,x_4,x_5)^⊤, we set

    f(x) = sin(π x_1 x_2) + 2
      (x_3-0.5)^2 +  x_4 + 1.5 x_1/|x_2|+|x_3| +2 x_1 (x_2+x_3).


Here,
X̆_2=(X̆_21,X̆_22,X_23,X_24)^⊤ where X̆_21 is generated from a Bernoulli distribution with 
success probability [ 1+exp{-(2X_11-2A_1-1)} ]^-1,  
X_22 = U_11, X_23 = U_12,
X_24 = 0.35 X_15 + U_13,
and X_25 = U_14,
where U_1l, l =1,…,3 are independent and uniformly
distributed on [0,1]; the noise variable ϵ is again generated from N(0,σ=0.5). 



In the above generative models, S^0 = (X^⊤_1,A_1,X̆^⊤_2)^⊤ and W^0 = X_1;
the vectors S=R (1,A_1, X̆_21, X̆_22)^⊤ and W = (1, X_11, X_12)^⊤ are respectively used 
to model the second and first stage Q-functions. 

The second stage blip functions are linear and correctly specified in both models, leading
to β_2^* = (0,^⊤_2,0)^⊤.

In both scenarios, for each subject i, we define the first-stage pseudo outcome as 

    Ỹ_i
    =Y_i + ϖ  R S_i^⊤β^*_2 {I(S_i^⊤β^*_2  >0) -A_2i},

its estimate Ỹ̂_i being calculated by substituting β̂_2n in for β^*_2.
The construction of the pseudo-outcome, specifically the projection S^⊤β^*_2, results in 
the violation of Assumption <ref> for both generative models.

In particular,  ϖ=0 corresponds to no second-stage effect modifier because
β^*_2 = 0; hence, we have P( |S^⊤^*_2| = 0) = 1.
Setting ϖ=1 instead implies that there is no second-stage treatment effect 
when R (A_1+ X̆_21) = 0 and a reasonably strong effect
when R (A_1+ X̆_21) = 1; in this case, we have β^*_2  > 0
and 0 < P( |S^⊤^*_2| = 0) < 1.
In addition to violating Assumption 6, both generative models also violate the key assumption of Corollary 1 
(i.e., E( A_1i - μ_1Ai|W_i, S_i,  I_n ) = 0, i ∈I_n) 
because the Stage 2 decision rule depends on  A_1. 




We compare the coverages of our estimators with those obtained for dWOLS using the m-out-of-n bootstrap (i.e., dWOLS_m,N^κ=0.05).  

Table <ref> presents the results based on 1000 datasets of size 2000. In the table, β^*_1,1, β^*_1,2 and β^*_1,3 
respectively the true parameter values corresponding to variables A_1 (i.e., we assume W has an intercept), A_1 X_11 and A_1 X_12.  

When ϖ=1, (β^*_1,1,β^*_1,2, β^*_1,3) ≈ (0.55,0.00,-0.57) and  (β^*_1,1,β^*_1,2, β^*_1,3) ≈ (0.77,0.00,-0.29) for the linear and non-linear outcome models, respectively (note: these are approximated by simulation). Also, when ϖ=0, (β^*_1,1,β^*_1,2, β^*_1,3) ≈ (-0.23,0.00,-0.29) and  (β^*_1,1,β^*_1,2, β^*_1,3) = (0,0,0) for the linear and non-linear outcome models, respectively. Despite violation
of the regularity assumption, the proposed method continues to provide valid confidence intervals for the parameters β_1,2, and β_1,3. However, the confidence intervals for β_1,1 exhibit less-than-nominal coverage in a majority of cases. As expected, when either of the outcome or the treatment assignment models are correctly specified the dWOLS_m,N^κ=0.05 provides valid, if typically conservative, confidence intervals for all parameters. However, when both of these models are misspecified, the coverage rates are more substantially compromised. 










 §.§ Value Function Estimates for Regular Case


In this section we plot the value functions corresponding to the estimated decision rules obtained
in the regular setting (see Section 5.1 of the main paper). 

All methods perform essentially identically as far as the underlying value function
when the outcome model is linear, even for N=250; these results are not included. 

In the case where the outcome model is not linear, we see that the proposed
method typically results in a value function that is closest to the true optimal value, 
followed by dWOLS and then standard Q-learning. However, the overall 
degree of discrepancy is generally small, particularly between dWOLS and
the proposed approach.


 

 



