A Chambaz

M J Van Der Laan

W Zheng

Targeted Covariate-Adjusted Response-Adaptive LASSO-Based Randomized Controlled Trials

Keywords:

We present a new covariate-adjusted response-adaptive randomized controlled trial design and inferential procedure built on top of it. The procedure is targeted in the sense that (i) the sequence of randomization schemes is group-sequentially determined by targeting a user-specified optimal randomization design based on accruing data and, (ii) our estimator of the user-specified parameter of interest, seen as the value of a functional evaluated at the true, unknown distribution of the data, is targeted toward it by following the paradigm of targeted minimum loss estimation. We focus for clarity on the case that the parameter of interest is the marginal effect of a binary treatment and that the targeted optimal design is the Neyman allocation, in an effort to produce an estimator with smaller asymptotic variance. For clarity too, we consider the case that the estimator of the conditional outcome given treatment and baseline covariates, a key element of the procedure, is obtained by LASSO regression. Under mild assumptions, the resulting sequence of randomization schemes converges to a limiting design, and the TMLE estimator is consistent and asymptotically Gaussian. Its asymptotic variance can be estimated too. Thus we can build valid confidence intervals of given asymptotic levels. A simulation study confirms our theoretical results.

Introduction 1.Overview

This technical report is devoted to the study of a so-called group-sequential CARA randomized controlled trial (RCT), with a particular focus on incorporating more flexible (i.e., data-adaptive) techniques to model the response. A CARA RCT is Covariate-Adjusted: the treatment randomization schemes are allowed to be a function of the patients' pre-treatment covariates. In addition, a CARA RCT is Response-Adaptive: the investigators have the opportunity to adjust these schemes during the course of the trial based on accruing information, including previous responses, in order to meet some pre-specified objectives. In a group-sequential CARA RCT, the latter adjustments are made at interim time points given by sequential inclusion of blocks of c patients, where c ≥ 1 is a pre-specified integer. We consider the case of c = 1 for simplicity of exposition, though the discussions generalize to any c > 1.

