Off-Policy Estimation of Long-Term Average Outcomes with Applications to Mobile Health

Introduction

Due to the recent advancement in mobile device and sensing technology, health scientists are more and more interested in developing mobile health (mHealth) interventions. In mHealth, mobile devices (e.g., wearables and smartphones) are used to deliver interventions to individuals as they go about their daily lives. In general, there are two types of mHealth treatments. Most are “pull” treatments that reside on the individual’s mobile device and allow the individual to access treatment content as needed. This work focuses on the second type, the “push” treatment, typically in the form of a notification or a text message that appears on a mobile device. There is a wide variety of possible treatment messages (e.g., behavioral, cognitive, and motivational message and reminders). These treatments are generally intended to impact a near time, proximal outcome, such as stress or behaviors such as physical activity over some subsequent minutes/hours. The mHealth intervention policies, often called just-in-time adaptive interventions in the mHealth literature (Nahum-Shani et al. 2018), are decision rules that map the individual’s current state (e.g., past behaviors as well as current observations of location, social activity, stress and urges to smoke) to a particular treatment at each of many time points. Many mHealth interventions are designed for long-term use in chronic disease management (Lee et al. 2018). The vast majority of current mHealth interventions deploy expert-derived policies with limited use of data evidence (for an example see Kizakevich et al. 2014), however the long-term efficacy of these policies on the health behavior is not well understood. An important first step toward developing data-based, effective mHealth interventions is to properly measure the long-term performance of these policies. In this work, we provide an approach for conducting inference about the optimality of one or more mHealth policies of interest. Our optimality criterion is the long-term average of the proximal outcomes should a particular mHealth policy be followed. We develop a flexible method to estimate the performance of an mHealth policy using a historical dataset in which the treatments are decided by a possibly different policy.

This work is motivated by HeartSteps (Klasnja et al. 2015), an mHealth physical activity intervention. To design this intervention, we are conducting a series of studies. The first, already completed, study was for 42 days. The last study will be for one year. Here we focus on the intervention component involving activity suggestions. These suggestions may be delivered at each of the five individual-specified times per day. While in the first study there were $42\times 5=210$ time points per individual, in the year-long study there will be about 2,000 time points per individual. The proximal outcome is the step count in the 30 minutes following each of the five times per day. Our goal is to use the data collected from the first 42-day study to predict and estimate the long-term average of proximal outcomes for a variety of policies that could be used to decide whether or not to send the activity suggestion at each time point in the year-long study. The 42-day study was a Micro-Randomized Trial (MRT) (Klasnja et al. 2015; Liao et al. 2016). In an MRT, a known stochastic policy, also called a behavior policy, is used to decide when and which type of treatment to provide at each time point. A partial list of MRTs in the field or completed can be found at the website http://people.seas.harvard.edu/~samurphy/JITAI_MRT/mrts4.html. From an experimental point of view, the stochastic behavior policy is used to conduct sequential randomizations within each individual. Here the adjective, “stochastic”, means that at each time point each individual is randomized between the possible treatments. In this work we focus on settings in which the randomization probabilities are known functions of the individual’s past data; this is the case with MRTs by design.

The rest of the article is organized as follows. Section 2 provides a review of Markov Decision Processes. In Section 3, we review related work. Section 4 develops an estimator for the long-term average proximal outcome; then Section 5 provides the asymptotic distribution of this estimator. As the estimation requires tuning parameters, in Section 6 we provide a procedure to select the tuning parameters. Simulations are used to assess the coverage probability of the proposed confidence intervals in various settings. A case study using data from the 42-day MRT of HeartSteps is presented in Section 7. We end with a discussion of future work in Section 8.

Distributional Assumptions and Goal

The data for each individual is of the form

Note that the MDP does not specify the distribution of the actions. And indeed the distribution of the actions may not satisfy the Markovian property. In an MRT, the actions, $\{A_{t}\}_{t=1}^{T}$, are randomized with probabilities that can depend on the entire history prior to time point $t$, $H_{t}=\{S_{1},A_{1},\dots,S_{t}\}$. Denote the distribution of $A_{t}\,|\,H_{t}$ by $\pi_{t}^{b}(\cdot\,|\,H_{t})$. We call $\pi_{b}=\{\pi_{1}^{b},\ldots,\pi_{T}^{b}\}$, a stochastic behavior policy. Throughout we assume that $\pi_{b}$ is known (as is the case in an MRT) and that the probabilities are strictly positive, i.e., $\pi_{t}^{b}(a\,|\,H_{t})\geq p_{\min}>0$ for all $a\in\mathcal{A}$, $H_{t}$ and $t\leq T$.

