RoME: A Robust Mixed-Effects Bandit Algorithm for Optimizing Mobile Health Interventions

Authors: easton  huch
jieru  shi
madeline r abbott
jessica r golbus
alexander  moreno
walter h dempsey
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-11 17:23:58

Abstract: Mobile health leverages personalized and contextually tailored interventions optimized through bandit and reinforcement learning algorithms. In practice, however, challenges such as participant heterogeneity, nonstationarity, and nonlinear relationships hinder algorithm performance. We propose RoME, a **Ro**bust **M**ixed-**E**ffects contextual bandit algorithm that simultaneously addresses these challenges via (1) modeling the differential reward with user- and time-specific random effects, (2) network cohesion penalties, and (3) debiased machine learning for flexible estimation of baseline rewards. We establish a high-probability regret bound that depends solely on the dimension of the differential-reward model, enabling us to achieve robust regret bounds even when the baseline reward is highly complex. We demonstrate the superior performance of the RoME algorithm in a simulation and two off-policy evaluation studies.

Main paper: 
 introduction: 
 Mobile health (mHealth) uses smart devices to deliver digital notification interventions to users. These notifications nudge users toward healthier attitudes and behaviors. Because mHealth applications can monitor and react to users and their environment, they offer the promise of personalized, contextually tailored interventions. In practice, achieving this promise requires bandit or reinforcement learning (RL) algorithms that can accurately learn mHealth intervention effects, including how they vary by user, over time, and based on context. In this regard, contextual bandit algorithms are appealing due to strong empirical performance in other settings, their ability to customize intervention decisions based on changing contexts, and their simplicity and extensibility relative to full RL algorithms.One unique aspect of mHealth (Greenewald et al., 2017) is the presence of a "do-nothing" action in which the mHealth app does not intervene. As mHealth rewards typically involve human decision making (e.g., whether the user chooses to exercise), the distribution of the the corresponding baseline reward-the reward under the do-nothing action-typically evolves in a complex, nonstationary fashion as time progresses and contexts change. In contrast, the treatment effects-the difference in the expected value of the reward relative to the do-nothing action-tend to be much more stable and can be adequately modeled using classical stationary models, such as linear models (Robinson, 1988) . Greenewald et al. (2017) introduced an action-centered contextual bandit algorithm with sublinear regret in this setting by replacing the observed reward, R t , with a "pseudo-reward," Rt , where t indexes the time points. Denoting the binary action as A t ∈ {0, 1} with A t = 0 corresponding to the do-nothing action, the pseudo-reward for A t can be written asEQUATIONwhere π t is the (known) probability of treatment assignment from the Thompson-sampling algorithm.We took an additional step to derive an equivalent expression (second equality), which reveals that Rt is proportional to an inverse-probability-weighted (IPW) estimator. Crucially, Rt provides an unbiased estimate of the treatment effect regardless of the nonstationarity of the baseline outcome.Although Greenewald et al. (2017) 's analysis shows sublinear regret in nonstationary settings, our simulations show that other methods can achieve lower regret over finite time horizons. The reasons include its failure to address treatment effect heterogeneity, pool information across users, and model the baseline outcome, which leads to highly variable pseudo-rewards and slow learning. We propose to generalize their method to overcome these limitations by 1. imposing hierarchical structure in the treatment effect model with shared fixed effects and random effects for users and time (i.e., a mixed-effects model), 2. efficiently pooling across users and time via nearest-neighbor regularization (NNR), and 3. denoising rewards with flexible supervised learning models via the debiased machine learning (DML) framework of Chernozhukov et al. (2018) .We name the resulting method RoME to highlight that it is a Robust Mixed-Effects algorithm. The primary contributions are (a) the RoME method, (b) a high-probability regret bound that relies solely on the dimension of the differential-reward model, which is typically much smaller than that of the complex baseline reward model, and (c) empirical comparisons demonstrating the superior performance of RoME in simulation and two real-world mHealth studies. 
related work: 
 Related work in statistics considers the estimation of treatment effects with longitudinal data. Cho et al. (2017) present a semiparametric random-effects method for estimating subject-specific static treatment effects. Qian et al. (2020) discusses the statistical challenges of applying linear mixedeffects models to mHealth data and shows that the resulting estimates may lack a causal interpretation.In the bandit literature, several works address a growing pool of users. Ma et al. (2011) introduced NNR in recommender systems using homophily principles-similar nodes are more likely to be connected than dissimilar ones (McPherson et al., 2001) . Cesa-Bianchi et al. (2013) adapted NNR for bandit settings, improving regret. Subsequent improvements focused on scalability and algorithm modifications for stronger regret bounds (Vaswani et al., 2017; Yang et al., 2020) .Other bandit approaches explicitly address the longitudinal setting in which treatment effects may evolve over time. The intelligentpooling method of Tomkins et al. (2021) employs a Gaussian linear mixed-effects model with both user-and time-specific random effects. Hu et al. (2021) provides a related approach based on generalized linear mixed effects models. Unlike our approach, both methods require the (generalized) linear conditional mean model to be correctly specified. Aouali et al. (2023) and Lee et al. (2024) also propose bandit algorithms with mixed effects; methodologically, however, their use of mixed effects is dissimilar from that of the works above and our own. Aouali et al. (2023) employs random effects to relate treatments within categories (e.g., movies within a genre), and Lee et al. (2024) considers the case in which a single, discrete action leads to a multivariate reward that is affected by a user-specific random effect. In contrast, the other approaches listed above (and our own) do not assume any relationships among the actions, and they associate each discrete action with a single scalar-valued reward in an iterative fashion.Other work develops semiparametric methods that lead to sublinear regret without requiring correct specification of the conditional mean reward. Krishnamurthy et al. (2018) and Kim & Paik (2019) propose improvements to the approach of Greenewald et al. (2017) that lead to stronger regret bounds. Kim et al. (2021 Kim et al. ( , 2023) ) develop a doubly robust approach for both linear and generalized linear contextual bandits. However, in contrast to Tomkins et al. (2021) , Hu et al. (2021) , and our approach, these semiparametric methods are not immediately applicable to longitudinal settings.3 Setting 
problem statement: 
 Recognizing the connection between mHealth policy learning and contextual bandit methods, we now consider a contextual multi-armed bandit environment with one control arm (the do-nothing action) denoted by a = 0 and q non-baseline arms corresponding to different actions or treatments. Individuals are indexed by i = 1, 2, . . . and decision points are indexed by t = 1, 2, . . .. For each individual i at time t, the algorithm observes a context vector S i,t ∈ S, chooses an action A i,t ∈ [q] := {0, . . . , q}, and receives a reward R i,t ∈ R. The differential reward describes the difference in expected rewards between choosing a non-baseline arm and the control arm, defined as ∆ i,t (s, a):= E [R i,t |S i,t = s, A i,t = a] -E [R i,t |S i,t = s, A i,t = 0].We denote the observation history up to time t for individual i as H i,t := (S i1 , A i1 , R i1 , . . . , S it ).Actions are selected according to stochastic policies, π i,t : H i,t × S → P([q]). We use π i,t (a|s) to denote the probability of action a ∈ [q] given current context s ∈ S for a fixed (implicit) history. As in Greenewald et al. (2017) , we assume that the probability of the control action is bounded: Assumption 1. There exists 0 < π min , π max < 1 s.t π min < π i,t (0|a), 1 -π i,t (0|a) < π max ∀i, t, a.As noted in Greenewald et al. (2017) , this assignment probability controls the number of messages that participants receive. It ensures that participants receive sufficient messages that they do not disengage, but not so many that they are overwhelmed or fatigued. We now defineA * = argmax a∈[q] E [R i,t |S i,t , A i,t = a] , Ā * = argmax a∈[q]\{0} E [R i,t |S i,t , A i,t = a] ,the optimal treatment arm (including control) and the alternative arm, respectively. Respecting the restrictions in Assumption 1, we can now express the optimal policy, π * i,t , in two cases:1. A * = 0: π * i,t (0|S i,t ) = π max and π * i,t ( Ā * |S i,t ) = 1 -π max . 2. A * > 0: π * i,t (A * |S i,t ) = 1 -π min and π * i,t (0|S i,t ) = π min .These cases represent the appropriate boundary condition given the optimal action; we maximize the probability of the do-nothing action if A ⋆ = 0 and minimize it if A ⋆ ̸ = 0, selecting the next-best action with the highest possible probability. In Section 5.2, we define regret relative to this optimal policy, π * i,t . Lastly, we consider a study design that progresses in stages, where individuals enter the study sequentially (Friedman et al., 2010) , as illustrated in Figure 1 . Staged recruitment serves two purposes in our analysis: it increases the fidelity of the theoretical setting by mimicking real-world recruitment strategies, and it allows us to estimate changes in the differential reward over time. Concretely, we first observe individual i = 1 at time t = 1. Then we observe individuals i = (1, 2) at times t = (2, 1). Subsequent stages k involve individuals i ≤ k, each having had k -i + 1 observations, respectively. Define O k = {(i, t) : i ≤ k & t ≤ k + 1 -i} as the set of all observations in stage k. A notation summary is presented in Appendix F.Following Greenewald et al. (2017) , we assume that the differential reward is linear:EQUATIONwhere x(s, a) ∈ R p×1 is a feature vector depending on the context and action. We allow the baseline reward, g t (s), to be a nonlinear and nonstationary function of context and time. Combined, these two who have yet to be recruited. For simplicity, we assume one participant is recruited at each stage.assumptions imply a partially linear model for the conditional mean reward:EQUATIONwhere δ a>0 is an indicator of a non-baseline action. By incorporating both the linear predictor x(s, a) ⊤ θ i,t δ a>0 and the nonlinear baseline g t (s), we can maintain interpretability of the action effects while also allowing flexibility in specifying the observed rewards. The individual rewards are then realized according to the modelR i,t = x(s, a) ⊤ θ i,t δ a>0 + g t (s) + ϵ i,t, where we make the following standard assumption (Abeille & Lazaric, 2017) regarding ϵ i,t : Assumption 2. ϵ i,t is conditionally mean zero (i.e., E[ϵ i,t |H i,t ] = 0) and sub-Gaussian with varianceσ 2 : E[exp(ηϵ i,t )|H i,t ] ≤ exp(η 2 σ 2 /2) for η > 0.Assumption 2 ensures that large errors are sufficiently rare and that Rt is proportional to an unbiased estimate of ∆ i,t (s, a). However, Rt suffers from high variance, which leads to suboptimal performance in practice. Our proposed algorithm RoME, however, relies on an improved estimator, Rf i,t (s, ā), that is also unbiased but achieves lower variance by denoising the rewards according to a working model, f i,t (s, a), for the conditional mean r i,t (s, a)EQUATIONDebiasing termEQUATION)where ā is a preselected non-baseline arm; i.e., we condition on A i,t ∈ {0, ā}. The first component is a model-based prediction of the differential reward, and the second component is a debiasing term that uses IPW to correct for potential errors in the specification of f . This improved estimator is a form of pseudo-outcome, a common technique in the causal inference literature (Bang & Robins, 2005; Kennedy, 2023; Nie & Wager, 2021) . We refer to it as a Doubly Robust Differential Reward because it produces an unbiased estimate of ∆ i,t (s, a) as long as either π i,t or f i,t is correctly specified (see proof in Appendix E.1); though, in practice, the benefit stems from (a) the robustness to misspecification of f (because π i,t is known) and (b) the variance reduction introduced by denoising the outcomes. Lemma 5 and Remark 2 in Appendix E.2 provide a proof and discussion of the variance reduction. In the following, we abbreviate Rf i,t (s, a) as Rf i,t , with the context and action implied. 