The trial protocol pre-specifies the observed data structure, scientific parameters of interest, analysis methods, and a criterion characterizing an optimal randomization scheme. Here, some baseline covariates and a primary outcome of interest are measured on each patient. We choose the marginal treatment effect of a binary treatment as our parameter of interest, ψ 0 . It is analyzed using targeted minimum loss estimation (TMLE) on top of the so-called LASSO (least absolute shrinkage and selection operator) methodology [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] that we choose to illustrate the application of data-adaptive techniques to model the response. The TMLE methodology was first introduced by (van der [START_REF] Van Der Laan | Targeted maximum likelihood learning[END_REF] in the independent identically distributed setting. Its extension to adaptive RCTs was considered in (van der [START_REF] Van Der Laan | The construction and analysis of adaptive group sequential designs[END_REF] and [START_REF] Chambaz | Inference in targeted goup sequential covariate-adjusted randomized clinical trials[END_REF], upon which this technical report relies. The extension based on LASSO that we present here encompasses the parametric approach of [START_REF] Chambaz | Inference in targeted goup sequential covariate-adjusted randomized clinical trials[END_REF] as a special case. For concreteness, we choose the so-called Neyman design as our optimal randomization scheme. The Neyman design minimizes the Cramér-Rao lower bound on the asymptotic variances of a large class of estimators of ψ 0 . The resulting Neyman allocation probabilities are evaluated conditionally on the baseline covariates. By targeting the Neyman design, we aim at improving the efficiency of the study, i.e., at reaching a valid result using as few blocks of patients as possible. We emphasize that the results and procedures presented here are generally applicable to other parameters and optimal randomization schemes.

We show that, under mild conditions, the resulting TMLE estimator of ψ 0 is consistent and asymptotically normal regardless of the consistency of the LASSO estimator of the conditional expectation of the response given treatment and baseline covariates. Furthermore, the resulting targeted CARA design converges to a fixed limiting design, which equals the Neyman design if the LASSO estimator is consistent and if the Neyman design belongs to a user-supplied set of randomization schemes. The general framework that combines CARA RCTs with machine-learning techniques is presented in a separate article. Before we delve into the main contents, let us motivate our discussion with a bird's eye view of the landscape of CARA designs.

Literature Review

Adaptive randomization has a long history that can be traced back to the 1930s. We refer to (Rosenberger, 1996, Rosenberger, Sverdlov, and[START_REF] Rosenberger | Adaptive randomization for clinical trials[END_REF], [START_REF] Hu | The theory of response-adaptive randomization in clinical trials[END_REF], Section 1.2) and (Jennison and Turnbull, 2000, Section 17.4) for a comprehensive historical perspective. Many articles are devoted to the study of response-adaptive randomizations, which select current treatment probabilities based on responses of previous patients, but not on the covariates of the current patients. We refer to [START_REF] Hu | The theory of response-adaptive randomization in clinical trials[END_REF][START_REF] Chambaz | Targeting the optimal design in randomized clinical trials with binary outcomes and no covariate: Simulation study[END_REF][START_REF] Rosenberger | Adaptive randomization for clinical trials[END_REF] for a bibliography on that topic. In a heterogeneous population, however, it is often desirable to take into account the patients' characteristics for treatment assignment. CARA randomization tackles the issue of heterogeneity by dynamically calculating the allocation probabilities based on previous responses and current and past values of certain covariates. Compared to the broader literature on response-adaptive randomization, the advances in CARA procedures are relatively recent, but growing steadily. Among the first approaches, [START_REF] Rosenberger | Covariate-adjusted response-adaptive designs for binary response[END_REF]Agarwal, 2001, Bandyopadhyay and[START_REF] Bandyopadhyay | Adaptive designs for normal responses with prognostic factors[END_REF] considered randomization procedures defined as explicit functions of the conditional responses, which are modeled by generalized linear models. Though these procedures are not defined based on formal optimality criteria, their general goal is to allocate more patients to their corresponding "better" treatment arm. [START_REF] Atkinson | Adaptive biased-coin designs for skewing the allocation proportion in clinical trials with normal responses[END_REF] presented a biased-coin design with skewed allocation, which is determined by sequentially maximizing a function that combines the variance of the parameter estimate, based on a Gaussian linear model for the conditional response, and the conditional treatment effect given covariates. Up till here, very little work had been devoted to asymptotic properties of CARA designs. Subsequently, [START_REF] Zhang | Asymptotic properties of covariate-adjusted responseadaptive designs[END_REF], [START_REF] Zhang | A new family of covariate-adjusted response-adaptive designs and their properties[END_REF] established the efficiency theory for CARA designs converging to any given target design, when the responses follow a generalized linear model, and proposed a covariate-adjusted doubly-adaptive biased coin design whose asymptotic variance achieves the efficiency bound. [START_REF] Chang | Sequential estimation for covariate-adjusted response-adaptive designs[END_REF] proposed a sequential estimation of CARA designs under generalized linear models for the response. This procedure allocates treatment based on the patients' baseline covariates, accruing information and sequential estimates of the treatment effect and uses a stopping rule that depends on the observed Fisher information. With regard to hypothesis testing, [START_REF] Shao | A theory for testing hypotheses under covariate-adaptive randomization[END_REF], [START_REF] Shao | Validity of tests under covariate-adaptive biased coin randomization and generalized linear models[END_REF] provided asymptotic results for valid tests under generalized linear models for the responses. Most recently, progress has also been made in CARA designs in the longitudinal settings, see for example [START_REF] Biswas | On a class of optimal covariate-adjusted response-adaptive designs for survival outcomes[END_REF][START_REF] Huang | Longitudinal covariate-adjusted response-adaptive randomized designs[END_REF][START_REF] Sverdlov | Utility of covariate-adjusted response-adaptive randomization in survival trials[END_REF].

To tackle the issue of restrictive modeling assumptions, Chambaz and van der Laan (2013) proposed a TMLE analysis of a CARA design where the treatment allocation is conditional on a summary measure of the covariates that takes only finitely many values. Under such a framework, the treatment effect is defined nonparametrically, and the consistency and asymptotic normality of its estimator is robust to misspecification of the parametric working model for the response. However, assigning treatment based on such summary measures is perhaps too restrictive in real-life RCTs where response to treatment may be correlated with a large number of a patient's baseline characteristics, some of which being continuous. Moreover, although a misspecified parametric working model for the response does not hinder the consistency of the treatment effect estimator, it may affect its efficiency and the convergence of the CARA design to the targeted optimal design.

In this technical report, we generalize the results of Chambaz and van der Laan (2013) to address the two issues mentioned above. We adopt a loss-based approach to the construction of more flexible CARA randomization schemes while exploiting data-adaptive estimators for the estimation of the response model, in search for greater efficiency through better variable adjustments and more accurate estimation of the variance of the estimator.

Organization

The remainder of this technical report is organized as follows. In Section 2 we introduce our LASSObased group-sequential CARA RCT design and the TMLE procedure built on top of it to infer the marginal treatment effect of a binary treatment. Section 3 is devoted to the presentation of results pertaining to the convergence of our targeted CARA design and to the asymptotics, consistency and central limit theorem, of the TMLE estimator. A simulation study is described and its results summarized in Section 4. The technical report closes on a discussion in Section 5.

2 Targeted CARA RCT using LASSO

In the introduction, we have outlined the motivation to use data-adaptive procedures to estimate the conditional response given treatment and covariates. For concreteness of the formal theoretical development, we consider here the LASSO estimator, which is a shrinkage and selection method for generalized regression models that optimizes a loss function of the regression coefficients subject to the constraint that the L 1 norm of the coefficient vector be upper-bounded by a given value. The parametric estimators considered in [START_REF] Chambaz | Inference in targeted goup sequential covariate-adjusted randomized clinical trials[END_REF]) are a special case of a LASSO estimator.

We begin by establishing the key features of the trial, namely, the parameter of interest, analysis method, and the optimal randomization scheme. Then, we describe the data generating process (including estimation of the response model using LASSO and adaptation of the randomization scheme) and the targeted maximum likelihood estimation procedure.

Observed Data Structure, Parameter of Interest and Optimal Design

Prior to data collection, the trial protocol notably specifies the observed data structure, parameter of interest, and the optimal randomization design to target, both expressed in terms of features of the true, unknown data-generating process in the population of interest. In this technical report we consider a simple situation, judging by the definition of the data and our choice of parameter of interest. The range of application of the methods presented here extends beyond this limited yet instructive framework.

Sections 2.1.1, 2.1.2 and 2.1.3 are respectively devoted to the presentation and discussion of the observed data structure, parameter of interest, and optimal randomization design.

Observed Data Structure

The data structure O writes as O ≡ (W, A,Y ), where W ∈ W consists of the baseline covariates (some of which may be continuous), A ∈ A ≡ {0, 1} is the binary treatment of interest, and Y ∈ Y is the primary outcome of interest. We assume that the outcome space O ≡ W × A × Y is bounded. Without loss of generality, we may then assume that Y ∈ Y ≡ (0, 1) is bounded away from 0 and 1.

Every distribution of O consists of three components. On one hand, the marginal distribution of W and the conditional distribution of Y given (A,W ) form a couple which is given by nature. On the other hand, the conditional distribution of A given W , also know as (a.k.a.) randomization scheme, is controlled by the investigators of the RCT. To reflect this dichotomy, we denote P Q,g the distribution of O whose couple formed by the marginal distribution of W and the conditional distribution of Y given (A,W ) equals Q and whose randomization scheme equals g ∈ G , with G the set of all randomization schemes. For a given Q, we denote Q W the related marginal distribution of W and Q Y the related conditional expectation of Y given (A,W ). Moreover, we denote Q 0 the true couple in our population of interest, which is unknown to us, and we assume that this Q 0 does not vary during the whole duration of the RCT. Thus, for any Q and g, P Q 0 ,g is the true, partially unknown distribution of O when one relies on g, and

E P Q,g (Y |A,W ) = Q Y (A,W ), P Q,g (A = 1|W ) = g(1|W ) = 1 -g(0|W ) P Q,g -almost surely.

Parameter of Interest

The parameter of interest under consideration in this technical report is the marginal treatment effect on an additive scale:

ψ 0 ≡ E P Q 0 ,g {Q Y,0 (1,W ) -Q Y,0 (0,W )} = (Q Y,0 (1, w) -Q Y,0 (0, w))dQ W,0 (w),
which evidently depends on P Q 0 ,g only through Q 0 . Of particular interest in medical, epidemiological and social sciences research, this parameter can be interpreted causally under additional assumptions on the data-generating process [START_REF] Pearl | Causality[END_REF]. Central to our approach is seeing ψ 0 as the value at any P Q 0 ,g of the mapping Ψ : M → [-1, 1] characterized over the set M of all possible distributions of O by

Ψ(P Q,g ) ≡ E P Q,g {Q Y (1,W ) -Q Y (0,W )} = (Q Y (1, w) -Q Y (0, w))dQ W (w).
(1)

The mapping Ψ enjoys a remarkable property: it is pathwise differentiable (think "smooth") with an efficient influence curve (think "gradient") which provides insight into the asymptotic properties of all regular and asymptotically linear (think "well-behaved") estimators of Ψ(P Q 0 ,g ). The following lemma makes the latter statement more formal-we refer the reader to [START_REF] Bickel | Efficient and adaptive estimation for semiparametric models[END_REF][START_REF] Van Der | Asymptotic statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics[END_REF][START_REF] Van Der Laan | Unified methods for censored longitudinal data and causality[END_REF] for definitions and proofs.

Lemma 1. The mapping Ψ : M → [-1, 1] is pathwise differentiable at every P Q,g ∈ M with respect to (wrt) the maximal tangent space. Its efficient influence curve at P Q,g is D * (P Q,g ) which satisfies D * (P Q,g )(O) = D * W (P Q,g )(W ) + D * Y (Q, g)(O) with D * W (P Q,g )(W ) ≡ Q Y (1,W ) -Q Y (0,W ) -Ψ(P Q,g ), D * Y (Q, g)(O) ≡ 2A -1 g(A|W ) (Y -Q Y (A,W )) .
The variance Var P Q,g D * (P)(O) is a generalized Cramér-Rao lower bound for the asymptotic variance of any regular and asymptotically linear estimator of Ψ(P Q,g ) when sampling independently from P

Q,g . Moreover, if either Q Y = Q ′ Y or g = g ′ then E P Q,g D * (P Q ′ ,g ′ )(O) = 0 implies Ψ(P Q,g ) = Ψ(P Q ′ ,g ′ ).
The last statement of Lemma 1, often referred to as a "double-robustness" property, shows that one can seek help from D * to protect oneself against model misspecifications when estimating ψ 0 . This is especially relevant in our setting where we know precisely what is the randomization scheme g at play when one samples an observation from P Q 0 ,g .

Optimal Design

Suppose our goal of adaptation is to reach a randomization scheme of higher efficiency, i.e., to obtain a valid estimate of ψ 0 using as few blocks of patients as possible. By Lemma 1, the asymptotic variance of a regular, asymptotically linear estimator is lower-bounded by min g∈G Var P Q 0 ,g D * (P Q 0 ,g ). In this light, the Neyman design [START_REF] Hu | The theory of response-adaptive randomization in clinical trials[END_REF]

g 0 ≡ arg min g∈G Var P Q 0 ,g D * (P Q 0 ,g ) = arg min g∈G E P Q 0 ,g (Y -Q Y,0 (A,W )) 2 g 2 (A|W ) (2)
can be considered as an optimal randomization design ("optimal design" for short). Since its definition involves the unknown Q 0 , the optimal design g 0 is unknown too. It is readily seen that g 0 is characterized by

g 0 (1|W ) = σ 0 (1,W )/(σ 0 (1,W ) + σ 0 (0,W )), where σ 2 0 (A,W ) is the conditional variance of Y given (A,W ) under Q 0 .
It therefore appears that, under this randomization scheme, the treatment arm with higher probability for a patient with baseline covariates W is the one for which the conditional variance of the outcome is higher.

If we knew the optimal design then we could undertake the covariate-adjusted trial consisting in drawing independently observations from P Q 0 ,g 0 . The next task would be to build a regular, asymptotically linear estimator with asymptotic variance Var P Q 0 ,g 0 D * (P Q 0 ,g 0 ) based on the resulting data. In the present situation, we are going to "target" g 0 at some pre-determined interim steps. By targeting g 0 we mean estimating g 0 based on past observations and relying on the resulting estimator to collect the next block of data. In addition to targeting g 0 , each interim analysis will also consist in building an adaptive, targeted, regular and asymptotically linear estimator of ψ 0 . The details of this procedure are presented in Section 2.2.

Data-Generating Mechanism and Estimation Procedures

Describing the data-generating mechanism amounts to presenting how we target the optimal design g 0 at each interim step, which involves the estimation of the conditional expectation Q Y,0 . We initiate the description in Section 2.2.1, describe a LASSO estimation procedure of Q Y,0 in Section 2.2.2 and the related targeting procedure of g 0 in Section 2.2.3. By then, the data-generating mechanism is fully characterized by recursion.

Initiating the Data-Generating Mechanism

In the sequel, we denote O i ≡ (W i , A i ,Y i ) the ith observation that we sample. The indexing reflects the time ordering of the data collection: j < i implies that O j was collected before or at the same time as O i . For convenience, we let O n ≡ (O 1 , . . . , O n ) be the ordered vector of first n observations, with convention O 0 ≡ / 0. In the adaptive trial, the treatment A i is drawn conditionally on W i from the Bernoulli law with parameter g i (1|W i ), where the randomization scheme g i : A → [0, 1] depends on past observations O i-1 . We set g n ≡ (g 1 , . . . , g n ), the ordered vector of first n randomization schemes. The data-generating distribution of O n is denoted P Q 0 ,g n . It is formally characterized by the following factorization of the density of O n wrt the product of the dominating measures: for any g ∈ G ,

P Q 0 ,g n (O n ) = n ∏ i=1 P Q 0 ,g i (O i ) = n ∏ i=1 Q W,0 (W i ) × g i (A i |W i ) × P Q 0 ,g (Y i |A i ,W i ).
Let g b be the balanced randomization scheme, for which each arm is assigned with probability 1/2 regardless of baseline covariates. For a pre-specified n 0 , we first draw n 0 independent observations O 1 , . . . , O n 0 from P Q 0 ,g b . At an interim point, suppose one has thus far drawn n observations O n ∼ P Q 0 ,g n . An estimator of Q Y,0 is obtained based on O n . The next randomization scheme g n+1 is defined using the latter estimator and (O n , g n ), then the (n + 1)th observation O n+1 is drawn from P Q 0 ,g n+1 . We will describe the estimation of Q Y,0 and construction of g n+1 in the two following sections.

LASSO Estimation of the Outcome's Conditional Expectation

Consider {b n } n≥1 and {d n } d≥1 two non-decreasing, possibly unbounded sequences over R + and, for some M > 0 and every n ≥ 1, introduce the subset

B M,n ≡ β ∈ ℓ 1 : β 1 ≤ min(b n , M) and ∀ j ≥ d n , β j = 0 (3) of ℓ 1 ≡ β ∈ R N : ∑ j∈N |β j | < ∞ . Let φ j : j ∈ N be a uniformly bounded set of functions from A × W to R.
Without loss of generality, we may assume that φ j ∞ = 1 for all j ∈ N, where • ∞ denotes the supremum norm. For all β ∈ ℓ 1 , we denote

Φ β : A × W → R the function such that Φ β (A,W ) ≡ ∑ j∈N β j φ j (A,W ).
The construction of our LASSO estimators of Q Y,0 relies on a working model Q 1 and on a loss function L for Q Y,0 , both specified by the investigators. This means that Q Y,0 is the minimizer of

Q Y → P Q 0 ,g L(Q Y ) over the set of all conditional expectations of Y given (A,W ), of which Q 1 is a user-specified subset (the value of g ∈ G plays no role in this statement). For instance, they can take Q 1 ≡ {Q Y,β ≡ Φ β : β ∈ B M,n } with M = 1,
and the least-square loss function L characterized over the latter by

L(Q Y,β )(O) ≡ (Y -Q Y,β (A,W )) 2 . ( 4
)
They can also take

Q 1 ≡ {Q Y,β ≡ expit(Φ β ) : β ∈ B M,n } with M a deterministic upper-bound on | logit(Y )|
(recall that Y is assumed bounded away from 0 and 1), and the quasi negative-log-likelihood loss function L characterized over the latter by

-L(Q Y,β )(O) ≡ Y log(Q Y,β (A,W )) + (1 -Y ) log 1 -Q Y,β (A,W ) . (5)
Note that in both cases, for all

β ∈ B M,n , Q Y,β ∞ is upper-bounded by a deterministic upper-bound on |Y |.
Recall that we have already drawn n observations O n ∼ P Q 0 ,g n . Given a user-specified reference g r ∈ G that is bounded away from 0 and 1, we estimate Q Y,0 with Q Y,β n , where

β n ∈ arg min β ∈B M,n 1 n n ∑ i=1 L(Q Y,β )(O i ) g r (A i |W i ) g i (A i |W i ) . ( 6
)
The above minimization with the constraint β 1 ≤ min(b n , M), see (3), can be rewritten as a minimization free of the latter constraint by adding a term of the form λ n β 1 to the empirical criterion, where λ n depends on b n . This is the so-called LASSO procedure introduced by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] for the sake of obtaining estimators with fewer nonzero parameter values, thus effectively reducing the number of variables upon which the given solution is dependent. Note that when d n = d is held constant by choice, (6) should be interpreted as a standard parametric procedure rather than as a LASSO.

Adapting Towards the Optimal Design

We now turn to the construction of the next randomization scheme g n+1 .

Our optimal design minimizes g → Var P Q 0 ,g D * (P Q 0 ,g ) over the class G of all randomization schemes, see (2). We adopt a loss-based approach, by defining g n+1 as the minimizer in g of an estimator of Var P Q 0 ,g D * (P Q 0 ,g ) over a user-specified class of randomization schemes. This approach is applicable in the largest generality. In the case that W is discrete, or if one is willing to assign treatment only based on a discrete summary measure V of W , g n+1 can be defined explicitly as an estimator of the Neyman design based on Q Y,β n and observations O n ; we refer the readers to [START_REF] Chambaz | Inference in targeted goup sequential covariate-adjusted randomized clinical trials[END_REF] for details.

