Post-Contextual-Bandit Inference
Year: 2021

Abstract:
Contextual bandit algorithms are increasingly replacing non-adaptive A/B tests in e-commerce, healthcare, and policymaking because they can both improve outcomes for study participants and increase the chance of identifying good or even best policies. To support credible inference on novel interventions at the end of the study, nonetheless, we still want to construct valid confidence intervals on average treatment effects, subgroup effects, or value of new policies. The adaptive nature of the data collected by contextual bandit algorithms, however, makes this difficult: standard estimators are no longer asymptotically normally distributed and classic confidence intervals fail to provide correct coverage. While this has been addressed in non-contextual settings by using stabilized estimators, variancestabilized estimators in the contextual setting pose unique challenges that we tackle for the first time in this paper. We propose the Contextual Adaptive Doubly Robust (CADR) estimator, a novel estimator for policy value that is asymptotically normal under contextual adaptive data collection. The main technical challenge in constructing CADR is designing adaptive and consistent conditional standard deviation estimators for stabilization. Extensive numerical experiments using 57 OpenML datasets demonstrate that confidence intervals based on CADR uniquely provide correct coverage.

heading: 1 Introduction
text: Contextual bandits, where personalized decisions are made sequentially and simultaneously with data collection, are increasingly used to address important decision-making problems where data is limited and/or expensive to collect, with applications in product recommendation [Li et al., 2010], revenue management [Kallus and Udell, 2020, Qiang and Bayati, 2016], and personalized medicine [Tewari and Murphy, 2017]. Adaptive experiments, whether based on bandit algorithms or Bayesian optimization, are increasingly being considered in place of classic randomized trials in order to improve both the outcomes for study participants and the chance of identifying the best treatment allocations [Athey et al., 2018, Quinn et al., 2019, Kasy and Sautmann, 2021, Bakshy et al., 2018].
But, at the end of the study, we still want to construct valid confidence intervals on average treatment effects, subgroup effects, or the value of new personalized interventions. Such confidence intervals are, for example, crucial for enabling credible inference on the presence or absence of improvement of novel policies. However, due to the adaptive nature of the data collection, unlike classic randomized trials, standard estimates and their confidence intervals actually fail to provide correct coverage, that is, contain the true parameter with the desired confidence probability (e.g., 95%). A variety of recent work has recognized this and offered remedies [Hadad et al., 2019, Luedtke and van der Laan, 2016], but only for the case of non-contextual adaptive data collection. Like classic
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
confidence intervals, when data comes from a contextual bandit – or any other context-dependent adaptive data collection – these intervals also fail to provide correct coverage. In this paper, we propose a novel asymptotically normal estimator for the value of a (possibly contextual) policy from context-dependent adaptively collected data. This asymptotic normality leads directly to the construction of valid confidence intervals.
Our estimator takes the form of a stabilized doubly robust estimator, that is, a weighted time average of an estimate of the so-called canonical gradient using plug in estimators for the outcome model, where each time point is inversely weighted by its estimated conditional standard deviation given the past. We term this the Contextual Adaptive Doubly Robust (CADR) estimator. We show that, given consistent conditional variance estimates which at each time point only depend on previous data, the CADR estimator is asymptotically normal, and as a result we can easily construct asymptotically valid confidence intervals. This normality is in fact robust to misspecifying the outcome model. A significant technical challenge is actually constructing such variance estimators. We resolve this using an adaptive variance estimator based on the importance-sampling ratio of current to past (adaptive) policies at each time point. We also show that we can reliably estimate outcome models from the adaptively-collected data so that we can plug them in. Extensive experiments using 57 OpenML datasets demonstrate the failure of previous approaches and the success of ours at constructing confidence intervals with correct coverage.