Suppose that a pre-specified time-invariant, Markovian policy, $\pi$, is being considered for use in future. Our goal is to conduct inference for the resulting average of the rewards over a large number of time points. In mHealth, the policy might be an expert-constructed policy. Considering a long time period makes most sense for individuals who are struggling with chronic problems or disorders for which, at this time, there is no general cure. Many health-behavior problems fall into this area including obesity, hypertension, adherence to medications for AIDs, mental illness and addictions. Let $\pi(a\,|\,s)$ be the probability of choosing the action, $a$, at the state, $s$. Given a dataset that consists of $n$ independent, identically distributed (i.i.d.) observations of $\mathcal{D}$, we aim to estimate the average reward of the policy, defined as

We propose to conduct inference about the long-term performance of each policy, $\pi$, via its average reward, $\eta^{\pi}$. This is because the average reward, $\eta^{\pi}$, is an asymptotic surrogate of the average of finite rewards over a long period of time. In fact, it can be shown that

where the leading constant depends on the mixing time of $P^{\pi}$ (see Theorem 7.5.10 in Hernández-Lerma & Lasserre 1999). In the case of HeartSteps, the goal is to use the data from the 42-day MRT study to estimate the average reward, $\eta^{\pi}$, for a variety of policies $\pi$. The average reward, $\eta^{\pi}$, provides a proxy for the average of the 30-min step counts when the policy, $\pi$, is used to determine whether to send the activity suggestions over a long time period (e.g., a year: $5\times 365$ time points).

Note that the data, $\cal D$, on each individual includes observations over $T$ time points and the actions are selected according to a behavior policy. However, as will be seen, the above assumptions including the Markovian and irreducibility assumptions will allow us to estimate the average reward over a long time period and under a different policy.

Related Work

The evaluation of a given target policy using data collected from a different policy (i.e., the behavior policy) is called off-policy evaluation. This has been widely studied in both the statistical and reinforcement learning (RL) literature. Many authors have evaluated and contrasted policies in terms of the expected sum of rewards over a finite number of time points (Murphy et al. 2001; Chakraborty & Moodie 2013; Jiang & Li 2015). However, because these methods often use products of weights with probabilities from the behavior policy in the denominator, the extension to problems with a large number of time points often suffers from a large variance (Thomas & Brunskill 2016; Jiang & Li 2015).

Closest to the setting of this work is the recent work by Murphy et al. 2016 and Liu et al. 2018. Murphy et al. 2016 considered the average reward setting. They assumed a linear model for the value function and constructed the estimating equations to estimate the average reward. However the linearity assumption of the value function is unlikely to hold in practice and difficult to validate (e.g., the value function involves the infinite sum of the rewards). Our method allows the use of a non-parametric model for the value function to increase robustness. Liu et al. 2018 also considered the average reward and proposed an estimator for the average reward based on estimating the ratio of the stationary distribution under the target policy divided by the stationary distribution under the behavior policy. However they did not provide confidence intervals or other inferential methods besides an estimator for the average reward. In addition, they restricted the behavior policy to be Markovian and time-stationary. In mHealth, the behavior policy can be determined by an algorithm based on the accruing data and thus violates this assumption (Liao et al. 2018; Dempsey et al. 2020).

Estimator for Off-Policy Evaluation

We assume that the dataset, $\mathcal{D}_{n}$, consists of $n$ trajectories:

Below we introduce the estimator for $\eta^{\pi}$. We follow the so-called “model-free” approach (i.e., does not require modeling the transition kernel, $P$) to estimate the average reward. Our estimator is based on the Bellman equation, also known as the Poisson equation (Puterman 1994); as will be discussed below this equation characterizes the average reward.