To proceed, we first note that, for all g ′ ∈ G ,

g 0 = arg min g∈G E P Q 0 ,g ′ (Y -Q Y,0 (A,W )) 2 g(A|W )g ′ (A|W ) .
This equality teaches us that for the sake of estimating g 0 using observations drawn from P Q 0 ,g ′ we may consider the loss function

L Q Y characterized over G by L Q Y (g)(O) ≡ (Y -Q Y (A,W )) 2 g(A|W ) ,
provided it is weighted by 1/g ′ (A|W ). Note that this loss function is indexed by a given Q Y .

Recall that we have already drawn n observations O n ∼ P Q 0 ,g n and estimated Q Y,0 with Q Y,β n . Now, given a class G 1 ⊂ G of randomization schemes uniformly bounded away from 0 and 1, we define the next randomization scheme as

g n+1 ∈ arg min g∈G 1 1 n n ∑ i=1 L Q Y,βn (g)(O i ) g i (A i |W i ) = arg min g∈G 1 1 n n ∑ i=1 Y -Q Y,β n (O i ) 2 g(A i |W i )g i (A i |W i ) (7)
This completes the description of our data-generating mechanism.

Targeted Maximum Likelihood Estimation

Given n observations O n ∼ P Q 0 ,g n and the estimator Q Y,β n of Q Y,0 defined in (6), we may carry out the estimation of the parameter of interest ψ 0 . We adopt the targeted minimum loss estimation methodology.

In the setting of a covariate-adjusted RCT with fixed design, a TMLE estimator is unbiased and asymptotically Gaussian regardless of the specification of the working model used for the estimation of Q Y,0 . It is known that unbiasedness and asymptotic normality still hold in the context of this technical report (CARA RCT for the estimation of ψ 0 based on copies of O), provided that the randomization schemes depend on W only through a summary measure taking finitely many values and that the working model used for the estimation of Q Y,0 be a simple linear model (this basically amounts to taking d n = d constant and b n = M) in Section 2.2.2), see [START_REF] Chambaz | Inference in targeted goup sequential covariate-adjusted randomized clinical trials[END_REF]). Yet by relying on more flexible randomization schemes and on more adaptive estimators of Q Y,0 we may achieve a greater efficiency through better estimation of the optimality criteria that may facilitate adaptation toward the optimal design, better adjustment of the variables that may directly improve on the estimation of the parameter of interest, and a more accurate estimation of the variance of the estimator. In a glimpse, the proposed strategy consists in targetedly fluctuating the initial estimator Q Y,β n by minimizing a pre-specified loss along a least favorable (wrt ψ 0 ) submodel through Q Y,β n , and then evaluating Ψ at the resulting updated estimator of Q 0 . Formally, consider the following one-dimensional parametric working model through Q Y,β n : for a given closed, bounded interval E ⊂ R containing 0 in its interior,

Q Y,β n (ε) ≡ expit logit(Q Y,β n ) + εH(g n ) : ε ∈ E , (8)
with notation H(g)(O) ≡ 2A-1 g(A|W ) for every g ∈ G 1 . This model passes through

Q Y,β n at ε = 0 in such a way that ∂ ∂ ε L(Q Y,β n (ε))| ε=0 = D * Y (Q Y,β n , g n ).
The optimal fluctuation parameter ε n minimizes the weighted empirical risk along the working model:

ε n ∈ arg min ε∈E 1 n n ∑ i=1 L(Q Y,β n (ε))(O i ) g n (A i |W i ) g i (A i |W i ) . (9) Set Q * Y,β n ≡ Q Y,β n (ε n ) then Q * β n ≡ (Q W,n , Q * Y,β n ) where Q W,n
is the empirical marginal distribution of the W . The TMLE estimator of ψ 0 is finally defined as

ψ * n ≡ 1 n n ∑ i=1 Q * Y,β n (1,W i ) -Q * Y,β n (0,W i ).
It satisfies ψ * n = Ψ(P Q * βn ,g ) for any g ∈ G .

Asymptotics

We first introduce further notation in Section 3.1 then we successively investigate the convergence of the targeted CARA design in Section 3.2 and the asymptotic behavior of the TMLE estimator in Section 3.3.

Notation

In general, given a known g ∈ G and an observation O drawn from P Q 0 ,g , Z ≡ g(A|W ) is a deterministic function of g and O. Note that Z should be interpreted as a weight associated with O and will be used as such. Therefore, we can augment O with Z, i.e., substitute (O, Z) for O, while still denoting (O, Z) ∼ P Q 0 ,g . In particular, during the course of our trial, conditionally on O i-1 , the randomization scheme g i is known and we can substitute

(O i , Z i ) = (O i , g i (A i |W i )) ∼ P Q 0 ,g i for O i drawn from P Q 0 ,g i .
f ≡ n -1 ∑ n i=1 f (O i , Z i ). Likewise, for any fixed P Q,g ∈ M , P Q,g f ≡ E P Q,g f (O, Z) and, for each i = 1, . . . , n, P Q 0 ,g i f ≡ E Q 0 ,g i [ f (O i , Z i )|O i-1 ], P Q 0 ,g n f ≡ n -1 ∑ n i=1 E Q 0 ,g i [ f (O i , Z i )|O i-1 ]. We endow the set {Q Y : P Q,g ∈ M } of all conditional expectations of Y given (A,W ) under P Q,g ∈ M with the norm • Y,0 characterized by Q Y -Q ′ Y 2 Y,0 ≡ E P Q 0 ,g r Q Y (A,W ) -Q ′ Y (A,W ) 2 .
Similarly, we endow the set G with the norm • A,0 characterized by

g -g ′ 2 A,0 ≡ E Q W,0 g(1|W ) -g ′ (1|W ) 2 .
For any class F of functions equipped with a norm

• and ε > 0, N(F , • , ε) is the ε-bracketing of F wrt • and J(1, F , • ) ≡ 1 0 log N(F , • ,
ε)dε is the corresponding bracketing entropy (evaluated at 1). Finally, the uniform norm of a real-valued operator

Π on F is Π F ≡ sup f ∈F |Π( f )|.

Convergence of the Targeted CARA Design

Our first concern is the convergence of Q Y,β n , see (6).

Proposition 1 (Convergence of Q Y,β n ). Consider the following assumptions:

A1. The conditional density under Q 0 of Y given (A,W ) wrt some dominating measure is bounded away from 0.

A2. There exists a unique β 0 ∈ n≥1 B M,n such that

β 0 ∈ arg min β ∈ n≥1 B n P Q 0 ,g r L(Q Y,β ).
A3. It holds that d n = O(n r ) for some r > 0 and

sup β ∈B M,n |(P n -P Q 0 ,g n )L(Q Y,β )| = o P (1). Under A1-A3, Q Y,β n -Q Y,β 0 Y,0 → 0 in probability.
In summary, Q Y,β n converges to Q Y,β 0 in probability for the norm • Y,0 if such a limit exists and if the dimension of the LASSO parameters grows polynomially wrt the sample size. Note that the limit Q Y,β 0 depends on the user-supplied reference design g r . Based on [START_REF] Dümbgen | Nemirovski's inequalities revisited[END_REF] we show that A3 holds for instance when L is given by ( 4) and

Q 1 = {Φ β : β ∈ B 1,n }.
The proof of Proposition 1 relies on empirical process theory for martingales, chaining arguments and tools developed by [START_REF] Van Handel | On the minimal penalty for markov order estimation[END_REF]. It requires that the sequence {g n } n≥1 of randomization schemes be uniformly bounded away from 0 and 1, which is guaranteed by specification of the user-supplied set G 1 .

We now turn to the convergence of the targeted CARA design {g n } n≥1 , see ( 7), toward a fixed, limiting design g * 0 ∈ G 1 . Proposition 2 (Convergence of the targeted CARA Design). Consider the setup of Proposition 1 and the following additional assumptions:

A4. There exists a unique g * 0 ∈ G 1 such that g * 0 ∈ arg min

g∈G 1 P Q 0 ,g r L Q Y,β 0 (g) g r . ( 10
) A5. The class 1/G 1 ≡ {1/g : g ∈ G 1 } satisfies the finite entropy condition J(1, 1/G 1 , • A,0 ) < ∞.
Under A1-A5, g n (1|W )g * 0 (1|W ) A,0 → 0 in probability. We have already emphasized that through the choice of G 1 , the investigators of the RCT benefit from a great flexibility in treatment allocation. The main constraint on G 1 is A4, a condition on the complexity/richness of the class. We refer the reader to (van der Vaart, 1998, Examples 19.7-19.11, Lemma 19.15) for typical examples. They notably include "well-behaved" parametric classes and VC classes. In particular, G 1 can consist of randomization schemes such that the allocation probabilities only depend on W through a discrete summary measure of it, as considered in [START_REF] Chambaz | Inference in targeted goup sequential covariate-adjusted randomized clinical trials[END_REF].

The limiting randomization scheme g * 0 depends on the user-supplied reference design g r only through Q Y,β 0 : replacing g r with any g ∈ G in (10) does not alter the definition of g 0 * . Furthermore, g * 0 can be interpreted as the most optimal design in G 1 given the limiting conditional outcome model Q Y,β 0 :

g * 0 ∈ arg min g∈G 1 Var P Q 0 ,g D * Y (Q Y,β 0 , g) = arg min g∈G 1 Var P Q 0 ,g D * Y (Q 0 , g) + P Q 0 ,g Q Y,0 -Q Y,β 0 2 g 2 .
Comparing the above equality with (2) yields that g * 0 = g 0 , the Neyman design, whenever Q Y,β 0 = Q Y,0 . In general, g * 0 minimizes an objective function writing as the sum of the Cramér-Rao lower bound and a second-order term residual. This underscores the motivation for using a flexible estimator in estimating Q Y,0 : by minimizing this second-order residual of the limiting conditional outcome model, we are closer to adapting toward the desired optimal design.

Consistency and Central Limit Theorem

As with the initial LASSO estimators of the conditional outcome, we are firstly concerned with the convergence of the updated estimators Q * Y,β n :

Proposition 3 (Consistency). Consider the setups of Propositions 1 and 2 and the follow additional assumption:

A6. There exists a unique ε 0 ∈ E such that

ε 0 ∈ arg min ε∈E P Q 0 ,g * 0 L(Q Y,β 0 (ε)).
Assume that A1-A6 are met and define

Q * Y,β 0 ≡ Q Y,β 0 (ε 0 ). It holds that Q * Y,β n -Q * Y,β 0 Y,0 → 0 in probability. Moreover, ψ * n consistently estimates ψ 0 .
If Q Y,β 0 = Q Y,0 then ε 0 = 0: the updating procedure preserves the consistency of the initial estimator

Ψ(P Q Y,βn ,g ) for any g ∈ G . More importantly, Proposition 3 guarantees that even if Q Y,β 0 = Q Y,0 then ψ *
n still consistently estimates ψ 0 , by double-robustness. Nonetheless, the convergence of the updated estimators Q * Y,β n (to the truth or otherwise) is crucial for studying the asymptotic behavior of ψ * n .

