Statistical Inference with M-Estimators on Adaptively Collected Data
Year: 2021

Abstract:
Bandit algorithms are increasingly used in real-world sequential decision-making problems. Associated with this is an increased desire to be able to use the resulting datasets to answer scientific questions like: Did one type of ad lead to more purchases? In which contexts is a mobile health intervention effective? However, classical statistical approaches fail to provide valid confidence intervals when used with data collected with bandit algorithms. Alternative methods have recently been developed for simple models (e.g., comparison of means). Yet there is a lack of general methods for conducting statistical inference using more complex models on data collected with (contextual) bandit algorithms; for example, current methods cannot be used for valid inference on parameters in a logistic regression model for a binary reward. In this work, we develop theory justifying the use of M-estimators—which includes estimators based on empirical risk minimization as well as maximum likelihood—on data collected with adaptive algorithms, including (contextual) bandit algorithms. Specifically, we show that M-estimators, modified with particular adaptive weights, can be used to construct asymptotically valid confidence regions for a variety of inferential targets.

heading: 1 Introduction
text: Due to the need for interventions that are personalized to users, (contextual) bandit algorithms are increasingly used to address sequential decision making problems in health-care [Yom-Tov et al., 2017, Liao et al., 2020], online education [Liu et al., 2014, Shaikh et al., 2019], and public policy [Kasy and Sautmann, 2021, Caria et al., 2020]. Contextual bandits personalize, that is, minimize regret, by learning to choose the best intervention in each context, i.e., the action that leads to the greatest expected reward. Besides the goal of regret minimization, another critical goal in these real-world problems is to be able to use the resulting data collected by bandit algorithms to advance scientific knowledge [Liu et al., 2014, Erraqabi et al., 2017]. By scientific knowledge, we mean information gained by using the data to conduct a variety of statistical analyses, including confidence interval construction and hypothesis testing. While regret minimization is a within-experiment learning objective, gaining scientific knowledge from the resulting adaptively collected data is a between-experiment learning objective, which ultimately helps with regret minimization between deployments of bandit algorithms. Note that the data collected by bandit algorithms are adaptively collected because previously observed contexts, actions, and rewards are used to inform what actions to select in future timesteps.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
There are a variety of between-experiment learning questions encountered in real-life applications of bandit algorithms. For example, in real-life sequential decision-making problems there are often a number of additional scientifically interesting outcomes besides the reward that are collected during the experiment. In the online advertising setting, the reward might be whether an ad is clicked on, but one may be interested in the outcome of amount of money spent or the subsequent time spent on the advertiser’s website. If it was found that an ad had high click-through rate, but low amounts of money was spent after clicking on the ad, one may redesign the reward used in the next bandit experiment. One type of statistical analysis would be to construct confidence intervals for the relative effect of the actions on multiple outcomes (in addition to the reward) conditional on the context. Furthermore, due to engineering and practical limitations, some of the variables that might be useful as context are often not accessible to the bandit algorithm online. If after-study analyses find some such contextual variables to have sufficiently strong influence on the relative usefulness of an action, this might lead investigators to ensure these variables are accessible to the bandit algorithm in the next experiment.
As discussed above, we can gain scientific knowledge from data collected with (contextual) bandit algorithms by constructing confidence intervals and performing hypothesis tests for unknown quantities such as the expected outcome for different actions in various contexts. Unfortunately, standard statistical methods developed for i.i.d. data fail to provide valid inference when applied to data collected with common bandit algorithms. For example, assuming the sample mean of rewards for an arm is approximately normal can lead to unreliable confidence intervals and inflated type-1 error; see Section 3.1 for an illustration. Recently statistical inference methods have been developed for data collected using bandit algorithms [Hadad et al., 2019, Deshpande et al., 2018, Zhang et al., 2020]; however, these methods are limited to inference for parameters of simple models. There is a lack of general statistical inference methods for data collected with (contextual) bandit algorithms in more complex data-analytic settings, including parameters in non-linear models for outcomes; for example, there are currently no methods for constructing valid confidence intervals for the parameters of a logistic regression model for binary outcomes or for constructing confidence intervals based on robust estimators like minimizers of the Huber loss function.
In this work we show that a wide variety of estimators which are frequently used both in science and industry on i.i.d. data, namely, M-estimators [Van der Vaart, 2000], can be used to conduct valid inference on data collected with (contextual) bandit algorithms when adjusted with particular adaptive weights, i.e., weights that are a function of previously collected data. Different forms of adaptive weights are used by existing methods for simple models [Deshpande et al., 2018, Hadad et al., 2019, Zhang et al., 2020]. Our work is a step towards developing a general framework for statistical inference on data collected with adaptive algorithms, including (contextual) bandit algorithms.