First consider the setting where the state space, $\mathcal{S}$, is finite and the induced Markov chain, $P^{\pi}$, is irreducible. Recall that in this setting the average reward, $\eta^{\pi}$, is a constant given in (2). Define the relative value function by

The average reward of the target policy, $\pi$, is independent of state and satisfies (2). $(\eta^{\pi},Q^{\pi})$ is the unique solution of the Bellman equation (4) up to a constant for $Q^{\pi}$. The stationary distribution of the induced transition kernel, $P^{\pi}$, exists.

The above argument motivates a coupled estimator in which we use the estimated Bellman error to form an objective function. In particular, for each $(\eta,Q)$, we replace the Bellman error, $\mathcal{T}_{\pi}(S_{t},A_{t};Q)-\eta-Q(S_{t},A_{t})$, in (6) by an estimate of the “projection” of the Bellman error into a second function class, $\mathcal{G}$:

where for each $(\eta,Q)$, $\hat{g}_{n,\pi}(\cdot,\cdot;\eta,Q)$ is an estimator for the projection of the Bellman error given by

We can see that for every $(\eta,Q)$, $\hat{g}_{n,\pi}(\cdot,\cdot;\eta,Q)$ is a penalized estimator for the projected Bellman error $g^{*}_{\pi}(\cdot,\cdot;\eta,Q)$ in (7). On the other hand, the objective function in (8) is a plug-in version of the objective function in (6) where we replace the Bellman error by $\hat{g}_{n,\pi}(\cdot,\cdot;\eta,Q)$. Compared to the classic empirical risk minimization, $(\hat{\eta}_{n}^{\pi},\hat{Q}^{\pi}_{n})$ solves a nested optimization problem in the sense that the objective function (8) depends on $\hat{g}_{n,\pi}(\cdot,\cdot;\eta,Q)$ which itself is the solution of another, lower-level optimization (9).

The penalty term, $\lambda_{n}J_{1}^{2}(Q)$, is used to balance between the model fitting (i.e., the squared estimated Bellman error) and the complexity of the relative value function measured by $J_{1}(Q)$. Similarly, $\mu_{n}J^{2}_{2}(g)$ is used to control the overfitting in estimating the projected Bellman error when the function class, $\mathcal{G}$, is complex. In the case where the function space is $k$-th order Sobolev space, the regularizer is typically defined by the $k$-th order derivative to capture the smoothness of function. In the case where the function space is Reproducing Kernel Hilbert Space (RKHS), the regularizer is the endowed norm. In Supplement D, we provide a closed-form solution of the estimator when both $\mathcal{Q}$ and $\mathcal{G}$ are RKHSs.

So far we have focused on evaluating a single target policy. In practice, one might want to compare the target policy to some reference policy or contrast multiple target policies of interest. Suppose we are interested in $K$ different target policies, $\left\{\pi_{j}\right\}_{j=1}^{K}$. The above procedure (8) can be applied to estimate $\left\{\eta^{\pi_{j}}\right\}_{j=1}^{K}$. In the next section, we will provide the result of the joint asymptotic distribution of $\left\{\hat{\eta}_{n}^{\pi_{j}}\right\}_{j=1}^{K}$ (see Corollary 1). This can be used, for example, to construct the confidence interval of the difference of the average rewards between two policies.

Theoretical Results

The assumption of a bounded reward is mainly to simplify the proof and can be relaxed to the sub-Gaussian case, that is, the error $R_{t+1}-r(S_{t},A_{t})$ is sub-Gaussian for all $t\leq T$. The boundedness assumption on the shifted relative value function can be ensured by assuming certain smoothness assumptions on the transition distribution (Ortner & Ryabko 2012) or assuming geometric convergence to the stationary distribution (Hernández-Lerma & Lasserre 1999). The boundedness assumption, (i), for members of the function class, $\mathcal{Q}$, is used to simplify the proof; a truncation argument can be used to avoid this assumption.

Recall that $g^{*}_{\pi}(\cdot,\cdot;\eta,Q)$ is a projected Bellman error in (7) into a function class, $\mathcal{G}$. We make the following assumptions about $\mathcal{G}$.