Proposition 4 (Central Limit Theorem for ψ * n ). Consider the setups of Propositions 1, 2 and 3 and the following additional assumption:

A7. For any deterministic function F, F(O) = 0 P Q 0 ,g * 0 -almost surely implies that F = 0. Assume that A1-A7 are met. For both β = β 0 and β = β n , introduce d * Y,β and q * Y,β characterized by d * Y,β (O, Z) ≡ 2A -1 Z Y -Q * Y,β (A,W ) , q * Y,β ≡ Q * Y,β (1,W ) -Q * Y,β (0,W )
and, for any g ∈ G ,

Σ n ≡ 1 n n ∑ i=1 d * Y,β n (O i , Z i ) + D * W (P * Q βn ,g (W i )) 2 . ( 11
)
Then (Σ n /n) -1/2 (ψ * n -ψ 0 ) converges in distribution to the standard normal distribution.

The asymptotic results in Proposition 4 underpin the statistical analysis of the proposed targeted CARA RCT. In particular, denoting ξ 1-α/2 the (1 -α/2)-quantile of the standard normal distribution,

ψ * n ± ξ 1-α/2 (Σ n /n) 1/2 is a confidence interval of asymptotic level (1 -α).

Simulation Study

We present here the results of a simulation study of the performances of the targeted procedure exposed in the previous sections.

Simulation Scheme

We rely on the same simulation scheme as in [START_REF] Chambaz | Inference in targeted goup sequential covariate-adjusted randomized clinical trials[END_REF]. For completeness, let us recall that Q 0 is such that:

• the baseline covariate W equals (U,V ), where U and V are independently drawn with U uniformly distributed on [0, 1] and

Q W,0 (V = 1) = 1/2, Q W,0 (V = 2) = 1/3, Q W,0 (V = 3) = 1/6;
• the conditional distribution of Y given (A,W ) is the Gamma distribution with conditional mean

Q Y,0 (A,W ) = 2U 2 + 2U + 1 + AV + 1 -A 1 +V
and conditional variance

σ 2 0 (Y |A,W ) = U + A(1 +V ) + 1 -A 1 +V 2 .
The marginal treatment effect on an additive scale satisfies ψ 0 = 91 72 ≃ 1.264. We target the optimal designs corresponding to eight parametric working models G 11 , . . . , G 18 that we present in Table 1.

working model parametric form dimension optimal variance G 11 θ 0 1 18.50

G 12 ∑ 3 v=1 θ v 1{V = v} 3 18.18 G 13 θ 0 + θ 1 U 2 18.37 G 14 ∑ 3 v=1 θ v 1{V = v} + θ 4 U 4 18.05 G 15 θ 0 + ∑ 3 v=1 θ v 1{V = v}U 4 18.12 G 16 ∑ 3 v=1 θ v 1{V = v} + θ 4 U + ∑ 3 v=2 θ 3+v 1{V = v}U 6 18.01 G 17 θ 0 + ∑ 3 v=1 θ v 1{V = v}U + ∑ 3 v=1 θ 4+v 1{V = v}U 2 7 18.36 G 18 ∑ 3 v=1 θ v 1{V = v} + θ 4 U + θ 5 U 2 + ∑ 3 v=2 θ 4+v 1{V = v}U + ∑ 3 v=2 θ 6+v 1{V = v}U 2 9 18.03
Table 1: Parametric working models G 1k (k = 1, . . . , 8). In the second column, we report the parametric forms of logit ((g θ (W ) -δ )/(1 -2δ )) for generic elements g θ ∈ G 1k (k = 1, . . . , 8). We set δ = 10 -2 . In the third column, we give the dimensions of the models. In the fourth column, we report the numerical values of arg min g∈G 1k Var P Q 0 ,g D * (P Q 0 ,g )(O) (k = 1, . . . , 8), with precision 10 -2 . Recall that Var P Q 0 ,g b D * (P Q 0 ,g )(O) = 23.87, with precision 10 -2 .

In addition to the latter parametric working models, we consider eight statistical procedures for the estimation of the conditional expectation Q Y,0 . Four of them consist in parametric estimation on smalldimensional models Q 11 , . . . , Q 14 . In contrast, the four others rely on moderate-dimensional parametric models, ℓ 1 -penalization and cross-validation to select the best regularization parameter. We denote Q 15 , . . . , Q 18 these "machine-learning", as opposed to "parametric", procedures/models, which embody the LASSO estimating procedure of Section 2.2.2. We summarize in Table 2 what are Q 11 , . . . , Q 18 . All procedures involve the logistic loss, even though the support of the marginal distribution of Y under P 0 is R + , not [0, 1]. In fact, given a sample O 1 , . . . , O n , we first scale Y 1 , . . . ,Y n to [0, 1], then regress the scaled outcomes on (A,W ) based on the logistic loss and one procedure among Q 11 , . . . , Q 18 , then scale back the resulting conditional expectation to the original range of the observed outcomes.

Set B = 1000 and let n = (250,500,750,1000,1250,1500,1750,2000,2250,2500) be a sequence of sample sizes. For each combination (k, l) ∈ {1, . . . , 8} 2 , we repeatedly simulate B = 1000 times a targeted CARA RCT based on G 1k and Q 1l , performing an update of the randomization scheme and the computation of the TMLE of ψ 0 at every intermediate sample size n i (1 ≤ i ≤ 10), which we denote ψ * n i ,klb . The simulations are mutually independent. The associated 95%-confidence intervals I n i ,klb rely on estimated working model parametric form dimension parametric