heading: 2 Problem Formulation
text: We assume that the data we have after running a contextual bandit algorithm is comprised of contexts {Xt}Tt=1, actions {At}Tt=1, and primary outcomes {Yt}Tt=1. T is deterministic and known. We assume that rewards are a deterministic function of the primary outcomes, i.e., Rt = f(Yt) for some known function f . We are interested in constructing confidence regions for the parameters of the conditional distribution of Yt given (Xt, At). Below we consider T ! 1 in order to derive the asymptotic distributions of estimators and construct asymptotically valid confidence intervals. We allow the action space A to be finite or infinite. We use potential outcome notation [Imbens and Rubin, 2015] and let {Yt(a) : a 2 A} denote the potential outcomes of the primary outcome and let Yt := Yt(At) be the observed outcome. We assume a stochastic contextual bandit environment in which {Xt, Yt(a) : a 2 A} i.i.d. ⇠ P 2 P for t 2 [1 : T ]; the contextual bandit environment distribution P is in a space of possible environment distributions P. We define the history Ht := {Xt0 , At0 , Yt0}tt0=1 for t 1 and H0 := ;. Actions At 2 A are selected according to policies ⇡ := {⇡t}t 1, which define action selection probabilities ⇡t(At, Xt,Ht 1) := P (At|Ht 1, Xt). Even though the potential outcomes are i.i.d., the observed data {Xt, At, Yt}Tt=1 are not because the actions are selected using policies ⇡t which are a function of past data, Ht 1. Non-independence of observations is a key property of adaptively collected data.
We are interested in constructing confidence regions for some unknown ✓⇤(P) 2 ⇥ ⇢ Rd, which is a parameter of the conditional distribution of Yt given (Xt, At). This work focuses on the setting in which we have a well-specified model for Yt. Specifically, we assume that ✓⇤(P) is a conditionally
maximizing value of criterion m✓, i.e., for all P 2 P, ✓⇤(P) 2 argmax
✓2⇥ EP [m✓(Yt, Xt, At)|Xt, At] w.p. 1. (1)
Note that ✓⇤(P) does not depend on (Xt, At) and it is an implicit modelling assumption that such a ✓⇤(P) exists for a given m✓. Note that this formulation includes semi-parametric models, e.g., the model could constrain the conditional mean of Yt to be linear in some function of the actions and context, but allow the residuals to follow any mean-zero distribution, including ones that depend on the actions and/or contexts.
To estimate ✓⇤(P), we build on M-estimation [Huber, 1992], which classically selects the estimator ✓̂ to be the ✓ 2 ⇥ that maximizes the empirical analogue of Equation (1):
✓̂T := argmax ✓2⇥
1
T
TX
t=1
m✓(Yt, Xt, At). (2)
For example, in a classical linear regression setting with |A| < 1 actions, a natural choice for m✓ is the negative of the squared loss function, m✓(Yt, Xt, At) = (Yt X>t ✓At)2. When Yt is binary, a natural choice is instead the negative log-likelihood function for a logistic regression model, i.e., m✓(Yt, Xt, At) = [YtX>t ✓At log(1 + exp(X >
t ✓At))]. More generally, m✓ is commonly chosen to be a log-likelihood function or the negative of a robust loss function such as the Huber loss. If the data, {Xt, At, Yt}Tt=1, were independent across time, classical approaches could be used to prove the consistency and asymptotic normality of M-estimators [Van der Vaart, 2000]. However, on data collected with bandit algorithms, standard M-estimators like the ordinary least-squares estimator fail to provide valid confidence intervals [Hadad et al., 2019, Deshpande et al., 2018, Zhang et al., 2020]. In this work, we show that M-estimators can still be used to provide valid statistical inference on adaptively collected data when adjusted with well-chosen adaptive weights.

heading: 3 Adaptively Weighted M-Estimators
text: We consider a weighted M-estimating criteria with adaptive weights Wt 2 (Ht 1, Xt, At) given
by Wt = q ⇡stat (At,Xt) ⇡t(At,Xt,Ht 1)
. Here {⇡stat }t 1 are pre-specified stabilizing policies that do not depend on data {Yt, Xt, At}t 1. A default choice for the stabilizing policy when the action space is of size |A| < 1 is just ⇡stat (a, x) = 1/|A| for all x, a, and t; we discuss considerations for the choice of {⇡stat } T t=1 in Section 3.3. We call these weights square-root importance weights because they are the square-root of the standard importance weights [Hammersley, 2013, Wang et al., 2017]. Our proposed estimator for ✓⇤(P), ✓̂T , is the maximizer of a weighted version of the M-estimation criterion of Equation (2):
✓̂T := argmax ✓2⇥
1
T
TX
t=1
Wtm✓(Yt, Xt, At) =: argmax ✓2⇥ MT (✓).
Note that MT (✓) defined above depends on both the data {Xt, At, Yt}Tt=1 and weights {Wt}Tt=1. We provide asymptotically valid confidence regions for ✓⇤(P) by deriving the asymptotic distribution of ✓̂T as T ! 1 and by proving that the convergence in distribution is uniform over P 2 P. Such convergence allows us to construct a uniformly asymptotically valid 1 ↵ level confidence region, CT (↵), for ✓⇤(P), which is a confidence region that satisfies
lim inf T!1 inf P2P PP,⇡ (✓⇤(P) 2 CT (↵)) 1 ↵. (3)
If CT (↵) were not uniformly valid, then there would exist an ✏ > 0 such that for every sample size T , CT (↵)’s coverage would be below 1 ↵ ✏ for some worst-case PT 2 P. Confidence regions which are asymptotically valid, but not uniformly asymptotically valid, fail to be reliable in practice [Leeb and Pötscher, 2005, Romano et al., 2012]. Note that on i.i.d. data it is generally straightforward to show that estimators that converge in distribution do so uniformly; however, as discussed in Zhang et al. [2020] and Appendix D, this is not the case on data collected with bandit algorithms.
To construct uniformly valid confidence regions for ✓⇤(P) we prove that ✓̂T is uniformly asymptotically normal in the following sense:
⌃T (P) 1/2M̈T (✓̂T )
p
T (✓̂T ✓ ⇤(P)) D ! N (0, Id) uniformly over P 2 P, (4)
where M̈T (✓) := @ 2 @2✓MT (✓) and ⌃T (P) := 1 T PT t=1 EP,⇡stat ⇥ ṁ✓⇤(P)(Yt, Xt, At) ⌦2 ⇤ . We define ṁ✓ := @ @✓m✓. Similarly we define respectively m̈✓ and ... m✓ as the second and third partial derivatives of m✓ with respect to ✓. For any vector z we define z⌦2 := zz>.