algorithm components: 
 We now construct Rf i,t , introduce the weighted least-squares loss function used by RoME to estimate parameters, and use nearest-neighbor regularization to pool information across users.To construct Rf i,t , we fit a working model f i,t for the conditional mean reward via the DML framework, which requires three main conditions: (a) Neyman orthogonality, which makes the pseudo-outcome Rf i,t robust to misspecification of f i,t , (b) sufficiently fast convergence of the working conditional mean model f i,t , and (c) independence of f i,t and the observed results. (a) is ensured by the structure of the doubly robust differential reward. (b) requires us to introduce the following assumption, which states that the mean squared error (MSE) of the working conditional mean is asymptotically smaller than the inverse square root of the number of stages: Assumption 3. The working model,f i,t , satisfies E p(s,ā) {r i,t (s, ā) -f i,t (s, ā)} 2 = o P (k -1/2 ) and E p(s) {r i,t (s, 0) -f i,t (s, 0)} 2 = o P (k -1/2 ). Further, |f i,t | ≤ 2M .Examples of supervised learning methods that converge at this rate (or faster) include parametric models of fixed dimension, local mean smoothers with small dimension (Wasserman, 2006, Chapter 5), and random forests and feedforward neural networks under certain regularity conditions (Biau, 2012; Farrell et al., 2021) . Assumption 3 guarantees that RoME achieves better asymptotic variance than action-centered (AC) bandits Greenewald et al. (2017) , resulting in the tighter regret bound in Section 5.2. Because our differential reward in ( 4) is robust to misspecification of f i,t , our regret bound can be adapted to settings in which Assumption 3 does not hold, provided f i,t converges to a bounded function at the rate specified in Assumption 3; this setting would result in a weaker bound of the same asymptotic order. The final requirement of independence, (c), can be accomplished via sample-splitting techniques. Below, we present two such options.Option 1: The first option exploits the fact that outcomes across different users are independent. Consequently, we can perform sample splitting as follows.Step 1: Randomly assign each user i to one of J folds. Let I j denote the set containing user indexes for the j-th fold and let I ∁ j denote its complement.Step 2: For each fold, use a supervised learning algorithm to estimate the working model for r i,t (s, a) denoted f (j) i,t (s, a) using I ∁ j .Step 3: Construct the pseudo-outcomes using (4). Option 2: Because we do not expect an adversarial environment in mHealth, we may be willing to assume that the errors, ϵ i,t , are independent across time. In this case, we can adapt the procedure given above by randomly assigning data to folds at the level of (i, t). In this case, I j , contains the (i, t) pairs assigned to fold j. This option is more powerful than option 1 in that it enables us to learn heterogeneity across users in the nuisance function, f i,t , leading to greater variance reduction.Lastly, we note two strategies that can reduce the computational demands of sample splitting. The first is to employ estimators that satisfy a leave-one-out stability condition, such as bagged estimators that use subsampling (Chen et al., 2022) ; these estimators eliminate the need to perform sample splitting. The second is to fit f i,t in an online fashion.Having formed Rf i,t (s, ā), we now explain how to estimate the parameters, θ i,t ∈ R p , determining the differential reward, ∆ i,t (s, a). We assume that θ i,t consists of three components:Assumption 4. θ i,t = θ shared + θ user i + θ time t for all i, t.The parameter θ shared is a fixed effect shared across all i, t. In contrast, θ user i and θ time t are random effects specific to a given user and time point, respectively. We place these parameters in a single column vector, θ, as vec(θ shared , θ user 1 , . . . , θ user K , θ time 1 , . . . , θ time K ), where K is the total number of stages. We then form a corresponding feature vector, ϕ i,t = ϕ(x i,t ), where ϕ inserts x i,t into the positions corresponding to θ shared , θ user i , and θ time t with zeros in all other locations; i.e., letting p := dim(x i,t ), the feature vector ϕ(x i,t ) ∈ R (2K+1)p is defined as follows:EQUATIONWe could then estimate θ i,t via a weighted least squares (WLS) regression that minimizes the following objective function:EQUATIONwhere σ2 i,t = π i,t (0|S i,t ) • (1 -π i,t (0|S i,t )). However, ℓ WLS (θ) is over-parameterized and does not take advantage of available network information. The next section modifies ℓ WLS (θ) by adding regularization, producing the final loss function used in RoME.We assume access to network information relating neighboring users and time points. These networks are defined by graphs G user = (V user , E user ) and G time = (V time , E time ), where V denotes vertices and E, edges. In practice, these networks can be formed via known clusters (e.g., shared schools, companies, geographic locations) or similar covariates. In the case of G time , a reasonable choice is to construct a graph of neighboring sequential time points; i.e., t = 1 ↔ t = 2 ↔ t = 3 . . . . We assume the following bounds on the random effects:Assumption 5. There exists D user , D time , B user , B time , B shared ∈ R + such that ∥θ user i -θ user j ∥ 2 2 ≤ D user , ∥θ time i -θ time j ∥ 2 2 ≤ D time , ∥θ user i ∥ ≤ B user ,and ∥θ time i ∥ ≤ B time and ∥θ shared ∥ ≤ B shared for all i, j.The assumption 5 suggests that we may be able to improve our estimates of θ by regularizing neighboring values of θ user i and θ time i toward each other. This can be accomplished by using an L 2 Laplacian penalty. We illustrate the idea using G user ; the application to G time is similar. First, we form the incidence matrix, Q ∈ R K×K , with the entry Q v,e corresponding to the v-th vertex (user) and the e-th edge. Denote the vertices of e as v i and v j with i > j. Q v,e is then equal to1 if v = v i , -1 if v = v j ,and 0 otherwise. The Laplacian matrix is then defined as L = QQ ⊤ . Similar to Yang et al. (2020) , we can then form a network cohesion penalty across users as follows:tr(Θ ⊤ user L user Θ user ) = (i,j)∈Euser ∥θ user i -θ user j ∥ 2 2 ,where Θ user := (θ user 1 , . . . , θ user K ) ⊤ ∈ R K×p and L user is the Laplacian matrix for users. The penalty is small when θ user i and θ user j are close for connected users. We employ a shared regularization hyperparameter λ, for G user and G time , and include a standard L 2 penalty on θ with hyperparameter γ. The full penalization matrix V 0 ∈ R (2K+1)p×(2K+1)p isEQUATIONwhere L user ⊗ := L user ⊗ I p , the Kronecker product of L user with I p . We then adapt ℓ WLS (θ) to include a corresponding penalty term:EQUATIONIn Section 5, we show that (8) leads to a Thompson sampling algorithm for action selection. 