Q 11 ∑ 3 v=1 θ v 1{V = v} + θ 4 U + θ 5 A 5 Q 12 θ 0 + A θ 1 U + ∑ 3 v=2 θ v 1{V = v} +(1 -A) θ 4 U + ∑ 3 v=2 θ 3+v 1{V = v} 7 Q 13 A ∑ 3 v=1 θ v 1{V = v} + θ 4 U +(1 -A) ∑ 3 v=1 θ 4+v 1{V = v} + θ 8 U 8 Q 14 A ∑ 3 v=1 θ v 1{V = v} + θ 4 U + θ 5 U 2 +(1 -A) ∑ 3 v=1 θ 5+v 1{V = v} + θ 9 U + θ 10 U 2 10 LASSO Q 15 A ∑ 3 v=1 θ v 1{V = v} + θ 4 U + θ 5 U 2 +(1 -A) ∑ 3 v=1 θ 5+v 1{V = v} + θ 9 U + θ 10 U 2 10 Q 16 A ∑ 3 v=1 θ v 1{V = v} + ∑ 5 l=1 θ 3+l U l +(1 -A) ∑ 3 v=1 θ 8+v 1{V = v} + ∑ 5 l=1 θ 11+l U l 16 Q 17 A ∑ 3 v=1 θ v 1{V = v} + ∑ 10 l=1 θ 3+l U l +(1 -A) ∑ 3 v=1 θ 13+v 1{V = v} + ∑ 10 l=1 θ 16+l U l 26 Q 18 A ∑ 3 v=1 θ v 1{V = v} + ∑ 20 l=1 θ 3+l U l +(1 -A) ∑ 3 v=1 θ 23+v 1{V = v} + ∑ 20 l=1 θ 26+l U l 46
Table 2: Working models Q 1k (k=1, . . . , 8) for the conditional expectation Q Y,0 . In the second column, we report the parametric form of logit ((q θ (A,W ) -δ )/(1 -2δ )) for generic elements q θ ∈ Q 1k (k = 1, . . . , 8). We set δ = 10 -2 . In the third column, we give the dimensions of the models. All working models are exploited in combination with the quasi negative-log-likelihood loss function (5

i ,kl = 1 B B ∑ b=1 ψ * 2 n i ,klb - 1 B B ∑ b=1 ψ n i ,klb 2
and the empirical coverage of the corresponding confidence interval

Ĉn i ,kl = 1 B B ∑ b=1 1{ψ 0 ∈ I n i ,klb }.
The simulation study is conducted using R (R Core Team, 2014) and the package glmnet [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF].

Discussion of the Results

Coverage

We propose an evaluation of the coverage performances based on testing. For every (k, l) ∈ {1, . . . , 8} 2 and n i (1 ≤ i ≤ 10), the statistic B × Ĉn i ,kl follows a Binomial distribution with parameter (B, π n i ,kl ) for some π n i ,kl ∈ [0, 1]. Denote p95 n i ,kl the exact p-value of the one-sided binomial test of H 95 n i ,kl : "π n i ,kl ≥ 95%" against "π n i ,kl < 95%". Under H 95 n i ,kl , p95 n i ,kl is drawn from the uniform distribution on [0, 1]. For every n i (1 ≤ i ≤ 10), we carry out one-sample Kolmogorov-Smirnov tests of the null stating that the common law of { p95 3. The first row corresponds to the choice L = {1, . . . , 8}. It teaches us that the expected 95%-coverage is generally not guaranteed. One may wonder if the same conclusion holds when focusing in turn on the parametric procedures (set L = {1, . . . , 4}) or on the LASSO procedures (set L = {5, . . . , 8}). Inspecting the second and third rows of Table 3 does not reveal an interesting pattern. One may now wonder to what extent the 95%-coverage is deficient. To answer this question, we proceed similarly. We denote p94 n i ,kl the exact p-value of the one-sided binomial test of H 94 n i ,kl : "π n i ,kl ≥ 94%" against "π n i ,kl < 94%". Under H 94 n i ,kl , p94 n i ,kl is drawn from the uniform distribution on [0, 1]. For every n i (1 ≤ i ≤ 10), we carry out a one-sample Kolmogorov-Smirnov test of the null stating that the common law of { p94 n i ,kl : 1 ≤ k ≤ 8, 1 ≤ l ≤ 8} is the uniform distribution on [0, 1] against the alternative that the common law is stochastically smaller than the uniform distribution on [0, 1]. The p-values of these tests are reported in the fourth row of Table 3. The conclusion is clear and satisfactory: even if the 95%-confidence intervals fail to guarantee the wished coverage, one can safely consider them as valid 94%-confidence intervals.

n i ,kl : 1 ≤ k ≤ 8, l ∈ L } (L ⊂ {1, . . . ,

Standard Deviation

Here we investigate how the targeted CARA RCT behaves in terms of standard deviation of the produced estimators. As in the previous subsection, the investigation relies on testing. For every (k, l) ∈ {1, . . . , 8} 2 and n i (1 ≤ i ≤ 10), we first compute the statistic

T n i ,kl = 1 B ∑ B b=1 (Σ n i ,klb ) 1/2 -Ŝn i ,kl 1/2 1 B ∑ B b=1 Σ n i ,klb -1 B ∑ B b=1 (Σ n i ,klb ) 1/2 2 1/2 ,
where Σ n i ,klb is the estimated variance of the TMLE produced at intermediate sample size n i by the bth simulated targeted CARA RCT based on G 1k and Q 1l , see (11). Thus, T n i ,kl sheds some light on the estimation of the standard deviation of the TMLE ψ * n i at sample size n i by (Σ n i /n) 1/2 for the targeted CARA RCT based on G 1k and Q 1l .

For every n i (1 ≤ i ≤ 10), we perform a Lilliefors test of normality based on the sample {T n i ,kl : 1 ≤ k ≤ 8, l ∈ L } with L = {1, . . . , 8}. The p-values of these tests are reported in Table 4. They teach us that there is no stark evidence of non-normality across the ten intermediate sample sizes. This first conclusion justifies the next step: for every n i (1 ≤ i ≤ 10), we perform a one-sided Student test of "µ n i ≥ 0" against "µ n i < 0", where µ n i denotes the mean of the common distribution of {T n i ,kl : 1 ≤ k ≤ 8, l ∈ L } with L = {1, . . . , 8}. The p-values of these tests are reported in the two first rows of Table 4. Adjusting for multiple testing in terms of the Benjamini and Yekutieli procedure for controlling the false discovery rate at the 5% level, we conclude that estimating the variance as in ( 11) is over-optimistic at least for intermediate sample sizes smaller than or equal to n 3 = 750. One may wonder if the same conclusions hold when focusing in turn on the procedures (set L = {1, . . . , 4}) or on the LASSO procedures (set L = {5, . . . , 8}). Inspecting separately the third and fourth rows of Table 4 on one hand then the fifth and sixth rows on the other hand leads to the conclusion that estimating the variance as in ( 11) is overoptimistic only for intermediate sample sizes smaller than or equal to n 2 = 500, still adjusting for multiple testing in terms of the Benjamini and Yekutieli procedure for controlling the false discovery rate at the 5% level.

The gap between the conclusions reached when considering all procedures or the parametric and LASSO ones separately may be simply explained by a loss of power due to the reduction of sample size (64 versus 32), or by subtle differences induced by the nature of Q 1l . In any case, in light of Section 4.2.1, the under-estimation of the true variance based on ( 11 In the first row we report the p-values of the Lilliefors tests of normality of the sample {T n i ,kl : 1 ≤ k, l ≤ 8} (1 ≤ i ≤ 10). In the second row, we report the p-values of the Student tests of "µ n i ≥ 0" against "µ n i < 0", where µ n i denotes the mean of the common distribution of {T n i ,kl : 1 ≤ k, l ≤ 8}. In the third and fourth rows (fifth and sixth rows, respectively), we report the p-values of the same Lilliefors and Student tests based on the samples {T n i ,kl : 1 ≤ k ≤ 8, 1 ≤ l ≤ 4} corresponding to parametric procedures (on the samples {T n i ,kl : 1 ≤ k ≤ 8, 5 ≤ l ≤ 8} corresponding to LASSO procedures, respectively).

Discussion

We have presented in this technical report a new group-sequential CARA RCT design and inferential procedure built on top of it. The procedure is targeted in the sense that (i) the sequence of randomization schemes is group-sequentially determined by targeting a user-specified optimal randomization design based on accruing data and, (ii) our estimator of the user-specified parameter of interest, seen as the value of a functional evaluated at the true, unknown distribution of the data, is targeted toward it by following the paradigm of targeted minimum loss estimation. We focused for clarity on the case that the parameter of interest is the marginal effect of a binary treatment and that the targeted optimal design is the Neyman allocation, in an effort to produce an estimator with smaller asymptotic variance, but our methodology extends beyond this instructive framework. For clarity too, we considered the case that the estimator of the conditional outcome given treatment and baseline covariates, a key element of the procedure, is obtained by LASSO regression, although our methodology can hinge on a wide class of data-adaptive estimators. Under mild assumptions, the resulting sequence of randomization schemes converges to a limiting design, and the TMLE estimator is consistent and asymptotically Gaussian, with an asymptotic variance that we can estimate too. Thus we can build valid confidence intervals of given asymptotic levels. A simulation study confirms our theoretical results. Across 64 different choices of pairs of working models and 10 intermediate sample sizes ranging from 250 to 2500, there is no empirical evidence that our 95%-confidence intervals do not provide at least 94%-coverage, based on 1000 independent replications. In addition, in the same framework, there is no empirical evidence that our estimators of the variances of the TMLE estimators are over-optimistic for sample sizes larger than 500, adjusting for multiple testing in terms of the Benjamini and Yekutieli procedure for controlling the false discovery rate at the 5% level.

For smaller sample sizes, the under-estimation is slight at most. We will soon make available a R package to allow interested readers to test the procedure. The proofs of our theoretical results will be presented in a forthcoming article. In this article, we will also describe and study a more general targeted CARA RCT design and related TMLE methodology involving possibly aggressive data-adaptive/machine-learning procedures and not only LASSO regression. In the future, we will also consider alternative strategies to randomly assign successive patients to the treatment arms in such a way that the overall empirical conditional distribution of treatment given baseline covariates be as close as possible to the current best estimator of the targeted optimal design. This will require both new theoretical developments and simulation studies.

  ). Models Q 11 , Q 12 , Q 13 , Q 14 are straightforwardly fitted by relying on the R function glm. Models Q 15 , Q 16 , Q 17 , Q 18 are LASSO-fitted by relying on the R function glmnet. variances of the TMLE as given in (11). For each combination (k, l) and intermediate sample size n i , we compute the empirical variance of the corresponding TMLE Ŝn

Table 3 :

 3 Evaluating the coverage performances based on testing. The first row gives p-values of Kolmogorov-Smirnov tests of the null consisting of the intersection of all H 95 n i ,kl . The second and third rows give p-values of Kolmogorov-Smirnov tests of the nulls consisting in the intersections of all H 95 n i ,kl based on parametric procedures (second row) and of all H 95 n i ,kl based on LASSO procedures (third row). The fourth row gives p-values of Kolmogorov-Smirnov tests of the null consisting of the intersection of all H 94 n i ,kl . the alternative that the common law is stochastically smaller than the uniform distribution on [0, 1]. Rejecting the null for its alternative indicates a defective coverage. The p-values of four such Kolmogorov-Smirnov tests are reported in Table

	8}) is the uniform distribution on [0, 1] against

Table 4 :

 4 ) is necessarily slight at most. Investigating the targeted CARA RCT in terms of standard deviation of the produced estimators.

	n i	250	500	750	1000 1250 1500 1750 2000 2250 2500
	Lilliefors	0.670	0.330	0.866 0.033 0.538 0.837 0.133 0.528 0.466 0.022
	Student	< 0.001 < 0.001 0.002 0.006 0.008 0.012 0.007 0.008 0.044 0.038
	Lilliefors	0.755	0.043	0.270 0.021 0.543 0.620 0.206 0.172 0.685 0.206
	Student	< 0.001 < 0.001 0.013 0.026 0.025 0.026 0.021 0.036 0.226 0.420
	Lilliefors	0.561	0.894	0.864 0.517 0.500 0.314 0.251 0.783 0.971 0.283
	Student	< 0.001 < 0.001 0.044 0.059 0.087 0.116 0.084 0.063 0.050 0.011

Acknowledgements.

Antoine Chambaz acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-13-BS01-0005 (project SPADRO). Mark van der Laan was supported by NIH grant R01 AI074345-06A1. Wenjing Zheng acknowledges the support of the Embassy of France in the United States through a Chateaubriand Fellowship-Science, Technology, Engineering, and Mathematics (STEM).


Xavier Guillonneau

Florian Sennlaub

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Letter to the Editors Comment on "Ccl2, Cx3cr1 and Ccl2/Cx3cr1 chemokine deficiencies are not sufficient to cause age-related retinal degeneration" by Luhmann et al. (Exp Eye Res. 2012 Dec 8. doi:pii: S0014-4835(12) 00342-9.) Christophe Combadière, William Raoul, Xavier Guillonneau, Florian Sennlaub We are writing to you concerning a recent publication in this journal by Luhmann et al. (Exp Eye Res. 2013;107:80). In this manuscript the authors describe the absence of age-related subretinal macrophage/microglial cells (Mφ/MC) and retinal degeneration in notably the Cx3cr1 -/-mice. The authors suggest that their study contradicts Chinnery et al. (Neurobiol Aging. 2012;33(8):1769) and our previous reports of age -and light-induced (Combadiere et al., J Clin Invest. 2007;117(10):2920 ;Raoul et al., J Neuroimmunol. 2008;198(1-2):56) subretinal Mφ/MC accumulation and associated photoreceptor degeneration. This phenotype was described in both, pigmented C57BL/6 and BALB albino background, in Cx3cr1 -/-knockout and Cx3cr1 GFP/GFP knockin mice (see references cited above). The title of the manuscript could suggest to the reader that our results are not reproducible and the authors state that their "findings are in clear contrast to several reports that describe an agerelated retinal degeneration in other chemokine knockout mouse lines", including the Cx3cr1 -/- mice. We would like to explain why we think that such conclusions cannot be drawn from the presented data.

A contamination with the rd8 mutation of the Crb1 gene has recently been recognized to cause early onset, severe retinal degeneration in Ccl2 -/-Cx3cr1 -/-mice, as described in several publications, independently of their Ccl2 and Cx3cr1 deletions. Luhman et al. question "whether the previously reported phenotypes of different chemokine single knockout mice might also have been affected by the rd8 mutation". As communicated to the authors, the Cx3cr1 -/-mice and Cx3cr1 GFP/GFP mice used in our experiments are not contaminated by the rd8 mutation.

We have previously reported that ambient light conditions are crucial for the Cx3cr1-dependent increase in subretinal Mφ/MC accumulation, as raising Cx3cr1 -/-BALB mice in the dark prevented (Combadiere et al. 2007, cited above), and a bright-light-challenge that does not induce subretinal inflammation in C57BL6/J mice, induced the phenotype in Cx3cr1 GFP/GFP C57BL6/J mice (Raoul et al. 2008, cited above). We observed the age-dependent presence of Mφ/MC in Cx3cr1 -/-and Cx3cr1 GFP/GFP C57BL6/J mice that were raised under 12-h light/12-h dark cycles (100-500 lux at the cage level, with no additional cover in the cage; Combadiere et al. 2007, cited above). Since our publication in 2007, we have observed this phenotype in over 30 Cx3cr1-deficient C57BL6/J mice with age (12 months and older; raised in the indicated light conditions) and in over 40 light-exposed Cx3cr1-deficient C57BL6/J mice, compared to a similar amount of wildtype controls. Chinnery et al. reproduced the increased subretinal Mφ/MC accumulation in Cx3cr1 GFP/GFP C57BL6/J and Cx3cr1 GFP/GFP BALB mice (also rd8-free, personal

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

communication) at 230-275 lux in young adults. Under their conditions the accumulation further increased with age and was not statistically different anymore in Cx3cr1 GFP/+ Cx3cr1 GFP/GFP at a massive 180 Mφ/MCs/mm 2 in 20-month-old BALB mice, suggesting that a maximal Mφ/MC cell accumulation can be reached in the subretinal space. Ng and Streilein (Invest Ophthalmol Vis Sci. 2001;42(13):3301) showed that subretinal Mφ/MC accumulate in albino mice and subretinal Mφ/MCs are cleared when transferred from light to dark conditions. In Cx3cr1 -/-mouse brains, macrophages clear less efficiently from the injection site compared to wildtype macrophages (Cardona et al., Nat Neurosci. 2006;9(7):917). In preliminary studies (Levy et al., ARVO 2011), we have shown that Cx3cr1 -/-macrophages present a similar impaired clearance from the subretinal space, which could explain the subretinal Mφ/MC accumulation in Cx3cr1-deficient animals we described.