heading: 3.1 Intuition for Square-Root Importance Weights
text: The critical role of the square-root importance weights Wt = q ⇡stat (At,Xt) ⇡t(At,Xt,Ht 1)
is to adjust for instability in the variance of M-estimators due to the bandit algorithm. These weights act akin to standard importance weights when squared and adjust a key term in the variance of M-estimators from depending on adaptive policies {⇡t}Tt=1, which can be ill-behaved, to depending on the prespecified stabilizing policies {⇡stat }Tt=1. See Zhang et al. [2020] and Deshpande et al. [2018] for more discussion of the ill-behavior of the action selection probabilities for common bandit algorithms, which occurs particularly when there is no unique optimal policy.
As an illustrative example, consider the least-squares estimators in a finite-arm linear contextual bandit setting. Assume that EP [Yt|Xt, At = a] = X>t ✓⇤a(P) w.p. 1. We focus on estimating ✓⇤a(P) for some a 2 A. The least-squares estimator corresponds to an M-estimator with m✓a(Yt, Xt, At) = At=a(Yt X > t ✓a) 2. The adaptively weighted least-squares (AW-LS) estimator is ✓̂AW-LST,a := argmax✓a{ PT t=1 Wt At=a(Yt X > t ✓a) 2 }. For simplicity, suppose that the stabilizing policy does not change with t and drop the index t to get ⇡sta. Taking the derivative of this criterion, we get 0 = PT t=1 Wt At=aXt Yt X>t ✓̂ AW-LS T,a , and rearranging terms gives
1 p T
TX
t=1
Wt At=aXtX > t ⇣ ✓̂AW-LST,a ✓ ⇤ a(P) ⌘ =
1 p T
TX
t=1
Wt At=aXt Yt X
> t ✓ ⇤ a(P) . (5)
Note that the right hand side of Equation (5) is a martingale difference sequence with respect to history {Ht}Tt=0 because EP,⇡[Wt At=a(Yt X>t ✓⇤a(P))|Ht 1] = 0 for all t; by law of iterated expectations and since Wt 2 (Ht 1, Xt, At), EP,⇡[Wt At=a(Yt X>t ✓⇤a(P))|Ht 1] equals
EP ⇥ Wt⇡t(a,Xt,Ht 1)EP ⇥ Yt X > t ✓ ⇤ a(P)|Ht 1, Xt, At = a ⇤ Ht 1 ⇤
= (i)
EP ⇥ Wt⇡t(a,Xt,Ht 1)EP ⇥ Yt X > t ✓ ⇤ a(P)|Xt, At = a ⇤ Ht 1 ⇤ = (ii) 0.
(i) holds by our i.i.d. potential outcomes assumption. (ii) holds since EP [Yt|Xt, At = a] = X>t ✓ ⇤
a(P). We prove that (5) is uniformly asymptotically normal by applying a martingale central limit theorem (Appendix B.4). The key condition in this theorem is that the conditional variance converges uniformly, for which it is sufficient to show that the conditional covariance of Wt At=a Yt X>t ✓ ⇤ a(P) given Ht 1 equals some positive-definite matrix ⌃(P) for every t, i.e.,
EP,⇡ h W 2t At=aXtX > t Yt X > t ✓ ⇤ a(P) 2 Ht 1 i = ⌃(P). (6)
By law of iterated expectations, EP,⇡[W 2t At=aXtX>t (Yt X>t ✓⇤a(P))2 Ht 1] equals
EP  EP,⇡  ⇡sta(At, Xt)
⇡t(At, Xt,Ht 1) At=aXtX
> t Yt X > t ✓ ⇤ a(P) 2 Ht 1, Xt Ht 1
(7)
= (a)
EP  EP,⇡sta  At=aXtX > t Yt X > t ✓ ⇤ a(P) 2 Ht 1, Xt Ht 1
= (b)
EP  EP,⇡sta  At=aXtX > t Yt X > t ✓ ⇤ a(P) 2 Xt Ht 1
= (c)
EP  EP,⇡sta  At=aXtX > t Yt X > t ✓ ⇤ a(P) 2 Xt
= (d)
EP,⇡sta [ At=aXtX>t (Yt X>t ✓⇤a(P))2] =: ⌃(P).
Above, (a) holds because the importance weights change the sampling measure from the adaptive policy ⇡t to the pre-specified stabilizing policy ⇡sta. (b) holds by our i.i.d. potential outcomes
assumption and because ⇡sta is a pre-specified policy. (c) holds because Xt does not depend on Ht 1 by our i.i.d. potential outcomes assumption. (d) holds by the law of iterated expectations. Note that ⌃(P) does not depend on t because ⇡sta is not time-varying. In contrast, without the adaptive weighting, i.e., when Wt = 1, the conditional covariance of At=a Yt X>t ✓ ⇤ a(P)
on Ht 1 is a random variable, due to the adaptive policy ⇡t.
In Figure 1 we plot the empirical distributions of the z-statistic for the least-squares estimator both with and without adaptive weighting. We consider a two-armed bandit with At 2 {0, 1}. Let ✓⇤1(P) := EP [Yt(1)] and m✓1(Yt, At) := At(Yt ✓1)2. The unweighted version, i.e., the ordinary least-squares (OLS) estimator, is ✓̂OLST,1 := argmax✓1 1 T PT t=1 m✓1(Yt, At). The adaptively weighted version is ✓̂AW-LST,1 := argmax✓1 1 T PT t=1 Wtm✓1(Yt, At). We collect data using Thompson Sampling and use a uniform stabilizing policy where ⇡sta(1) = ⇡sta(0) = 0.5. It is clear that the least-squares estimator with adaptive weighting has a z-statistic that is much closer to a normal distribution.
t=1 At(✓̂
OLS
T,1 ✓
⇤
1(P)) on the left and
⇣
1
p T
PT
t=1
q
0.5
⇡t(1)
At
⌘
(✓̂AW-LST,1 ✓
⇤
1(P)) on the right.
The square-root importance weights are a form of variance stabilizing weights, akin to those introduced in Hadad et al. [2019] for estimating means and differences in means on data collected with multi-armed bandits. In fact, in the special case that |A| < 1 and (Xt, At) = [ At=1, At=2, ..., At=|A|]
>, the adaptively weighted least-squares estimator is equivalent to the weighted average estimator of Hadad et al. [2019]. See Section 4 for more on Hadad et al. [2019].