thompson sampling algorithm: 
 We now describe the Thompson sampling procedure and provide a high-probability regret bound.The minimizer of the weighted regularized least squares loss in (8) can be expressed as θ = V -1 b, whereV = V 0 + n i=1 K-i+1 t=1 σ2 i,t ϕ i,t ϕ ⊤ i,t , b = n i=1 K-i+1 t=1 Rf i,t σ2 i,t ϕ i,t .In analogy to Bayesian linear regression, the penalized and weighted Gram matrix, V , plays the role of a (scaled) posterior precision matrix under a multivariate Gaussian prior for θ. This motivates an algorithm in which we randomly perturb the point estimate, θ, scaling by an inverse square-root matrix, V -1/2 , times a mean-zero random deviate, η. The perturbed estimate, θ, can be written as θ = θ + β(δ)V -1/2 η, where letting W > 0 be the max number of neighbors a time point or user has, β(δ) is defined asEQUATIONEQUATION)The value v is defined in Corollary 4 and can be set to a sufficiently large constant. The hyperparameter δ determines the probability with which the regret bound holds. The random deviate, η, is drawn from a distribution, D T S , satisfying the technical conditions given in Definition 1 in Appendix D.In practice, we recommend setting D T S to a multivariate Gaussian distribution or t-distribution to simplify calculations. Given θ, we then compute Āi,t = argmax a∈[q]\{0} ϕ ⊤ i,t θ and randomly select A i,t = 0 with probabilityEQUATIONor A i,t = Āi,t with probability 1 -π i,t (0| Āi,t ). The action selection is summarized in Algorithm 1.Algorithm 1 Action selection input V , b, π min , π max , ϕ i,t , β(δ) θ = V -1 b Sample η ∼ D T S θ = θ + β(δ)V -1/2 η, with β(δ) given in (9) Āi,t = argmax a∈[q]\{0} ϕ ⊤ i,t θCalculate π i,t (0| Āi,t ) according to (11) With probability π i,t (0| Āi,t ), play action 0; otherwise, play action Āi,tHaving detailed the action selection algorithm, we now summarize RoME in Algorithm 2. The algorithm includes an outer loop that iterates through the stages from k = 1 to k = K. Within each stage, we handle users sequentially, iteratively selecting actions according to Algorithm 1, observing rewards, and updating V and b accordingly.Algorithm 2 RoME Initialize V 0 according to (7) V ← V 0 b ← 0 for k = 1, . . . , K do for (i, t) = (1, k),(2, k -1), . . . , (k, 1) do Choose an arm according to Algorithm 1 using V and b Observe rewards, R i,tConstruct pseudo-rewards, Rf (m) i,t , as explained in Section 4.1 Set ϕ = ϕ(x i,t ) as described in (5) Update weighted Gram matrix:V ← V + σ2 i,t ϕ i,t ϕ ⊤ i,t Update feature-outcome products: b ← b + σ2 i,tRf i,t ϕ i,t end for end forWe first define regret, Regret K , in terms of the following average across stages:K k=1 1 k (i,t)∈O k \O k-1 π ⋆ i,t ( Ā * i,t |S i,t ) • x(S i,t , Ā * i,t ) ⊤ θ * i,t -π i,t ( Āi,t |S i,t ) • x(S i,t , Āi,t ) ⊤ θ * i,t .The inner sum calculates regret (the difference between the expected reward of the optimal action and the chosen action) weighted by the probability of the action, following Greenewald et al. (2017) , for all individuals at stage k. We then average these values to account for the growing number of pairs (i, t) as the stages progress. The cumulative regret is the sum of these averaged values over all stages.The rewards associated with the baseline action are excluded since ∆ i,t (S i,t , 0) = 0. Regret K is a form of pseudo-regret because the randomness attributed to {ϵ i,t } K t=1 has been eliminated (Audibert et al., 2003) . Before presenting the regret bound, we introduce one more assumption: Assumption 6. ∥x(s, a)∥ ≤ 1 for all s and a. Also, there exists a known value M ∈ R + bounding both ∥θ i,t ∥ and |g t (s)| for all i, t, and s. Further, |f i,t | ≤ 2M .Theorem 1. Under Assumptions 1-6, with probability at least 1 -δ, Regret K satisfiesO β K (δ ′ ) + γ K (δ ′ ) 1 + 4 p • log(K) cKpwhere δ ′ = δ/4K and min{π(0|s), 1 -π(0|s)} > 1/c. The exact regret bound is given and both β K and γ K are defined in Appendix E.Theorem 1's proof is in Appendix E. It resembles Abeille & Lazaric (2017) , albeit with adjustments to address (a) the nonlinear baseline and (b) the increasing pool of users and time points. The nonlinear baseline requires deriving the new sub-Gaussian variance factor for the centered pseudo-reward and adapting classic contextual bandit theory (Abbasi-Yadkori et al., 2011) to our weighted least squares setting. The increasing pool of users and time points requires modifying the Abeille & Lazaric (2017) proof with a careful use of stages and the Cauchy-Schwarz inequality to decouple the weighted feature norms from inverse stage scaling, which allows us to apply standard techniques to the sum of norms and bound the sum of inverse stages by the harmonic number.DML is particularly advantageous for β K (δ ′ ) and γ K (δ ′ ), which rely on the convergence rate of the model f towards the true mean differential reward r, as discussed further in Appendix E.2. Another interesting aspect of the regret bound is its dependence solely on the dimension, p, of the differential-reward model, rather than the dimension of the baseline reward. This characteristic allows us to achieve strong regret bounds even in cases where the baseline reward is highly complex.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: The primary simulation results can be entirely reproduced using the provided code on Github. However, the Valentine study data, used for our real-world analysis, is not publicly accessible. Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 This section presents results from applying RoME in a simulated mHealth study with staggered recruitment and an off-policy comparison study for cardiac rehabilitation mHealth interventions. The simulations were implemented using Python and the results were generated using individual compute nodes with two 3.0 GHz Intel Xeon Gold 6154 processors and 180 GB of RAM. Case studies were implemented using R 4.2.2 and results were generated on a cluster composed of individual compute nodes with 2.10 GHz Intel Xeon Gold 6230 processors and 192 GB of RAM.In this section, we compare RoME to four competing methods in simulation. We implemented RoME using Option 2 with a bagged ensemble of stochastic gradient trees (Gouk et al., 2019; Mastelini et al., 2021) trained online via the River library (Montiel et al., 2021) . We programmed the linear algebra operations using SuiteSparse (Davis & Hu, 2011) to take advantage of the sparse nature of the Laplacian matrices and ϕ i,t from Algorithm 2.Our competing baselines are (a) Standard: Standard Thompson sampling for linear contextual bandits, (b) AC: the Action-Centered contextual bandit algorithm (Greenewald et al., 2017) , (c) IntelPooling:The intelligentpooling method of Tomkins et al. (2021) fixing the variance parameters close to their true values, and (d) Neural-Linear: a method that uses a pre-trained neural network to transform the feature space for the baseline reward (similar to the Neural Linear method of Riquelme et al. (2018) ).We compare these methods under three settings: Homogeneous Users, Heterogeneous Users, and Nonlinear. The first two involve a linear baseline model and time-homogeneous parameters, with the second setting having distinct user parameters. The third setting, designed to mirror the challenges of mHealth studies, includes a nonlinear baseline and both user-and time-specific parameters. For each setting, we simulate 200 stages following the staged recruitment regime depicted in Figure 1 , and we repeat the full 200-stage simulation 50 times. Appendix A.1 provides details on the setup and a link to our implementation. Appendices A.2-A.5 offer more results, including additional method comparisons, a rectangular data array simulation, hyperparameter sensitivity analyses, and pairwise statistical comparisons. The latter show that RoME outperformed each competitor in at least 48 of 50 repetitions in the Nonlinear setting, indicating even higher statistical confidence than Figure 2 suggests. RoME substantially outperforms the competing algorithms in the additional comparisons and is several orders of magnitude faster than would be required in a standard mHealth study, producing over 20,000 decisions in as little as 45 seconds (see Table 2 in Appendix A.2).In this section, we compare RoME to the above algorithms via off-policy evaluation using data from the Valentine Study (Jeganathan et al., 2022) , a prospective, randomized-controlled, remotely administered trial designed to evaluate an mHealth intervention to supplement cardiac rehabilitation for low-and moderate-risk patients. In the analyzed dataset, participants were randomized to receive or not receive contextually tailored notifications promoting low-level physical activity and exercise throughout the day. The left of Figure 3 shows the estimated improvement in the average reward over the original constant randomization, averaged over stages (K = 120) and participants (N=108).We see that RoME achieved the highest average reward; its average performance (across bootstrap replications) is higher than the 75th percentile of all other methods.To further analyze the improvement of RoME relative to the other methods, we test whether the proposed algorithm significantly improves cumulative rewards using paired t-tests with one-sided alternative hypotheses. The null hypothesis (H 0 ) for these tests is that the two algorithms being compared achieve the same average reward. The alternative hypothesis (H 1 ) is that the algorithm listed in the column achieves higher average rewards than the algorithm listed in the row. The right of Figure 3 displays the p-values obtained from these pairwise t-tests. The dark shade of the last column indicates that the proposed RoME algorithm achieves significantly higher rewards than the other five competing algorithms (p ≤ 0.01 for all pairwise comparisons). The results imply that RoME would have achieved a 3.5% increase in step count relative to the constant randomization policy that was actually implemented in this study. This effect translates to an increase of 140 steps per day, given a baseline of 1,000 steps per hour and four one-hour measurement windows. See Appendix B for additional implementation details. Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: Algorithm 2 requires data splits, and we devoted two paragraphs on Page 5 to discussing options for this split. Additionally, we clearly listed all of the hyperparameters involved in forming objective functions and provided guidance on how to make suitable choices for them. Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