The mice used in the studies described by Luhmann et al. were raised under a 12h/12h dark-light cycle with a mean luminescence during the light period at the bottom of the cage of 33 ± 28 lx. Furthermore, the animals had access to protection from light exposure (e.g., paper roll and excess of bedding) inside the cage, which allowed them to burrow. The light-conditions used in the described experiments were, therefore, significantly dimmer than the light conditions used in the previous publications.

We never suggested that Cx3cr1 -/-or Cx3cr1 GFP/GFP C57BL6/J mice raised under such dim-light conditions accumulate subretinal Mφ/MCs or display age-related photoreceptor degeneration. In fact we have emphasized in all our related publications that light conditions are crucial for obtaining the phenotype. Therefore, we believe that the only conclusion that can be drawn safely from the presented data is that none of the mouse strains used develop subretinal Mφ/MC accumulation or retinal degeneration under very dim light conditions. The comparison of theresults described by Luhmann et al. to those obtained in our studies and their interpretation as a "clear contrast" to our published findings (which suggests our results were not reproducible in their laboratory) are not justified, as the conditions used in the two studies were not comparable.


Olivier Nicol

Jérémie Mary

Jeremie Mary

Philippe Preux
email: philippe.preux@univ-lille3.

Improving offline evaluation of contextual bandit algorithms via bootstrapping techniques

In many recommendation applications such as news recommendation, the items that can be recommended come and go at a very fast pace. This is a challenge for recommender systems (RS) to face this setting. Online learning algorithms seem to be the most straight forward solution. The contextual bandit framework was introduced for that very purpose. In general the evaluation of a RS is a critical issue. Live evaluation is often avoided due to the potential loss of revenue, hence the need for offline evaluation methods. Two options are available. Model based methods are biased by nature and are thus difficult to trust when used alone. Data driven methods are therefore what we consider here. Evaluating online learning algorithms with past data is not simple but some methods exist in the literature. Nonetheless their accuracy is not satisfactory mainly due to their mechanism of data rejection that only allow the exploitation of a small fraction of the data. We precisely address this issue in this paper. After highlighting the limitations of the previous methods, we present a new method, based on bootstrapping techniques. This new method comes with two important improvements: it is much more accurate and it provides a measure of quality of its estimation. The latter is a highly desirable property in order to minimize the risks entailed by putting online a RS for the first time. We provide both theoretical and experimental proofs of its superiority compared to state-of-the-art methods, as well as an analysis of the convergence of the measure of quality.

Introduction

Under various forms, and under various names, recommendation has become a very common activity over the web. One can think of movie recommendation (Netflix), e-commerce (Amazon), online advertising (everywhere), news recommendation (Digg), personalized radio stations (Pandora) or even job recommendation (LinkedIn)... All these applications have their own characteristics. Yet the common key idea is to take advantage of some information we may have on a user (profile, demographics, time of the day etc.) in order to identify the most attractive content to serve him/her in a given context. Note that an element of content is generally referred to as an item. To perform recommendation, a piece of software called a recommender system (RS) can use past user/item interactions such as clicks or ratings. In particular, the typical approach to recommendation is to train a predictor of ratings and/or clicks of users on items on past data and use the resulting predictions to make personalized recommendations. This approach is based on the implicit assumption that past behavior can be used to predict future behavior.

In this paper we consider a recommendation applications in which the aforementioned assumption is not reasonable. We refer to this setting as dynamic recommendation. In dynamic recommendation, only a few tens of different items are available for recommendation at any given moment. These items have a limited lifetime and are continuously replaced by new ones, with different characteristics. We also consider that the tastes of users change, sometimes dramatically due to external parameters that we do not control. Many examples of dynamic recommendation exist on the web. The most popular one is news recommendation that can be found on specialized websites such as Digg or on general web portals (Yahoo!) and websites of various media (newspapers, TV channels...). Other examples can be mentioned such as private auctions in which the user can buy a limited set of items that changes everyday. Another example is a RS that can only recommend the K most recent items (this may apply to movies, videos, songs...). This problem has begun to be addressed recently with online learning solutions, by considering the contextual bandit framework. Nonetheless this is not the case in most of the recommendation literature. In all the textbooks, dynamic recommendation is handled with content-based recommendation. The idea is to consider an item as a set of features and to try to predict the taste of a user with regards to these features, using an offline predictor as before. Yet we argue that although this can be a good idea in some special cases, this is not the way to go in general:

1. It requires a continuous labeling effort of new items.

2. We are limited by what the expert labels: things can be hard to label such as the appeal for a picture, the quality of a textual summary, etc.

3. Tastes are not static: the appeal of a user to some kind of news can be greatly impacted by the political context. Similarly the appeal towards clothing can be impacted by fashion, movie stars etc.

For such systems, the best way to compare the performance of two algorithms is to perform A/B testing on a subset of the web audience [START_REF] Kohavi | Controlled experiments on the web: Survey and practical guide[END_REF]). Yet there is almost no e-commerce website that would let a new RS go live just for testing, even on a small portion of the audience for fear of hurting the user experience and loosing money.

The entailed engineering effort can also be discouraging. Therefore being able to evaluate offline a RS is crucial. In classical recommendation, the measures of prediction accuracy and other metrics that can be computed on past data are well accepted and trusted in the community. Nevertheless for the reasons we gave above, they are irrelevant for online learning algorithms designed for dynamic recommendation. This paper is about the offline evaluation of such algorithms. Some fairly simple replay methodologies do exist in the literature. Nonetheless they have a well known, yet very little studied drawback which is that they use a very small portion of the data at hand. One may argue that the tremendous amount of data available with web applications makes this a marginal problem. Yet in this paper we will exhibit that this is a major issue that generates a huge bias when evaluating online learning algorithms for dynamic recommendation. Furthermore we will explain that acquiring more data does not solve the problem. Then we will propose a solution to this issue, that builds on the previously introduced methods and on different elements of bootstrapping theory. This solution is backed by a theoretical analysis as well as an empirical study. They both clearly exhibit the improvement in terms of evaluation accuracy brought by our new method. Furthermore the use of bootstrapping allows us to estimate the distribution of our estimation and therefore to get a sense of its quality. This is a highly desirable property, especially considering that such an evaluation method is designed in order to decide whether we should risk putting a new algorithm online. The fast theoretical convergence of this estimation is also proved in our analysis. Note that the experiments are run on synthetic data, for reasons that we will detail and also on a large publicly available dataset.

Background on bandits and notations

We motivated the need for online learning solutions in order to deal with dynamic recommendation. A natural way to model this situation is as a reinforcement learning problem [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF], and more precisely using the contextual bandit framework [START_REF] Lu | Contextual multi-armed bandits[END_REF] that was introduced for the very purpose of news recommendation.

Contextual Bandits

The bandit problem is also known in the literature as the multi-arm bandit problem and other variations. This problem can be traced back to [START_REF] Robbins | Some aspects of the sequential design of experiments[END_REF]Munro in 1952 (Robbins, 1952) and even [START_REF] Thompson | On the likelihood that one unknown probability exceeds another in view of the evidence of two samples[END_REF][START_REF] Thompson | On the likelihood that one unknown probability exceeds another in view of the evidence of two samples[END_REF]. There are many variations in the definition of the problem; the contextual bandit framework as it is defined in [START_REF] Langford | The epoch-greedy algorithm for multi-armed bandits with side information[END_REF]) is as follows:

Let X be an arbitrary input space and A = {1..K} be a set of K actions. Let D be a distribution over tuples (x, r) with x ∈ X and r ∈ {0, 1} K a vector of rewards: in the (x, r) pair, the j th component of r is the reward associated to action a j , performed in the context x. In recommendation, the context is the information we have on the user (session, demographics, profile etc.) and an action is the recommendation of an item.

A contextual bandit game is an iterated game in which at each round t:

• (x t , r t ) is drawn from D.

• x t is provided to the player.

• The player chooses an action a t ∈ A based on x t and on the knowledge it gathered from the previous rounds of the game.

• The reward r t [a t ] is revealed to the player whose score is updated.

• The player updates his knowledge based on this new experience.

It is important to note the partial observability of this game: the reward is known only for the action performed by the player. This is what makes offline evaluation a complicated problem. A typical goal is to maximize the sum of reward obtained after T steps. To succeed a player has to learn about D and try to exploit that information. Therefore at time t, a player faces a classical exploration/exploitation dilemma: either perform an action he is uncertain about in order to improve his model of D (explore), either perform the action he believe to be optimal, although it may not be (exploit).

A simpler variant of this problem in which no contextual information is given, called the multi armed bandit problem (MAB) was extensively studied in the literature. One can for instance mention the UCB algorithm [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF] that optimistically deals with the dilemma by performing the action with higher upper confidence bound on the estimated reward expectation. The contextual bandit problem is less studied due to the additional complexity and additional assumptions entailed by the context space.

The most popular algorithm is without doubt LinUCB [START_REF] Li | A contextual-bandit approach to personalized news article recommendation[END_REF], although a few others exist such as epochgreedy [START_REF] Langford | The epoch-greedy algorithm for multi-armed bandits with side information[END_REF]. LinUCB is basically an extension of the classical UCB that uses the contexts under the assumption of normality and that X = R d . The reward expectation of an action is estimated via a linear regression on the context vectors and the confidence bound is computed using the dispersion matrix of the context vectors. These two state-of-the-art algorithms are the ones we will evaluate when running experiments.

Evaluation

We define a contextual algorithm A as taking as input an ordered list of (x, a, r) triplets (history) and outputting a policy π. A policy π maps X to A, that is chooses an action given a context. Note that we are also interested in evaluating policies. In our setting which is the most popular one, an algorithm is said optimal when maximizing the expectation of the sum of rewards after T steps:

G A (T, D) def = E D T t=1 r t [A(x t )].
For convenience, we define the per-trial payoff as the average click rate after T steps:

g A (T, D) def = G A (T ) T
Note that for a static algorithm (ie. that always outputs the same policy π), we have that:

∀T, g A (T, D) = g A (1, D) def = g A (D).
Note that from this point on, we will simplify the notations by systematically dropping D. g A (T ) is thus the quantity we wish to estimate as the measure of performance of a bandit algorithm.

In order to minimize the risks entailed by playing live a new algorithm, we are also interested in the quality of the estimation. Bootstrapping will enable us to estimate it. To do so we need additional notations. CT R A (T ) denotes the distribution of the per-trial payoff of A after T steps (so g A (T ) is its mean). Besides estimating g A (T ), our second goal is the computation of an estimator quality assessment ξ (CT R A (T )). Note that typically, ξ can be a quantile, a standard error, a bias or what we will consider here for simplicity: a confidence region around the mean of CT R A (T ) (aka g A ).

The time acceleration issue with replay methodologies

This section describes the replay methodology, that we call replay and that was introduced by [START_REF] Langford | Exploration scavenging[END_REF] and analyzed for the setting we consider by [START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF]. This section also highlights the method's limitations that we overcome in this paper and is crucial to understand the significance of our contribution.