heading: 3.2 Asymptotic Normality and Confidence Regions
text: We now discuss conditions under which the adaptively weighted M-estimators are asymptotically normal in the sense of Equation (4). In general, our conditions differ from those made for standard M-estimators on i.i.d. data because (i) the data is adaptively collected, i.e., ⇡t can depend on Ht 1 and (ii) we ensure uniform convergence over P 2 P, which is stronger than guaranteeing convergence pointwise for each P 2 P.
Condition 1 (Stochastic Bandit Environment). Potential outcomes {Xt, Yt(a) : a 2 A} i.i.d. ⇠ P 2 P over t 2 [1 : T ].
Condition 1 implies that Yt is independent of Ht 1 given Xt and At, and the conditional distribution Yt | Xt, At is invariant over time. Also note that action space A can be finite or infinite. Condition 2 (Differentiable). The first three derivatives of m✓(y, x, a) with respect to ✓ exist for every ✓ 2 ⇥, every a 2 A, and every (x, y) in the joint support of {P : P 2 P}. Condition 3 (Bounded Parameter Space). For all P 2 P, ✓⇤(P) 2 ⇥, a bounded open subset of Rd. Condition 4 (Lipschitz). There exists some real-valued function g such that (i) sup
P2P,t 1 EP,⇡stat [g(Yt, Xt, At) 2] is bounded and (ii) for all ✓, ✓0 2 ⇥,
|m✓(Yt, Xt, At) m✓0(Yt, Xt, At)|  g(Yt, Xt, At)k✓ ✓ 0 k2.
Conditions 3 and 4 together restrict the complexity of the function m in order to ensure a martingale law of large numbers result holds uniformly over functions {m✓ : ✓ 2 ⇥}; this is used to prove the consistency of ✓̂T . Similar conditions are commonly used to prove consistency of M-estimators based on i.i.d. data, although the boundedness of the parameter space can be dropped when m✓ is a concave function of ✓ for all Yt, At, Xt (as it is in many canonical examples such as least squares) [Van der Vaart, 2000, Engle, 1994, Bura et al., 2018]; we expect that a similar result would hold for adaptively weighted M-estimators. Condition 5 (Moments). The fourth moments of m✓⇤(P)(Yt, Xt, At), ṁ✓⇤(P)(Yt, Xt, At), and m̈✓⇤(P)(Yt, Xt, At) with respect to P and policy ⇡ sta t are bounded uniformly over P 2 P and t 1. For all sufficiently large T , the minimum eigenvalue of ⌃T,P := 1 T PT t=1 EP,⇡stat ⇥ ṁ✓⇤(P)(Yt, Xt, At) ⌦2 ⇤ is bounded above ṁ2 > 0 for all P 2 P.
Condition 5 is similar to those of Van der Vaart [2000, Theorem 5.41]. However, to guarantee uniform convergence we assume that moment bounds hold uniformly over P 2 P and t 1. Condition 6 (Third Derivative Domination). For B 2 Rd⇥d⇥d, we define kBk1 :=Pd
i=1 Pd j=1 Pd k=1 |Bi,j,k|. There exists a function ... m(Yt, Xt, At) 2 Rd⇥d⇥d such that (i)
sup P2P,t 1 EP,⇡stat
⇥ k ... m(Yt, Xt, At)k21 ⇤ is bounded and (ii) for all P 2 P there exists some ✏...m > 0
such that the following holds with probability 1,
sup ✓2⇥:k✓ ✓⇤(P)k✏...m
k ... m✓(Yt, Xt, At)k1  k ... m(Yt, Xt, At)k1.
Condition 6 is again similar to those in classical M-estimator asymptotic normality proofs [Van der Vaart, 2000, Theorem 5.41]. Condition 7 (Maximizing Solution). (i) For all P 2 P, there exists a ✓⇤(P) 2 ⇥ such that (a) ✓⇤(P) 2 argmax✓2⇥ EP ⇥ m✓(Yt, Xt, At) Xt, At ⇤ w.p. 1, (b) EP ⇥ ṁ✓⇤(P)(Yt, Xt, At) Xt, At ⇤ = 0 w.p. 1,
and (c) EP ⇥ m̈✓⇤(P)(Yt, Xt, At) Xt, At ⇤ 0 w.p. 1. (ii) There exists some positive definite matrix H such that 1T PT t=1 EP,⇡stat ⇥ m̈✓⇤(P)(Yt, Xt, At) ⇤ ⌫