discussion: 
 This paper introduces RoME, a robust mixed-effects contextual bandit algorithm for mHealth studies. RoME adapts DML and network cohesion penalties to dynamic settings, enabling researchers to efficiently pool information across users and over time. In addition to the methodological contribution, we also prove a high-probability regret bound in a challenging asymptotic regime that involves a growing pool of users and time points. Our implementation of RoME is several orders of magnitude faster than would be required in practice and could easily be adapted for use in future mHealth studies.We see several promising directions for improving RoME in future work. One such direction might consider data-adaptive strategies for setting the hyperparameters or constructing a network in settings where the network is not known a priori. Another interesting direction is to address the long-term impacts of mHealth treatments, such as accumulating treatment fatigue. Future work could also consider computational improvements that would enable large-scale deployment of RoME, which would be especially relevant in nonclinical settings. This appendix provides details regarding the setup of our simulation study, additional results not presented in the main paper, and results from additional simulations employing different competitors and/or simulation designs. We summarize the differences between the competing algorithms in Table 1 .User × √ √ N/A × × IntelPooling √ √ √ √ × × × Neural-Linear √ × × N/A N/A √ √ Standard √ × × N/A N/A × × AC √ × × N/A N/A √ ×Table 1 : This table compares the methods tested in the simulation study based on the design components listed in the columns, with √ indicating that the method includes the component listed in the column, × indicating that it does not, and N/A indicating that the component is not applicable. This table can be used to investigate the effect of certain model components on cumulative regret; for example, comparing RoME and RoME-BLM helps us understand the impact of modeling the nonlinear baseline-the only column that differs between these two methods.The code for the simulation study is fully containerized and publicly available at https://github. com/eastonhuch/RoME. The simulation study involves a generative model of the following form for user i at time t:R it = g(S it ) + x(S it , A it ) ⊤ θ it + ϵ it , ϵ it ∼ N (0, 1)Here S it = (s 1 , s 2 ) ∈ R 2 is a context vector, with both dimensions iid ∼ U (-1, 1). We set x(s, a) = a (1, s 1 , s 2 ). For simplicity, we set g to a time-homogeneous function. The specific nature of the function varies across the following three settings mentioned in Section 6.1:• Homogeneous Users: Standard contextual bandit assumptions with a linear baseline and no user-or time-specific parameters. The linear baseline is g(S it ) = 2 -2s 1 + 3s 2 , and the causal parameter is θ it = (1, 0.5, -4) such that the optimal action varies across the context space. • Heterogeneous Users: Same as the above but each user's causal parameter has iid N (0, 1) noise added to it. • Nonlinear: The general setting discussed in the paper with a nonlinear baseline, user-specific parameters, and time-specific parameters. The base causal parameter and user-specific parameters are the same as in the previous two settings. The nonlinear baseline and timespecific parameter are shown in Figure 4 .We assume that the data are observed via a staged recruitment scheme, as illustrated in Figure 1 . For computational convenience, we update parameters and select actions in batches. If, for instance, we observe twenty users at a given stage, we update our estimates of the relevant causal parameters and select actions for all twenty users simultaneously. This strategy offers a slight computational advantage with limited implications in terms of statistical performance.For simplicity, we assume that the nearest-neighbor network is known and set the relevant hyperparameters accordingly. We took care to set the hyperparameters such that RoME performs a similar amount of shrinkage compared to other methods, especially IntelPooling, which effectively uses a separate penalty matrix for users and time. To accomplish this, we used the same penalty matrices for RoME and IntelPooling, effectively generalizing the scalar γ penalty discussed in the main paper.For the purposes of calculating the constant B given in Algorithm 2, we used a separate value of γ for both users and time and set them equal to the maximum eigenvalue of the corresponding penalty matrix. We use 5 neighbors within the DML methods and set the other hyperparameters as follows: λ = 1, and δ = 0.01, and v k = 1 (i.e., v k is constant).For the Neural-Linear method, we generate a 200 × 200 array of baseline rewards to train the neural network prior to running the bandit algorithm. Consequently, the results shown in the paper for Neural-Linear are better than would be observed in practice because we allowed the Neural-Linear method to leverage data that we did not make available to the other methods. This setup offers the computational benefit of not needing to update the neural network within bandit replications, which substantially reduces the necessary computation time.Aside from the input features used, the Neural-Linear method has the same implementation as the Standard method. The Neural-Linear method uses the output from the last hidden layer of a neural network to model the baseline reward. However, we use the original features (the state vectors) to model the advantage function because the true advantage function is, in fact, linear in these features.Our neural networks consisted of four hidden layers with 10, 20, 20, and 10 nodes, respectively. The first two employ the ReLU activation function (Nair & Hinton, 2010) while the latter two employ the hyperbolic tangent. We chose to use the hyperbolic tangent for the last two layers because Snoek et al. (2015) found that smooth activation functions such as the hyperbolic tangent were advantageous in their neural bandit algorithm. The loss function was the mean squared error between the neural networks' output and the baseline reward on a simulated data set. We trained our networks using the Adam optimizer (Kingma & Ba, 2014) with batch sizes of 200 for between 20 and 50 epochs. We simulated a separate validation data set to check that our model had converged and was generating accurate predictions.Figure 5 compares the true baseline reward function in the nonlinear setting (left) to that estimated by the neural network (right). We see that the neural network produced an accurate approximation of the baseline reward, which helps explain the good performance of the Neural-Linear method relative to other baseline approaches, such as Standard.We include the Neural-Linear method primarily to demonstrate that correctly modeling the baseline is not sufficient to ensure good performance. In the mHealth, algorithms should also be able to (1) efficiently pool data across users and time and (2) leverage network information. The Neural-Linear method satisfies neither of these criteria. Note that a neural network could be used to model the baseline rewards as part of our algorithm. Future work could consider allowing the differential rewards themselves to also be complex nonlinear functions, which could be accomplished by combining our method with Neural-Linear. We leave the details to future work.Our results in the main paper indicate that RoME outperforms existing methods, especially in the complex longitudinal setting that we are targeting. These positive results beg the question: What aspects of RoME contribute most to its superior performance?In this section, we describe an additional simulation study designed to answer this question. In it, we compare RoME (as implemented in the main paper) to three additional comparison methods: RoME-BLM, RoME-SU, and NNR-Linear. Each comparison isolates particular aspects of our algorithm to understand its contribution to RoME's performance.RoME-BLM is a version of RoME that employs Bagged Linear Models as the baseline supervised learning algorithm. Unlike the method implemented in the main paper (which uses bagged stochastic gradient trees), this baseline model is misspecified in the nonlinear setting; i.e., it cannot accurately model the nonlinear baseline reward function. Consequently, this comparison will allow us to gauge the impact of the supervised learning method used in our method. In particular, it should allow us to validate whether the algorithm is robust to baseline misspecification as the theory suggests.RoME-SU is a Single-User version of our algorithm. Rather than having a separate θ i for each user, this method estimates a single advantage function that is shared across all users. It does, however, employ DML, time effects, and nearest-neighbor regularization (for the time effects). Consequently, this comparison will help us understand the impact of the user-specific parameters in our simulation.NNR-Linear is similar to the existing network-cohesion bandit algorithms described in the main paper. It includes a distinct θ i for each user and regularizes them toward each other using a Laplacian penalty. It does not, however, use DML or time effects. As a result, this comparison isolates the impact of these two factors and should provide evidence that RoME is more broadly applicable than existing network-cohesion bandit algorithms.Figure 6 displays the cumulative regret as a function of the stage for these methods across the three settings used in the main paper. We see that all methods perform similarly in panel (a), the setting with homogeneous users. However, in panel (b) we see that the regret achieved by RoME-SU begins to resemble a straight upward-sloping line because it cannot appropriately model the user-specific advantage functions.In panel (c), we observe a much wider spread of performance, with RoME and RoME-BLM performing similarly and substantially outperforming the other methods. The fact that RoME and RoME-BLM perform similarly indicates that our algorithm is, in fact, robust to misspecification as expected. The comparison to NNR-Linear reveals that the DML and time effects-the new components of our method compared to existing network-cohesion bandit algorithms-do contribute meaningfully to RoME's superior performance. We found in additional simulations studies not reported here that the gap between NNR-Linear and our method widens as the magnitude of the time effects increases. Table 3 shows pairwise comparisons between methods across the three settings. The individual cells indicate the percentage of repetitions (out of 50) in which the method listed in the row outperformed the method listed in the column. The asterisks indicate p-values below 0.05 from paired two-sided ttests on the differences in final regret. The Avg column indicates the average pairwise win percentage.RoME and RoME-BLM perform well across all three settings and the difference between them is statistically indistinguishable. These methods perform about equally well compared to IntelPooling and NNR-Linear in the first setting, but they perform better in the other two. In the Nonlinear setting in particular, RoME and RoME-BLM substantially outperform all other methods, and the pairwise differences are statistically significant at the 5% level. Table 3 : Pairwise comparisons between methods in the three settings of the main simulation. Each cell indicates the percent of repetitions (out of 50) in which the method listed in the row outperformed the method listed in the column in term of final regret. Asterisks indicate p-values below 0.05 from paired two-sided t-tests on the differences in final regret. RoME and RoME-BLM perform well in all three settings, and their final regret is statistically indistinguishable. They substantially outperform all other methods in the Nonlinear setting.The main simulation involves simulating data from a triangular data array. At the 200-th (final) stage, the algorithm has observed 200 rewards for user 1, 199 rewards for user 2, and so on.In this section, we simulate actions and rewards under a rectangular array with 100 users and 100 time points. Although we still follow the staged recruitment regime depicted in Figure 1 , at stage 100 we stop sampling actions and rewards for user 1; at stage 101 we stop sampling for user 2; and so on until we have sampled 100 time points for all 100 users. Aside from the shape of the data array, the setup is the same for this simulation as for the main simulation. Nonlinear settings using a rectangular array of data in which we observe 100 time points for 100 users in a stagewise fashion as depicted in Figure 1 . We observe the same performance ordering as in Figure 6 , but here the relative differences are even larger.The cumulative regret for these methods as a function of time-not stage-is shown in Figure 7 .Qualitatively, the results are nearly identical to those from the main simulation (compare to Figure 2 ). The differences in final regret are even larger in this simulation than they were in the main simulation, presumably because this study exhibits greater variation in time effects due to the balance across time.Figure 8 displays comparisons similar to those shown in Figure 6 but for the rectangular simulation. We again see that the results are qualitatively similar to the main simulation. The methods perform similarly in the Homogeneous Users setting. RoME-SU performs poorly in the Heterogeneous Users setting, but the other methods perform similarly. Finally, RoME and RoME-BLM substantially outperform the other methods in the nonlinear setting. Again, the differences in performance are even larger in this simulation than in the main simulation. Table 4 : Pairwise comparisons between methods in the three settings of the simulation with a rectangular array of data. As in Table 3 , each cell indicates the percent of repetitions (out of 50) in which the method listed in the row outperformed the method listed in the column in term of final regret. Asterisks indicate p-values below 0.05 from paired two-sided t-tests on the differences in final regret. The qualitative results are similar to those of Table 3 . RoME and RoME-BLM perform well across all three settings and substantially outperform all other methods in the Nonlinear setting.simulation. In particular, in this simulation study, RoME and RoME-BLM outperform all other methods in 100% of repetitions in the nonlinear setting.In the simulation study, we chose hyperparameters for the fairest comparison possible, employing the same regularization parameters across methods. To assess robustness, we performed eight sensitivity analyses that alter hyperparameters for our methods while fixing those of the other methods; this approach gives the competing algorithms an advantage. Under these alternative hyperparameters, RoME and RoME-BLM still dramatically outperform the other methods in the Nonlinear Setting. We provide detailed results below.Three analyses led to no meaningful changes to Figure 2 or Table 3 . These involved rescaling (1) γ,(2) λ, and (3) and the bounds B, D by a factor of 10. The remaining simulations resulted in minor performance differences, especially in the Heterogeneous setting. In the following, we summarize these five analyses.• Adding (low and medium) noise to the network: NNR-Linear slightly outperforms our methods in the Heterogeneous Setting, winning 50-60% of simulations (compared to 40-42% without noise) because NNR-Linear has access to a higher quality network.• Increasing to 10 neighbors: RoME, RoME-BLM perform much better than RoME-SU (the Single-User ablation) in the Heterogeneous Setting, presumably because this change enforced stronger network cohesion among highly similar users.• Rescaling σ by 10: NNR-Linear outperforms our methods in the Heterogeneous setting. NNR-Linear and IntelPooling outperform RoME-SU (but not RoME) in the Nonlinear Setting. Both occur because NNR-Linear and IntelPooling use the true value of σ.• Setting δ to 0.05: NNR-Linear outperforms RoME-SU in the Nonlinear Setting, likely because this change led to insufficient exploration.The results for the medium-noise sensitivity analysis are displayed in Table 5 . Results for the remaining sensitivity analyses are available with our code.Personalizing treatment delivery in mobile health is a common application for online learning algorithms. We focus here on the Valentine study, a prospective, randomized-controlled, remotely administered trial designed to evaluate an mHealth intervention to supplement cardiac rehabilitation for low-and moderate-risk patients (Jeganathan et al., 2022; Golbus et al., 2023) . We aim to use smartwatch data (Apple Watch and Fitbit) obtained from the Valentine study to learn the optimal timing of notification delivery given the users' current context.Prior to the start of the trial, baseline data was collected from each of the participants (e.g., age, gender, baseline activity level, and health information). During the study, participants are randomized to either receive a notification (A t = 1) or not (A t = 0) at each of 4 daily time points (morning, lunchtime, mid-afternoon, evening), with probability 0.25. Contextual information was collected frequently (e.g., number of messages sent in prior week, step count variability in prior week, and pre-decision point step-counts).Since the goal of the Valentine study is to increase participants' activity levels, we define the reward, R t , as the step count for the 60 minutes following a decision point (log-transformed to eliminate skew). Our application also uses a subset of the baseline and contextual data; this subset contains the variables with the strongest association with the reward. Table 6 shows the features available to the bandit in the Valentine study data set.For baseline variables, we use the participant's device model (Z 1 , Fitbit coded as 1), the participant's step count variability in the prior week (Z 2 ), and a measure of the participant's pre-trial activity level based on an intake survey (Z 3 , with larger values corresponding to higher activity levels). 7 : ANOVA analysis of the pseudo-outcomes in the Valentine study. The small p-values constitute strong evidence that the treatment effects differ by participant and over time.At every decision point, before selecting an action, the learner sees two state variables: the participant's previous 30-minute step count (S 1 , log-transformed) and the participant's phase of cardiac rehabilitation (S 2 , dummy coded). The cardiac rehabilitation phase is defined based on a participant's time in the study: month 1 represents Phase I, month 2-4 represents Phase II, and month 5-6 represents Phase III.The motivation for Assumption 4 arises from an exploratory analysis we performed using data from the Valentine study. We constructed the pseudo-reward for each observation as suggested by Equation ( 4) and then conducted an ANOVA test. The results show clear heterogeneity between individuals and across time, motivating the adoption of user-and time-specific random effects. The results of the ANOVA test are displayed in Table 7 . Figure 9 shows the shape of the heterogeneity in the treatment effects over time.Figure 9 : The time heterogeneity in the pseudo-outcomes. We calculated the pseudo-outcomes using 4, then averaged them across participants and plotted them over time. This exploratory analysis shows evidence that the causal effects vary substantially over time.The Valentine study collected the sensor-based features at 4 decision points per day for each study participant. The reward for each message was defined to be log(0.5 + x), where x is the step count of the participant in the 60 minutes following the notification. As noted in the introduction, the baseline reward, i.e. the step count of a subject when no message is sent, not only depends on the state in a complex way but is likely dependent on a large number of time-varying observed variables. Both of these characteristics (complex, time-varying baseline reward function) suggest using our proposed approach.We generated 100 bootstrap samples and ran our contextual bandit on them, considering the binary action of whether or not to send a message at a given decision point based on the contextual variables S 1 and S 2 . Each user is considered independently and with a cohesion network, for maximum personalization and independence of results. To guarantee that messages have a positive probability of being sent, we only sample the observations with notification randomization probability between 0.01 and 0.99. In the case of the algorithm utilizing NNR, we chose four baseline characteristics (gender, age, device, and baseline average daily steps) to establish a measure of "distance" between users. For this analysis, the value of k representing the number of nearest neighbors was set to 5. To utilize bootstrap sampling, we train the Neural-Linear method's neural network using out-of-bag samples. The neural network architecture comprises a single hidden layer with two hidden nodes. The input contains both the baseline characteristics and the contextual variables and the activation function applied here is the softplus function, defined as softplus(x) = log(1 + exp (x)).We performed an offline evaluation of the contextual bandit algorithms using an inverse propensity score (IPS) version of the method from Li et al. (2010) , where the sequence of states, actions, and rewards in the data are used to form a near-unbiased estimate of the average expected reward achieved by each algorithm, averaging over all users. We further investigate the primary contributors to the algorithm's performance, we conducted a parallel ablation study using real-world data analysis. In this study, we compared RoME (as detailed in the main paper) with two additional methods for reference: RoME-SU and NNR-Linear. Each comparison focuses on specific elements of our algorithm to discern their impact on RoME's overall performance.The differences between these methods are illustrated in Figure . 10. In this particular dataset, NNR plays a more significant role in driving the overall superior performance of our algorithm.In the implemented Valentine study, the treatment was randomized with a constant probability p t = 0.25 at each time t. To conduct off-policy evaluation using our proposed algorithm and the competing variations of the TS algorithm, we outline the IPS estimator for an unbiased estimate of the per-trial expected reward based on what has been studied in Li et al. (2010) .Given the logged data D = {s t = s t , A t = a t , R t = r t } T t=1 collected under the policy p = {p t } T t=1 , and the treatment policy being evaluated π = {π t } T t=1 , the objective of this offline estimator is to reweight the observed reward sequence {R t } T t=1 to assign varying importance to actions based on the propensities of both the original and new policies in selecting them. Lemma 1 (Unbiasedness of the IPS estimator). Assuming the positivity assumption in logging, which states that for any given s and a, if p t (a|s) > 0, then we also have π t (a|s) > 0, we can obtain an unbiased per-trial expected reward using the following IPS estimator:RIPS = 1 T T t=1 π t (a t |s t ) p t (a t |s t ) r t (12)As mentioned in the previous section, we restrict our sampling to observations with notification randomization probabilities ranging from 0.01 to 0.99. This selection criterion ensures the satisfaction of the positivity assumption. The proof essentially follows from definition, we have:Proof.E[R IPS ] = E p 1 T T t=1 π t (a t |s t ) p t (a t |s t ) R t (a t , s t ) = 1 T T t=1 π t (a t |s t ) p t (a t |s t ) R t (a t , s t ) × p t (a t |s t ) = 1 T T t=1 π t (a t |s t )R t (a t , s t ) = E π 1 T T t=1 R t (a t , s t )To address the instability issue caused by reweighting in some cases, we use a Self-Normalized Inverse Propensity Score (SNIPS) estimator. This estimator scales the results by the empirical mean of the importance weights, and still maintains the property of unbiasedness. To further enhance the competitive performance of our proposed RoME algorithm, we performed an additional comparative analysis using a real-world data set from the Intern Health Study (IHS) (NeCamp et al., 2020) . This micro-randomized trial investigated the use of mHealth interventions aimed at improving the behavior and mental health of individuals in stressful work environments.The estimates obtained represent the improvement in average reward relative to the original constant randomization, averaging across stages (K = 30) and participants (N = 1553). The available IHS data consist of 20 multiple-imputed data sets. We apply the algorithms to each imputed data set and perform a comparative analysis of the competing algorithms. The results presented in Figure 11 shows our proposed RoME algorithm achieved significantly higher rewards than the other three competing ones and demonstrated performance comparable to the AC algorithm. These findings further support the advantages of our proposed algorithm.Prior to the start of the trial, baseline data was collected on each of the participants (e.g., institution, specialty, gender, baseline activity level, and health information). During the study, participants are randomized to either receive a notification (A t = 1) or not (A t = 0) every day, with probability 3/8. Contextual information was collected frequently (e.g., step count in prior five days, and current day in study).We define the reward, R t , as the step count on the following day (cubic root). Our application also uses a subset of the baseline and contextual data; this subset contains the variables with the strongest association to the reward. Table 8 shows the features available to the bandit in the IHS data set.At every decision point, before selecting an action, the learner sees two state variables: the participant's previous 5-day average daily step count (S 1 , cubic root) and the participant's day in the study (S 2 , an integer from 1 to 30).We run our contextual bandit on the IHS data, considering the binary action of whether or not to send a message at a given decision point based on the contextual variables S 1 and S 2 . Each user is considered independently and with a cohesion network, for maximum personalization and independence of results. To guarantee that messages have a positive probability of being sent, we only sample the observations with notification randomization probability between 0.01 and 0.99. For the algorithm employing NNR, we defined participants in the same institution as their own "neighbors". This definition enables the flexibility for the value of k, representing the number of nearest neighbors, to vary for each participant based on their specific institutional context. Furthermore, in our study setting, we assume that individuals from the same institution enter the study simultaneously as a group. Due to the limited access to prior data, we are unable to build the neural linear models as in the Valentine Study.We utilized 20 multiple-imputed data sets and performed an offline evaluation of the contextual bandit algorithms on each data set. The result is presented below in Figure 11 . Similar to Section B.4, here we also compared RoME (as detailed in the main paper) with RoME-SU and NNR-Linear. The differences between these methods are illustrated in Figure . 12. In this specific dataset, both NNR and DML exhibit comparable performance individually, but their combined effect significantly enhances the overall performance of the algorithm.This appendix briefly discusses two details of Algorithm 1. The first is efficient computation of V -1 . Because V is a fairly large matrix, a full matrix inversion is expensive. Fortunately, however, we can dramatically reduce the necessary computational requirements because each additional (i, t) involves a rank-one perturbation to V . Consequently, we can apply the Sherman-Morrison formula to speed up computations. We leveraged this trick in our implementation of RoME for the simulation study by using efficient rank-one updates available in the SuiteSparse library (Davis & Hu, 2011) .The second detail is the requirements for the distribution, D T S : Definition 1. D T S is a multivariate distribution on R d absolutely continuous with respect to Lebesgue measure which satisfies: 1. (anti-concentration) that there exists a strictly positive probability p such that for any u ∈ R d with ∥u∥ = 1, P(u ⊤ η ≥ 1) ≥ p; and 2. (concentration) there exists c, c ′ positive constants such that ∀δ ∈ (0, 1), P (∥η∥≤ cd log(c ′ d/δ)) ≥ 1 -δ.While a Gaussian prior satisfies Definition 1, this approach allows us to move beyond Bayesian posteriors to generic randomized policies. As discussed in the main paper, we recommend setting Going forward, we use the notation ∆ f i,t (s, ā) := f i,t (s, ā) -f i,t (s, 0) to denote the prediction of the differential reward.Lemma 2. If either p i,t = π i,t or f i,t = r i,t , thenE Rf i,t |s, ā = ∆ i,t (s, ā).That is, the pseudo-reward is an unbiased estimator of the true differential reward.Rf i,t = R it -f i,t (s, A i,t ) δ Ai,t=ā -π i,t (0|s) + ∆ f i,t (s, ā)Case I: π's are correctly specifiedThen E R it δ Ai,t=ā -π i,t (0|s) s, ā = r i,t (s, ā) -r i,t (s, 0) = ∆ i,t (s, ā)andE f i,t (s, A i,t ) δ Ai,t=ā -π i,t (0|s) s, ā = f i,t (s, ā) -f i,t (s, 0) = ∆ f i,t (s, ā) so that E R it -f i,t (s, A i,t ) δ Ai,t=ā -π i,t (0|s) + ∆ f i,t (s, ā) s, ā = ∆ i,t (s, ā) -∆ f i,t (s, ā) + ∆ f i,t (s, ā) = ∆ i,t (s, ā) Case II: f correctly specified E R it δ Ai,t=ā -π i,t (0|s) s, ā = 1 -p i,t (0|s) 1 -π i,t (0|s) r i,t (s, ā) - p i,t (0|s) π i,t (0|s) r i,t (s, 0)andE f i,t (s, A i,t ) δ Ai,t=ā -π i,t (0|s) s, ā = 1 -p i,t (0|s) 1 -π i,t (0|s) f i,t (s, ā) - p i,t (0|s) π i,t (0|s) f i,t (s, 0) = 1 -p i,t (0|s) 1 -π i,t (0|s) r i,t (s, ā) - p i,t (0|s) π i,t (0|s) r i,t (s, 0)andEQUATIONE.2 Preliminaries Lemma 3. Let X be a mean-zero sub-Gaussian random variable with variance factor v 2 and Y be a bounded random variable such that |Y | ≤ B for some 0 ≤ B < ∞. Then XY is sub-Gaussian with variance factor v 2 B 2 .Proof. Recall that X being mean-zero sub-Gaussian means thatP (|X| ≥ t) ≤ 2 exp - t 2 2v 2 . Now note that |XY | ≤ |X|B so that if |XY | > t, then |X|B > t.Thus by monotonicityP (|XY | ≥ t) ≤ P (|X|B ≥ t) = P |X| ≥ t B ≤ 2 exp - t 2 2B 2 v 2 as desired. Lemma 4. If X, Y are sub-Gaussian with variance factors v 2 x , v 2 y , respectively, then αX + βY is sub-Gaussian with variance factor α 2 v 2 x + β 2 v 2 y ∀α, β ∈ R.Proof. Recall the equivalent definition of sub-Gaussianity that X, Y are sub-Gaussian iff for some b, a > 0 and all λ > 0E exp (λ(X -EX)) ≤ exp(λ 2 v 2 x /2) E exp (λ(Y -EY )) ≤ exp(λ 2 v 2 y /2) Then E exp (λ(αX + βY -αEX -βEY )) ≤ E exp (2αλ(X -EX)) E exp (2βλ(Y -EY )) ≤ exp(2α 2 λ 2 v 2 x ) exp(2β 2 λ 2 v 2 y ) = exp((α 2 v 2 x + β 2 v 2 y )λ 2 )The following Lemma gives the sub-Gaussianity and variance of the difference between the pseudoreward and its expectation. We see that in the variance, all terms except those involving the inverse propensity weighted noise variance vanish as f i,t becomes a better estimate of r i,t . Note that means and variances may be implicitly conditioned on the history. Lemma 5. If π i,t is correctly specified and σ2 i,t ≥ 1 c , the difference between the pseudo-reward and its expectation (taken wrt the action and noise) is mean zero sub-Gaussian with varianceEQUATION+ (r i,t (s, 0) -f i,t (s, 0)) 2 [1 -π i,t (0|s)] + Var(ϵ i,t ) π i,t (0|s) -2(r i,t (s, ā) -f i,t (s, ā))(r i,t (s, 0) -f i,t (s, 0))Proof. We need to show that it is sub-Gaussian and upper bound its variance. We write the difference asEQUATIONis bounded and thus (not necessarily mean zero) sub-Gaussian. Since ϵ i,t is sub-Gaussian, its denominator is bounded, and the remaining terms are deterministic, the entire difference between the pseudo-reward and its mean is sub-Gaussian. NowEQUATIONsince E[ Rf i,t |s, ā] is not random. Now we expand the first term on the rhs.EQUATIONFor the first term on the rhs of Eqn. 15,EQUATIONEQUATION2+ E ϵ i,t δ Ai,t=ā -π i,t (0|s) 2 = (r i,t (s, ā) -f i,t (s, ā)) 2 + E[ϵ 2 i,t ] 1 -π i,t (0|s) + (r i,t (s, 0) -f i,t (s, 0)) 2 + E[ϵ 2 i,t ] π i,t (0|s)so that plugging this into Eqn. 15, we haveE Rf i,t (s, ā) 2 = (r i,t (s, ā) -f i,t (s, ā)) 2 + E[ϵ 2 i,t ] 1 -π i,t (0|s) + (r i,t (s, 0) -f i,t (s, 0)) 2 + E[ϵ 2 i,t ] π i,t (0|s) + 2(∆ i,t (s, ā) -∆ f i,t (s, ā))∆ f i,t (s, ā) + ∆ f i,t (s, ā) 2 and plugging this into Eqn. 14 we obtain the variance.EQUATION-(∆ i,t (s, ā) -∆ f i,t (s, ā)) 2as desired. Now note thatEQUATIONEQUATIONEQUATION+ 2 E p(s,ā) [(r i,t (s, ā) -f i,t (s, ā)) 2 ] E p(s) [(r i,t (s, 0) -f i,t (s, 0)) 2 ]Proof. It follows immediately by taking expectations and applying Cauchy Schwartz.Corollary 2. (r i,t (•, 0) -f i,t (•, 0)) 2 = o P (k -1/2 ) and (r i,t (•, ā) -f i,t (•, ā)) 2 = o P (k -1/2 )Proof. Note thatP (r i,t (•, 0) -f i,t (•, 0)) 2 k 1/2 > C ≤ E p(H k-1 ) E p(s) (r i,t (s, 0) -f i,t (s, 0)) 2 k 1/2 C = E p(H k-1 ) o P (k -1/2 )k 1/2 C Assumption 3 = E p(H k-1 ) o P (1) C ≤ δ for sufficiently large k so that (r i,t (•, 0) -f i,t (•, 0)) 2 = o P (k -1/2). In the second to last line we used that since r i,t , f i,t are bounded, E p(s) (r i,t (s, 0) -f i,t (s, 0)) 2 is bounded. For bounded random variables, convergence in probability implies convergence in mean. A similar result holds when using ā instead of 0.EQUATION+ 2o P (k -1/2 )Proof. This follows immediately from the previous two Corollaries.Corollary 4. For all δ > 0, there exists C > 0 s.t. w.p. at least 1 -δ,Var Rf i,t (s, ā) ≤ Ck -1/2 π i,t (0|s) + Var(ϵ i,t ) 1 -π i,t (0|s) + Ck -1/2 [1 -π i,t (0|s)] + Var(ϵ i,t ) π i,t (0|s) + 2Ck -1/2 ≡ v 2 k Proof. By the definition of o P , if (f i,t -r i,t ) 2 = o P (k -1/2 ) for any A > 0, δ > 0, there exists K s.t. if K > k, then w.p. at least 1 -δ, (f i,t -r i,t ) 2 ≤ Ak -1/2 .Then for all k ∈ N, since f, r are bounded(f i,t -r i,t ) 2 ≤ 16B 2 I(k ≤ K) + I(k > K)Ak -1/2 .C ≡ max A, 16B 2 K 1/2 . Remark 1. In the limit as k → ∞ σ2 i,t Var Rf i,t (s, ā) ≤ π i,t (0|s)Var(ϵ i,t ) + (1 -π i,t (0|s))Var(ϵ i,t ) = Var(ϵ i,t )In the next remark, we show what the variance would be if we did not use DML and estimate f i,t ≈ r i,t , but used only the inverse propensity weighted observed reward as the pseudo-reward. This was done in Greenewald et al. (2017) . In this case, there are terms dependent on the mean reward that do not vanish as the number of stages goes to infinity.Remark 2. If we instead used as our pseudo-reward the inverse propensity weighted observed rewardEQUATIONthis would be unbiased with variancev 2 k ≡ r i,t (s, ā) 2 + Var(ϵ i,t ) 1 -π i,t (0|s) + r i,t (s, 0) 2 + Var(ϵ i,t ) π i,t (0|s) -∆ i,t (s, ā) 2Proof. The unbiasedness is clear from our proof of Lemma 2. For the variance,Var R it δ Ai,t=ā -π i,t (0|s) = E (r i,t (s, ā) + ϵ i,t ) 2 1 -π i,t (0|s) + (r i,t (s, 0) + ϵ i,t ) 2 π i,t (0|s) -∆ i,t (s, ā) 2 = r i,t (s, ā) 2 + Var(ϵ i,t ) 1 -π i,t (0|s) + r i,t (s, 0) 2 + Var(ϵ i,t ) π i,t (0|s) -∆ i,t (s, ā) 2Here we collect some important results. We first adapt an important concentration inequality for regularized least-squares estimates. Our proof follows the same basic strategy as Abbasi-Yadkori et al. ( 2011), but with modifications due to (a) the use of weighted least squares (b) the use of pseudo-rewards to estimate differential rewards (c) replacing the scaled diagonal regularization with Laplacian regularization.Lemma 6 (Adapted from Theorem 2 in Abbasi-Yadkori et al. ( 2011)). Let θk be the stage k regularized least squares (RLS) estimate from Algorithm 1 and θ * the ground truth where we assume ∥θ * ∥ ≤ S. For any δ > 0, w.p. at least 1 -δ the estimatesEQUATIONwhere v 2 k is the variance factor for the difference between the pseudo-reward and its mean at stage k. In particular, settingx k = V k ( θk -θ * ) implies ∥ θk -θ * ∥ V k ≤ v k 2 log det(V k ) 1/2 det(V 0 ) -1/2 δ + ∥θ * ∥ V0holds w.p. at least 1 -δ for all k ≥ 1.Proof. Let m i,t = σi,t ϕ(x i,t ) andρ i,t = σi,t [ Rf i,t (s, ā) -E[ Rf i,t |s, ā]]. Further let ξ k ≡ (i,t)∈O k-1 σ2 i,t [ Rf i,t (s, ā) -E[ Rf i,t |s, ā]]ϕ(x i,t ) = (i,t)∈O k-1 σi,t [ Rf i,t (s, ā) -E[ Rf i,t |s, ā]]m i,t = (i,t)∈O k-1 m i,t ρ i,tThen noting thatV k = (i,t)∈O k-1 m i,t m ⊤ i,t + V 0 and b k = (i,t)∈O k-1 σ 2 i,tRf i,t ϕ i,t and letting W k be the diagonal matrix of weights σ2 i,t , we haveθk = V -1 k b k = V -1 k (ξ k + Φ ⊤ k W k E[R f k k |ā, s]) = V -1 k ξ k + V -1 k Φ ⊤ k W k ∆ k by Lemma 2 = V -1 k ξ k + V -1 k Φ ⊤ k W k Φ k θ * = V -1 k ξ k + V -1 k (Φ ⊤ k W k Φ k + V 0 )θ * -V -1 k V 0 θ * = V -1 k ξ k + θ * -V -1 k V 0 θ * k and thus θk -θ * = V -1 k (ξ k -V 0 θ * ) which gives |x ⊤ k θk -x ⊤ k θ * | ≤ ∥x k ∥ V -1 k (∥ξ k ∥ V -1 k + ∥V 0 θ * ∥ V -1 k )Now since ξ k is sub-Gaussian with variance factor v 2 k , by Theorem 1 in Abbasi-Yadkori et al. ( 2011), w.p. 1 -δ,∥ξ k ∥ 2 V -1 k ≤ 2v 2 k log det(V k ) 1/2 det(V 0 ) -1/2 δ Further note that since V 0 ⪯ V k , then V -1 k ⪯ V -1 0 . Thus ∥V 0 θ * ∥ 2 V -1 k = θ * ⊤ V ⊤ 0 V -1 k V 0 θ * ≤ θ * ⊤ V ⊤ 0 V -1 0 V 0 θ * = θ * ⊤ V 0 θ * . Finally, setting x k = V k ( θk -θ * ) implies ∥ θk -θ * ∥ 2 V k = ( θk -θ * ) ⊤ V k ( θk -θ * ) = |x ⊤ k θk -x ⊤ k θ * | ≤ ∥x k ∥ V -1 k   v k 2 log det(V -1 k ) 1/2 det(V 0 ) -1/2 δ + θ * ⊤ V 0 θ *   = ∥V k ( θk -θ * )∥ V -1 k   v k 2 log det(V -1 k ) 1/2 det(V 0 ) -1/2 δ + θ * ⊤ V 0 θ *   = ∥ θk -θ * ∥ V k   v k 2 log det(V -1 k ) 1/2 det(V 0 ) -1/2 δ + θ * ⊤ V 0 θ *   and dividing both sides by ∥ θk -θ * ∥ V k implies ∥ θk -θ * ∥ V k ≤ v k 2 log det(V k ) 1/2 det(V 0 ) -1/2 δ + θ * ⊤ V 0 θ *In the above inequality, we will need to bound ∥θ * ∥ V0 in order to obtain useful confidence sets. Lemma 7.∥θ * ∥ V0 ≤ √ γB shared + √ λW K( D user + D time ) + γW K(B user + B time ).Proof. Recall that V 0 = diag(γI p , λL user ⊗ + γI Kp , λL time ⊗ + γI Kp ), (18) Thenθ * ⊤ V 0 θ * = γ∥θ shared ∥ 2 + λtr Θ ⊤ user L user Θ user + γtr Θ ⊤ user I Kp Θ user + λtr Θ ⊤ time L time Θ time + γtr Θ ⊤ time I Kp Θ time = γ∥θ shared ∥ 2 + λ (i,j)∈Euser ∥θ user i -θ user j ∥ 2 2 + λ (i,j)∈Etime ∥θ time i -θ time j ∥ 2 2 + γ∥Θ user ∥ 2 F + γ∥Θ time ∥ 2 F ≤ γB 2 shared + λW K(D user + D time ) + γW K(B 2 user + B 2 time )by Assumption 5 where we assume that each user and time have at most m neighbors. Now using√ a + b ≤ √ a + √ b for a, b > 0 we have θ * ⊤ V 0 θ * ≤ √ γB shared + √ λW K( D user + D time ) + γW K(B user + B time ).We will also want to bound the log determinant terms. The next two Lemmas show us how to do so. Lemma 8.det(V K+1 ) ≤ 3K(K+1) 8 + γ(2K + 1)p + 2λM Kp (2K + 1)p (2K+1)pProof. By AM-GM inequality,det(V K+1 ) = det   (i,t)∈O K σ2 i,t ϕ(x i,t )ϕ(x i,t ) ⊤ + V 0   ≤   tr (i,t)∈O K σ2 i,t ϕ(x i,t )ϕ(x i,t ) ⊤ + V 0 (2K + 1)p   (2K+1)p Now tr   (i,t)∈O K σ2 i,t ϕ(x i,t )ϕ(x i,t ) ⊤   ≤ (i,t)∈O K σ2 i,t ∥ϕ(x i,t )∥ 2 2 ≤ 3 4 K(K + 1) 2 .Recall that V 0 = diag(γI p , λL user ⊗ + γI Kp , λL time ⊗ + γI Kp ), so that tr(V 0 + γ(2K + 1)p + 2λM Kp (2K + 1)p(2K+1)p Lemma 9. log det(V K ) det(V 0 ) ≤ (2K + 1)p log 3K(K + 1) γ8(2K + 1)p + 1 + 2λM K γ(2K + 1) Proof. Note that det(V 0 ) ≥ γ (2K+1)p , so that det(V K ) det(V 0 ) ≤ 3K(K+1) 8 + γ(2K + 1)p + 2λM Kp γ(2K + 1)p (2K+1)pTaking the log, we havelog det(V K ) det(V 0 ) ≤ log 3K(K+1) 8 + γ(2K + 1)p + 2λM Kp γ(2K + 1)p (2K+1)p ≤ (2K + 1)p log 3K(K + 1) γ8(2K + 1)p + 1 + 2λM K γ(2K + 1)Corollary 5. If Assumption 3, then for any δ > 0, there exists C > 0 s.t. w.p. at least 1 -δ the estimates { θk } ∞ k=0 in Algorithm 2 satisfies for anyEQUATIONIn particular, settingx k = V k-1 ( θk-1 -θ * ) implies ∥ θk -θ * ∥ V k ≤ C k 1/2 + σ 2 c 2 2 log det(V k ) 1/2 det(V 0 ) -1/2 δ/2 + ∥θ * ∥ V0holds w.p. at least 1 -δ for all k ≥ 1.Proof. Use Corollary 3 and Lemma 6, each with δ/2. Then w.p. at least 1 -δ the result holds.We next state a slightly modified form of a standard result of RLS (Lemma 11 in Abbasi-Yadkori et al. ( 2011)) that helps to guarantee that the prediction error is cumulatively small. This bounds the sum of quadratic forms where the matrix is the inverse Gram matrix and the arguments are the feature vectors. We use such terms to construct a martingale in the regret bound so that we can bound such terms and the martingale.Proposition 1. Let λ ≥ 1 and γ ≥ 1. For any arbitrary sequence (x i,t ) (i,t)∈O k , letV k+1 ≡ (i,t)∈O k σ2 i,t ϕ(x i,t )ϕ(x i,t ) ⊤ + V 0 ,be the regularized Gram matrix. ThenK k=1 (i,t)∈O k \O k-1 ∥ϕ(x i,t )∥ 2 V -1 k ≤ 2c log det(V K+1 ) det(V 0 ) .where c is a constant such that 0 < 1 c < σ2 i,t ∀i, t ∈ N.Proof. By Lemma 11 in Abbasi-Yadkori et al. ( 2011), we haveK k=1 (i,t)∈O k \O k-1 σ2 i,t ∥ϕ(x )∥ 2 V -1 k ≤ 2 log det(V K+1 ) det(V 0 ) .The lower bound on the weights impliesK k=1 (i,t)∈O k \O k-1 ∥ϕ(x i,t )∥ 2 V -1 k ≤ 2c log det(V K+1 ) det(V 0 )as desired.Finally, we state Azuma's concentration inequality which describes concentration of super-martingales with bounded differences and is useful in controlling the regret due to the randomization of Thompson sampling.Proposition 2 (Azuma's concentration inequality). If a super-martingale (Y t ) t≥0 corresponding to a filtration F t satisfies |Y t -Y t-1 | < c t some constant c t for all t = 1, . . . , T then for any α > 0:P (Y T -Y 0 ≥ α) ≤ exp - α 2 2 T t=1 c 2 t .We first decompose the regret boundK k=1 1 k (i,t)∈O k \O k-1 π ⋆ i,t x(S i,t , a * i,t ) ⊤ θ * i,t -π i,t x(S i,t , A i,t ) ⊤ θ * i,t = K k=1 1 k (i,t)∈O k \O k-1 (π ⋆ i,t -π i,t )x(S i,t , A i,t ) ⊤ θ * i,t + π ⋆ i,t x(S i,t , a * i,t ) ⊤ θ * i,t -x(S i,t , A i,t ) ⊤ θ * i,t ≤ K k=1 1 k (i,t)∈O k \O k-1 (π ⋆ i,t -π i,t )x(S i,t , A i,t ) ⊤ θ * i,t + K k=1 1 k (i,t)∈O k \O k-1 x(S i,t , a * i,t ) ⊤ θ * i,t -x(S i,t , A i,t ) ⊤ θ * i,tFor the first term,K k=1 1 k (i,t)∈O k \O k-1 (π ⋆ i,t -π i,t )x(S i,t , A i,t ) ⊤ θ * i,t, one can apply a nearly identical proof to that of Greenewald et al. (2017) . For the second term we need a novel strategy, and thus focus on that term. The proof follows closely from Abeille & Lazaric (2017) with several adjustments. Assumption 6 implies that we only need to consider the unit ball X = {∥x∥ ≤ 1}. Then the second term of the regret can be decomposed intoK k=1 1 k (i,t)∈O k \O k-1 (ϕ(x ⋆ i,t )) ⊤ θ ⋆ -ϕ(x i,t ) ⊤ θk R T S (K) + K k=1 1 k (i,t)∈O k \O k-1 ϕ(x i,t ) ⊤ θk -ϕ(x i,t ) ⊤ θ ⋆ R RLS (K)where ϕ(x ⋆ i,t ) is the context vector under the optimal action and θ ⋆ is the true parameter value. The first term is the regret due to the random deviations caused by sampling θk and whether it provides sufficient useful information about the true parameter θ ⋆ k . The second term is the concentration of the sampled term around the true linear model for the advantage function. Definition 2. We define the filtration F k as the information accumulated up to stage k before the sampling procedure, that is, F k = (F 1 , σ(x 1 , r 2 , x 2 , . . . , x k-1 , r k-1 )), and filtration F x k as the information accumulated up to stage k and including the sampled context, that is,F t = (F 1 , σ(x 1 , r 2 , x 2 , . . . , x k-1 , r k-1 , x k )).Bounding R RLS (T ). We decompose the second term into the variation of the point estimate and variation of the random sample around the point estimate:K k=1 1 k (i,t)∈O k \O k-1 ϕ(x i,t ) ⊤ θk -ϕ(x i,t ) ⊤ θk + K k=1 1 k (i,t)∈O k \O k-1 ϕ(x i,t ) ⊤ θk -ϕ(x i,t ) ⊤ θ ⋆ kThe first term describes the deviation of the TS linear predictor from the RLS one, while the second term describes the deviation of the RLS linear predictor from the true linear predictor. The first term is controlled by the construction of the sampling distribution D T S , while the second term is controlled by the RLS estimate being a minimizer of the regularized cumulative squared error in (8).In particular, the first term will be small when the TS estimate concentrates around the RLS one, while the second will be small when the RLS estimate concentrates around the true parameter vector. The next proposition gives a lower bound on the probability that, for all stages, both the RLS parameter vector concentrates around the true parameter vector and the TS parameter vector concentrates around the RLS one.Recall thatβ k (δ) = v k 2 log det(V k ) 1/2 det(V 0 ) 1/2 δ/2 1/2 + B where B = √ γB shared + √ λM K( √ D user + √ D time ) + √ γM K(B user + B time ). This is, from Lemma 6 and Corollary 7, an upper bound on a 1 -δ confidence set on the RLS estimator.Proposition 3. Let Êk denote the event that θk concentrates around the true parameter for all l ≤ k,i.e., Êk = {∀l ≤ k, ∥ θl -θ ⋆ l ∥ V l ≤ β l (δ ′ )}. Let γ k (δ) ≡ β k (δ ′ ) cd log c ′ d δ Let Ẽk denote the event that θl concentrates around the estimated parameter for all l ≤ k, i.e., Ẽk = {∀l ≤ k, ∥ θl -θl ∥ V l ≤ γ l (δ ′ )}. Let E k = Êk ∩ Ẽk . Then P (E k ) ≥ 1 -δ/2.Proof. Let δ ′ = δ/4K, then Lemma 6 and a union bound give usP ( ÊK ) = P (∩ K k=1 {∥ θk -θ ⋆ ∥ V k ≤ β k (δ ′ )}) = 1 - K k=1 P (∥ θk -θ ⋆ ∥ V k > β k (δ ′ )) = 1 - K k=1 δ ′ = 1 -δ ′ K = 1 -δ/4.Applying the TS sampling distribution and θk= θk + β k (δ ′ )V -1/2 k η k where η t is drawn i.i.d. from D T S we have P ∥ θk -θk ∥ V k ≤ β k (δ ′ ) cd log c ′ d δ ′ = P ∥η k ∥ ≤ cd log c ′ d δ ′ ≥ 1 -δ ′ .by Definition 1. A union-bound argument yields the conclusion.We can then bound R RLS (K) by leveraging Lemma 6 and decomposing the error viaR RLS (K) ≤ K k=1 1[E K ] k   (i,t)∈O k \O k-1 |ϕ(x i,t ) ⊤ ( θk -θk )|   + K k=1 1[E K ] k   (i,t)∈O k \O k-1 |ϕ(x i,t ) ⊤ ( θk -θ * )|  By definition of the event E K , we have|ϕ(x i,t ) ⊤ ( θk -θk )|1[E k ] ≤ ∥ϕ(x i,t )∥ V -1 k γ k (δ ′ ), |ϕ(x i,t ) ⊤ ( θk -θ ⋆ )|1[E k ] ≤ ∥ϕ(x i,t )∥ V -1 k β k (δ ′ ) so from Proposition 1, we have K k=1 1[E K ] k   (i,t)∈O k \O k-1 |ϕ(x i,t ⊤ ( θk -θk )|   ≤ γ K (δ ′ ) K k=1   (i,t)∈O k \O k-1 1 k ∥ϕ(x i,t )∥ V -1 k   ≤ γ K (δ ′ ) K k=1 (i,t)∈O k \O k-1 1 k 2 K k=1 (i,t)∈O k \O k-1 ∥ϕ(x i,t )∥ 2 V -1 k ≤ γ K (δ ′ ) K k=1 1 k K k=1 (i,t)∈O k \O k-1 ∥ϕ(x i,t )∥ 2 V -1 k ≤ γ K (δ ′ ) H K (i,t)∈O k ∥ϕ(x i,t )∥ 2 V -1 k ≤ γ K (δ ′ ) H K 2c log det(V K+1 ) det(V 0 ) .Using a similar derivation for the β k (δ ′ ) case, we obtainR RLS (K) ≤ (β K (δ ′ ) + γ K (δ ′ )) K k=1 (i,t)∈O k \O k-1 1 k 2 2c log det (V K+1 ) det (V 0 ) ≤ (β K (δ ′ ) + γ K (δ ′ )) K k=1 1 k 2c log det (V K+1 ) det (V 0 ) ≤ (β K (δ ′ ) + γ K (δ ′ )) H K 2c (2K + 1)p log 3K(K + 1) γ8(2K + 1)p + 1 + 2λM K γ(2K + 1)with probability at least 1 -δ/2 by Proposition 3, where H K is the harmonic number. Note that H K ∼ log(K) for large K.Bounding R T S (T ). Leveraging Abeille & Lazaric (2017) , Definition 1 lets us bound R T S (K)EQUATIONWe re-write the sum in (20) as:K k=1 1 k (i,t)∈O k \O k-1 ∥ϕ(x i,t )∥ V -1 k + K k=1 1 k (i,t)∈O k \O k-1 E ∥ϕ(x ⋆ i,t )( θ)∥ V -1 k |F k -∥ϕ(x i,t )∥ V -1 k R T S 2The first term is bounded by Proposition 1:K k=1 1 k (i,t)∈O k \O k-1 ∥ϕ(x i,t )∥ V -1 k ≤ 2cH K log det(V K+1 ) det(V 0 )The second term is a martingale by construction and so we can apply Azuma's inequality. Under Assumption 6, so since V k ≤ 1 λ I we haveE ∥ϕ(x i,t ) ⋆ ( θ)∥ V -1 k |F t -∥ϕ(x i,t )∥ V -1 k ≤ 2 √ λ , a.s.This provides the upper-boundR T S (K) ≤ 4γ K (δ ′ ) d 8K λ log 4 δ + 2cH K (2K + 1)p log 3K(K + 1) γ8(2K + 1)p + 1 + 2λM K γ(2K + 1).Overall bound. Putting together the two bounds under a union bound argument yields the upper bound in Theorem 1; specifically, we haveβ K (δ ′ ) + γ K (δ ′ ) 1 + 4 d 2cH K (2K + 1)p log 3K(K + 1) γ8(2K + 1)p + 1 + 2λM K γ(2K + 1) + 4γ K (δ ′ ) p 8K λ log 4 δFor convenience, we summarize below some key notation used in the main paper.• a = 0: control action• q: number of non-baseline treatment arms• i = 1, 2, . . .: index for individuals (later in the paper, we consider N individuals)• t = 1, 2, . . .: index for decision points (later in the paper, we consider T decision points)• S i,t ∈ S: context vector observed for individual i at decision point t 2k+1) : the set of all parameters, including the individual-specific time-invariant parameters and the shared time-specific parameters, at stage k of the sequential recruitment setting • D user , D time , B user , B time , M : terms for bound on parameters and mean estimate] ∈ R p(In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: In Section 6.2, the experiment results include hypothesis testing, and we have clearly indicated the p-values in the Figure legend to denote statistical significance.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: This paper is driven by the current challenges in mobile health policy learning, aiming to enhance policy learning to maximize the benefits of mobile interventions. As a result, we emphasize the societal impacts, focusing on how the algorithm has the potential to improve public health.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: The paper proposes a new contextual bandit algorithm, namely DML-TS-NNR, for mHealth which can effectively incorporate the nonstationarity, nonlinearity and user heterogeneity of the health outcomes.