Below we make assumptions on the complexity of the function classes, $\mathcal{Q}$ and $\mathcal{G}$. These assumptions are satisfied for common function classes, for example RKHS and Sobolev spaces (Van de Geer 2000; Zhao et al. 2016; Steinwart & Christmann 2008; Györfi et al. 2006). We denote by $\cal N(\epsilon,\mathcal{F},\|\cdot\|)$ the $\epsilon$-covering number of a set of functions, $\mathcal{F}$, with respect to the norm, $\|\cdot\|$.

The upper bound on $J_{2}(g^{*}_{\pi}(\cdot,\cdot;\eta,Q))$ in (i) is realistic when the transition kernel is sufficiently smooth (see Farahmand et al. 2016 for an example of MDP satisfying this condition). We use a common $\alpha\in(0,1)$ for both $\mathcal{Q}$ and $\mathcal{G}$ in (ii) to simply the proof.

Now we are ready to state the theorem about the convergence rate for $(\hat{\eta}_{n}^{\pi},\hat{Q}_{n}^{\pi})$ in terms of the Bellman error.

Let $(\hat{\eta}_{n}^{\pi},\hat{Q}_{n}^{\pi})$ be the estimator defined in (8). Suppose Assumptions 1-5 hold and the tuning parameters, $(\lambda_{n},\mu_{n})$, satisfy $\tau^{-1}n^{-\frac{1}{1+\alpha}}\leq\mu_{n}\leq\tau\lambda_{n}$ for some constant, $\tau>0$. Then the following bounds hold with probability at least $1-\delta$,

where the leading constants depend only on $(\tau,R_{\max},Q_{\max},G_{\max},C_{1},C_{2},C_{3},\alpha)$.

In Lemma in Supplement B, we show that up to a constant, $|\hat{\eta}_{n}^{\pi}-\eta^{\pi}|\lesssim\|\mathcal{T}_{\pi}(\cdot,\cdot;\hat{Q}_{n}^{\pi})-\hat{\eta}_{n}^{\pi}-\hat{Q}_{n}^{\pi}(\cdot,\cdot)\|^{2}$ and thus $\hat{\eta}_{n}^{\pi}$ is a consistent estimator for $\eta^{\pi}$ when $\lambda_{n}=o_{P}(1)$. When the tuning parameters are chosen such that $\lambda_{n}\asymp\mu_{n}$ and $\lambda_{n}\asymp n^{-1/(1+\alpha)}$, the Bellman error at $(\hat{\eta}_{n}^{\pi},\hat{Q}_{n}^{\pi})$ has the optimal rate of convergence, i.e., $\|\mathcal{T}_{\pi}(\cdot,\cdot;\hat{\eta}_{n}^{\pi},\hat{Q}_{n}^{\pi})-\hat{\eta}_{n}^{\pi}-\hat{Q}_{n}^{\pi}(\cdot,\cdot)\|^{2}=O_{P}(n^{-1/(1+\alpha)})$. The proof of Theorem 1 is provided in Supplement B.

In the following, we provide the asymptotic distribution of the estimated average reward. This requires additional notation as follows. Define $d^{\pi}(s,a)=\pi(a\,|\,s)d^{\pi}(s)$; $d^{\pi}$ is the density of the stationary distribution of the state-action under the target policy, $\pi$. For each $t\geq 1$, denote by $d_{t}(s,a)$ the density of the state-action pair in the trajectory, $\mathcal{D}$, under the behavior policy. Let $\bar{d}_{T}(s,a)$ be the average density over $T$ decision times. Motivated by the least favorable direction in partial linear regression problems (Van de Geer 2000; Zhao et al. 2016), we define the direction function, $e^{\pi}(s,a)$, by

To see this, note that $\int Q(s,a)d^{\pi}(s,a)dsda=\int\sum_{a^{\prime}}\pi(a^{\prime}\,|\,s^{\prime})Q(s^{\prime},a^{\prime})P(s^{\prime}\,|\,s,a)d^{\pi}(s,a)dsdads^{\prime}$. The numerator in (10) is a ratio between the stationary distribution of state-action pair under target policy, $\pi$, and the average distribution of state-action pair in the trajectory, $\mathcal{D}$, under the behavior policy. The denominator is the expectation of the ratio under the stationary distribution. As a result of the denominator, we have $\int e^{\pi}(s,a)d^{\pi}(s,a)dsda=1$.