heading: H for all P 2 P and all sufficiently large T .
text: For matrices A,B, we define A ⌫ B to mean that A B is positive semi-definite, as used above. Condition 7 (i) ensures that ✓⇤(P) is a conditionally maximizing solution for all contexts Xt and actions At; this ensures that {ṁ✓⇤(P)(Yt, Xt, At)}Tt=1 is a martingale difference sequence with respect to {Ht}Tt=1. Note it does not require ✓⇤(P) to always be a conditionally unique optimal solution. Condition 7 (ii) is related to the local curvature at the maximizing solution and the analogous condition in the i.i.d. setting is trivially satisfied; we specifically use this condition to ensure we can replace M̈(✓⇤(P)) with M̈(✓̂T ) in our asymptotic normality result, i.e., that M̈(✓⇤(P)) 1M̈(✓̂T ) P ! Id uniformly over P 2 P.
Condition 8 (Well-Separated Solution). For all sufficiently large T , for any ✏ > 0, there exists some > 0 such that for all P 2 P,
inf ✓2⇥:k✓ ✓⇤(P)k2>✏
⇢ 1
T
TX
t=1
EP,⇡stat ⇥ m✓⇤(P)(Yt, Xt, At) m✓(Yt, Xt, At) ⇤ .
A well-separated solution condition akin to Condition 8 is commonly assumed in order to prove consistency of M-estimators, e.g., see Van der Vaart [2000, Theorem 5.7]. Note that the difference between Condition 7 (i) and Condition 8 is that the former is a conditional statement (conditional on Xt, At) and the latter is a marginal statement (marginal over Xt, At, where At is chosen according to stabilizing policies ⇡stat ). Condition 7 (i) means there is a ✓⇤(P) solution for all contexts Xt and actions At that does not need to be unique, however Condition 8 assumes that marginally over Xt, At there is a well-separated solution. Condition 9 (Bounded Importance Ratios). {⇡stat }Tt=1 do not depend on data {Yt, Xt, At}Tt=1. For all t 1, ⇢min  ⇡stat (At,Xt) ⇡t(At,Xt,Ht 1)  ⇢max w.p. 1 for some constants 0 < ⇢min  ⇢max < 1.
Note that Condition 9 implies that for a stabilizing policy that is not time-varying, the action selection probabilities of the bandit algorithm ⇡t(At, Xt,Ht 1) must be bounded away from zero w.p. 1. Similar boundedness assumptions are also made in the off-policy evaluation literature [Thomas and Brunskill, 2016, Kallus and Uehara, 2020]. We discuss this condition further in Sections 3.3 and 6. Theorem 1 (Uniform Asymptotic Normality of Adaptively Weighted M-Estimators). Under Conditions 1-9 we have that ✓̂T P ! ✓⇤(P) uniformly over P 2 P. Additionally,
⌃T (P) 1/2M̈T (✓̂T )
p
T (✓̂T ✓ ⇤(P)) D ! N (0, Id) uniformly over P 2 P. (8)
The asymptotic normality result of equation (8) guarantees that for d-dimensional ✓⇤(P),
lim inf T!1 inf P2P
PP,⇡ ✓h ⌃T (P) 1/2M̈T (✓̂T ) p T (✓̂T ✓ ⇤(P)) i⌦2  2d,(1 ↵) ◆ = 1 ↵.
Above 2d,(1 ↵) is the 1 ↵ quantile of the 2 distribution with d degrees of freedom. Note that the region CT (↵) := ✓ 2 ⇥ : [⌃T (P) 1/2M̈T (✓̂T ) p T (✓̂T ✓⇤(P))]⌦2  2d,(1 ↵) defines a
d-dimensional hyper-ellipsoid confidence region for ✓⇤(P). Also note that since M̈T (✓̂T ) does not concentrate under standard bandit algorithms, we cannot use standard arguments to justify treating ✓̂T as multivariate normal with covariance M̈T (✓̂T ) 1⌃T (P)M̈T (✓̂T ) 1. Nevertheless, Theorem 1 can be used to guarantee valid confidence regions for subset of entries in ✓⇤(P) by using projected confidence regions [Nickerson, 1994]. Projected confidence regions take a confidence region for all parameters ✓⇤(P) and project it onto the lower dimensional space on which the subset of target parameters lie (Appendix A.2).