heading: 1.1 Problem Statement and Notation
text: The data. Our data consists of a sequence of observations indexed t = 1, . . . , T comprising of context X(t) ∈ X , action A(t) ∈ A, and outcome Y (t) ∈ Y ⊂ R generated by an adaptive experiment, such as a contextual bandit algorithm. Roughly, at each round t = 1, 2, . . . , T , an agent formed a contextual policy gt(a | x) based on all past observations, then observed an independently drawn context vector X(t) ∼ Q0,X , carried out an action A(t) drawn from its current policy gt(· | X(t)), and observed an outcome Y (t) ∼ Q0,Y (· | A(t), X(t)) depending only on the present context and action. The action and context measurable spaces X ,A are arbitrary, e.g., finite or continuous. We assume outcomes are bounded and, without loss of generality, take Y = [−1, 1]. More formally, we let O(t) := (X(t), A(t), Y (t)) and make the following assumptions about the sequence O(1), . . . , O(T ) comprising our dataset. First, we assume X(t) is independent of all else given A(t) and has a time-independent marginal distribution that we denote by Q0,X . Second, we assume A(t) is independent of all else given O(1), . . . , O(t − 1), X(t) and we set gt(· | X(t)) to its (random) conditional distribution given O(1), . . . , O(t − 1), X(t). Third, we assume Y (t) is independent of all else given X(t), A(t) and has a time-independent conditional distribution given X(t) = x,A(t) = a that is denoted by Q0,Y (· | A,X). The distributions Q0,X and Q0,Y are unknown, while the policies gt(a | x) are known, as would be the case when running an adaptive experiment. To simplify presentation we endow A with a base measure µA (e.g., counting for finite actions or Lebesgue for continuous actions) and identify policies gt with conditional densities with respect to (w.r.t.) µA. In the case of K <∞ actions, policies are maps from X to the K-simplex. Note that, as the agent updates its policy based on already collected observations, gt is a random O(1), . . . , O(t − 1)-measurable object. This is the major departure from the setting considered in other literature on off-policy evaluation, which only consider a fixed logging policy, gt = g, that is independent of the data. See Section 1.2.
The target parameter. We are interested in inference on a generalized average causal effect expressed as a functional of the unknown distributions above, Ψ0 = Ψ(Q0,X , Q0,Y ), where for any distributions QX , QY , we define
Ψ(QX , QY ) := ∫ yQX(dx)g ∗(a | x)dµA(a)QY (dy | a, x),
where g∗(a | x) : A × X → [−G,G] is a given fixed, bounded function. Two examples are: (a) when g∗ is a policy (conditional density), then Ψ0 is its value; (b) when g∗ is the difference between two policies then Ψ0 is the difference between their values. A prominent example of the latter is when A = {+1,−1} and g∗(a | x) = a, which is known as the average treatment effect. If we include an indicator for x being in some set, then we get the subgroup effect. Defining the conditional mean outcome,
Q̄0(a, x) := EQ0,Y (·|x,a)[Y ] = ∫ yQ0,Y (dy | a, x),
we note that the target parameter only depends on Q0,Y via Q̄0, so we also overload notation and write Ψ(QX , Q̄) = ∫ Q̄(a, x)QX(dx)g(a | x)dµA(a) for any function Q̄ : A×X → Y . Note that when |A| <∞ and µA is the counting measure, the integral over a is a simple sum.
Canonical gradient. We will make repeated use of the following function: for any conditional density (a, x) 7→ g(a | x), any probability distribution QX over the context space X , and any function Q̄ : A×X , we define the function D′(g, Q̄) : O → R by
D′(g, Q̄)(x, a, y) := g∗(a | x) g(a | x)
(y − Q̄(a, x)) + ∫ Q̄(a′, x)g∗(a′ | x)dµA(a′).
Further, define D(g,QX , Q̄) = D′(g,QX , Q̄) − Ψ(QX , Q̄), which coincides with the so-called canonical gradient of the target parameter Ψ w.r.t. the usual nonparametric statistical model comprising all joint distributions over O [van der Vaart, 2000, van der Laan and Robins, 2003].
Integration operator notation. For any policy g and distributions QX , QY , denote by PQ,g the induced distribution on O := X × A × Y . For any function f : O → R, we use the integration operator notation
PQ,gf = ∫ f(x, a, y)QX(dx)g(a | x)dµA(a)QY (dy | a, x),
that is, the expectation w.r.t. PQ,g alone. Then, for example, for anyO(1), . . . , O(s−1)-measurable random function f : O → R, we have that PQ0,gsf = EQ0,gs [f(O(s)) | O(1), . . . , O(s− 1)].