DML-TS-NNR allows the baseline outcome (outcome under "no action") to be highly complex and imposes a parametric structure only to the differential reward. The differential reward model parameter includes user specific and time specific parameters, with minimum assumptions that parameters of proximal users (according to a known network) and proximal time are similar. Laplacian regularization is employed to efficiently pool information across users and time. Finally, debiased machine learning methods are employed in the construction of a doubly robust estimate of the differential reward, which reduces the variance of existing IPW estimates. The paper also provides a regret upper bound which solely depends on the dimension of the differential reward model. 

The paper is well-motivated from mHealth applications and provides a novel algorithm that is robust to the frequently occurring challenges such as nonstationarity, heterogeneity, and complexity of baseline rewards. These points are well supported by theory and extensive experiments.

Some concerns raised by reviewers were as follows,
- robustness to hyperparameters
- scalability
- the fact that the network structure should be known.

As for the first point, authors have provided in their rebuttal an extensive sensitivity analysis to show the robustness of the algorithm to the hyperparameter values, and also suggested to use the Empirical Bayes method to set the hyperparameters in practice.

As for the second point, authors provided the exact computation time in their rebuttal, showing that the speed is competitive to other existing methods and is fast enough in mHealth scenarios.

As for the third point, as authors claimed, it is true that there exist cases where the network structure is known, and most existing works exploiting the network structure assume the knowledge of the network. Authors also provided some experiments for the case where the network is not known perfectly and showed robustness of their method.