heading: 3.3 Choice of Stabilizing Policy
text: When the action space is bounded, using weights Wt = 1/ p ⇡t(At, Xt,Ht 1) is equivalent to using square-root importance weights with a stabilizing policy that selects actions uniformly over A; this is because weighted M-estimators are invariant to all weights being scaled by the same constant. It can make sense to choose a non-uniform stabilizing policy in order to prevent the square-root importance weights from growing too large and to ensure Condition 9 holds; disproportionately up-weighting a few observations can lead to unstable estimators. Note that an analogue of our stabilizing policy exists in the causal inference literature, namely, “stabilized weights" use a probability density in the numerator of the weights to prevent them from becoming too large [Robins et al., 2000].
We now discuss how to choose stabilizing policies {⇡stat }t 1 in order to minimize the asymptotic variance of adaptively weighted M-estimators. We focus on the adaptively weighted least-squares estimator when we have a linear outcome model EP [Yt|Xt, At] = X>t ✓At :
✓̂AW-LS := argmax ✓2⇥
⇢ 1
T
TX
t=1
Wt Yt X
> t ✓At 2 . (9)
Recall that our use of adaptive weights is to adjust for instability in the variance of M-estimators induced by the bandit algorithm in order to construct valid confidence regions; note that weighted estimators are not typically used for this reason. On i.i.d. data, the least-squares criterion is weighted like in Equation (9) in order to minimize the variance of estimators under noise heteroskedasticity; in this setting, the best linear unbiased estimator has weights Wt = 1/ 2(At, Xt) where 2(At, Xt) := EP [(Yt X>t ✓⇤At(P)) 2 |Xt, At]; this up-weights the importance of observations with low noise variance. Intuitively, if we do not need to variance stabilize, {Wt}t 1 should be determined by the relative importance of minimizing the errors for different observations, i.e., their noise variance.
In light of this observation, we expect that under homoskedastic noise there is no reason to up-weight some observations over others. This would recommend choosing the stabilizing policy to make Wt =p ⇡stat (At, Xt)/⇡t(At, Xt,Ht 1) as close to 1 as possible, subject to the constraint that the stabilizing policies are pre-specified, i.e., {⇡stat }t 1 do not depend on data {Yt, Xt, At}t 1 (see Appendix C for details). Since adjusting for heteroskedasticity and variance stabilization are distinct uses of weights, under heteroskedasticity, we recommend that the weights are combined in the following sense: Wt = 1/ 2(At, Xt) p ⇡stat (At, Xt)/⇡t(At, Xt,Ht 1). This would mean that to minimize variance, we still want to choose the stabilizing policies to make ⇡stat (At, Xt)/⇡t(At, Xt,Ht 1) as close to 1 possible, subject to the pre-specified constraint.