heading: 1.2 Related Literature and Challenges for Post-Contextual-Bandit Inference
text: Off-policy evaluation. In non-adaptive settings, where gt = g is fixed and does not depend on previous observations, common off-the shelf estimators for the mean outcome under g∗ include the Inverse Propensity Scoring (IPS) estimator [Beygelzimer and Langford, 2009, Li et al., 2011] and and the Doubly Robust (DR) estimator [Dudı́k et al., 2011, Robins et al., 1994]:
Ψ̂IPS := 1
T T∑ t=1 D′(g, 0), Ψ̂DR := 1 T T∑ t=1 D′(g, ̂̄Q) where ̂̄Q is an estimator of the outcome model Q̄0(a, x). If we use cross-fitting to estimate ̂̄Q [Chernozhukov et al., 2018], then both the IPS and DR estimators are unbiased and asymptotically normal, permitting straightforward inference using Wald confidence intervals (i.e., ±1.96 of the estimated standard error). There also exist many variants of the IPS and DR estimators that, rather than plugging in the importance sampling (IS) ratios (g∗/gt)(A(t) | X(t)) and/or outcome-model estimators, instead choose them directly with the aim to minimize error [e.g. Kallus, 2018, Farajtabar et al., 2018, Thomas and Brunskill, 2016, Wang et al., 2017, Kallus and Uehara, 2019].
Inference challenges in adaptive settings. In the adaptive setting, it is easy to see that, if in the tth term for DR we use an outcome model ̂̄Qt−1 fit using only the observations O(1), . . . , O(t− 1), then both the IPS and DR estimators both remain unbiased. However, neither generally converges to a normal distribution. One key difference between the non-adaptive and adaptive settings is that the IS ratios (g∗/gt)(A(t) | X(t)) can both diverge to infinity or converge to zero. As a result of this, the above two estimators may either be dominated by their first terms or their last terms. At a more theoretical level, this violates the classical condition of martingale central limit theorems that the conditional variance of the terms given previous observations stabilizes asymptotically.
Stabilized estimators. The issue for inference due to instability of the DR estimator terms was recognized by Luedtke and van der Laan [2016] in another setting. They work in the non-adaptive setting but consider the problem of inferring the maximum mean outcome over all policies when the optimal policy is non-unique. Their proposal is a so-called stabilized estimator, in which each term is inversely weighted by an estimate of its conditional standard deviation given the previous terms.
This stabilization trick has been also been reused for off-policy inference from non-contextual bandit data by Hadad et al. [2019], as the stabilized estimator remains asymptotically normal, permitting inference. In their non-contextual setting, an estimate of the conditional standard deviation of the terms can easily be obtained by the inverse square root propensities. In contrast, in our contextual setting, obtaining valid stabilization weights is more challenging and requires a construction involving adaptive training on past data. Zhang et al. [2021] use similar weights to Hadad et al. [2019] in order to stabilize an M -estimator fit to adaptively collected data for the purpose of inference on a true model for Q̄0. While they do consider a contextual setting, their method requires correct specification, else the inferential target depends on the logging policies. We focus only on off-policy evaluation but do not require well-specified models.
Concurrent work of Zhan et al. [2021]. The paper of Zhan et al. [2021], which was published after the submission of our paper, is also motivated by establishing an unbiased and asymptotically normal estimator from data collected by a contextual bandit. Both our method and theirs rely on adaptively weighted DR scores and yield, under certain conditions, asymptotically normal estimators. However, the way they construct their adaptive weights is quite different from ours. In particular, their adaptive weights at time t (a) do not incorporate an outcome model and do not rely on it and (b) depend only on the known propensity score, which is plugged-in. As a result of these different constructions, our analysis and theirs require different sets of conditions. While we require slightly more stringent conditions on the exploration rate and the complexity of the logging policy class than they do (our assumptions 5 and 6), they require some form of asymptotic stabilization of the logging policies (their conditions (8) and (12)). An instance where the logging policy does not stabilize is, for example, when various actions have the same expected reward on certain subsets of the context space. Our approach in this paper is to use the variance-stabilization weights for the contextual bandit setting, just like Hadad et al. [2019] do the in the non-contextual bandit setting and Luedtke and van der Laan [2016] do in an iid setting with multiple optimal policies. Although Zhan et al. [2021] also tackle the contextual bandit setting, they do not use the variance-stabilization weights; correspondingly their conditions and results apply in different, complementary settings. The key challenge we tackle in applying variance-stabilization weights in the contextual bandit setting is dealing with the need to use sequentially estimated weights.

heading: 1.3 Contributions
text: In this paper, we construct and analyze a stabilized estimator for policy evaluation from contextdependent adaptively collected data, such as the result of running a contextual bandit algorithm. This then immediately enables inference. After constructing a generic extension of the stabilization trick, the main technical challenge is to construct a sequence of estimators σ̂1, . . . , σ̂T of the conditional standard deviations that are both consistent and such that for each t, σ̂t only uses the previous data points O(1), . . . , O(t−1). We show in extensive experiments across a large set of contextual bandit environments that our confidence intervals uniquely achieve close to nominal coverage.

heading: 2 Construction and Analysis of the Generic Contextual Stabilized Estimator
text: In this section, we give a generic construction of a stabilized estimator in our contextual and adaptive setting. That is, given generic plug-ins for outcome model and conditional standard deviation. We then provide conditions under which the estimator is asymptotically normal, as desired. To develop CADR, we will then proceed to construct appropriate plug in estimators in the proceeding sections.

heading: 2.1 Construction of the Estimator
text: Outcome and variance estimators. Our estimator uses a sequence ( ̂̄Qt)t≥1 of estimators of the outcome model Q̄0, such that, for every t, ̂̄Qt is O(1), . . . , O(t)-measurable, that is, is trained using only the data up to time t. A key part of our estimator is the conditional variance estimators. Also, we require estimates of the conditional standard deviation of the canonical gradient. Define:
σ0,t := σ0,t(gt), where σ20,t(g) := VarQ0,g ( D′(g, ̂̄Qt−1)(O(t)) | O(1), . . . , O(t− 1)) .
Let (σ̂t)t≥1 be a given sequence of estimates of σ0,t such that σ̂t is O(1), . . . , O(t− 1)-measurable, that is, is estimated using only the data up to time t.
The generic form of the estimator. The generic contextual stabilized estimator is then defined as:
Ψ̂T :=
( 1
T T∑ t=1 σ̂−1t
)−1 1
T T∑ t=1 σ̂−1t D ′(g, ̂̄Qt−1). (1)