First of all, as [START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF], we assume that we have a dataset S that was acquired online using an random uniform allocation of the actions for T steps. This data collection phase can be referred as exploration policy and is our unique information on D. This random decision making implies that any point in X × A has a non null probability to belong to S; this allows the evaluation of any policy. In a nutshell, the replay methodology on such a dataset works as follows: for each record (x t , a t , r t ) in S, the algorithm A is asked to choose an action given x t . If this action is a t , r t is revealed to A and taken into account to evaluate its performance. If the action is different, the record is discarded. This method is proved to be unbiased in some cases [START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF]. Note that the fact it needs the data to be acquired uniformly at random is quite restrictive. This problem is well studied and replay can be extended to allow the use of data acquired by a different but known logging policy at a cost of increased variance [START_REF] Langford | Exploration scavenging[END_REF]. Some work has been done to reduce this variance and even allow the use of data for which the logging policy is unknown [START_REF] Dudík | Doubly robust policy evaluation and learning[END_REF][START_REF] Strehl | Learning from logged implicit exploration data[END_REF]. Note also that if the evaluated bandit algorithm is close from the logging policy, we may even further reduce the variance [START_REF] Bottou | Counterfactual reasoning and learning systems[END_REF]. Finally there exist ways to adapt this method to the case where a list of items can be recommended [START_REF] Langford | Exploration scavenging[END_REF]. Although we do not take into account these considerations and keep the simplest assumption in this paper for clarity, our method is based on the same ideas as replay and could therefore be extended similarly as what is presented in the works we just cited.

Another issue with replay is well-known but not studied at all up to our knowledge. In average, only T K records are used. Therefore replay outputs an estimate ĝA 1 for a visual example of this problem. Note that in any situation except K = 1, CT R A T K = CT R A (T ), and the same thing goes for the confidence region ξ.

T K which follows the distribution CT R A T K of mean g A T K . It is important to have in mind that g A T K = g A (T ) if and only if T = +∞ or A is a static policy. See figure
One may argue that when evaluating a RS, plenty of data is available and therefore that T is almost infinite. Consequently one may also consider replay to be almost unbiased. This is true with the classical contextual bandit framework considered in the literature. With dynamic recommendation, the main application for this method, this could not be more wrong. Indeed, we argued that in this context, everything changes, especially the available items. For instance, in news recommendation a news remains active from a few hours to two days and its CTR systematically decreases over time [START_REF] Agarwal | Spatio-temporal models for estimating clickthrough rate[END_REF]. Moreover we also mentioned reasons to believe that the user tastes may change as well. Therefore when evaluating a contextual bandit algorithm, we want to evaluate its behavior against a very dynamic environment and in particular its ability to adapt to it. The use of replay in such a context is often justified by the fact that the environment can be considered static for small periods of time. This is not necessary but makes the understanding of our point easier. When an algorithm faces a "static" region of a dataset of size T i , when being replayed, it only has Ti K instead of T i steps to learn and exploit that knowledge. It is impossible to solve this problem by considering more data since new data would concern the next region, where different news with different CTR are available for recommendation, and potentially users with different tastes. In fact whatever assumptions we use to characterize how things evolve, using replay is equivalent to playing an algorithm with time going K times faster than in reality. This generates a huge bias. Note that it is most likely because of time acceleration that a non-contextual algorithm which looks a lot like UCB won a challenge evaluated by replay on the Yahoo! R6B dataset (Yahoo! Research, 2012) (news recommendation). See chapter 4 of Nicol [START_REF] Nicol | Data-driven evaluation of Contextual Bandit algorithms and applications to Dynamic Recommendation[END_REF] for more details.

As a conclusion, we consider in this work a classical contextual bandit framework with a fixed number of steps T . We assume that no more than T records can be acquired. Yet it is clear that if we manage to deal with this problem without adding data, we would also be able to deal as well with the problem of evaluating dynamic recommendation for which using more data may not be possible.

Bootstrapped Replay on Expanded Data

Now that the shortcoming of the replay method has been understood, we look for an other offline evaluation protocol that does not suffer from the time acceleration issue. The idea we propose stems from the idea of bootstrapping, introduced by [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF]. Thus let us remind the standard bootstrap approach. Basically, the idea is to compute the empirical distribution DT of an estimator θT computed on T observations. To do so, one only has access to a dataset S of size T . Therefore B new datasets S 1 , ...S B of size T are generated by making bootstrap resamples from S. A bootstrap resample is generated via drawing T samples with replacement. Note that this bootstrap procedure is a way to approximate D, the underlying distribution of the data. Computing θT on all the S i yields CT R(T ), an estimation of CT R(T ). From a theoretical point of view and under mild assumptions, ξ CT R(T ) converges with no bias at a speed in O(1/T ). This means that under a assumption of the concentration speed of a statistic we are able to estimate the confidence interval of the mean of the statistic much faster than its mean. Recall that ξ can be any measure of accuracy (defined in terms of bias, variance, confidence intervals, . . . ) over the statistic we want to study. Here we are interested in confidence intervals over CT R A (T ). The core idea of the evaluation protocol we propose in this paper is inspired by bootstrapping and inherits its theoretical properties. Using our notations, here is the description of this new method. From a dataset of size T with K possible choices of action at each step -we do not require K to be constant over time -,we generate B datasets of size K.T by sampling with replacement, following the non-parametric bootstrap procedure. Then for each dataset S b we use the classical replay method to compute an estimate ĝA (b) (T ). Therefore A is evaluated on T records on average. This step can be seen as a subsampling step that allows to return in the classical bootstrap setting. Thus note that it would not work for a purely deterministic policy, that for obvious reason would not take advantage of the data ex-Algorithm 1 Bootstrapped Replay on Expanded Data (BRED). We sketch this algorithm so that it looks very much like replay in [START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF]. Thus an algorithm is a function that maps an history of events h (b) to a policy π which itself maps contexts to actions. This makes the learning process of the algorithm appear clearly. A computationally efficient implementation would be slightly different. Notice also that for clarity, we only compute ĝA (T ) and omit CT R A (T ). Input

• A (contextual) bandit algorithm A • A dataset S of T triplets (x, a, r) • An integer B Output: An estimate of g A h (b) ← ∅, ∀b ∈ {1..B} /*empty history*/ G (b) A ← 0, ∀b ∈ {1..B} T (b) ← 0, ∀b ∈ {1..B} /* Bootstrap loop*/ for b ∈ {1..B} do /* estimation of CT R (b) A (T )*/ for i ∈ {1..T × K} do Sample with replacement (x, a, r) of S x ← JITTER(x) /*optional*/ π ← A(h (b) ) if π(x) = a then add (x, a, r) to h (b) G (b) A ← G (b) A + r T (b) ← T (b) + 1 end if end for end for return 1 B B b=1 G (b) A T (b)
pansion (an assumption in the formal analysis will reflect this fact). ĝA (T ) is given by averaging the ĝA (b) (T ). Together, the ĝA (b) (T ) are also an estimation of CT R A (T ), the distribution of the CTR of A after T interactions with D on which we can compute our estimator quality assessment ξ. More formally, the bootstrap estimator of the density of CT R A (T ) is estimated as follows:

CT R A (T )(x) = 1 B B b=1 I √ T ĝ(b) A (T ) -ĝA (T ) σA (T ) ≤ x
where σA (T ) is the empirical deviation obtained when computing the bootstrap estimates ĝA (b) (T ). The complete procedure, called Bootstrapped Replay on Expanded Data (BRED), is implemented in algorithm 1.

To complete the BRED procedure, one last detail is necessary. Each record of the original dataset S is contained K times in expectation in each expanded dataset S b . Therefore a learning algorithm may tend to overfit which would bias the estimator. To prevent this from happening, we introduce a small amount of Gaussian noise on the contexts. This technique is known as jittering and is well studied in the neural network field (Koistinen & Holmström, 1992).

The goal is the same, that is avoiding overfitting. In practice however it is slightly different as neural network are generally not learning online but on batches of data, each data being used several times during learning. In bootstrapping theory this technique is known as the smoothed bootstrap and was introduced by [START_REF] Silverman | The bootstrap: To smooth or not to smooth?[END_REF]. We mentioned that the bootstrap resampling is a way to approximate D. The smoothed bootstrap goes further by performing a kernel density estimation (KDE) of the data and sampling from it. Sampling from a KDE of the data where the kernel is Gaussian of bandwidth h is equivalent to sampling a record uniformly from S and applying a random noise sampled from N (0, h 2 ), which is what jittering does. The usual purpose of doing so in statistics is to get a less noisy empirical distribution for the estimator. Note that here we perform a partially smoothed bootstrap as we only apply a KDE on the part of D that generates the contexts.

Theoretical analysis

In this section, we make a theoretical analysis of our evaluation method BRED. The core loop in BRED is a bootstrap loop; henceforth, to complete this analysis, we first restate the theorem 1 which is a standard result of the bootstrap asymptotic analysis [START_REF] Kleiner | The big data bootstrap[END_REF]. Notice a small detail: each bootstrap step estimates a realization of CT R A (T ). The number of evaluations -which is also the number of non rejects -is a random variable denoted T (b) . Theorem 1. Suppose that:

• A is a recommendation algorithm which generates a fixed policy over time (this hypothesis can be weakened as discussed in remark 2),

• K items may be recommended at each time step,

• ξ (CT R A (T )) admits an expansion as an asymptotic series

ξ (CT R A (T )) = z + p 1 √ T + . . . + p α T α/2 + o 1 T α/2
where z is a constant independent of the distribution D (as defined in Sec. 2.1), and the p i are polynomials in the moments of CT R A (T ) under D (this hypothesis is discussed and explained in remark 1),

• The empirical estimator ξ CT R A (T ) admits a similar expansion:

z + p 1 √ T + . . . + p α T α/2 + O 1 T α/2 .
(1)

Then, for T ≤ T (b) × K and assuming finite first and second moments of CT R A (T ), with high probability:

ξ (CT R A (T )) -ξ CT R A (T ) = O Var( p (1) α -p α |D T ) √ T • B + O 1 T + O 1 T √ T (2)
where D T is the resampled distribution of D using T realizations.

Proof. it is actually a straightforward adaptation of the proof of theorem 3 of [START_REF] Kleiner | The big data bootstrap[END_REF]. Also note that this theorem is a reformulation of the bootstrap main convergence result as introduced by [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF]. Now, we use theorem 1 to bound the error made by BRED in the theorem 2.

Theorem 2. Assuming that

ξ (CT R A (T )) = z + p 1 √ T + . . . + p α T α/2 + o 1 T α/2
Then for algorithm A producing a fixed policy over time, BRED applied on a dataset of size T evaluates the expectation of the CT R A with no bias and with high probability for B and T large enough:

ξ (CT R A (T )) -ξ CT R A (T ) = O 1 T
This means that the convergence of the estimator of