heading: 4 Related Work
text: Villar et al. [2015] and Rafferty et al. [2019] empirically illustrate that classical ordinary least squares (OLS) inference methods have inflated Type-1 error when used on data collected with a variety of regret-minimizing multi-armed bandit algorithms. Chen et al. [2020] prove that the OLS estimator is asymptotically normal on data collected with an ✏-greedy algorithm, but their results do not cover settings in which there is no unique optimal policy, e.g., a multi-arm bandit with two identical arms (Appendix E). Recent work has discussed the non-normality of OLS on data collected with bandit algorithms when there is no unique optimal policy and proposed alternative methods for statistical inference. A common thread between these methods is that they all utilize a form of adaptive weighting. Deshpande et al. [2018] introduced the W-decorrelated estimator, which adjusts the OLS estimator with a sum of adaptively weighted residuals. In the multi-armed bandit setting, the W-decorrelated estimator up-weights observations from early in the study and down-weights observations from later in the study [Zhang et al., 2020]. In the batched bandit setting, Zhang et al. [2020] show that the Z-statistics for the OLS estimators computed separately on each batch are jointly asymptotically normal. Standardizing the OLS statistic for each batch effectively adaptively re-weights the observations in each batch.
Hadad et al. [2019] introduce adaptively weighted versions of both the standard augmented-inverse propensity weighted estimator (AW-AIPW) and the sample mean (AWA) for estimating parameters of simple models on data collected with bandit algorithms. They introduce a class of adaptive “variance stabilizing” weights, for which the variance of a normalized version of their estimators converges in probability to a constant. In their discussion section they note open questions, two of which this work addresses: 1) “What additional estimators can be used for normal inference with adaptively collected data?” and 2) How do their results generalize to more complex sampling designs, like data collected with contextual bandit algorithms? We demonstrate that variance stabilizing adaptive weights can be used to modify a large class of M-estimators to guarantee valid inference. This generalization allows us to perform valid inference for a large class of important inferential targets: parameters of models for expected outcomes that are context dependent.
Recently, adaptive weighting has also been used in off-policy evaluation methods for when the behavior policy (policy used to collect the data) is a contextual bandit algorithm [Bibaut et al., 2021, Zhan et al., 2021]. In this literature the estimand is the value, or average expected reward, of a prespecified policy (note this is a scalar value). In contrast, in our work we are interested in constructing confidence regions for parameters of a model for an outcome (that could be the reward)—for example, this could be parameters of a logistic regression model for a binary outcome. We believe in the future there could be theory that could unify these adaptive weighting methods for these different estimands.
An alternative to using asymptotic approximations to construct confidence intervals is to use highprobability confidence bounds. These bounds provide stronger guarantees than those based on asymptotic approximations, as they are guaranteed to hold for finite samples. The downside is that these bounds are typically much wider, which is why much of classical statistics uses asymptotic approximations. Here we do the same. In Section 5, we empirically compare our to the selfnormalized martingale bound [Abbasi-Yadkori et al., 2011], a high-probability bound commonly used in the bandit literature.

heading: 5 Simulation Results
text: In this section, Rt = Yt. We consider two settings: a continuous reward setting and a binary reward setting. In the continuous reward setting, the rewards are generated with mean EP [Rt|Xt, At] = X̃>t ✓ ⇤ 0(P) + AtX̃ > t ✓ ⇤
1(P) and noise drawn from a student’s t distribution with five degrees of freedom; here X̃t = [1, Xt] 2 R3 (Xt with intercept term), actions At 2 {0, 1}, and parameters ✓⇤0(P), ✓ ⇤
1(P) 2 R3. In the binary reward setting, the reward Rt is generated as a Bernoulli with success probability EP [Rt|Xt, At] = [1 + exp( X̃>t ✓⇤0(P) AtX̃>t ✓⇤1(P))] 1. Furthermore, in both simulation settings we set ✓⇤0(P) = [0.1, 0.1, 0.1] and ✓⇤1(P) = [0, 0, 0], so there is no unique optimal arm; we call vector parameter ✓⇤1(P) the advantage of selecting At = 1 over At = 0. Also in both settings, the contexts Xt are drawn i.i.d. from a uniform distribution.
In both simulation settings we collect data using Thompson Sampling with a linear model for the expected reward and normal priors [Agrawal and Goyal, 2013] (so even when the reward is
binary). We constrain the action selection probabilities with clipping at a rate of 0.05; this means that while typical Thompson Sampling produces action selection probabilities ⇡TSt (At, Xt,Ht 1), we instead use action selection probabilities ⇡t(At, Xt,Ht 1) = 0.05 _ 0.95 ^ ⇡TSt (At, Xt,Ht 1) to select actions. We constrain the action selection probabilities in order to ensure weights Wt are bounded when using a uniform stabilizing policy; see Sections 3.2 and 6 for more discussion on this boundedness assumption. Also note that increasing the amount the algorithm explores (clipping) decreases the expected width of confidence intervals constructed on the resulting data (see Section 6).
To analyze the data, in the continuous reward setting, we use least-squares estimators with a correctly specified model for the expected reward, i.e., M-estimators with m✓(Rt, Xt, At) = (Rt X̃>t ✓0 AtX̃>t ✓1)
2. We consider both the unweighted and adaptively weighted versions. We also compare to the self-normalized martingale bound [Abbasi-Yadkori et al., 2011] and the W-decorrelated estimator [Deshpande et al., 2018], as they were both developed for the linear expected reward setting. For the self-normalized martingale bound, which requires explicit bounds on the parameter space, we set ⇥ = {✓ 2 R6 : k✓k2  6}. In the binary reward setting, we also assume a correctly specified model for the expected reward. We use both unweighted and adaptively weighted maximum likelihood estimators (MLEs), which correspond to an M-estimators with m✓(Rt, Xt, At) set to the negative log-likelihood of Rt given Xt, At. We solve for these estimators using Newton–Raphson optimization and do not put explicit bounds on the parameter space ⇥ (note in this case m✓ is concave in ✓ [Agresti, 2015, Chapter 5.4.2]). See Appendix A for additional details and simulation results.
In Figure 4 we plot the empirical coverage probabilities and volumes of 90% confidence regions for ✓⇤(P) := [✓⇤0(P), ✓⇤1(P)] and ✓⇤1(P) in both the continuous and binary reward settings. While the confidence regions based on the unweighted least-squares estimator (OLS) and the unweighted MLE have significant undercoverage that does not improve as T increases, the confidence regions based on the adaptively weighted versions, AW-LS and AW-MLE, have very reliable coverage. For the confidence regions for ✓⇤1(P) based on the AW-LS and AW-MLE, we include both projected confidence regions (for which we have theoretical guarantees) and non-projected confidence regions. The confidence regions based on projections are conservative but nevertheless have comparable volume to those based on OLS and MLE respectively. We do not prove theoretical guarantees for the non-projection confidence regions for AW-LS and AW-MLE, however they perform well across in our simulations. Both types of confidence regions based on AW-LS have significantly smaller volumes than those constructed using the self-normalized martingale bound and W-decorrelated estimator. Note that the W-decorrelated estimator and self-normalized martingale bounds are designed for linear contextual bandits and are thus not applicable for the logistic regression model setting. The confidence regions constructed using the self-normalized martingale bound have reliable coverage as well, but are very conservative. Empirically, we found that the coverage probabilities of the confidence regions based on the W-decorrelated estimator were very sensitive to the choice of tuning parameters. We use 5, 000 Monte-Carlo repetitions and the error bars plotted are standard errors.