Reviews: review: paper_summary: paper_summary: The paper introduces a novel robust mixed-effects bandit algorithm, named "DML-TS-NNR", designed to optimize mobile health (mHealth) interventions. mHealth aims to deliver personalized and contextually tailored notifications to promote healthier behaviors. The proposed algorithm addresses key challenges in mHealth, such as participant heterogeneity, nonstationarity, and nonlinearity in rewards. The main contributions of the paper are:
- Modeling Differential Rewards: Incorporates user- and time-specific parameters.
- Network Cohesion Penalties: Uses penalties to pool information across users and time.
- Debiased Machine Learning: Employs this technique for flexible baseline reward estimation.

The algorithm's high-probability regret bound is solely dependent on the differential reward model's dimension. The effectiveness of DML-TS-NNR is demonstrated through simulations and two off-policy evaluation studies.
strength_weakness: strength_weakness: - The algorithm's reliance on complex calculations, such as log-determinants and matrix inverses, may limit its scalability for large datasets.
- The assumption of a known network with binary edges may not always be practical. Real-world networks can be more complex, and this assumption might limit the algorithm's applicability.
questions: questions: - How robust is the algorithm to violations of the assumptions regarding the network structure and hyperparameters? Can the authors provide guidance on tuning these parameters in practice?
- How can the algorithm be extended to consider long-term effects and treatment fatigue? Are there plans to address these aspects in future work?
limitations: limitations: While the authors have made strides in addressing some key limitations, such as computational demand and network assumptions, there are still areas that may need further attention:
- The high computational demand due to log-determinants and matrix inverses is a limitation. Future work should explore more efficient computational methods or approximations.
- The algorithm focuses on immediate rewards and does not consider long-term effects or treatment fatigue. Addressing these aspects is crucial for real-world applicability and effectiveness.
- The need for correctly specified hyperparameters can be challenging in practice. Providing more robust methods for hyperparameter tuning and addressing potential misspecifications would improve the algorithm's robustness and ease of use.
paper_summary: paper_summary: The authors propose a novel contextual bandit algorithm that addresses individual heterogeneity, nonstationarity, and nonlinearity of the reward function. This algorithm involves three distinct steps to manage these challenges.
strength_weakness: strength_weakness: The selection of tuning parameters plays a crucial role but is not discussed in the proposed method. It is unclear whether the promising performance of the proposed algorithm is due to the choice of a specific tuning parameter. Providing some robustness checks would enhance the soundness of the simulation experiments.