We make the following smoothness assumption about $e^{\pi}$ and $q^{\pi}$, akin to the assumptions used in partially linear regression literature (Van de Geer 2000; Zhao et al. 2016).

The last assumption is a contraction-type property. This assumption will be used to control the variance of a remainder term caused by the estimation of $Q^{\pi}$.

The parameter, $\beta$, in Assumption 7 is akin to the discount factor, $\gamma$, in the discounted reward setting. Intuitively, this is related to the “mixing rate” of the Markov chain induced by the target policy $\pi$. A similar assumption was imposed in Van Roy 1998 (Assumption 7.2 on p. 99). Now we are ready to present our main result, the asymptotic normality of the estimated average reward, $\hat{\eta}_{n}^{\pi}$.

Suppose the conditions in Theorem 1 hold. In addition, suppose Assumption 6 and 7 hold and $\lambda_{n}=a_{n}n^{-1/2}$ with $a_{n}\rightarrow 0$. The estimator, $\hat{\eta}_{n}^{\pi}$, in (9) is $\sqrt{n}$-consistent and asymptotically normal: $\sqrt{n}(\hat{\eta}_{n}^{\pi}-\eta^{\pi})\Rightarrow\textbf{N}(0,\sigma^{2})$, where

From Theorem 2, the variance in estimating the average reward parameter, $\eta^{\pi}$, depends on the length of trajectory and the ratio between the stationary distribution of the state-action pair induced by the target policy (i.e., $d^{\pi}$) and the average state-action distribution in the training data (i.e., $\bar{d}_{T}$). To gain intuition of how these impact the asymptotic variance of $\hat{\eta}^{\pi}_{n}$, consider a simplified setting where the conditional variance of $R_{t+1}+\sum_{a^{\prime}}Q^{\pi}(S_{t+1},a^{\prime})-\eta^{\pi}-Q^{\pi}(S_{t},A_{t})$ given $(S_{t},A_{t})$ is a constant, denoted by $\sigma_{0}^{2}$. It can be shown that the asymptotic variance becomes $\sigma^{2}=\frac{\sigma_{0}^{2}}{T}(1+\|(d^{\pi}/\bar{d}_{T})-1\|^{2})$. Thus the smaller $\|(d^{\pi}/\bar{d}_{T})-1\|^{2}$ (i.e., the ratio, $d^{\pi}/\bar{d}_{T}$, close to one), the smaller the asymptotic variance of the estimated average reward. Although here we focus only on the asymptotic properties of $\hat{\eta}_{n}^{\pi}$ for large $n$ (recall $n$ is the number of i.i.d. trajectories), one can see that increasing length of the trajectory, $T$, reduces the asymptotic variance.

Now we present the result for evaluating a class of policies, $\Pi=\{\pi_{1},\dots,\pi_{K}\}$. Denote by $\hat{\eta}_{n}^{\pi_{j}}$ the estimated average reward of the policy, $\pi_{j}$, using (8).

Suppose the conditions in Theorem 1 and 2 hold for each $\pi\in\Pi$. Let $\epsilon_{t}^{\pi}=\frac{d^{\pi}(S_{t},A_{t})}{\bar{d}_{T}(S_{t},A_{t})}[R_{t+1}+\sum_{a^{\prime}}\pi(a^{\prime}\,|\,S_{t+1})Q^{\pi}(S_{t+1},a^{\prime})-\eta^{\pi}-Q^{\pi}(S_{t},A_{t})]$ for each $\pi\in\Pi$. Then the estimated average rewards, $\{\hat{\eta}_{n}^{\pi_{1}},\dots,\hat{\eta}_{n}^{\pi_{K}}\}$, jointly converge in distribution to a multivariate Gaussian distribution:

Simulation