ξ (CT R A (T )) is much faster than the convergence of the estimator of g A (T ) (which is in O(1/ √ T ).
This will allow a nice control of the risk that ĝA (T ) may be badly evaluated.

The sketch of the proof of theorem 2 is the following: first we prove that the replay strategy is able to estimate the moments of the distribution of CT R A fast enough with respect to T . The second step consists in using classical results from bootstrap theory to guarantee the unbiased convergence of the aggregation CT R A (T ) to the true distribution with an O( 1T ) speed. The rational behind this is that the gap introduced by the subsampling will be of the order of

O( 1 √ T B ).
Proof. At each iteration of the bootstrap loop (indexed by b), BRED is estimating the CTR using the replay method on a dataset of size T ′ = K × T . As the actions in S were chosen uniformly at random, we have

E(T (b) ) = T ′ /K = T .
As the policy is fixed, we can use the multiplicative Chernoff's bound as in [START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF] to obtain for all bootstrap step b:

Pr T (b) -T ≥ γ 1 T ≤ exp - T γ 2 1 3
for any γ 1 > 0 (where Pr(e) denotes the probability of event e). A similar inequality can be obtained with

E( G A ) = T g A : Pr G A -T g A ≥ γ 2 T g A ≤ exp - T g A γ 2 2 3 Thus with γ 1 = 3 T ln 4 δ and γ 2 = 3 T g A ln 4
δ using a union bound over probabilities, we have with probability at least 1δ:

1 -γ 1 ≤ T (b) T ≤ 1 + γ 1 g A 1 -γ 2 ≤ G A T ≤ g A 1 + γ 2 which implies G A T (b) -g A ≤ (γ 1 + γ 2 )g A 1 -γ 1 = O g a T ln 1 δ
So with high probability the first moment of CT R A (T (b) ) as estimated by the replay method admits an asymptotic expansion in 1/ E(T (b) ) = 1/T . Now we need to focus on higher order terms. All the moments are finite because the reward distribution over S is bounded. Recall that by hypothesis ξ (CT R A (T )) admits a α th order term:

p α = E D CT R A (T (b) ) α
The Chernoff's bound can be applied to

|(T (b) ) α -T α | and | G α A -T α g α A | leading to G α A (T (b) ) α -g α A = O g a T α 2 ln 1 δ With probability at least 1 -δ. So for a large enough T (b) , ξ CT R A (T (b)
) admits a expansion in polynomials of 1/ √ T . Thus theorem 1 applies and the aggregation of all the ĝ(b) A (T (b) ) allows an estimation of CT R A (T ). For a large enough number of bootstrap iterations (the value of B in BRED), we obtain a convergence speed in O(1/T ) with high probability, which concludes the proof.

After this analysis, we make two remarks about the assumptions that were needed to establish the theorems.

Remark 1:

The key point of the theorems is the existence of an asymptotic expansion of CT R A (T ) and CT R A (T ) in polynomials of 1/ √ T . This is a natural hypothesis for CT R A (T ) because the CTR is an average of bounded variables (probabilities of click). Note that the proof of theorem 2 shows that although T (b) is random the expansion remains valid anyway. For a contextual bandit algorithm A producing a fixed policy, the mean is going to concentrate according to the central limit theorem (CLT). Furthermore this hypothesis, omnipresent in bootstrap theory [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF], is for instance justified in econometrics by the fact that all the common estimators respect it [START_REF] Horowitz | The bootstrap[END_REF]. Yet this assumption is not verified for a static deterministic policy.

Remark 2 Let us consider algorithms that produce a policy which changes over time (a learning algorithm in particular). After a sufficient amount of recommendations, a reasonable enough algorithm will produce a policy that will not change any longer (if the world is stationary). Thus again, the CLT will apply and we will observe a convergence of ĝA (T ) to its limit in 1/ √ T . Nevertheless noth- here) is actually efficient to avoid overfitting issues. The rather small error rates for very small datasets are due to the fact that on too small datasets all the recommender algorithm tend to make random choices which are not very hard to evaluate.

ing holds true here when the algorithm is actually learning. This is due to the fact that the Chernoff bound no longer applies as the steps are not independent. However the behavior of classical learning algorithms are smooth, especially when randomized (see [START_REF] Auer | The nonstochastic multiarmed bandit problem[END_REF] for an example of a randomized version of UCB). [START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF] argue that in this case convergence bounds may be derived for replay (which then could be applied to BRED) at the cost of a much more complicated analysis including smoothness assumptions. For non reasonable algorithms and thus in the general case, no guarantees can be provided. By the way note that a very intuitive way to justify Jittering is to consider that it helps the Chernoff bound being "more true" in the case of a learning algorithm.

Experiments in realistic settings

As we proved that BRED has promising theoretical guarantees in the setting introduced in [START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF], let us now compare its empirical performance to that of the replay method

Synthetic data and discussing Jittering

The first set of experiments was run on synthetic data. Indeed, we needed to be able to compare the errors of estimation of the two methods on various fixed size datasets relatively to the ground truth: an evaluation against the model itself. Before going any further, let us describe the model we used. It is a linear model with Gaussian noise (as in [START_REF] Li | A contextual-bandit approach to personalized news article recommendation[END_REF]) and was built as follows:

• a fixed action set (or news set) of size K = 10.

• The context space is X = R 15 . Each context x is generated as a sum c + n where c ∼ N (0, 1) and n ∼ N (0, 1 2 ).

• The CTR of a news i displayed in a context x is given linearly by q i +w T i c. Note the non-contextual element q i and that the noise n is ignored.

• Finally there are two kinds of news: (i) 4 "universal" news that are interesting in general like Obama is reelected and for which q i ∼ U (0.4, 0.5) is high and w i = 0 15 . (ii) 6 specific news like New Linux distribution released for which q i ∼ U (0.1, 0.2) is low and w i consists of zeros except for a number m of relevant weights sampled from N (0, 1 5 ).

A non contextual approach would perform decently by quickly learning the q i values. Yet LinUCB [START_REF] Li | A contextual-bandit approach to personalized news article recommendation[END_REF], a contextual bandit algorithm will do better by learning when to recommend the specific news. Figure 2 displays the results and interpretation of an experiment which consists in evaluating LinUCB(α = 1) using the different methods. It is clear that BRED converges much faster than the replay method.

Remark: As it can be seen on Figure 2, jittering is very important to obtain good performance when evaluating a learning algorithm. Empirically, a good choice for the level of jitter seems to be a function in O 1 √ T , with T the size of the dataset. Note that this is proportional to the standard distribution of the posterior of the data. The results confirm our intuition: jittering is very important when the dataset is small but gets less and less necessary as the dataset grows.

Real data

Adapting replay to a real life dataset, corresponding to dynamic recommendation is straightforward although it leads to biased estimations. BRED really needs the assumption of a static world in order to perform the bootstrap resamples. Therefore BRED needs to be run on successive windows of the dataset on which a static assumption can be made. This creates a bias/variance trade-off: if the windows are too big, some dynamic properties of the world may be erased (bias). On the contrary, too smalls window will lead to very variate bootstrap estimates. To simplify things, we ran experiments assuming a static world on small portions of the Yahoo! R6B dataset. We actually took the smallest number of portions such that a given portion has a fixed pool of news (≈ 630 portions). This experiment is similar to what is done in [START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF]: the authors measured the error of the estimated CTR of UCB (α = 1) by the replay method on datasets of various sizes relatively to what they call the ground truth: an evaluation of the same algorithm on a real fraction of the audience. As we obviously cannot do that, we used a simple trick. For each portion i of size T i with K i news, we computed an estimation of the ground truth g A Ti Ki

by averaging the estimated CTR of UCB using the replay method on 100 random permutations of the data. For each portion the experiment consists in subsampling T i /K i records and evaluating UCB using replay and BRED on this smaller dataset to estimate the ground truth using less data, faking time acceleration. The results and interpretation are shown on Figure 3: the better accuracy of BRED is very clearly illustrated.

Conclusion

In this paper, we studied the problem of recommendation system evaluation, sticking to a realistic setting: we focused on obtaining a methodology for practical offline evaluation, providing a good estimate using a reasonable amount of data. Previous methods are proved to be asymptotically unbiased with a low speed of convergence on a . The difference between the estimated CTR and the actual one on some batches extracted from the Yahoo! R6B dataset for a UCB. Batches are build as explained in section 6.2. The closer to 0, the better. Please note that the replay method tends to under-estimate the true CTR for small batches. This is due to the fact that UCB does not have enough time to reach its actual CTR. static dataset, but yield counter-intuitive estimates of performance on real datasets. Here, we introduce BRED, a method with a much faster speed of convergence on static datasets (at the cost of loosing unbiasedness) which allows it to be much more accurate on dynamic data. Experiments demonstrated our point; they were performed on a publicly available dataset made from Yahoo! server logs and on synthetic data presenting the time acceleration issue. This paper was also meant to highlight the time acceleration issue and the misleading results given by a careless evaluation of an algorithm. Finally our method comes with a very desirable property in a context of minimizing the risks entailed by putting online a new RS: an extremely accurate estimation of the variability of the estimator it provides.

An interesting line of future work is the automatic selection of the Jittering bandwidth. Note that this problem is extensively studied in the context of KDE [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF].

A possible extension of this work is to use BRED to build a "safe controller". Indeed, when a company uses a recommendation system that behaves according to a certain policy π that reaches a certain level of performance, the hope is that when changing the recommendation algorithm, the performance will not drop. As an extension of the work presented here, it is possible to collect some data using the current policy π, compute small variations of π with tight confidence intervals over their CTR and then replace the current policy π with the improved one. This may be seen as a kind of "gradient" ascent of the CTR in the space of policies.

Figure 2 .

 2 Figure 2. Mean absolute error of the CTR of LinUCB estimated by different methodologies. Conducted on artificial datasets as described by section 6.1. The lower, the better. Jittering (h = 50 √ T

  Figure3. The difference between the estimated CTR and the actual one on some batches extracted from the Yahoo! R6B dataset for a UCB. Batches are build as explained in section 6.2. The closer to 0, the better. Please note that the replay method tends to under-estimate the true CTR for small batches. This is due to the fact that UCB does not have enough time to reach its actual CTR.

Supplementary material

The detailed implementation of replay using our notations is given in algorithm 2. Note that apart from notations, no modification are made with regards to the original presentation in [START_REF] Langford | Exploration scavenging[END_REF][START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF]. Figure 4 is the same experiment as in section 6.1 but with a noncontextual algorithm UCB: this plot exhibits both the importance of Jittering and the improvement brought by our method compared to the state of the art. Algorithm 2 Replay method [START_REF] Langford | Exploration scavenging[END_REF][START_REF] Li | Unbiased offline evaluation of contextual-bandit-based news article recommendation algorithms[END_REF]. Remark: for the sake of the precision of the specification of the algorithm, we use a history h which is the list of triplets (x, a, r) that have yet been used to estimate the performance of the algorithm A. The goal is to avoid hiding internal information maintenance in A; a real implementation may be significantly different for the sake of efficiency, by learning incrementally. Input:

• A contextual bandit algorithm A

• A set S of L triplets (x, a, r)

Output: An estimate of g A h ← ∅ G A ← 0 T ← 0 for t ∈ {1..L} do Get the t-th element (x, a, r) of S π ← A(h) if π(x) = a then add (x, a, r) to h G A ← G A + r T ← T + 1 end if end for return G A T