I also need some motivation for Assumption 4. For example, how well does the dataset under investigation support this assumption? In its current form, Assumption 4 appears to be arbitrary.

I also found the assumption of a known network structure to be quite restrictive, as network structures are typically unknown in practice. This limitation hinders the practicality of the proposed algorithm.
questions: questions: I found the comments on page 4 about the doubly robustness concerning. As the methods discussed in Section 4.1 are all nonparametric estimators (hence, the model misspecification is not an issue), the benefit of incorporating DML is to improve estimation efficiency.
limitations: limitations: not noted.
paper_summary: paper_summary: The paper "A Robust Mixed-Effects Bandit Algorithm for Assessing Mobile Health Interventions" introduces the DML-TS-NNR (Debiased Machine Learning Thompson Sampling with Nearest-Neighbor Regularization) algorithm. This novel contextual bandit algorithm is designed to address challenges in mobile health (mHealth) interventions, such as participant heterogeneity, nonstationarity, and nonlinearity in rewards. The algorithm incorporates user- and time-specific incidental parameters, network cohesion penalties, and debiased machine learning to flexibly estimate baseline rewards.
strength_weakness: strength_weakness: ### Novelty:

While the combination of debiased machine learning, network cohesion penalties, and Thompson sampling is innovative, a more detailed comparison with existing methods, particularly in terms of theoretical and practical advantages, would further highlight the unique contributions of the proposed DML-TS-NNR algorithm.