In this section, we conduct a simulation study to evaluate the performance of the proposed method. The generative model is given as follows. We follow the state generative model in Luckett et al. 2020. Specifically, the state, $S_{t}=(S_{t,1},S_{t,2})$, is a two-dimensional vector and the action, $A_{t}\in\{0,1\}$, is binary. Given the current state, $S_{t}$, and action, $A_{t}$, the next state, $S_{t+1}=(S_{t+1,1},S_{t+1,2})$, is generated by $S_{t+1,1}=(3/4)(2A_{t}-1)S_{t,1}+(1/4)S_{t,1}S_{t,2}+\boldsymbol{N}(0,0.5^{2})$ and $S_{t+1,2}=(3/4)(1-2A_{t})S_{t,2}+(1/4)S_{t,1}S_{t,2}+\boldsymbol{N}(0,0.5^{2})$. Note that receiving a treatment $(A_{t}=1)$ increases the value of $S_{t,1}$ while decreases $S_{t,2}$. The reward is generated by $R_{t+1}=S_{t+1,1}+(1/2)S_{t+1,2}+(1/4)(2A_{t}-1)$. For each trajectory in the training data, the state variables are generated as independent standard normal random variables and the behavior policy is to choose $A_{t}=1$ with a fixed probability 0.5. We evaluate and compare two natural policies: the “always treat” policy, $\pi_{1}(a\,|\,s)=1$, and “no treatment” policy, $\pi_{2}(a\,|\,s)=0$.

We consider different scenarios of the number of the trajectories, $n\in\{25,40\}$, and the length of each trajectory, $T\in\{25,50,75\}$. In each scenario, we generate 500 simulated dataset and for each dataset we construct the $95\%$ confidence intervals of $\eta^{\pi_{1}},\eta^{\pi_{2}}$ and $\eta^{\pi_{1}}-\eta^{\pi_{2}}$. The coverage probability of each confidence interval is calculated over 500 repetitions. The simulation result is reported in Table 1. When the number of trajectories is small (i.e., $n=25$), the simulated coverage probability is slightly smaller than the claimed value, 0.95, especially when the length of the trajectory, $T$, is small. It can be seen that the coverage probability slightly improves when $T$ increases. When $n=40$, the coverage probability becomes closer to 0.95 as desired. Overall, the simulation result demonstrates the validity of the inference and the selection procedure for the tuning parameters. It suggests that it is necessary to perform a small-sample correction when both $n$ and $T$ are small. This is left for future work.

Case Study: HeartSteps

We apply the method to the data collected in the first study in HeartSteps (Klasnja et al. 2015; Liao et al. 2016; Klasnja et al. 2019). Below we refer to this study by HS1 for simplicity. HS1 was a 42-day MRT with 44 healthy sedentary adults. We focus on the activity suggestion intervention component. There were five individual-specified times in a day which were roughly separated by 2.5 hours and corresponded to the individual’s morning commute, mid-day, mid-afternoon, evening commute, and post-dinner times. At each decision time, an activity suggestion was sent with a fixed probability 0.6 only if the participants were considered to be available for treatment. For example, the participants were considered unavailable when they were currently physically active (e.g., walking or running) or driving a vehicle. The activity suggestions were intended to motivate near-time walking. Each participant wore a Jawbone wrist tracker and the minute-level step count data was recorded.

We construct the state based on the participant’s step count data (e.g., the 30-min step count prior to the decision time and the total step count from yesterday), location, temperature and number of the notifications received over the last seven days. We also include in the state the time slot index in the day (1 to 5) and the indicator measuring how the step count varies at the current time slot over the last seven days. The reward is formed by the log transformation of the total step count collected in 30-min window after the decision time. The log transformation is performed as the step count data is positively skewed (Klasnja et al. 2019). The step count data might be missing because the Jawbone tracker recorded data only when there were steps occurred. We use the same imputation procedure as in Klasnja et al. 2019. The state related to the step count are constructed based on the imputed step counts. The variables in the state are chosen to be predictive of the reward. In particular, each variable is selected, at the significance level of 0.05, based on a marginal Generalized Estimating Equation (GEE) analysis. In the analysis, we exclude seven participants’ data as in the primary analysis in Klasnja et al. 2019 (three due to technical issues and four due to early dropout). In addition, from the 37 participants’ data we exclude the decision times when participants were traveling abroad or experiencing technical issues or when the reward (i.e., post 30-min step count) is considered as missing (see Klasnja et al. 2019 for details).