heading: 2.2 Asymptotic normality guarantees
text: We next characterize the asymptotic distribution of Ψ̂T under some assumptions. Assumption 1 (Non degenerate efficiency bound). infg PQ0,gD2(g, Q̄0, Q0,X) > 0.
Assumption 1 states that there is no fixed logging policy g such that the efficiency bound for estimation of Ψ(Q̄0, Q0,X) in the nonparametric model, from i.i.d. draws of PQ0,g , is zero. If assumption 1 does not hold, there exists a logging policy g such that, if O = (X,A, Y ) ∼ PQ0,g , then (g∗(A | X)/g(A | X))Y equals Ψ(Q̄0, Q0,X) with probability 1. In other words, if assumption 1 does not hold, there exists a logging policy g such that Ψ(Q̄0, Q0,X) can be estimated with no error with probability 1 from a single draw of PQ0,g . Thus, it is very lax. An easy sufficient condition for Assumption 1 is that the outcome model has nontrivial variance in that VarQ0,X ( ∫ Q̄(a,X)g∗(a | x)dµA(a)) > 0.
Assumption 2 (Consistent standard deviation estimators.). σ̂t − σ0,t t→∞−−−→ 0 almost surely.
In the next section we will proceed to construct specific estimators σ̂t that satisfy Assumption 2, leading to our proposed CADR estimator and confidence intervals.
Assumption 3 (Exploration rate). ‖ supa∈A,x∈X g−1t (a | x)‖∞ = ω(t−1/2).
Assumption 3 requires that the exploration rate of the adaptive experiment does not decay too quickly. Crucially, Assumption 3 requires that we have nontrivial propensities of pulling any one arm in any one context, conditioned on the history. While this may be satisfied by inherently randomized algorithms such as -greedy and Thompson sampling, this is not satisfied by algorithms such as UCB, which are deterministic once given the history. This limitation is the same as for Hadad et al. [2019], Zhan et al. [2021]. Remark 1. One can weaken Assumption 3 assumption to just requiring ω(t−1) decay by instead slightly strengthening Assumption 1 to require that σ20,t & ‖ supa∈A,x∈X g−1t (a | x)‖∞, which would be implied by g∗ placing some mass on actions with decaying propensity and Y having everywhere-nontrivial conditional variance. Here at & bt means that for some constant c > 0, we have at ≥ cbt for all t ≥ 1. However, as we will in fact strengthen Assumption 3 anyway in Assumption 6 below, which we will use in constructing estimators σ̂t that satisfy Assumption 2, we do not discuss this potential weakening of Assumption 3.
Based on these assumptions, we have the following asymptotic normality result: Theorem 1. Denote ΓT := ( T−1 ∑T t=1 σ̂ −1 t )−1 . Under Assumptions 1 to 3, it holds that
Γ−1T √ T ( Ψ̂T −Ψ0 ) d−→ N (0, 1).
Remark 2. Theorem 1 does not require the outcome model estimator to converge at all. As we will see in Section 3, our conditional variance estimator does require that the outcome model converges to a fixed limit Q̄1, but this limit does not have to be the true outcome model Q̄0. In other words, consistency of the outcome model is not required at any point of our analysis.