### Experimental Validation:

The experimental validation primarily relies on controlled datasets. Including more diverse real-world testing scenarios, such as various health conditions or treatment modalities, would provide a more comprehensive assessment of the DML-TS-NNR algorithm's effectiveness and generalizability.


### Technical Details:

Some aspects of the algorithm, such as the optimization process for debiased machine learning and the derivation of the network cohesion penalties, could be explained in greater detail to enhance clarity and understanding. The choice of evaluation metrics and their suitability for different types of mHealth interventions could be discussed more extensively.
questions: questions: 1. Can the authors elaborate on why specific baseline methods (e.g., Standard Thompson Sampling, Action-Centered contextual bandit) were chosen for comparison? What advantages do these techniques offer for evaluating the effectiveness of the DML-TS-NNR algorithm in mHealth interventions?

2. How does the DML-TS-NNR algorithm perform in more diverse and dynamic real-world settings, such as different health conditions or treatment modalities? Are there plans to test the approach in more varied environments, including chronic diseases or mental health applications?

3. What potential challenges could arise regarding the scalability of the DML-TS-NNR algorithm for very large datasets or real-time applications? How does the algorithm handle computational and communication overhead in such scenarios?

4. Can the authors provide more details on the optimization process for debiased machine learning and its impact on overall performance and accuracy of reward estimation?

5. What are the potential ethical considerations or privacy concerns associated with the proposed method, particularly in using sensitive health data for mHealth interventions? How does the DML-TS-NNR algorithm address these concerns?
limitations: limitations: 1. Potential challenges and mitigation strategies for deploying the DML-TS-NNR algorithm in real-world settings, particularly regarding data diversity and environmental variability.
2. Discussing the ethical implications and privacy concerns related to continuous monitoring and intervention in high-stakes healthcare applications. This includes addressing the implications of using sensitive patient data and the potential impact of intervention decisions on different demographic groups.
3. Further exploring the trade-offs between computational efficiency and intervention accuracy, particularly in scenarios where rapid decision-making is crucial for clinical outcomes.
paper_summary: paper_summary: The paper introduces a robust contextual bandit algorithm to optimize personalized mobile health interventions. The proposed algorithm leverages Thompson sampling, mixed-effect models, debiased machine learning, and nearest-neighbor regularization techniques to address the problems of user and time heterogeneity, information pooling, and complex baseline rewards. The paper establishes a high-probability regret bound and demonstrates the algorithm’s effectiveness compared to existing methods through a simulation and two off-policy evaluation studies.
strength_weakness: strength_weakness: The proposed method combines several existing techniques, including Thompson sampling, mixed-effect models, debiased machine learning, and nearest-neighbor regularization. Could the authors elaborate on the specific challenges associated with the proposed method, whether in terms of methodology, computation, or theoretical proof?
questions: questions: -	It is mentioned in the appendix that the nearest neighbor network is assumed to be known in the simulation study, and that the number of nearest neighbors k is set to 5 in the Valentine Study. Could the authors provide details on how the nearest neighbor network is defined in the simulation study and how the number of nearest neighbors is chosen in the Valentine Study? Is the algorithm's performance sensitive to the construction of the nearest neighbor network?
-	The paper lists several feasible choices for the working model f_{i,t}. However, its construction in the simulation and real data analysis was not detailed. It would be helpful if the authors could discuss how the working model is chosen in practice.
-	Assumption 1 ensures that the randomization probability for action 0 is always positive, but it does not specify this for other actions. Is a positivity assumption required for all actions in the proposed method? If so, the authors should discuss how this is ensured for other actions.
-	Computational efficiency is crucial for online algorithms. It would be beneficial if the authors could provide the computation time of the proposed method, e.g., in the simulation study or the Valentine Study. How does it compare with the baseline methods? 
-	In the single-column vector \theta, user-specific parameters are numbered from 1 to K. Should they instead be numbered from 1 to N, as in the \Theta_{user} matrix? Besides, it seems that N has not been defined in the main text.
limitations: limitations: The authors have discussed the limitations of the proposed method.

Comments: comment: Thank you once again for your insightful comments—we really appreciate your support! We will be sure to incorporate these discussions into our revised manuscript.
Thank you for your thoughtful and detailed response. I particularly appreciate the in-depth discussion of the technical challenges, which I found highly informative. Accordingly, I have adjusted my rating upward.
Thank you very much for your feedback; we will be sure to incorporate your comments into the revised manuscript.
Thank you for your responses. I currently don't have any additional questions. I will make my final decision after further discussions with the other reviewers and the AC.