We consider three target policies. The first policy, $\pi_{\text{nothing}}$, is “do nothing”. The second policy, $\pi_{\text{always}}$, is the “always treat” policy. Recall that in HeartSteps the activity suggestion can be sent only when the participant is available. So here the “always treat” policy refers to always send the suggestion whenever the participant is available. The third policy, $\pi_{\text{location}}$, is based on the location. Specifically, we consider the policy that sends the activity suggestion when the participant is at either home or work location and available. This policy is of interest because people at home or work are in a more structured environment and thus might be able to better respond to an activity suggestion as compared with at other locations. In HS1, about $44\%$ of the available decision times were at times that the participants were at their home or work location. Thus the policy, $\pi_{\text{location}}$, is different from the “always treat” policy, $\pi_{\text{always}}$.

In the implementation, we use the RKHS with the radial basis function kernel to form the function classes, $\mathcal{Q}$ and $\mathcal{G}$. The tuning parameters are selected based on the procedure described in Section 6. The estimated average reward of the location-based policy, $\pi_{\text{location}}$, is $3.155$ with the 95% confidence interval , $[2.893,3.417]$, which is slightly better than the “do nothing” policy. Specifically, the estimated average reward of $\pi_{\text{nothing}}$ is $2.962$ and the 95% confidence interval of the difference, $\eta^{\pi_{\text{location}}}-\eta^{\pi_{\text{nothing}}}$, is $[-0.016,0.402]$. Translating back to the raw step count as in Klasnja et al. 2019, the location-based policy is able to increase the average 30-min step count roughly by $22\%$ (i.e., $\exp(3.16-2.96)-1=1.22$), corresponding to $55$ steps (the mean post-decision time step count is 248 across all decision times in the dataset). However if we compare the “always treat” policy ($\hat{\eta}^{\pi_{\text{always}}}=3.127$, 95% confidence interval is $[2.840,3.413]$) with the location-based policy, $\pi_{\text{location}}$, we see no indication that providing treatment only at home or work is better than always providing treatment (the 95% confidence interval of $\eta^{\pi_{\text{location}}}-\eta^{\pi_{\text{always}}}$ is $[-0.161,0.217]$). Recall that the sample size for this study is $n=37$ thus this non-significant finding may be due to the small sample.

Discussion

In this work we developed a flexible method to conduct inference about the the long-term average outcomes for given target policies using data collected from a possibly different behavior policy. We believe that this is an important first step towards developing data-based just-in-time adaptive interventions. Below we discuss some directions for future research.

In many MRT studies, a natural choice of the proximal outcome to assess the effectiveness of the intervention is binary. For example, in the Substance Abuse Research Assistance study (Rabbi et al. 2018), the proximal outcome was whether the individual completed a daily survey. An interesting open question is how to extend the method to the binary reward setting, which would require carefully choosing the model to represent the relative value function and/or the loss functions used in estimating the Bellman error and solving the Bellman equation.

Non-stationarity occurs mainly because of the unobserved aspects of the current state (e.g., the engagement and/or burden) in many mHealth applications. It will be interesting to generalize the average reward framework to incorporate the non-stationarity detected in the observed trajectory. Alternatively, one can consider evaluating the treatment policy in the indefinite horizon setting where there is an absorbing state (akin to the individual disengaging from the mobile app) and thus we aim to conduct inference about the expected total rewards until the absorbing state is reached.

We focused on evaluating and contrasting multiple pre-specified treatment policies. An important next step is to extend the method to learn the optimal policy that would lead to the largest long-term average reward and to develop the inferential methods to assess the usefulness of certain variables in the policy.

Acknowledgment

This work was supported by National Institute on Alcohol Abuse and Alcoholism (NIAAA) of the National Institutes of Health under award number R01AA23187, National Institute on Drug Abuse (NIDA) of the National Institutes of Health under award numbers P50DA039838 and R01DA039901, National Institute of Biomedical Imaging and Bioengineering (NIBIB) of the National Institutes of Health under award number U54EB020404, National Cancer Institute (NCI) of the National Institutes of Health under award number U01CA229437, and National Heart, Lung, and Blood Institute (NHLBI) of the National Institutes of Health under award number R01HL125440. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

References