heading: 3 Construction of the Conditional Variance Estimator and CADR
text: We now tackle the construction of σ̂t satisfying our assumptions; namely, they must be adaptively trained only on past data at each t and they must be consistent. Observe that σ20,t = σ 2 0(gt, ̂̄Qt−1),
Algorithm 1 The CADR Estimator and Confidence Interval Input: Data O(1), . . . , O(T ), policies g1, . . . , gT , target g∗, outcome regression estimator for t = 1, 2, . . . , T do
Train ̂̄Qt−1 on O(1), . . . , O(t− 1) using the outcome regression estimator Set D′t,s = D(gs,
̂̄Qt−1)(O(t)) for s = t, . . . , T // (note index order compared to next line) Set σ̂2t = 1 t−1 ∑t−1 s=1 gt(A(s)|X(s)) gs(A(s)|X(s)) (D ′ s,t) 2 − ( 1 t−1 ∑t−1 s=1 gt(A(s)|X(s)) gs(A(s)|X(s))D ′ s,t
)2 end for Set ΓT = ( 1 T ∑T t=1 σ̂ −1 t
)−1 Return estimate Ψ̂T = ΓTT ∑T t=1 σ̂ −1 t D ′ t,t and confidence intervals CIα = [Ψ̂T ±ζ1−α/2ΓT / √ T ]
where we define σ20(g, Q̄) := Φ0,2(g, Q̄)− (Φ0,1(g, Q̄))2, Φ0,i(g, Q̄) := PQ0,g(D′)i(g, Q̄), i = 1, 2.
Designing an O(1), . . . , O(t − 1)-measurable estimator of σ20,t presents several challenges. First, while we can only use observations O(1), . . . , O(t − 1) to estimate it, σ20,t is defined as a function of integrals w.r.t. PQ0,gt , from which we have only one observation, namely O(t). Second, our estimation target σ20,t = σ0(gt,
̂̄Qt−1) is random as it depends on gt and ̂̄Qt. Third, gt, ̂̄Qt depend on the same observations O(1), . . . , O(t− 1) that we have at our disposal to estimate σ20,t.
Representation via importance sampling. We can overcome the first difficulty via importance sampling, which allows us to write Φ0,i(g, Q̄), i = 1, 2 as integrals w.r.t. PQ0,gs , s = 1, . . . , t − 1, i.e., the conditional distributions of observations O(s), s = 1, . . . , t− 1 given their respective past. Namely, for any s ≥ 1, i = 1, 2, we have that
Φ0,i(g, Q̄) = PQ0,gs g
gs (D′)i(g, Q̄). (2)
Dealing with the randomness of the estimation target. We now turn to second challenge. Since σ20,t can be written in terms of Φ0,i(gt, ̂̄Qt−1) for i = 1, 2, Eq. (2) suggests perhaps an approach based on sample averages of (gt/gs)(D′)i(gt, ̂̄Qt−1) over s. However, whenever s < t, the latter is an O(1), . . . , O(t − 1)-measurable function due to the dependence on gt and ̂̄Qt. Namely, PQ0,gs{(gt/gs)(D′)i(gt, ̂̄Qt−1)} does not coincide in general with the conditional expectation EQ0,gs [((gt/gs)(D ′)i(gt, ̂̄Qt−1))(O(s)) | Ō(s − 1)], as would arise from a sample average. We
now look at solutions to overcome this difficulty, considering first ̂̄Qt−1 and then gt. Dealing with the randomness of ̂̄Qt−1. We propose an estimator of σ20(g, ̂̄Qt−1) for any fixed g. While requiring that ̂̄Qt−1 converges to the true outcome regression function Q̄0 is a strong requirement, most reasonable estimators will at least converge to some fixed limit Q̄1. As a result, under an appropriate stochastic convergence condition on ( ̂̄Qt−1)t≥1, Φ0,i(g, ̂̄Qt−1) can be reasonably approximated by the corresponding Cesaro averages, defined for i = 1, 2 as
Φ̄0,i,t(g) := 1
t− 1 t−1∑ s=1 Φ0,i(g, ̂̄Qs−1) = 1t− 1 t∑ s=1 EQ0,gs [ ((g/gs)(D ′)i(g, ̂̄Qs−1))(O(s)) | Ō(s− 1)] . These are easy to estimate from the corresponding sample averages, defined for i = 1, 2 as
Φ̂i,t(g) := 1
t− 1 t∑ s=1 ((g/gs)(D ′)i(g, ̂̄Qs−1))(O(s)),
since for each i = 1, 2, the difference Φ̂i,t(g) − Φ̄0,i,t(g) is the average of a martingale difference sequence (MDS). We then define our estimator of σ20(g, ̂̄Qt−1) as σ̂t(g) := Φ̂2,t(g)− (Φ̂1,t(g))2. (3)
From fixed g to random gt. So far, we have proposed and justified the construction of σ̂t(g) as an estimator of σ0,t(g, ̂̄Qt−1) for a fixed g. We now discuss conditions under which σ̂t(gt) is valid estimator of σ0,t(gt, ̂̄Qt−1). When g is fixed, for each i = 1, 2, the error Φ̂i,t(g) − Φ0,i,t(g, ̂̄Qt−1) decomposes as the sum of the MDS average Φ̂i,t(g) − Φ̄0,i,t(g) and of the Cesaro approximation error Φ̄0,i,t(g)−Φ0,i(g, ̂̄Qt−1). Both differences are straightforward to bound. For a random gt, the term Φ̂i,t(gt)− Φ̄0,i,t(gt) is no longer an MDS average. Fortunately, under a complexity condition on the logging policy class G, we can bound the supremum of the martingale empirical processes {|Φ̂i,t(g)− Φ̄0,i,t(g)| : g ∈ G}, which in turn gives us a bound on |Φ̂i,t(gt)− Φ̄0,i,t(gt)|.
Consistency guarantee for σ̂2t . Our formal consistency result relies on the following assumptions. Assumption 4 (Outcome regression estimator convergence). There exists β > 0, and a fixed function Q̄1 : A×X → R such that ‖ ̂̄Qt − Q̄1‖1,Q0,X ,g∗ = O(t−β) almost surely. The next assumption is a bound on the bracketing entropy (see, e.g., [van der Vaart and Wellner, 1996] for definition) of the logging policy class. Assumption 5 (Complexity of the logging policy class). There exists a fixed reference function gref(a | x), a class of conditional densities G, and some p > 0 such that gt ∈ G ∀t ≥ 1 almost surely, supg∈G ‖g/gref‖∞ ≤ G, ‖g∗/gref‖∞ ≤ G, and
logN[ ]( ,G/gref , ‖ · ‖2,Q0,X ,gref ) . −p,
where G/gref := {g/gref : g ∈ G}.
Next, we require a condition on the exploration rate that is stronger than Assumption 3.
Assumption 6 (Exploration rate (stronger)). ‖ supa∈A,x∈X gref(a | x)/gt(a | x)‖∞ . tα(β,p), where α(β, p) := min(1/(3 + p), 1/(1 + 2p), β).
Theorem 2. Suppose that Assumptions 4 to 6 hold. Then, σ̂2t − σ20,t = o(1) almost surely. Remark 3. While we theoretically require the existence of a logging policy class G with controlled complexity, we do not actually need to know G to construct our estimator. Moreover, while we require a bound on the bracketing entropy of the logging policy class G, we impose no restriction on the outcome regression model complexity, permitting us to use flexible black-box regression methods. Usually we simply take gref(a | x) = 1 so the boundedness in Assumption 5 is simply boundedness of densities, but the above allows for slightly more flexibility in cases where we do not necessarily explore all actions equally well in Assumption 6.
Remark 4. Assumption 4 requires ( ̂̄Qt) to be a sequence of regression estimator, such that for every t ≥ 1, ̂̄Qt is fitted on O(1), . . . , O(t) and for which we can guarantee a rate of convergence to some fixed limit Q̄1. Note that this can at first glance pose a challenge since observations O(1), . . . , O(t) are adaptively collected. In the appendix, we give guarantees for outcome regression estimation over a nonparametric model using an importance sampling weighted empirical risk minimization.
Remark 5. Assumptions 5 and 6 put restrictions on the logging policy, that is, the contextual bandit algorithm that was run. Assumption 5 is not very restrictive, especially for large p. It is generally satisfied for p arbitrarily close to 0 for -greedy with outcome regressions in a VC-subgraph class such as trees or for Thompson sampling with parametric models. For larger p it can accommodate very general policy classes, even non-Donsker classes (p > 2). Assumption 6, on the other hand, may be restrictive as it can require a lot of exploration. While it appears that at least ω(1/t) exploration might be necessary for some normality as in Theorem 1 (see Remark 1), the rates needed in Assumption 6 may be too strong than actually needed. We therefore view it as a high-level assumption, as it does not change the method itself and it may not need to be satisfied in practice, and we leave it as an open question how to tighten our analysis to be able to rely on less exploration.
CADR asymptotics. Our proposed CADR estimator is now given by plugging our estimates σ̂t from Eq. (3) into Eq. (1), as summarized in Algorithm 1 As an immediate corollary of Theorems 1 and 2 we have our main guarantee for this final estimator, showing CADR is asymptotically normal, whence we immediately obtain asymptotically valid confidence intervals.
Corollary 1 (CADR Asymptotics and Inference). Suppose that Assumptions 1 and 4 to 6 hold. Let σ̂t be given as in Eq. (3). Denote ΓT := ( T−1 ∑T t=1 σ̂ −1 t )−1 . Then,
Γ−1T √ T ( Ψ̂T −Ψ0 )
d−→ N (0, 1). Moreover, letting ζα denote the α-quantile of the standard normal distribution,
Pr [ Ψ(Q0,X , Q̄0) ∈ [ Ψ̂T ± ζ1−α/2ΓT / √ T ]] T→∞−−−−→ 1− α.