heading: 6 Discussion
text: Immediate questions We assume that ratios ⇡stat (At, Xt)/⇡t(At, Xt,Ht 1) are bounded for our theoretical results; this precludes ⇡t(At, Xt,Ht 1) from going to zero for a fixed stabilizing policy. For simple models, e.g., the AW-LS estimator, we can let these ratios grow at a certain rate and still guarantee asymptotic normality (Appendix B.5); we conjecture similar results hold more generally.
Generality and robustness This work assumes that we have a well-specified model for the outcome Yt, i.e., that ✓⇤(P) 2 argmax✓2⇥ EP [m✓(Yt, Xt, At)|Xt, At] w.p. 1. Our theorems use this assumption to ensure that {Wtṁ✓(Yt, Xt, At)}t 1 is a martingale difference sequence with respect to {Ht}t 0. On i.i.d. data it is common to define ✓⇤(P) to be the best projected solution, i.e., ✓0(P) 2 argmax✓2⇥ EP,⇡ [m✓(Yt, Xt, At)]. Note that the best projected solution, ✓⇤(P), depends on the distribution of the action selection policy ⇡. It would be ideal to also be able to perform inference for a projected solution on adaptively collected data.
Another natural question is whether adaptive weighting methods work in Markov Decision Processes (MDP) environments. Taking the AW-LS estimator introduced in Section 3.1 as an example, our conditional variance derivation in Equation (7) fails to hold in an MDP setting, specifically equality (c). However, the conditional variance condition can be satisfied if we instead use weights Wt = {[⇡stat (At, Xt)p sta(Xt)]/[⇡t(At, Xt,Ht 1)PP(Xt|Xt 1, At 1)]}1/2 where PP are the state
transition probabilities and psta is a pre-specified distribution over states. In general though we do not expect to know the transition probabilities PP and if we tried to estimate them, our theory would require the estimator to have error op(1/ p T ), below the parametric rate.
Trading-off regret minimization and statistical inference objectives In sequential decisionmaking problems there is a fundamental trade-off between minimizing regret and minimizing estimation error for parameters of the environment using the resulting data [Bubeck et al., 2009, Dean et al., 2018]. Given this trade-off there are many open problems regarding how to minimize regret while still guaranteeing a certain amount of power or expected confidence interval width, e.g., developing sample size calculators for use in justifying the number of users in a mobile health trial, and developing new adaptive algorithms [Liu et al., 2014, Erraqabi et al., 2017, Yao et al., 2020].
Acknowledgements and Disclosure of Funding
We thank Yash Nair for feedback on early drafts of this work.
Research reported in this paper was supported by National Institute on Alcohol Abuse and Alcoholism (NIAAA) of the National Institutes of Health under award number R01AA23187, National Institute on Drug Abuse (NIDA) of the National Institutes of Health under award number P50DA039838, National Cancer Institute (NCI) of the National Institutes of Health under award number U01CA229437, and by NIH/NIBIB and OD award number P41EB028242. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.
This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE1745303. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Review summaries:
This paper studies how to conduct statistical inference with M-estimators based on data collected from adaptive algorithms. It proposes an adaptively weighted M-estimator with square-root importance weights and proves that under fairly standard assumptions, the estimator (uniformly) has asymptotic normal distribution and it can be used to construct uniformly valid confidence intervals. By studying M-estimation problems, this paper substantially broadens the scope of estimation problems in this recent literature of statistical inference with adaptively collected data. 

The authors analyze square root importance sampling weighted M-estimators over parametric classes for adaptively collected data.   This work is motivated by the fact that the usual unweigthed M-estimators that one would use for i.i.d. data fail to be asymptotically normally distributed under adaptively collected data.  Under appropriate conditions, the proposed weighting scheme makes the resulting M-estimator asymptotically normal, uniformly over the parameter space.  The authors demonstrate their claims through simulations.