heading: 4 Empirical Evaluation
text: In this section and in sections E.2-E.5 of the appendix, we present an extensive and reproducible1 empirical evaluation on public datasets that demonstrates the robustness of CADR confidence intervals using contextual bandit data with comparison to several baselines.
To produce the 8 figures in the paper (Figure 1 in the main paper and Figures 2-8 in the appendix) we worked with: 57-72 real-world classification datasets (OpenML Curated Classification benchmarking suite 2018 [Bischl et al., 2017]) converted to contextual bandit datasets (wide-spread evaluation approach in the offline and online contextual bandit literature); 4 different target policies; 6 estimators (CADR, and 5 other popular baseline estimators in the literature); 8 different training procedures (sequential cross-fitting vs. cross-time cross-fitting in Figures 3 and 3; misspecified vs. well-specified outcome model family in Figures 4 and 5; weighted vs. unweighted outcome model fitting in Figures 6 and 7; large data vs. small data in Figures 3 and 8). For each (dataset, target policy, estimator, training procedure)-combination, 64 simulations were run to obtain and report tight standard errors of the estimators’ performance.
In the subsections that follow, we present our approach and the main results.

heading: 4.1 Baseline Estimators
text: We compare CADR to several benchmarks (our experiments focus on the case of finitely-many actions, A = {1, . . . ,K}). All baselines take the following form for a choice of wt, ωt, ̂̄Qt:
Ψ̂T =
( 1
T wt
)−1 1
T T∑ t=1 wtD̃ ′ t, CIα = Ψ̂T ± ζ1−α/2 √√√√√∑Tt=1 w2t (D̃′t − Ψ̂)2(∑T
t=1 wt
)2  ,
where D̃′t = ωt(Y (t)− ̂̄Qt−1(A(t), X(t))) + K∑ a=1 ̂̄Qt−1(a,X(t))g∗(a | X(t)). The Direct Method (DM) sets wt = 1, ωt = 0 and fits ˆ̄Qt−1(a, ·) by running some regression method for each a on the data {(X(s), Y (s)) : 1 ≤ s ≤ t − 1, A(s) = a}. We will use either linear regression or decision-tree regression, both using default sklearn parameters. Note that even in non-contextual settings, where ˆ̄Qt−1 is a simple per-arm sample average, ˆ̄Qt−1 may be biased due to adaptive data collection [Xu et al., 2013, Luedtke and van der Laan, 2016, Bowden and Trippa, 2017, Nie et al., 2018, Hadad et al., 2019, Shin et al., 2019]. Inverse Propensity Score Weighting (IPW) sets wt = 1, ωt = (g∗/gt)(A(t) | X(t)), ˆ̄Qt = 0. Doubly Robust (DR) sets wt = 1, ωt = (g
∗/gt)(A(t) | X(t)) and fits ˆ̄Qt−1 as in DM. More Robust Doubly Robust (MRDR) [Farajtabar et al., 2018] is the same as DR but when fitting ˆ̄Qt−1 we reweight each data point by g
∗(A(s)|X(s))(1−gs(A(s)|X(s))) gs(A(s)|X(s))2 . None of the above are generally asymptotically normal under adaptive data collection [Hadad et al., 2019]. Adaptive Doubly Robust (ADR; a.k.a. stabilized onestep estimator for multi-armed bandit data) [Luedtke and van der Laan, 2016, Hadad et al., 2019] is the same as DR but sets wt = g −1/2 t (A(t)|X(t)). ADR is unbiased and asymptotically normal for multi-armed bandit logging policies but is biased for context-measurable adaptive logging policies, which is the focus of this paper. Finally, note that our proposal CADR takes the same form as DR but with wt = σ̂−1t using our adaptive conditional standard deviation estimators σ̂t in Eq. (3).
1The code can be found at https://github.com/mdimakopoulou/post-contextual-bandit-inference.

heading: 4.2 Contextual Bandit Data from Multiclass Classification Data
text: To construct our data, we turn K-class classification tasks into a K-armed contextual bandit problems [Dudı́k et al., 2014, Dimakopoulou et al., 2017, Su et al., 2019], which has the benefits of reproducibility using public datasets and being able to make uncontroversial comparisons using actual ground truth data with counterfactuals. We use the public OpenML Curated Classification benchmarking suite 2018 (OpenML-CC18; BSD 3-Clause license) [Bischl et al., 2017], which has datasets that vary in domain, number of observations, number of classes and number of features. Among these, we select the classification datasets which have less than 100 features. This results in 57 classification datasets from OpenML-CC18 used for evaluation and Table 1 in section E.1 of the appendix summarizes the characteristics of these datasets.
Each dataset is a collection of pairs of covariates X and labels L ∈ {1, . . . ,K}. We transform each dataset to the contextual bandit problem as follows. At each round, we draw X(t), L(t) uniformly at random with replacement from the dataset. We reveal the context X(t) to the agent, and given an arm pull A(t), we draw and return the reward Y (t) ∼ N (1{A(t) = L(t)}, 1). To generate our data, we set T = 10000 and use the following -greedy procedure. We pull arms uniformly at random until each arm has been pulled at least once. Then at each subsequent round t, we fit ̂̄Qt−1 using the data up to that time in the same fashion as used for the DM estimator above using decision-tree regressions. We set Ãx(t) = arg maxa=1,...,K
̂̄Qt−1(a,X(t)) and t = 0.01 · t−1/3. We then let gt(a | x) = t/K for a 6= Ãx(t) and gt(Ãx(t) | x) = 1 − t + t/K. That is, with probability t we pull a random arm, and otherwise we pull ÃX(t)(t).
We then consider four candidate policies to evaluate: (1) “arm 1 non-contextual”: g∗(1 | x) = 1 and otherwise g∗(a | x) = 0 (note that the meaning of label “1” changes by dataset), (2) “arm 2 noncontextual”: g∗(2 | x) = 1 and otherwise g∗(a | x) = 0, (3) “linear contextual”: we sample a new
dataset of size T using a uniform exploration policy, then fit ̂̄QT as above using linear regression, fix a∗ = arg maxa∈{1,...,K}
̂̄QT (a, x), and set g∗(a∗ | x) = 1 and otherwise g∗(a | x) = 0, (4) “tree contextual”: same as “linear contextual” but fit ̂̄QT using decision-tree regression.

heading: 4.3 Results
text: Figure 1 shows the comparison of CADR estimator against DM, IPW, DR, ADR, and MRDR w.r.t. coverage, that is, the frequency over 64 replications of the 95% confidence interval covering the true Ψ0, for each of the 57 OpenML-CC18 datasets and 4 target policies. In each subfigure, each dot represents a dataset, the y-axis corresponds to the coverage of the CADR estimator and the x-axis corresponds to the coverage of one of the baseline estimators. The lines represent one standard error over the 64 replications. The dot is depicted in blue if for that dataset CADR has significantly (i.e., by more than one standard error) better coverage than the baseline estimator, in red if it has significantly worse coverage, and in black if the difference in coverage of both estimators is within one standard error. Results are averaged over 64 replications that we ran for each dataset, estimator, and target policy. The standard errors (computed as the sample standard deviation over the 64 replications divided by √ 64) are depicted and the 95% coverage targets are shown in dashed lines. In Figure 1, outcome models for CADR, DM, DR, ADR, and MRDR are fit using linear regression (with default sklearn parameters). In the appendix, we provide additional empirical results where we use decision-tree regressions, or where we use the MRDR outcome model for CADR, or where we use cross-fold estimation across time.
Across all of our experiments, we observe that the confidence interval of CADR has better coverage of the ground truth than any other baseline, which can be attributed to its asymptotic normality. The second best estimator in terms of coverage is DR. The advantages of CADR over DR are most pronounced when either (a) there is a mismatch between the logging policy and the target policy (e.g., compare the 1st and 2nd rows in Figure 1; the tree target policy is most similar to the logging policy, which also uses trees) or (b) when the outcome model is bad (either due to model misspecification such as with a linear model on real data or due to small sample size).

heading: 5 Conclusions
text: Adaptive experiments hold great promise for better, more efficient, and even more ethical experiments. However, they complicate post-experiment inference, which is a cornerstone of drawing credible conclusions from controlled experiments. We provided here a novel asymptotically normal estimator for policy value and causal effects when data were generated from a contextual adaptive experiment, such as a contextual bandit algorithm. This led to simple and effective confidence intervals given by adding and subtracting multiples of the standard error, making contextual adaptive experiments a more viable option for experimentation in practice.

heading: 6 Societal Impact and Limitations
text: Adaptive experiments hold particular promise in settings where experimentation is costly and/or dangerous, such as in medicine and policymaking. By adapting treatment allocation, harmful interventions can be avoided, outcomes for study participants improved, and smaller studies enabled. Being able to draw credible conclusions from such experiments make them viable replacements for classic randomized trials. Our confidence intervals offer one way to do so. At the same time, and especially subject to our assumption of vanishing but nonzero exploration, these experiments must be subject to the same ethical guidelines as classic randomized experiments. Additionally, the usual caveats of frequentist confidence intervals hold here, such as its interpretation only as a guarantee over data collection, this guarantee only being approximate in finite samples when we rely on asymptotic normality, and the risks of multiple comparisons and of p-hacking. Finally, we note that our inference focused on an average quantity, as such it focuses on social welfare and need not capture the risk to individuals or groups. Subgroup analyses may therefore be helpful in complementing the analysis; these can be conducted by setting g∗(a | x) to zero for some x’s. Future work may be necessary to further extend our results to conducting inference on risk metrics such as quantiles of outcomes.

Review summaries:
The current paper concerns inference of quantities such as treatment effects or policy values from data that is collected from contextual bandit experiments, Motivated by applications in healthcare, revenue management or product recommendations where randomized experimentation is not often the norm. They isolate a key issue (variance stabilization), identified in previous works on analysis of adaptively collected data (see paper for references), and propose estimators for average causal effects of the  `doubly robust` form to overcome this issue. The main result is a central limit theorem for this estimator, allowing for frequentist inference (i.e. confidence intervals, p values) in the standard way. The estimator relies heavily on a sequence of conditional variance estimators, which  they construct using a judicious use of averaging and a base regression model for the outcome ( which may be misspecified even in the limit). 

This study proposes a variance-stabilized estimator for policy evaluation from data obtained from contextual bandit algorithms, which is inspired by [Luedtke and van der Laan 2016]  and [Hadad et al. 2019].  In this study, the authors estimate the variance of the score function (canonical gradient) from samples, which standardizes the MDS for applying martingale CLT. [Luedtke, and van der Laan 2016] Statistical inference for the mean outcome under a possibly non-unique optimal treatment strategy   [Hadad et al 2019] Confidence Intervals for Policy Evaluation in Adaptive Experiments.

This paper is motivated by the fact that no asymptotically normal estimators exist when the data is collected context-dependent adaptively. This work provides such estimators (named CADR, Algo. 1) exhibiting asymptotic normality (Thm. 1).