
\section{Introduction}

There is growing interest in using online Reinforcement Learning (RL) to optimize the delivery of messages or other forms of prompts in digital  interventions. In mobile health, RL algorithms have been used to increase the effectiveness of the content and timing of intervention messages designed to promote physical activity \citep{DBLP:journals/corr/abs-1909-03539,yom2017encouraging,hochberg2016reinforcement} or to manage weight loss \citep{forman2019can}.
In other areas including the social sciences and education, RL algorithms are used to provide pre-trial nudges to encourage court hearing attendance \citep{stanford:pre-trialnudge}, to personalize math explanations \citep{cai2021bandit} or quiz questions during lecture videos \citep{NEURIPS2018_d8c9d05e}. Unlike areas such as games and some subareas of robotics, it can be extremely costly to run a digital intervention study. Further, when the study is a pre-registered clinical trial, once initiated, the trial protocol (including any online algorithms) cannot be altered, without jeopardizing trial validity. Thus design decisions are a "one-way door" \citep{bezos:1997};
namely, once we commit to a set of design decisions, they are irreversible for the duration of the trial. To prevent poor design decisions that could be detrimental to the effectiveness and the validity of study results, RL algorithms must undergo a thorough design and testing process before deployment. 

The development of an RL algorithm in digital interventions requires a multitude of design decisions. These decisions include how best to accommodate the lack of mechanistic models for dynamic human responses to digital interventions and how to ensure robustness of the algorithm to potentially non-stationary/non-Markovian outcome distributions. Further, one must not only ensure that the RL algorithm learns and quickly optimizes interventions, but also ensure that the algorithm runs stably and autonomously online within constrained amounts of time. These include decisions to ensure that the RL algorithm can obtain data in a timely manner. Time and budgetary considerations may restrict the ability to implement more complex RL algorithms. Further, it is important to ensure the data collected by the RL algorithm can be used to inform future studies and address scientific, causal inference questions. Addressing these challenges in a reproducible, replicable manner is critical if RL algorithms are to play a role in optimizing digital interventions. Therefore we need a framework for making
design decisions for RL algorithms intended to optimize digital interventions.
\\

\noindent The primary contributions of this work are two-fold: 
\begin{enumerate}
  \item \textbf{Framework for Guiding Design Decisions in RL Algorithms for Digital Interventions:} {\it We provide a framework for evaluating the design of an online RL algorithm to increase confidence that the inclusion of an online RL algorithm as part of the digital intervention will improve the effectiveness and maintain reproducibility and replicability of the intervention, in real-life implementation. }  
 
  Specifically, we extend the PCS (Predictability, Computability, Stability) data science framework of Yu \cite{yu2020veridical} to address the specific challenges in the development and evaluation of an online RL algorithm for personalizing digital interventions. 
    
  \item \textbf{Case Study:} 
 This case study concerns the development of an RL algorithm for Oralytics, a mobile health intervention study designed to encourage oral self-care behaviors. The study is planned to go into the field in late 2022.  {\it This case study provides a concrete implementation of the PCS framework in informing the design of an online RL algorithm.}  

 One tool is the development of multiple  simulation test-bed environments and the use of these environments along with PCS principles to evaluate the candidate RL algorithms.  
\end{enumerate}


\section{Review of Online Reinforcement Learning Algorithms}
\label{section:rl}
\begin{figure}[t]
    \centering
    \includegraphics[width=1.0\textwidth, trim={0 5cm 0 0cm},clip]{figures/rl_prob_def.pdf}
\caption{\textbf{Decision times and actions for the RL algorithm in the Oralytics study.} For each day $d$, there is a morning decision time and an evening decision time for user $i$. At every decision time $t$, the RL algorithm selects action $A_{i, t}$ given the current state $S_{i, t}$. After an action is taken, we observe reward $R_{i, t}$ from the user.}
\label{figs/rl_problem_def}
\end{figure}

Reinforcement learning (RL) \cite{sutton2018reinforcement,den2020reinforcement} is the area of machine learning which is concerned with learning how to best make a sequence of decisions. In digital intervention settings, the sequence of decisions concerns which treatment (e.g., motivational messages, reminders, types of feedback, etc.) to provide given the user's current state and history. Definitions of decision times, state, action, reward, and update times are provided below. The decision times and actions for the RL algorithm for Oralytics are provided in Figure~\ref{figs/rl_problem_def}.

\subsection{Decision Times} These are the times, indexed by $t$ at which the RL algorithm may deliver a treatment (via a smart device such as a desktop computer, smartwatch, smartphone, smart speaker, wearable, etc.). The cadence of the decision times (minute level, hourly, daily, etc.) depends on the type of digital intervention. For example in Oralytics, we have two decision times per day, namely, one hour prior to the set morning and evening brushing windows specified by the user.

\subsection{State} 
$S_{i,t} \in \mathbb{R}^d$ represents the $i$th user's state at time $t$.  $d$  is the number of features describing the user's state (e.g. current location, recent adherence to medication, current social setting, recent engagement with the intervention application, etc.). See Section~\ref{case_study_prob_def} for the state definition for Oralytics.

\subsection{Action}
$A_{i, t} \in \mathcal{A}$ represents the decision made by the RL algorithm for the $i$th user at decision time $t$. Treatment actions in digital interventions frequently include the action of not delivering any treatment at time $t$. For Oralytics, the action space is $\mathcal{A} := \{ 0, 1 \}$, where $A_{i, t} = 1$ represents sending the user an engagement message and $A_{i, t} = 0$ represents not sending an engagement  message. See Section~\ref{oralytics_context} for descriptions of the types of messages that could be sent in Oralytics.

\subsection{Reward} 
$R_{i, t} \in \mathbb{R}$ is the reward for the $i$th user at decision time $t$, recorded after taking action $A_{i, t}$. The definition of the reward depends on the type of digital intervention. Examples include successful completion of a math problem, taking a medication, level of physical activity, etc. The reward function is the mean of $R_{i, t}$ conditional on the current state $S_{i, t}$ and action $A_{i, t}$. In Oralytics the reward is the subsequent brushing duration; see Section~\ref{case_study_prob_def} for further discussion of the reward in Oralytics.   

\subsection{Update Times}
These are the times at which the RL algorithm is updated, which typically includes updating a model of the user (e.g., a model of the user's reward function). The RL algorithm updates using user $i$'s current history of past states, actions, and rewards up to time $t$, denoted $H_{i, t - 1} = \left\{ S_{i,s}, A_{i,s}, R_{i,s} \right\}_{s=1}^{t-1}$. Also, if the algorithm pools data across users, then the history of other users in the study, $H_{j, t - 1}$ for $i \not= j$, are used to update the model for user $i$. For Oralytics, the update cadence is once a week.

\paragraph{}
Generally, online RL algorithms are composed of two parts: (a) fitting a model of the user and (b) an action selection strategy. The simplest type of user model is a model for the reward function, $\mathbb{E}[R_{i, t} | S_{i, t}, A_{i, t}]$. In more general cases, a model for the sum of future rewards, conditional on current state, $S_{i, t}$ and action, $A_{i, t}$ is also learned. The action selection strategy of the RL algorithm uses the user's current state $S_{i, t}$, along with the learned user model, and outputs the treatment action $A_{i, t}$ at each decision time $t$. The user model is periodically updated using newly collected data $H_{i, t - 1}$; these updates can occur after each decision time or at longer time scales. For example in \cite{DBLP:journals/corr/abs-1909-03539}, the decision times are 5 times per day, but the update times are only nightly. 

An online RL algorithm should quickly learn which action to deliver in which states for each user. One of the most widely used and simplest RL algorithms is a contextual bandit algorithm \cite{wang2005bandit,langford2007epoch,tewari2017ads}. A contextual bandit algorithm, incrementally, as data accrues,  learns the action that will lead to maximal reward in each state.  The online algorithm sequentially updates an estimate of the reward function (mean of the reward conditional on state and action) and selects actions.  The performance of the algorithm is often measured by the sum of rewards; the faster the algorithm learns the higher the sum. 
\section{ PCS Framework For Designing RL Algorithms for Digital Intervention Development}
\label{section:pcs}

The PCS Framework \citep{yu2020veridical} incorporates best practices from machine learning and applied statistics to provide a set of guidelines for assessing the quality of a prediction algorithm when used to address problems in real life. The goal is to enhance the trust of the scientific community in the application of the prediction algorithm in terms of predictability of results, computability in the implementation of the algorithm, and stability in the performance results of the learning algorithm across perturbations. PCS has been adopted and extended to other domains; see Section~\ref{related_works} for a further discussion.

Similar to the prediction setting, there are a variety of design decisions one needs to make before deploying an RL algorithm, e.g., choosing the model class to use to approximate the reward function. While many of the original PCS principles can be used in the development and evaluation of RL algorithms, RL algorithm development also introduces new challenges for the PCS framework, particularly in the online digital intervention setting. First, the main task is not prediction, but rather in RL, the main goal is to select intervention actions so that average rewards across time are maximized for each user. We call this the goal of Personalization \cite{fan2006personalization}. We generalize the PCS framework to include an evaluation of the ability of an online RL algorithm to personalize. Second, in digital intervention settings, it is important to evaluate the ability of the online RL algorithm to maximize rewards under real-world constraints. For example, there are often time constraints on computations, budgetary constraints on software engineering development, and constraints on the algorithm in terms of obtaining data in real-time. Further, the algorithm must run stably online without constant human monitoring and adjustment.  The current PCS framework does not provide evaluation tools that deal with the above needs.
We extend the PCS framework to the context of designing and evaluating online RL algorithms. This  framework focuses on providing confidence that the online RL algorithm will lead to greater effectiveness under real-world constraints and with stability.   

\subsection{Personalization (P)} 
\label{pcs:P}
The PCS Framework uses (P) for predictability with the goal of ensuring that models used in data science have good predictive accuracy on both seen and unseen data. Predictive accuracy is a simple, commonly used metric for evaluating the model, but in some cases, multiple evaluation metrics or domain-specific metrics are appropriate. 
In our setting, the main task is Personalization. Namely, the online RL algorithm should learn to select actions to maximize {\it each} user's average rewards. Instead of defining a predictive accuracy metric, we want a metric to validate the extent of personalization. For example, when choosing a metric to evaluate RL algorithms for multiple users, one may be interested in not just the average over the users' sums of rewards, but other metrics that capture the variation in the sum of rewards across users. Let $N$ be the total number of users with $T$ total decision times. We suggest the following metrics:

\begin{itemize}
    \item \textbf{Average of Users' Average (Across Time) Rewards:} First average users' rewards across time, and then average across users. This metric is defined as $\frac{1}{N} \sum_{i=1}^N \big ( \frac{1}{T} \sum_{t=1}^T R_{i, t} \big)$. The metric serves as a global measure of the RL algorithm's performance.
    \item \textbf{The 25th Percentile of Users' Average (Across Time) Rewards:} The metric shows how well an RL algorithm performs for users that don't benefit as much, namely users in the lower quartile of average rewards across time.
    \item \textbf{Average Reward For Multiple Time-Points:} Average users' rewards across time for multiple time-points $t_0 = 1, 2,..., T$. Then average across users. This metric is defined as $\frac{\sum_{i=1}^N\sum_{t=1}^{t_0} R_{i, t}}{N \cdot t_0}$. The metric assesses the speed at which the RL algorithm learns across weeks in the trial.
\end{itemize}

\subsection{Computability (C)}
\label{pcs:C}
Computability focuses on the efficiency and scalability of algorithms, decisions, and processes. While the original PCS framework focused on the computability of training and evaluating models,
we also consider the ability to implement the algorithm within the constraints of the study. In the online RL setting, computability encompasses all issues related to ensuring that the RL algorithm can select actions and update in a timely manner while running online. 
The performance of the online RL algorithm must be evaluated under the constraints of the study;  key RL design constraints that could arise include:
\begin{itemize}
    \item \textbf{Timely Access to Reward and State Information: } Since RL algorithms for digital interventions must make decisions in an online fashion, the development team must ensure that the state feature data is available at each decision time and that both the reward and state feature data is available to the algorithm at the update times. For example, due to delays in communication between sensors, the digital application, and the cloud storage, the algorithm may not have timely access to the investigators' first choice of reward, necessitating the use of an alternate.
    \item \textbf{Engineering Budget: } One should consider the engineering budget, supporting software needed, and time available to deliver a production-ready algorithm. In this case, a simpler algorithm may be preferred over a sophisticated one because it is easier to implement, test, and set up monitoring systems for.
    \item \textbf{Off-Policy Evaluation and Causal Inference Considerations: } The investigative team often not only cares about the RL algorithm's ability to learn but also the ability to use the data collected by the RL algorithm to answer scientific questions after the study is over. These scientific questions can include topics such as off-policy evaluation \cite{thomas2016data,levine2020offline} and causal inference \cite{boruvka2018assessing,hadad2021confidence}. Thus, the algorithm may be constrained to select actions probabilistically with probabilities that are bounded away from zero and one. This enhances the ability of investigators to use the resulting data to address scientific questions with sufficient power \cite{pmlr-v149-yao21a}.
\end{itemize}

\subsection{Stability (S)}
\label{pcs:S}
Stability measures how the results change with minor perturbations and also emphasizes the need for documentation and reproducibility of results.  In online RL, stability plays two roles.  
First, the RL algorithm must run stably and automatically  without the need for constant human monitoring and adjustment. This is particularly critical as users abandon digital interventions that have inconsistent functionality (unstable RL algorithm) \cite{10.1145/2800835.2800943,info:doi/10.2196/jmir.2583}.
Second, the RL algorithm should perform well across a variety of potential real-world environments. 
A critical tool in assessing stability to perturbations of the environment is the use of simulation test-beds. Test-beds include a variety of plausible environmental variants, each of which encodes different concerns of the investigative team. Often important state features of the user's environment are not observed and we lack an accurate mechanistic model of user behavior. The following challenges in digital intervention problems are concerns that one could design testbeds for:

\begin{itemize}
    \item \textbf{User Heterogeneity:} There is likely some amount of user heterogeneity in response to actions, even when users are in the same context. User heterogeneity can be partially due to unobserved user traits (e.g. factors that are stable or change very slowly over time like family composition or personality type). 
    The amount of between user heterogeneity impacts whether an RL algorithm that pools data (partially or using clusters) across users to select actions will lead to improved rewards. 
    \item \textbf{Non-Stationarity:} Unobserved factors common to all users such as societal changes (e.g. new wave of the pandemic), and time-varying unobserved treatment burden (e.g. user's response to intervention may depend on how many days they have experienced the digital intervention) may make the distribution of the reward appear to vary with time, i.e., non-stationary.
    \item \textbf{High Noise Environments} 
    Digital interventions typically deliver treatments to users in highly noisy environments. This is in part because digital interventions deliver treatments to users in daily life, where many unobserved factors (e.g. social context, mood, or stress) can affect a user's responsiveness to an intervention. If unobserved, these factors produce noise. Moreover, the effect of digital prompts on a near-term reward tends to be small due to the nature of the intervention.  
    Therefore it is important to evaluate the algorithm's ability to personalize even in highly noisy, low signal-to-noise ratio environments.
\end{itemize}

\subsection{Simulation Environments for PCS Evaluation}
\label{section:sim_guidelines}
To utilize the PCS framework, we advocate using a simulation environment for designing and evaluating RL algorithms. We aim to compare RL algorithm candidates under real-world constraints (Computability). Thus, we build multiple variants of the simulation environment, each reflecting plausible user dynamics (i.e. state transitions and reward distributions) (Stability). This is followed by simulating a digital intervention study for each simulation environment variant and RL algorithm candidate pairing. Finally, we use multiple metrics to evaluate the performance of the RL algorithm candidates (Personalization). 

In the case of digital interventions, there is often no mechanistic model or physical process for user behavioral dynamics, which makes it difficult to accurately model transitions (e.g. modeling a user's future level of physical activity as a function of their past physical activity, location, local weather). 
Note that the goal of developing the simulators is not to conduct model-based RL \cite{agarwal2021persim}.
Rather, here the simulators represent a variety of plausible environments to facilitate the evaluation of the performance of potential RL algorithms in terms of Personalization, Computability, and Stability across these environments. Existing data and domain expertise is most naturally used to construct the simulation environments. However, as is the case for Oralytics, the previously collected data may be scarce, i.e., we have very few data points per user. Moreover, the data may only be partially informative, e.g., the data was collected under only a subset of the actions.
Next, we provide guidelines for how to build an environment simulator in such challenging settings.

\underline{Base Environment Simulator:}
In order to have the best chance possible of accurately evaluating how well different RL candidates will perform, we recommend first building a base environment simulator that mimics the existing data to the greatest extent possible.
This involves carefully choosing the set of time-varying features and reward generating model class that will be expressive enough to model the true reward distribution well.
To check how well simulated data generated by the model of the environment mimics the observed data, we recommend a variety of ways to compare distributions. 
This includes visual comparisons such as plotting histograms, comparing moments such as mean reward, between user variance, and within user variance, of the real data with the simulated data, and measuring how well the base model captured the variance in the data.
Examples of these checks done for Oralytics are in Appendix~\ref{sim_env_evaluation}.

\underline{Variant Environment Simulators:}
We recommend considering many variants or perturbed simulation environments to evaluate the stability of RL algorithms across multiple plausible environments. These variants can be used to evaluate the concerns of the investigative team.
For example, if the base simulator generates stationary rewards and the investigative team is concerned that the real reward distribution could not be stationary, a variant could incorporate non-stationarity into the environment dynamics. 

In the case that the previously collected data does not include particular actions, as was the case for Oralytics, we recommend consulting domain experts to recommend a range of potential realistic effect sizes (differences in mean reward under the new action versus a baseline action). For example, in Oralytics we only have data under no intervention and did not have data on rewards under the intervention. Thus, using the input of the domain experts on the team, we imputed several plausible treatment effects (varying by certain state features and the amount of heterogeneity in treatment effects across users). 
\section{Related Works}
\label{related_works}

\subsection{Digital Intervention Case Studies}

Liao \cite{DBLP:journals/corr/abs-1909-03539} describe the development of an online RL algorithm for HeartSteps V2, a physical activity mobile health application. They highlight how their design decisions address specific challenges in designing RL algorithms, such as adjusting for longer-term effects of the current action and accommodating noisy data. However, they do not provide general guidelines for how to make design decisions for RL algorithms in digital intervention development. 

Another related work is that of Figueroa \cite{figueroa2021adaptive} who provide an in-depth case study of the design decisions, and the associated challenges and considerations, for an RL algorithm for text messaging in a physical activity mobile application serving patients with diabetes and depression. Through this case study, they provide guidelines to others developing RL algorithms for mobile health applications.
Specifically, they first categorize the challenges they faced into 3 major themes: 1) choosing a model for decision making, 2) data handling, and 3) weighing algorithm performance versus effectiveness in the real world. Then for each of these challenges, they describe how they dealt with the challenges in the design process for their physical activity intervention. In contrast, by expanding the PCS framework, this work introduces general guidelines for comprehensively evaluating RL algorithms. Moreover, we make recommendations for how to design a variety of simulation testbeds even using only sparse and partially informative existing data in service of PCS. The generality of the PCS framework allows it to be more widely applicable. For example, Figueroa \cite{figueroa2021adaptive} have an existing dataset for all actions, which makes their recommendations less applicable to those designing algorithms with existing data only under a subset of actions. The PCS framework provides more holistic guidelines for facing challenges in developing simulation environments and evaluating algorithms, beyond those faced in any particular case study.

\subsection{Simulation Environments in Reinforcement Learning}
In RL, simulators  (generative models) may be used to derive a policy from the generative model underlying the simulator (model-based learning). Agarwal \cite{agarwal2021persim} use simulation as an intermediate step to learn personalized policies in a data-sparse regime with heterogeneous users, where they only observe a single trajectory per user. Wei \cite{wei2020learning} propose a framework for simulating in a data-sparse setting by using imitation learning to better interpolate traffic trajectories in an autonomous driving setting.   
In contrast, in PCS, the simulator is used as a crucial tool for using the framework to design, compare, and evaluate RL algorithm candidates for use in a particular problem setting.

There exist many resources aimed at improving the design and evaluation of RL algorithms through simulation; however, in constrast to this work, they do not provide guidelines for designing plausible simulation environments using existing data. RecSim \cite{ie2019recsim} gives a general framework but does not advise on the quality of the environment nor how to make critical design decisions such as reward construction, defining the state space, simulating unobserved actions, etc. MARS-Gym \cite{santana2020marsgym} provides a full end-to-end pipeline process (data processing, model design, optimization, evaluation) and open-source code for a simulation environment for marketplace recommender systems. OpenAI Gym \cite{brockman2016openai} is a collection of benchmark environments in classical control, games, and robotics where the public can run and compare the performance of RL algorithms.

There are also a handful of papers that build simulation environment testbeds using real data. Wang \cite{wang2021optimizing} evaluate their algorithm for promoting running activity with a simulation environment built using two datasets. Singh \cite{singh2020building} develop a simulation environment using movie recommendations to evaluate their safe RL approach. Korzepa \cite{10.1145/3386392.3399291} use a simulation environment to guide the design of personalized algorithms that optimize hearing aid settings. Hassouni \cite{10.1007/978-3-030-03098-8_31,https://doi.org/10.48550/arxiv.1804.03592} fit a realistic simulation environment using the US timekeeping research project data. Their simulation environment creates daily schedules of activities for each user (i.e. sleep, work, workout, etc.) where each user is one of many different user profiles (i.e. workaholic, athlete, retiree) for the task of improving physical activity. 


\subsection{PCS Framework Extensions}
The PCS framework has been extended to other learning domains such as causal inference \cite{https://doi.org/10.1111/insr.12427}, network analysis \cite{ward2021next}, and evaluating the interpretability of algorithms \cite{margot2021new}. Despite the variety of these tasks, they can all be framed as supervised learning problems in batch data settings that can be evaluated in terms of prediction accuracy on a hold-out dataset. PCS has not been extended to provide guidelines for developing an online decision-making algorithm. This extension is needed because of the additional considerations, as discussed above, present in a real-world RL setting. 
Additionally, while these papers focus on evaluating how well a model accurately predicts the outcome on training and hold-out datasets, our extension includes Personalization for evaluating how well an algorithm personalizes to each user. Dwivedi and Ward \cite{https://doi.org/10.1111/insr.12427,ward2021next} implement the original Computability principle by considering algorithm and process efficiency and scalability. Margot \cite{margot2021new} provides a new principle Simplicity which is based on the sum of the lengths of generated rules. In our case, we extend Computability to include the constraints of the study. Finally, these papers consider the stability of results across different changes to the data (e.g. bootstrapping or cross-validation) or design decisions (e.g. choice of representation space or the embedding dimension). Our framework focuses on how stable an algorithm is in plausible real world environments which could be complex (e.g. user-heterogeneity, non-stationary, high noise).
\section{Case Study: Oral Health}
\label{section:case_study}
In this case study, we illustrate the use of PCS principles in designing an RL algorithm for Oralytics. Two main challenges are 1) we do not have timely access to a lot of features and the reward is relatively noisy and 2) we have sparse, partially informative data to inform the construction of our simulation environment testbed. In addition, there are a number of study constraints. 
\begin{enumerate}
    \item Once the study is initiated, the trial protocol and algorithm cannot be altered without jeopardizing trial validity.
    \item We are using an online algorithm so we may not have timely access to certain desirable state features or rewards.
    \item We have a limited engineering budget.
    \item We must answer post-study scientific questions that require causal inference or off-policy evaluation.
\end{enumerate}

We highlight how we handle these challenges by using the PCS framework, despite being in a highly constrained setting. The case study is organized as follows. In Section~\ref{oralytics_context} we give background context and motivation for the Oralytics study. In Section~\ref{case_study_prob_def}, we explain the Oralytics sequential decision-making problem such as defining the state and reward. In Section~\ref{design_rl_alg_cands} we describe our process for designing RL algorithm candidates that can stably learn despite having a severely constrained features space and noisy rewards. Finally, in Section~\ref{design_sim_env}, we describe how we designed the simulation environment variants to evaluate the RL algorithm candidates; furthermore, we offer recommendations for designing realistic environment variants and for constructing such environments using data for only a subset of actions.


\subsection{Oralytics}
\label{oralytics_context}
Oralytics is a digital intervention for improving oral health. Each user is provided a commercially-available electric toothbrush with integrated sensors and Bluetooth connectivity as well as the Oralytics mobile application for their smartphone. There are two decision times per day (prior to the user's morning and evening brushing windows) when a message may or may not be delivered to the user via their smartphone. The types of messages focus on winning a gift for oneself, winning a gift for one's favorite charity, feedback on prior brushing, and educational information. Also, once a message is delivered to the user, the app records the messages so that a user is highly unlikely to ever receive the same message twice. Oralytics will be implemented with approximately 70 users in a clinical trial where the participant duration is 10 weeks; this means each user has $T=140$ decision times. The study duration is 2 years and the expected weekly incremental recruitment rate is around 4 users. The Oralytics mobile app will use an online RL algorithm to optimize message delivery (i.e. treatment actions) to maximize an oral health-related reward (see below). To inform the RL algorithm design, we have access to data from a prior oral health study, ROBAS 2 \citep{info:doi/10.2196/17347}, and input from experts in oral and behavioral health. The ROBAS 2 study used earlier versions of both the electric toothbrush and the Oralytics application to track the brushing behaviors of 32 users over 28 days. Importantly, in ROBAS 2, no intervention messages were sent to the users.

\subsection{The Oralytics Sequential Decision Making Problem}
\label{case_study_prob_def}

We now discuss how we designed the state space and rewards for our RL problem in collaboration with domain experts and the software team while considering various constraints. These decisions must be communicated and agreed upon with the software development team because they provide the RL algorithm with the necessary data at decision and update times, and execute actions selected by the RL algorithm.

1. \underline{Choice of Decision Times:}
We chose the decision times to be prior to each user's specified morning and evening brushing windows, as the scientific team thought this would be the best time to influence users' brushing behavior.

2. \underline{Choice of Reward:} The research team's first choice of reward was a measure of brushing quality derived from the toothbrush sensor data from each brushing episode. However, the brushing quality outcome is often not reliably obtainable because it requires (1) that the toothbrush dock is plugged in and (2) that the user is standing within a few feet of the toothbrush dock when brushing their teeth. Users could fail to meet these two requirements for a variety of reasons, e.g., the user brushes their teeth in a shared bathroom where they are unable to conveniently leave the dock plugged in. Thus, we selected brushing duration in seconds as the reward (Personalization) since 120 seconds is the dentist-recommended brushing duration and brushing duration is a necessary factor in calculating the brushing quality score. Additionally, brushing duration is expected to be reliably obtainable even when the user is far from the toothbrush dock when brushing (Computability).
Note that in Figure~\ref{figs/robas_2_dataset}, a small number of user-brushing episodes have durations over the recommended 120 seconds. Hence we truncate the brushing time to avoid optimizing for over-brushing. Let $D_{i,t}$ denote the user's brushing duration. The reward is defined as $R_{i,t} := \min(D_{i,t}, 180)$.

3. \underline{Choice of State Features At Decision Time:} 
To provide the best personalization, it is ideal that an RL algorithm has as many relevant state features as possible to make a decision, e.g., recent brushing, location, user's schedule, etc. However, our choice of the state space is constrained by the need to get features reliably before decision and update times, as well as our limited engineering budget.
For example, we originally wanted a feature for the evening decision time to be the morning's brushing outcome, however, this feature may not be accessible in a timely manner. This is because in order for the algorithm to receive the morning brushing data the Oralytics smartphone app requires the user to open the app and we do not expect most users to reliably open the app after every morning brush time before the evening brushing window. A further discussion of our choice of decision time state features can be found in Appendix~\ref{app:rl_alg_features}.

4. \underline{Choice of Algorithm Update Times}
In our simulations, we update the algorithm weekly. In terms of speed of learning (at least in idealized settings), it is best to update the algorithm after each decision time. However, due to Computability considerations, we chose a slower update cadence. Specifically, for the Oralytics app, the consideration was due to the fact that we are only able to update the policy used to select actions when the user opens the app. This means that if the user does not open the app for many days, even if we try to update the app after each decision time, we would be unable to do so. Since we believe that it could be common that users could fail to open the app for a few days at a time, we chose weekly updates. In the future, we will explore other update cadences as well, e.g., once a day.

\subsection{Designing the RL Algorithm Candidates}
\label{design_rl_alg_cands}
Here we discuss our use of the PCS framework to guide and evaluate the following design decisions for the RL algorithm candidates. There are some decisions that we've already made and other decisions that we encode as axes for our algorithm candidates to test in the simulation environment. See Appendix \ref{appendix:RLalgs} for further details regarding the RL algorithm candidates.

1. \underline{Choice of using a Contextual Bandit Algorithm Framework}
We understand that actions will likely affect a user's future states and rewards, e.g., sending an intervention message the previous day may affect how receptive a user is to an intervention message today.
This suggests that an RL algorithm that models a full Markov-Decision process (MDP) may be more suitable than a contextual bandit algorithm. However, the highly noisy environment and the limited data to learn from (140 decision times per user total), makes it difficult for the RL algorithm to accurately model state transitions. Due to errors in the state transition model, the estimates of the delayed effects of actions used in MDP-based RL algorithms can often be highly noisy or inaccurate. This issue is exacerbated by our severely constrained state space (i.e. we have few features and the features we get are relatively noisy). As a result, an RL algorithm that fits a full MDP model may not learn very much during the study, which could compromise Personalization and offer a poor user experience. To mitigate these issues, we use contextual bandit algorithms, which fit a simpler model of the environment. Using a lower discount factor (a form of regularization) has been shown to lead to learning a better policy than using the true discount factor, especially in data-scarce settings \cite{jiang2015dependence}. Thus, a contextual bandit algorithm can be interpreted as an extreme form of this regularization where the discount factor is zero. Finally, contextual bandits are the simplest algorithm for sequential decision making (Computability) and have been used to personalize digital interventions in a variety of areas
\cite{info:doi/10.2196/jmir.7994,yom2017encouraging,DBLP:journals/corr/abs-1909-03539,figueroa2021adaptive,cai2021bandit}.

2. \underline{Choice of a Bayesian Framework:} We consider contextual bandit algorithms that use a Bayesian framework, specifically Posterior (Thompson) Sampling algorithms \citep{DBLP:journals/corr/0001RKO17}. Posterior Sampling involves placing a prior on the parameters of the reward approximating function and updating the posterior distribution of the reward function parameters at each algorithm update time. This allows us to incorporate prior data and domain expertise into the initialization of the algorithm parameters. However, these naturally lead to stochastic algorithms (action selections are a not deterministic function of the data), which better facilitates causal inference analyses later on using the data collected in the study.

3. \underline{Choice of Constrained Action Selection Probabilities:} We constrain the action selection probabilities to be bounded away from zero and one, in order to facilitate off-policy and causal inference analyses once the study is over (Computability). With help from the domain experts, we decided to constrain the action selection probabilities of the algorithm to be in the interval $[0.35, 0.75]$. \\

\begin{figure}[t]
    \centering
    \includegraphics[width=0.3\textwidth,height=40mm]{figures/reward_histogram_robas_2.pdf}
\caption{Histogram of brushing durations in seconds for all user brushing sessions (twice a day) in ROBAS 2.
}
\label{figs/robas_2_dataset}
\end{figure}

The following are decisions we will test using the simulation environment.

4. \underline{Choice of the Reward Approximating Function: } An important decision in designing the contextual bandit algorithm is how to approximate the reward function. We consider two types of approximations, a Bayesian linear regression model (BLR) and a Bayesian zero-inflated Poisson regression model (ZIP) which are both relatively simple, well-studied, and well-understood. Formal specifications for BLR and ZIP as reward functions can be found in Appendix~\ref{reward_approx_functions}. We consider the ZIP because of the zero-inflated nature of brushing durations on our existing dataset ROBAS 2; see Figure~\ref{figs/robas_2_dataset}. We expected the ZIP to provide a better fit to the reward function by the contextual bandit and thus lead to increased average rewards. Moreover, the linear model for the reward function is easily interpretable by domain experts and allows them to critique and inform the model. To perform posterior sampling, both models are Bayesian with uninformative priors to both models. We discuss using informative priors further in Section \ref{sec:experiments}. From the perspective of Computability and Stability, the posterior for the BLR is of closed-form, which makes it easier to write software leading to efficient and stable updates than the case for the ZIP (the posterior distribution must be approximated and the posterior approximation scheme is another axis to algorithm design that the scientific team needs to consider); see Appendix \ref{app:algUpdatesActions} for further discussion on how to update the RL algorithm candidates. We design both types of algorithms with weekly updates to the posterior; this decision will be re-examined in the future.

5. \underline{Choice of Cluster Size: } We consider clustering users with cluster sizes $K$ = 1 (no pooling), $K$ = 4 (partial pooling), and $K = N$ (full pooling) to determine whether clustering in our setting will lead to higher sums of rewards (Personalization). Clustering-based algorithms pool data from multiple users to learn an algorithm per cluster (i.e. at update times, the algorithm uses $H_{i, t - 1}$ for all users $i$ in the same cluster, and at decision times, the same algorithm is used to select actions for all users in the cluster).
Clustering-based algorithms have been empirically shown to perform well when users within a cluster are similar \cite{zhu2018group,tomkins2021intelligentpooling}. In addition, we believe that clustering will facilitate learning within environments that have noisy within user rewards \cite{deshmukh2017multi,vaswani2017horde}. There is a trade-off between no pooling and full pooling. No pooling may learn a policy more specific to the user later on in the study, but may not learn as well earlier in the study when there is not a lot of data for that user. Full pooling may learn well earlier in the study because it can take advantage of all users' data, but may not personalize as well as a no pooling algorithm, especially if users are very heterogeneous. In addition, we chose $K = 4$ because that is the expected weekly recruitment rate for the study and the update cadence is also weekly. We consider the two extremes and partial pooling as a way to explore this trade-off. A further discussion on choices of cluster size can be found in Appendix~\ref{app:cluster_sizes}.

\subsection{Designing the Simulation Environment}
\label{design_sim_env}

We build a simulator that considers multiple variants for the environment, each encoding a concern by the research team, and allows us to evaluate the stability of results for each RL algorithm across the environmental variants (Stability). 

\underline{Fitting Base Models:}
Recall that the ROBAS 2 study did not involve intervention messages. However, we are still able to use the ROBAS 2 dataset to fit the base model for the simulation environment, i.e., a model for the reward (brushing duration in seconds) under no action. Two main approaches for fitting zero-inflated data are the zero-inflated model and the hurdle model \cite{feng2021comparison}. Throughout the model fitting process, we performed various checks on the quality of the model to determine whether the fitted model was sufficient (Appendix~\ref{sim_env_evaluation}). This included checks regarding whether the percentage number of zero brush times simulated by our model was comparable to that of the original ROBAS 2 dataset. Additionally, we checked whether the model accurately captured the mean and variance of the brushing durations for non-zero brush durations across users.

The first approach we took was to choose one model class (zero-inflated Poisson) and fit a single population-level model for all users in the ROBAS 2 study. However, a single population-level model was insufficient for fitting all users due to the high level of user heterogeneity (i.e. the between and within user variance of the simulated brushing durations from the fitted model was smaller than the between user and within user variance of brushing durations in the ROBAS 2 data). Thus, next, we decided to maintain one model class, but fit one model per user for all users. However, when we fit a zero-inflated Poisson to each user, we noticed that the model provided an adequate fit for some users, but other users exhibited more variability in their brushing durations. Namely, the within user variance simulated rewards from the model fit on those users was still lower than the within user variance of the ROBAS 2 user data used to fit the model.
Therefore, we considered a hurdle model \cite{feng2021comparison} because it is more flexible than the zero-inflated Poisson. For Poisson distributions, the mean and variance are equal, whereas in the hurdle model the mean and variance are not conflated. 

Ultimately, for each user we considered three model classes: 1) a zero-inflated Poisson, 2) a hurdle model with a square root transform, and 3) a hurdle model with a log transform, and chose one out of these three model classes for each user (Appendix~\ref{env_base_model}). Specifically, to select the model class for user $i$, we fit all three model classes using each user's data from ROBAS 2. Then we then chose the model class that had the lowest root mean squared error (Appendix~\ref{fitting_env_base_models}). Additionally, along with the base model that generates stationary rewards, we include an environmental variant with a non-stationary reward function; here day in study is used as a feature in the environment's reward generating model (Appendix~\ref{baseline_features}). 

\underline{Imputing Treatment Effect Sizes:}
To construct a model of rewards for when an intervention message is sent (for which we have no data on), we impute potential treatment effects with the interdisciplinary team and modify the fitted base model with these effects.
Specifically, we impute treatment effects on both the user's intent to brush, as well as treatment effects on the brushing duration when the user intends to brush. We impute both types of treatment effects because the investigative team's intervention messages were developed to both encourage users to brush more frequently and to brush for the recommended duration. Further, because the research team believes that the users may respond differently to the engagement messages depending on the context and depending on the user, we included context-aware, population-level, and user-heterogeneous effects of the engagement messages as environmental variants (Appendix~\ref{effect_sizes}). 

We use the following guidelines to guide the design of the effect sizes:
\begin{enumerate}
    \label{effect_size_rationale}
    \item In general for mobile health digital interventions, we expect the effect (magnitude of weight) of actions to be smaller than (or on the order of) the effect for baseline features, which include time of day and the user's previous day brushing duration (all features are specified in Appendix~\ref{baseline_features}).
    \item The variance in treatment effects (weights representing the effect of actions) across users should be on the order of the variance in the effect of features across users, i.e., looking at variance in parameters of fitted user-specific models.
\end{enumerate}

In accordance with guideline 1 above, to set the population level effect size, we take the absolute value of the weights (excluding that for the intercept term) of the base models fitted for each ROBAS 2 user and the average across users and features (e.g., the average absolute value of weight for time of day and previous day brushing duration). 
For the heterogeneous (user-specific) effect sizes, for each user, we draw a value from a normal centered at the population effect sizes. In accordance with guideline 2, the variance of the normal distributions is found by again taking the absolute value of the weights of the base models fitted for each user, averaging the weights across features, and taking the empirical variance across users.
In total there are eight environment variants as summarized in Table \ref{tab:variants}.
See Appendix \ref{appendix:sim_envs} for further details regarding the development of the simulation environments.

\begin{table}[h]
    \centering
    \begin{tabular}{|p{0.45\linewidth}|p{0.45\linewidth}|}
        \hline
        \textbf{S\_Pop:} Stationary Base Model, Population Effect Size & \textbf{NS\_Pop:} Non-Stationary Base Model, Population Effect Sizes \\ 
        \hline
        \textbf{S\_Het:} Stationary Base Model, Heterogeneous Effect Size & \textbf{NS\_Het:} Non-Stationary Base Model, Heterogeneous Effect Sizes \\ 
        \hline
    \end{tabular}
    \caption{\textbf{Four Environment Variants}  We consider 2 environment base models (stationary and non-stationary) and 2 effect sizes (population effect size, heterogeneous effect size).}
    \label{tab:variants}
\end{table}


\begin{comment}
\begin{table}[h]
    \centering
    \begin{tabular}{|p{0.5\linewidth}|p{0.5\linewidth}|}
        \hline
        \textbf{stat\_pop\_eff\_env:} Stationary Base Model and Small Population Effect Size & \textbf{non\_stat\_pop\_eff\_env:} Non-Stationary Base Model and Small Population Effect Sizes \\ 
        \hline
        \textbf{stat\_het\_eff\_env:} Stationary Base Model and Small Heterogeneous Effect Size & \textbf{non\_stat\_het\_eff\_env:} Non-Stationary Base Model and Small Heterogeneous Effect Sizes \\ 
        \hline
    \end{tabular}
    \caption{\textbf{Four Environment Variants}  We consider 2 environment base models (stationary and non-stationary) and 2 effect sizes (population effect size, heterogeneous effect size).}
    \label{tab:variants}
\end{table}
\end{comment}
\section{Experiment and Results}
\label{sec:experiments}

\begin{table}
\centering
\begin{tabular}{ |p{2.3cm}||p{2.3cm}|p{2.3cm}|p{2.3cm}|p{2.3cm}|}
 \hline
 \multicolumn{5}{|c|}{RL Algorithm Candidates} \\
 \hline
 \hline
 \multicolumn{1}{|c}{} & \multicolumn{4}{c|}{Average Rewards} \\
 \hline
 RL Algorithm & S\_Het & NS\_Het & S\_Pop & NS\_Pop \\
 \hline
      ZIP $k=1$ &  100.462 (0.789) &  103.072 (0.751) &  107.745 (0.854) &  109.892 (0.764) \\
      ZIP $k=4$ &  101.096 (0.776) &  103.665 (0.781) &  108.975 (0.819) &  110.921 (0.783) \\
   ZIP $k=N$ &  101.373 (0.799) &  104.042 (0.769) &  109.137 (0.822) &  111.295 (0.778) \\
      BLR $k=1$ &   97.910 (0.779) &  100.187 (0.750) &  104.279 (0.815) &  106.109 (0.752) \\
      BLR $k=4$ &  100.376 (0.764) &  102.868 (0.761) &  108.313 (0.837) &  110.100 (0.732) \\
   BLR $k=N$ &  \textbf{101.837 (0.806)} &  \textbf{104.655 (0.775)} &  \textbf{109.730 (0.829)} &  \textbf{112.033 (0.755)} \\
    \hline
 \multicolumn{1}{|c}{} & \multicolumn{4}{c|}{25th Percentile Rewards} \\
    \hline RL Algorithm & S\_Het & NS\_Het & S\_Pop & NS\_Pop \\
    \hline
      ZIP $k=1$ &  69.790 (1.293) &  74.240 (0.510) &  76.809 (1.193) &  79.311 (0.725) \\
      ZIP $k=4$ &  70.470 (1.151) &  74.199 (0.616) &  77.642 (1.415) &  81.491 (0.801) \\
   ZIP $k=N$ &  \textbf{71.688 (1.291)} &  75.218 (0.648) &  78.916 (1.323) &  81.194 (0.865) \\
      BLR $k=1$ &  67.975 (1.348) &  71.263 (0.680) &  73.413 (1.105) &  75.517 (0.804) \\
      BLR $k=4$ &  69.954 (1.226) &  73.871 (0.597) &  78.316 (1.326) &  80.975 (0.755) \\
   BLR $k=N$ &  71.656 (1.276) &  \textbf{75.479 (0.520)} &  \textbf{79.509 (1.332)} &  \textbf{82.732 (0.732)} \\
    \hline 
\end{tabular}
\caption{\textbf{Average and 25th Percentile Rewards.} Average and 25th percentile rewards are defined in Section~\ref{pcs:P}. The naming convention for environment variants is found in Table~\ref{tab:variants}. "BLR" denotes a Bayesian Linear Regression reward approximating function and "ZIP" denotes a Bayesian Zero-Inflated Poisson reward approximating function. "k" refers to the cluster size. Average rewards are averaged across time, users, and 50 trials. For the 25th percentile rewards, we average rewards across time, find the lower 25th percentile across users, and then averaged that across 50 trials. The value in the parenthesis is the standard error of the mean. Best performing algorithm candidates in each environment variant are bolded. 
}
\label{table/experiment_rewards}
\end{table}

\begin{figure}
    \centering
    \subfloat(a){\includegraphics[width=0.45\textwidth]{figures/across_time_reward_graphs/large_stat_het_eff_env_reward_across_time_graph.pdf}}
    \subfloat(b){\includegraphics[width=0.45\textwidth]{figures/across_time_reward_graphs/large_non_stat_het_eff_env_reward_across_time_graph.pdf}}
    \subfloat(c){\includegraphics[width=0.45\textwidth]{figures/across_time_reward_graphs/large_stat_pop_eff_env_reward_across_time_graph.pdf}}
    \subfloat(d){\includegraphics[width=0.45\textwidth]{figures/across_time_reward_graphs/large_non_stat_pop_eff_env_reward_across_time_graph.pdf}}
    \caption{\textbf{Average User Reward.} For each time step $t_0 \in [20, 40, 60, 80, 100, 120, 140]$ figure shows the mean and 1.96 * SE of average user reward ($\frac{1}{n} \sum_{i=1}^n \frac{1}{t_0} \sum_{s=1}^{t_0} R_{i,s}$) on average across trials. For \textbf{(a)} Stationary Base Model and Heterogeneous Effect Size, \textbf{(b)} Non-Stationary Base Model and Heterogeneous Effect Size, \textbf{(c)} Stationary Base Model and Population Effect Size, \textbf{(d)} Non-Stationary Base Model and Population Effect Size
    }
    \label{figs/rewards_across_time}
\end{figure}

We evaluate RL candidates in each of the environment variants (Stability). Following Personalization, we compare the RL candidates by their average sum of rewards, by the 25th percentile reward across all users and trials, and by their across time average reward. For each trial, we redraw 72 simulated users (approximately the expected sample size for the Oralytics study) with replacement in groups of 4. To simulate incremental recruitment of 4 users per week, the order in which we randomly select these groups of 4 is their simulated entry date. Every week in a trial, we check whether a group has completed the 10-week study duration and whether we should add another group of 4 users into the study. We cluster users by the day they enter the study (e.g. for cluster size $4$ the first four users are in the first cluster and the next four users are in the second cluster, etc). We have one RL algorithm instantiation per cluster (no data shared across clusters) and the RL algorithm has a weekly update cadence with the first update starting after one week (at time $t = 14$).

We then run 50 trials for each environmental variant and algorithm candidate pairing. The results are shown in Table~\ref{table/experiment_rewards}. Notice that the average rewards in Table~\ref{table/experiment_rewards} are lower than the 120 seconds of dentist-recommended brushing duration. This is because of the zero-inflated nature of our setting (i.e. the user does not brush). We see in the table that BLR and ZIP with cluster size $k = N$ perform better than other RL algorithm candidates in all environments for average reward and lower 25th percentile reward. This indicates that even though the users are heterogeneous, the heterogeneity is not so large that pooling users' data is ever detrimental.

Next, consider Figure \ref{figs/rewards_across_time}. Again, BLR and ZIP with cluster size $k = N$ perform better than other RL algorithm candidates in all environments across all time steps. We expected algorithms' average reward to increase and then decrease for non-stationary environments. This is because the RL algorithms do not incorporate non-stationarity in the reward approximating function and the modeled reward functions become increasingly biased with large time-steps $t0$. However, it seems that the effect size we considered was large enough to overcome this non-stationarity. 

From Table~\ref{table/experiment_rewards} and Figure~\ref{figs/rewards_across_time}, BLR with cluster sizes $k = 4, N$ is comparable to ZIP. Our hypothesis was that ZIP would perform better than BLR because ZIP has more complexity than BLR. However, there is a bias-variance tradeoff between ZIP and BLR. Although ZIP has more complexity than BLR, ZIP needs to fit more parameters and therefore requires more data to learn effectively; but we are in a data-sparse setting. In addition, ZIP is constrained by the performance of the approximate posterior sampling scheme, which may be unstable and difficult to debug. On the other hand, BLR has a stable closed-form posterior update and BLR with action centering specifically does not require the knowledge of the baseline features at decision time (See Appendix~\ref{update:blr}). This means that baseline features only need to be available at update time and we can incorporate more features that were not available in real-time at decision time.
Also, BLR with action centering is very robust, namely, it is guaranteed to be unbiased even when the baseline reward model is incorrect \cite{DBLP:journals/corr/abs-1909-03539}.

We highlight the following takeaways from our experiments:
\begin{enumerate}
    \item \textbf{BLR vs ZIP:} We prefer BLR over ZIP. BLR with pooling performs similarly or better than ZIP and is more reliable/computable. 
    \item \textbf{Cluster Size:} Large cluster sizes perform better especially when there are population treatment effects. They also perform well not only for the average user but also for users who do not benefit as much. This performance could be due to faster learning because these algorithms can leverage other users' data to make decisions. Even under heterogeneous treatment effects, the benefit of reducing variance by learning across users is not outweighed by user heterogeneity.
\end{enumerate}

There are other limitations to these experiments. 

\underline{Fixed Noise Variance for BLR:} The experiments were run using a fixed noise variance term fit using the ROBAS 2 data set. Assuming a fixed noise variance is unrealistic; we plan to learn the noise variance along with other parameters in the real study. The fixed noise variance term could be a reason that BLR performed comparably to ZIP. 

\underline{More Distinct and Complex Simulation Environments:}
We may not be looking widely enough across variants to find settings where these algorithms perform differently. With sufficient data per user in a highly heterogeneous user environment, then we expect cluster size $k=1$ to do the best. In future work, we aim to add in simulation environments with greater heterogeneity and less noise to see if large cluster sizes still perform well, and we aim to create more complex simulation environment variants that are more distinct (e.g. environments where users may differ by heterogeneous demographic features like age and gender). Also, we want to impute state features of interest in the real study that were not present in the data set, such as phone engagement. 

\underline{Additional RL Algorithm Candidate Considerations.}
We also aim to consider other axes for algorithm candidates such as algorithms with other update cadences (e.g. every night or biweekly) and algorithms with an informative prior. In initial simulations using algorithms with informative priors, we found that since the same (limited amount) of ROBAS 2 data was used to build both the simulation environment and the prior, the algorithms did not need to learn much to perform well. An open question is how to develop both simulation environments and informative priors in a realistic way using a very limited amount of data. Finally, we will also explore additional design decisions such as how to carefully design the baseline and advantage feature spaces for the RL algorithm.

These investigations will determine the final algorithm that will go into the actual study.
\section{Discussion and Future Work}
\label{sec:discussion}
In this paper, we present the first extension of the PCS framework for designing RL algorithms in digital intervention settings. The case study demonstrates how to use the PCS framework to make design decisions and only highlights of our ongoing work in designing the Oralytics RL algorithm. Our illustration helps fellow researchers balance and understand the benefits and drawbacks of certain aspects of their RL algorithm for their digital intervention study.   

\section{Simulation Environments}
\label{appendix:sim_envs}

\subsection{Baseline Feature Space of the Environment Base Models}
\label{baseline_features}
The ROBAS 2 dataset has a variety of features that we anticipate to be associated with brushing duration. These include the time of day (morning vs evening), weekday vs weekend, and summaries of the user's past brushing behavior. Together with domain experts in behavioral health and dentistry, we chose the following features to use to fit a model of the reward. Recall that the ROBAS 2 dataset only includes data under no brushing, so for now we are only fitting a model for the baseline reward model (i.e., the brushing duration under action $A_{i,t}$). In Appendix~\ref{effect_sizes} we discuss how to model brushing duration under action $1$.

\begin{enumerate}
    \item \textbf{Bias / Intercept Term} $\in \mathbb{R}$
    \item \textbf{Time of Day (Morning/Evening) $\in \{0, 1\}$}
    \item \textbf{Prior Day Total Brushing Duration (Normalized)} $\in \mathbb{R}$ 
    \item \textbf{Weekend Indicator (Weekday/Weekend) $\in \{0, 1\}$}
    \item \textbf{Proportion of Non-zero Brushing Sessions Over Past 7 Days} $\in [0, 1]$
    \item \textbf{Day in Study (Normalized) $\in [-1,1]$}
\end{enumerate}
We use these features to generate two types of base reward environments (Stationary and Non-Stationary). The \textbf{Stationary model of the base environment} uses the state function $g(S_{i,t}) \in \mathbb{R}^5$ that \textit{only includes the first five features above}. The \textbf{Non-Stationary model of the base environment} uses state $g(S_{i,t}) \in \mathbb{R}^6$ that corresponds to all of the above features.


\paragraph{\bf{Normalization of State Features}}
\label{normalizations}
We normalize features to ensure that state features are all in a similar range. The Prior Day Total Brushing Duration feature is normalized using z-score normalization (subtract mean and divide by standard deviation) and the Day in Study feature (originally in the range $[1:28]$ since the study length of ROBAS 2 is 28) is normalized to be between $[-1, 1]$. 
Note that when generating rewards, Day in Study was normalized based on Oralytic's anticipated 10 week study length (range is still $[-1,1]$).
$$\text{Normalized Total Brushing Duration in Seconds} = (\text{Brushing Duration} - 172) / 118 $$
$$\text{Normalized Day in Study When Fitting Model} = (\text{Day} - 14.5) / 13.5$$
$$\text{Normalized Day in Study When Generating Rewards} = (\text{Day} - 35.5) / 34.5$$



\subsection{Environment Base Model}
\label{env_base_model}
We consider three model classes: 1) a zero-inflated Poisson, 2) hurdle model with a square root transform, and 3) hurdle model with a log transform, and choose one out of these three model classes for each user. We define these three model classes below. Additionally, below $g(S)$ is the baseline feature vector of the current state defined in Appendix~\ref{baseline_features}, $w_{i,b}, w_{i,p}, w_{i,\mu}$ are user-specific weight vectors, $\sigma_{i,u}^2$ is the user specific variance for the normal component, and $\mathrm{sigmoid}(x) = \frac{1}{1 + e^{-x}}$ is the sigmoid function.

\paragraph{\bf{1) Zero-Inflated Poisson Model for Brushing Duration}}
\label{0_infl_poisson}
$$
Z \sim \text{Bernoulli}\left(1 - \mathrm{sigmoid}(g(S)^T w_{i,b}) \right)
$$
$$
Y \sim \text{Poisson} \left( \exp \left( g(S)^T w_{i,p} \right) \right)
$$
$$
\text{Brushing Duration in Seconds} : D = ZY
$$

\paragraph{\bf{2) Hurdle Model with Square Root Transform for Brushing Duration}}
$$
Z \sim \text{Bernoulli} \left(1 - \mathrm{sigmoid} \left( g(S)^T w_{i,b} \right) \right)
$$
$$
Y \sim \mathcal{N} \left( g(S)^T w_{i,\mu}, \sigma_{i,u}^2 \right)
$$
$$
\text{Brushing Duration} : D = ZY^2
$$
Note that the non-zero component of this model, $Y^2$, can also be represented as a constant times a non-central chi-squared, where the non-centrality parameter is the square of the mean of the normal distribution. 


\paragraph{\bf{3) Hurdle Model with Log Transform for Brushing Duration}}
$$
Z \sim \text{Bernoulli} \left(1 - \sigma \left( g(S)^T w_{i,b} \right) \right)
$$
$$
Y \sim \text{Lognormal} \left( g(S)^T w_{i,\mu}, \sigma_{i,u}^2 \right)
$$
$$
\text{Brushing Duration} : D = ZY
$$

Since we want to simulate brushing duration in seconds, we also round outputs of the hurdle models to the nearest whole integer. Notice that the zero-inflated Poisson model is a mixture model with a latent state. The Bernoulli draw $Z$ is latent and represents the user's intention to brush, and the Poisson models the user's brushing duration when they intend to brush (this is because the brush time can still be zero when the user intends to brush). On the other hand, the hurdle model provides a model for brushing duration conditional on whether the user brushed or not. The Bernoulli draw $Z$ in the hurdle model is observed.

Note that the Hurdle model is used for the simulation environment only and \textit{not} the RL algorithm. The Hurdle model conditions on a collider (e.g. whether the person brushes their teeth), thus potentially leading to causal bias \cite{cole2010illustrating,luque2019educational}. For example, consider an unobserved cause $U$, intervention $A$, whether the user brushed or not $Z$, brushing duration $D$, and a directed acyclic graph with $A \rightarrow Z, U \rightarrow Z, U \rightarrow D$, and $Z \rightarrow D$. Then conditioning on collider $Z$ of treatment opens a pathway from $A$ to $D$ through $U$ \cite{elwert2014endogenous}. Suppose in reality $A$ only impacts whether the user brushes their teeth but not the duration. Then if we condition on $Z$ to evaluate the impact of $A$ on $D$, we may erroneously learn that $A$ impacts the duration of brushing. This makes the Hurdle unsuitable as a model for an RL algorithm that aims to learn causal effects. 

\subsection{Fitting the Environment Base Models} 
\label{fitting_env_base_models}
We use ROBAS 2 data to fit the brushing duration model under action $0$ (no message). For all model classes, we fit one model per user. All models were fit using MAP with a prior $w_{i,b}, w_{i,p}, w_{i,\mu} \sim \mathcal{N}(0, I)$ as a form of regularization because we have sparse data for each user. Weights were chosen by running random restarts and selecting the weights with the highest log posterior density.

\underline{Fitting Hurdle Models:} For fitting hurdle models for user $i$, we fit the Bernoulli component and the non-zero brushing duration component separately. We use $D_{i,t}$ to denote the $i$th ROBAS 2 user's brushing duration in seconds at the $t$ time point. 
Set $Z_{i,t} = 1$ if the original observation $D_{i, t} > 0$ and $0$ otherwise.  We fit a model for this Bernoulli component.
We then fit a model for the normal component to either the square root transform $Y_{i, t} = \sqrt{D_{i, t}}$ or to inverse-log transform $Y_{i, t} = \exp(D_{i, t})$ of the $i$th user's non-zero brushing duration. 

\underline{Fitting Zero-Inflated Poisson Models:} For the zero-inflated Poisson model, we jointly fit parameters for both the Bernoulli and the Poisson components. Since the brushing durations in the ROBAS 2 data were integer values, we did not have to transform the observation to fit the zero-inflated Poisson model. 

\paragraph{\bf{Selecting the Model Class For Each User}}
To select the model class for user $i$, we fit all three model classes using user $i$'s data from ROBAS 2. We then chose the model class that had the lowest root mean squared error (RMSE). Namely, we choose the model class with the lowest $L_i$, where:

$$
L_i := \sqrt{\sum_{t = 1}^T (D_{i, t} - \hat{\mathbb{E}}[D_{i, t} | S_{i, t}])^2}
$$
Recall that $D_{i,t}$ is the brush time in seconds for user $i$ at decision time $t$.
Definitions of $\hat{\mathbb{E}}[D_{i, t} | S_{i, t}]$ for each model class are specified below in  Table~\ref{table/mean_for_model_classes}.

\begin{table}[h]
    \centering
\begin{tabular}{c|c}
 \hline
 \textbf{Model Class} & $\widehat{\mathbb{E}}[D_{i, t} | S_{i, t}]$ \\
 \hline
    Zero-Inflated Poisson & $\left[1 - \mathrm{sigmoid} \left(g(S_{i, t})^Tw_{i,b} \right) \right] \cdot \exp \left( g(S_{i, t})^Tw_{i,p} \right)$ \\
    Hurdle (Square Root) & $\left[1 - \mathrm{sigmoid}\left( g(S_{i, t})^T w_{i,b} \right) \right] \cdot \left[\sigma_{i,u}^2 + (g(S_{i, t})^T w_{i,\mu})^2\right]$  \\
    Hurdle (Log) & $\left[1 - \mathrm{sigmoid}(g(S_{i, t})^Tw_{i,b}) \right] \cdot \exp \left(g(S_{i, t})^Tw_{i,\mu} + \frac{\sigma_{i,u}^2}{2} \right)$ \\
    \hline
\end{tabular}
\caption{\textbf{Definitions of $\widehat{\mathbb{E}}[D_{i, t} | S_{i, t}]$ for each model class.} $\widehat{\mathbb{E}}[D_{i, t} | S_{i, t}]$ is the mean of user model $i$ fitted using data $\{(S_{i, t}, D_{i, t})\}_{t = 1}^T$. }
\label{table/mean_for_model_classes}
\end{table}

\begin{table}[h]
    \centering
\begin{tabular}{c|c}
 \hline
 \textbf{Model Class} &  $\widehat{\mathbb{E}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]$ \\
 \hline
    Hurdle (Square Root) & $\sigma_{i,u}^2 + (g(S_{i, t})^Tw_{i, \mu})^2$ \\
    Hurdle (Log) & $\exp(g(S_{i, t})^Tw_{i,\mu} + \frac{\sigma_{i,u}^2}{2})$ \\
    Zero-Inflated Poisson & $\frac{\exp(g(S_{i, t})^Tw_{i,p})\exp(\exp(g(S_{i, t})^Tw_{i,p}))}{\exp(\exp(g(S_{i, t})^Tw_{i,p})) - 1}$ \\
    \hline & $\widehat{\text{Var}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]$ \\
    \hline
    Hurdle (Square Root) &  $g(S_{i, t})^Tw_{i,\mu}^4 + 3\sigma_{i,u}^4 + 6\sigma_{i,u}^2(g(S_{i, t})^Tw_{i,\mu})^2 - \widehat{\mathbb{E}}[R_{i, t} | S_{i, t}, R_{i, t} > 0]^2$ \\
    Hurdle (Log) & $(\exp(\sigma_{i,u}^2) - 1) \cdot \exp(2g(S_{i, t})^T w_{i,\mu} + \sigma_{i,u}^2)$ \\
    Zero-Inflated Poisson & $\widehat{\mathbb{E}}[D_{i, t} | S_{i, t}, D_{i, t} > 0] \cdot (1 + \exp(g(S_{i, t})^Tw_{i,p}) - \widehat{\mathbb{E}}[D_{i, t} | S_{i, t}, D_{i, t} > 0])$ \\
    \hline
\end{tabular}
\caption{\textbf{Definitions of $\widehat{\mathbb{E}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]$} and $\widehat{\text{Var}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]$ \textbf{for each model class.} $\widehat{\mathbb{E}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]$ and $\widehat{\text{Var}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]$ is the mean and variance of the non-zero component of user model $i$ fitted using data $\{(S_{i, t}, D_{i, t})\}_{t = 1}^T$.}
\label{table/non_zero_mean_var_for_model_classes}
\end{table}

After the procedure was run, we obtained the following number of model classes for all users in the ROBAS 2 study in Table~\ref{table/num_model_class}.

\begin{table*}[!h]
    \centering
\begin{tabular}{c|c|c}
 \hline
 \textbf{Model Class} & \textbf{Stationary} & \textbf{Non-Stationary} \\
 \hline
    Hurdle with Square Root Transform & 9 & 7\\
    Hurdle with Log Transform & 9 & 8 \\
    Zero-Inflated Poisson & 14 & 17
\end{tabular}
\caption{\textbf{Number of selected model classes for the Stationary and Non-Stationary environments.}}
\label{table/num_model_class}
\end{table*}

\subsection{Checking the Quality of the Simulation Environment Base Model}
\label{sim_env_evaluation}

\paragraph{\textbf{Checking Moments}}
Using the chosen user-specific models, we simulate $100$ ``trials''. In each trial, for each of users in ROBAS 2, we use their respective model to generate a data trajectory $(S_{i,t}, R_{i,t})_{t=1}^{56}$ (note that the ROBAS 2 study had $2$ brushing windows per day for $28$ days for a total of $56$ brushing windows). We then compute the following metrics for each of the trials and averaged across trials:
\begin{enumerate}
    \item \textbf{Proportion of Missed Brushing Windows:} $$\frac{1}{N} \sum_{i=1}^N \frac{1}{T} \sum_{t=1}^T \mathbb{I}[D_{i,t} = 0]$$
    \item \textbf{Average Non-Zero Brushing Duration:} $$\frac{1}{N} \sum_{i=1}^N \frac{1}{\sum_{t=1}^T \mathbb{I}[D_{i,t} > 0]} \sum_{t=1}^T \mathbb{I}[D_{i,t} > 0] D_{i,t}$$
    \item \textbf{Variance of Non-Zero Brushing Durations:} Let $\widehat{\mathrm{Var}}(\{X_k \}_{k=1}^K)$ represent the empirical variance of $X_1, X_2, ..., X_K$. 
    $$\widehat{\mathrm{Var}} \big( \big\{ D_{i,t} : t \in [1 \colon T], D_{i,t} > 0 \big\}_{i=1}^N \big)$$
   
    \item \textbf{Variance of Average User Brushing Durations:} This metric measures the degree of between user variance in average brushing.
    $$\widehat{\mathrm{Var}} \bigg( \bigg\{ \frac{1}{T} \sum_{t=1}^T D_{i,t} \bigg\}_{i=1}^N \bigg)$$
    \item \textbf{Average of Variances of Within User Brushing Durations:} This metric measures the average amount of within user variance.
    $$\frac{1}{N} \sum_{i=1}^N \widehat{\mathrm{Var}} \big( \{ D_{i,t} \}_{t=1}^T \big)$$
\end{enumerate}

The base models slightly overestimate the proportion of missed brushing windows in the ROBAS 2 data set. Our base models also slightly underestimate the average brushing duration. Our base models also for the most part slightly overestimate the between-user and within-user variance of rewards.

\begin{table*}[h]
    \begin{tabular}{l|r|r|r}
    \toprule
            \textbf{Metrics} & \textbf{ROBAS 2} & \textbf{Stationary} & \textbf{Non-Stationary} \\
    \midrule
    Proportion of Missed Brushing Windows &    0.376674 &    0.403114 &    0.397812 \\
    Average Non-Zero BDs &  137.768129 &  131.308445 &  134.676955 \\
    Variance of Non-Zero BDs & 2326.518304 & 2392.955018 & 2253.177853 \\
    Variance of Average User BDs & 1415.920148 & 1699.126897 & 1399.615330 \\
    Average of Variances of Within User BDs & 1160.723506 & 1405.944459 & 1473.239769 \\
    \bottomrule
    \end{tabular}
\caption{\textbf{Comparing Moments Between Base Models and ROBAS 2 Data set.} Above we use BDs to abbreviate Brushing Durations. Values for the Stationary and Non-Stationary base models are averaged across 100 trials.}
\label{table/moments}
\end{table*}

\paragraph{\textbf{Measuring If A Base Model Captures the Variance in the Data}}
We measure how well the fitted base models captured (1) whether or not the user brushed, and (2) the variance of the brush time when the users did brush.
To measure point (1), for each user model $i$, we calculate the statistic:
\begin{equation}
\label{eqn:point1}
   U_i := \frac{1}{T} \sum_{t = 1}^T \bigg(\frac{\mathbb{I}[D_{i,t} > 0] - \widehat{\mathbb{E}}[\mathbb{I}[D_{i,t} > 0] | S_{i, t}]}{\widehat{\text{Var}}[\mathbb{I}[D_{i,t} > 0] | S_{i, t}]} \bigg)^2 
\end{equation}
where $\widehat{\mathbb{E}}[\mathbb{I}[D_{i,t} > 0] | S_{i, t}] = 1 - \mathrm{sigmoid}(S_{i, t}^Tw_{i,b})$ and $\widehat{\text{Var}}[\mathbb{I}[D_{i,t} > 0] | S_{i, t}] = \widehat{\mathbb{E}}[\mathbb{I}[D_{i,t} > 0] | S_{i, t}] \cdot \mathrm{sigmoid}(S_{i, t}^Tw_{i,b})$.

To measure point (2), for each user model $i$, we calculate the statistic:
\begin{equation}
\label{eqn:point2}
    U_i := \frac{1}{\sum_{t=1}^T \mathbb{I}[D_{i,t} > 0] } \sum_{t = 1}^T \mathbb{I}[D_{i,t} > 0] \bigg(\frac{D_{i, t}  - \widehat{\mathbb{E}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]}{\widehat{\text{Var}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]} \bigg)^2
\end{equation}

Definitions of $\widehat{\mathbb{E}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]$ and $\widehat{\text{Var}}[D_{i, t} | S_{i, t}, D_{i, t} > 0]$ for the non-zero component of each model class are specified in  Table~\ref{table/non_zero_mean_var_for_model_classes}.
For a user model to capture the variance in the data, $U_i$ should be close to 1. We calculate the empirical mean $\overline{U} = \frac{1}{N} \sum_{i = 1}^N U_i$ and standard deviation $\overline{\sigma_U} = \text{std}(U_i)$, and the approximate 95\% confidence interval is $\overline{U} \pm 1.96 * \frac{\overline{\sigma_U}}{\sqrt{N}}$.

\begin{table*}[h]
    \begin{tabular}{l|r|r}
    \toprule
            \textbf{Metric} & \textbf{Stationary} & \textbf{Non-Stationary} \\
    \midrule
    Eqn~\eqref{eqn:point1} $\overline{U}$ & 0.811 & 0.792 \\
    Eqn~\eqref{eqn:point1} $\overline{\sigma_U}$ & 0.146 & 0.150 \\
    Eqn~\eqref{eqn:point1} Confidence Interval &  (0.760, 0.861) & (0.739, 0.844) \\
    Eqn~\eqref{eqn:point2} $\overline{U}$ & 3.579 & 3.493 \\
    Eqn~\eqref{eqn:point2} $\overline{\sigma_U}$ & 4.861 & 4.876 \\
    Eqn~\eqref{eqn:point2} Confidence Interval &  (1.895, 5.263) & (1.803, 5.182) \\
    \midrule
    \bottomrule
    \end{tabular}
\caption{\textbf{Statistic $U_i$ for Capturing Variance in the Data.} Values are rounded to the nearest 3 decimal places.}
\label{table/capturing_variance}
\end{table*}

Results are in Table~\ref{table/capturing_variance}. We can see that after computing the statistic for each user, the confidence interval is close to 1. 
We understand that the confidence intervals do not contain 1 which implies that we are overestimating the amount of variance in $\mathbb{I}[D > 0]$, and underestimating the amount of variance in $D | D > 0$. In the future, we hope to improve upon this statistic by considering non-linear components in our base models.

\subsection{Imputing Treatment Effect Sizes for Simulation Environments}
\label{effect_sizes}

Recall the ROBAS 2 dataset does not have any data under action $1$ (send a message). Thus to model the reward under action $1$ we must impute. 

\paragraph{\bf{Treatment Effect Feature Space}} 
The following choice of the treatment effect (advantage) feature space was made after discussion with domain experts on which features are most likely to interact with the intervention (action).
The Stationary model uses state $h(S) \in \mathbb{R}^4$ corresponding to the following features:

\begin{enumerate}
    \item \textbf{Bias / Intercept Term $\in \mathbb{R}$}
    \item \textbf{Time of Day (Morning/Evening) $\in \{0, 1\}$}
    \item \textbf{Prior Day Total Brushing Duration (Normalized)} $\in \mathbb{R}$ 
    \item \textbf{Weekend Indicator (Weekday/Weekend) $\in \{0, 1\}$}
\end{enumerate}
The Non-Stationary model uses state $h(S) \in \mathbb{R}^5$ corresponding to all of the above features as well as the following:
\begin{enumerate}[resume]
    \item \textbf{Day in Study (Normalized) $\in \mathbb{R}$} 
\end{enumerate}

\paragraph{\bf{Imputation Approach}}

For the Zero-Inflated Poisson model, we impute treatment effects on both the user's intent to brush (Bernoulli component) and the user's brushing duration when they intend to brush (Poisson component). Similarly, for the Hurdle models, we impute treatment effects on both whether the user's brushing duration is zero (Bernoulli component) and the user's brushing duration when the duration is nonzero.

After incorporating effect sizes, brushing duration under action $A$ in state $S$ is $D$ where:
$$
Z \sim \text{Bernoulli}(1 - \mathrm{sigmoid}(g(S)^Tw_{i,b} + A * h(S)^T\Delta_{i,B}))
$$
$$D =
\begin{cases}
ZY^2, Y \sim \mathcal{N}(g(x)^Tw_{i,\mu} + A * h(x)^T\Delta_{i,N}, \sigma_u^2) & \text{for Hurdle square root model} \\
Z\exp(Y), Y \sim \mathcal{N}(g(x)^Tw_{i,\mu} + A * h(x)^T\Delta_{i,N}, \sigma_{i,u}^2) & \text{for Hurdle log normal model} \\
ZY, Y \sim \text{Poisson}(\exp(g(x)^Tw_{i,p} + A * h(x)^T\Delta_{i,N})) & \text{for Zero-Inflated Poisson model}
\end{cases}
$$
$\Delta_{i,B}, \Delta_{i,N}$ are user-specific effect sizes; we will also consider population-level effect sizes (same across all users) which we denote as $\Delta_{B}, \Delta_{N}$. $g(S)$ is the baseline feature vector as described in Appendix~\ref{baseline_features}, and $h(S)$ is the feature vector that interacts with the effect size as specified above.

\textbf{\paragraph{Heterogeneous versus Population-Level Effect Size}}
We consider realistic heterogeneous effect sizes (each user has a unique effect size) and a realistic population-level effect size (all users who share the same base model class also share the same effect size). 
\\

\underline{Population Level Effect Sizes:}
Recall that for the stationary base model we fit models for $Y$ and $Z$, and get user-specific parameters $w_{i,b}, w_{i,p} \in \mathbb{R}^5$ (Zero-Inflated Poisson model) and parameters $w_{i,b}, w_{i,\mu} \in \mathbb{R}^5$ (Hurdle models). We use these parameters to form the population effect sizes as follows: \\

\noindent Zero-Inflated Models Effect Sizes:
\begin{itemize}
    \item $\Delta_B = \mu_{B,\text{avg}}$ where $\mu_{B,\text{avg}} = \frac{1}{4} \sum_{d \in [2 \colon 5]} \frac{1}{N} \sum_{i=1}^N |w_{i,b}^{(d)}|$.
    \item $\Delta_N = \mu_{N,\text{avg}}$ where $\mu_{N,\text{avg}} = \frac{1}{4} \sum_{d \in [2 \colon 5]} \frac{1}{N} \sum_{i=1}^N |w_{i,p}^{(d)}|$.
\end{itemize}

\noindent Hurdle Models Effect Sizes:
\begin{itemize}
    \item $\Delta_B = \mu_{B,\text{avg}}$ where $\mu_{B,\text{avg}} = \frac{1}{4} \sum_{d \in [2 \colon 5]} \frac{1}{N} \sum_{i=1}^N |w_{i,b}^{(d)}|$.
    \item $\Delta_N = \mu_{N,\text{avg}}$ where $\mu_{N,\text{avg}} = \frac{1}{4} \sum_{d \in [2 \colon 5]} \frac{1}{N} \sum_{i=1}^N |w_{i,\mu}^{(d)}|$.
\end{itemize}

\underline{Heterogeneous Effect Sizes:}
To calculate the heterogeneous effect sizes, we again group users by their chosen base model (Zero-Inflated, Hurdle Square-Root, Hurdle Log). We then draw effect sizes for each user from a normal distribution specific to their base model:
$$
\Delta_{i,B} \sim \mathcal{N}(\mu_{B}, \sigma^2_{B})
$$
$$
\Delta_{i,B} \sim \mathcal{N}(\mu_{N}, \sigma^2_{N})
$$
$\mu_{B}, \mu_{N}$ are set to the population level effect sizes for that base model class as described above. To set $\sigma^2_{B}, \sigma^2_{N}$, we do the following: \\

\noindent Zero-Inflated Models
\begin{itemize}
    \item $\sigma_B$ is the empirical standard deviation over $\{ \mu_{i,B} \}_{i=1}^{N}$ where $\mu_{i,B} = \frac{1}{4} \sum_{d \in [2 \colon 5]} |w_{i,b}^{(d)}|$. We use $w_{i,b}^{(d)}$ to denote the $d^{\mathrm{th}}$ dimension of the vector $w_{i,b}$; we take the minimum over all dimensions excluding $d=1$ which represents the weight for the bias/intercept term.
    \item $\sigma_N$ is the empirical standard deviation over $\{ \mu_{i,N} \}_{i=1}^{N}$ where $\mu_{i,N} = \frac{1}{4} \sum_{d \in [2 \colon 5]} |w_{i,p}^{(d)}|$.
\end{itemize}

\noindent Hurdle Models:
\begin{itemize}
    \item $\sigma_B$ is the empirical standard deviation over $\{ \mu_{i,B} \}_{i=1}^{N}$ where $\mu_{i,B} = \frac{1}{4} \sum_{d \in [2 \colon 5]} |w_{i,b}^{(d)}|$.
    \item $\sigma_N$ is the empirical standard deviation over $\{ \mu_{i,N} \}_{i=1}^{N}$ where $\mu_{i,N} = \frac{1}{4} \sum_{d \in [2 \colon 5]} |w_{i,\mu}^{(d)}|$.
\end{itemize}
Histograms of $\Delta_{i,B}, \Delta_{i,N}$ and values of $\mu_{B}, \mu_{N}$ for each base model class are specified in Figure~\ref{figs/effect_sizes_hist}. 

\begin{figure*}[!h]
    \centering
    \subfloat(a){\includegraphics[width=0.3\textwidth]{figures/effect_size_histograms/large_eff_size_hist_bern_hurdle.pdf}}
    \subfloat(b){\includegraphics[width=0.3\textwidth]{figures/effect_size_histograms/large_eff_size_hist_sqrt_hurdle.pdf}}
    \subfloat(c){\includegraphics[width=0.3\textwidth]{figures/effect_size_histograms/large_eff_size_hist_log_hurdle.pdf}}
    \subfloat(d){\includegraphics[width=0.3\textwidth]{figures/effect_size_histograms/large_eff_size_hist_bern_zip.pdf}}
    \subfloat(e){\includegraphics[width=0.3\textwidth]{figures/effect_size_histograms/large_eff_size_hist_poisson_zip.pdf}}
    \caption{\textbf{$\Delta_{i,B}$'s, $\Delta_{i,N}$'s, $\mu_{B}, \mu_{N}$ for each base model class.} For  \textbf{(a)} Effect Sizes for Bernoulli Component of the Hurdle Models, \textbf{(b)} Effect Sizes for Non-Zero Component of Hurdle (Square Root) Model,  \textbf{(c)} Effect Sizes for Non-Zero Component of Hurdle (Log) Model, \textbf{(d)} Effect Sizes for Bernoulli Component of the ZIP Model, \textbf{(e)} Effect Sizes for Poisson Component of the ZIP Model 
    }
    \label{figs/effect_sizes_hist}
\end{figure*}





\section{Reinforcement Learning Algorithm Candidates}
\label{appendix:RLalgs}

\subsection{Feature Space for the RL Algorithm Candidates}
\label{app:rl_alg_features}
\paragraph{\bf{Advantage Feature Space}}
We use $f(S) \in \mathbb{R}^3$ to denote the feature space used by our RL algorithm candidates to predict the advantage, i.e., the immediate treatment effect, is the following:
\begin{enumerate}
    \item \bf{Bias / Intercept Term} $\in \mathbb{R}$
    \item \bf{Time of Day (Morning/Evening)} $\in \{0, 1\}$
    \item \bf{Prior Day Total Brushing Duration (Normalized)} $\in \mathbb{R}$ 
\end{enumerate}
The normalization procedure for Prior Day Brushing Duration is the same as described in Appendix~\ref{normalizations}.

\paragraph{\bf{Baseline Feature Space}}
We use $m(S) \in \mathbb{R}^4$ to denote the feature space used by the RL algorithm candidates to approximate the baseline reward function containing all the above covariates as well as the following:

\begin{enumerate}[resume]
    \item \bf{Weekend Indicator (Weekday/Weekend)} $\in \{0, 1\}$
\end{enumerate}

Note, {\it the feature space used by the RL algorithm candidates is different than the feature space used to model the reward in the simulation environments, specified in Appendix~\ref{baseline_features}} and Appendix~\ref{effect_sizes}; this means that the RL algorithms will have a misspecified reward model. Namely, the baseline feature space for the simulation environment has an additional \textbf{Proportion of Non-zero Brushing Sessions Over Past 7 Days} feature and the non-stationary variant has the \textbf{Day in Study (Normalized)} feature. The treatment effect feature space for the simulation environment has an additional \textbf{Weekend Indicator (Weekday/Weekend)} and the non-stationary variant has the \textbf{Day in Study (Normalized)} feature. 

The rationale for not including the \textbf{Day in Study (Normalized)} feature is although we wanted to capture potential non-stationarity in brushing outcomes in order to create a realistic simulation environment, our RL algorithm candidates do not have reward functions that vary arbitrarily over time. 
We do not include \textbf{Proportion of Non-zero Brushing Sessions Over Past 7 Days} and \textbf{Weekend Indicator (Weekday/Weekend)} to detect the robustness of RL algorithm candidates to a misspecified reward model.

\subsection{Decision 1: Reward Approximating Function}
\label{reward_approx_functions}

The first decision in designing the RL algorithm is the choice between using a linear model or a 0-Inflated Poisson model as the reward approximating function used by the posterior sampler (note this is separate from the reward model used to generate the environment). More information on how the posterior sampling algorithm performs action selection can be found in Appendix~\ref{action_selection}, and how the algorithm updates at update times can be found in Appendix~\ref{inference_methods}. We describe the two candidates below.

Note that the function $m$ for the RL algorithm's baseline reward model is only used at update times. The function $f$ for the RL algorithm's advantage model is used at both decision and update times.

\subsubsection{Bayesian Linear Regression Model}
The first candidate is to use the following reward generating model with action centering used in \cite{DBLP:journals/corr/abs-1909-03539} for the posterior sampler:
\begin{equation}
\label{eqn:blr}
    R_{i, t} = m(S_{i, t})^T \alpha_{i,0} + \pi_{i,t} f(S_{i, t})^T \alpha_{i,1} + (A_{i, t} - \pi_{i, t}) f(S_{i, t})^T \beta_i + \epsilon_{i,t}
\end{equation}
where $\alpha_{i,0} \in \mathbb{R}^4$ and $\alpha_{i,1}, \beta_i \in \mathbb{R}^3$. $\pi_{i,t}$ is the probability that action $A_{i,t}=1$ is selected by the RL algorithm for user $i$ in state $S_{i,t}$; we discuss how to compute these in Appendix~\ref{action_selection}. 
The RL algorithm models $\epsilon_{i,t}$ as being drawn from $\mathcal{N}(0, \eta^2)$ (the choice of $\eta^2$ is informed by the ROBAS 2 dataset). 
Additionally, we put uninformative normal priors on the parameters:
$\alpha_{i,0} \sim \mathcal{N}(0, \sigma_{prior}I_4)$, $\alpha_{i,1} \sim \mathcal{N}(0, \sigma_{prior}I_3)$, $\beta_i \sim \mathcal{N}(0, \sigma_{prior}I_3)$, where $\sigma_{prior} = 5$.

\subsubsection{Zero-Inflated Poisson Regression Model}
The second candidate is to use the zero-inflated reward generating model for the posterior sampler:
\begin{equation}
\label{eqn:zip_alg}
\begin{split}
Z_{i, t} &\sim \text{Bernoulli} \left(1 - \mathrm{sigmoid}(m(S_{i, t})^T \alpha_{i,p} + A_{i, t} \cdot f(S_{i, t})^T \beta_{i,b}) \right) \\
Y_{i, t} &\sim \text{Poisson} \left( \exp \left(m(S_{i, t})^T \alpha_{i,p} + A_{i, t} \cdot f(S_{i, t})^T \beta_{i,p} \right) \right) \\
R_{i, t} &= Z_{i, t}Y_{i, t}
\end{split}
\end{equation}
The above model closely resembles the zero-inflated Poisson model class used to develop the simulation environment in Appendix~\ref{env_base_model}; however, recall that the feature space used by the RL algorithm and the model to generate the environment are different. Additionally, here we model the reward directly, rather than the raw brushing duration.

Additionally, the posterior sampling algorithm will put the following uninformative normal priors on the parameters: 
$\alpha_{i,b}, \alpha_{i,p} \sim \mathcal{N}(0, \sigma_{prior}I_4)$ and $\beta_{i,b}, \beta_{i,p} \sim \mathcal{N}(0, \sigma_{prior}I_3)$, where $\sigma_{prior} = 5$.

\subsection{Decision 2: Cluster Size}
\label{app:cluster_sizes}
Clustering involves grouping $k$ users together and pooling all of their data together for the RL algorithm. This means that we have one RL algorithm instantiation per cluster (no data shared across clusters). For our experiments, we draw 72 simulated users (the expected sample size for the Oralytics study) with replacement and cluster these users at random (every possible cluster is equally likely). We then keep these cluster assignments fixed across the trials.

For simplicity in running our experiments, we consider randomly formed clusters, but we are thinking of clustering by entry date in the real study. We cannot predict how many users who share the same baseline feature will join within a relatively short period of time (e.g. we cannot depend on there being 4 females within the first two weeks). Entry date is reasonable because we are guaranteed to form a cluster for those users and the domain experts believe that users who enter the study around the same time will be similar. Users who enter near the end of the study may be very different from users who enter near the beginning because of societal factors (e.g. pandemic restrictions being lifted), seasonal influences (e.g. changes in user's mood in spring versus mid-winter), and fidelity (e.g. quality of on-boarding procedures and staff experience may improve over time). One natural approach is to cluster by baseline features, however, that is not feasible for a study where the recruitment rate is slow such as in Oralytics.
\section{RL Algorithm Posterior Updates and Posterior Sampling Action Selection}
\label{app:algUpdatesActions}

\subsection{Posterior Updates to the RL Algorithm at Update Time}
\label{inference_methods}
During the update step, the reward approximating function will update the posterior with newly collected data. Also, we make $M$ draws of the parameters from the updated posterior and use them for all decision times until the next update time. Here are the procedures for how the posterior updates for the Bayesian Linear Regression model and the 0-Inflated Poisson model. 
\subsubsection{Bayesian Linear Regression Model}
\label{update:blr}
Suppose we are selecting actions for decision time $t$.
Let $\phi(S_{i, t}, A_{i, t}) = [m(S_{i, t})$, $\pi_{i, t}f(S_{i, t})$, $(A_{i, t} - \pi_{i, t})f(S_{i, t})]$ be the joint feature vector and $\theta_i = [\alpha_{0,i}, \alpha_{i,1}, \beta_i]$ be the joint weight vector. Notice that Equation~\ref{eqn:blr} can be vectorized of the form: $R_{i, t} = \phi(S_{i, t}, A_{i, t})^T\theta_i + \epsilon_{i,t}$. 
Now let $\Phi_{i,1:t-1}$ be the matrix of all stacked vectors $\{ \phi(S_{i, s}, A_{i, s}) \}_{s=1}^{t-1}$, and $\mathbf{R}_{i,1:t-1}$ be a vector of stacked rewards $\{ R_{i, s} \}_{s=1}^{t-1}$, where we have batch data of the $t-1$ decision times before the current update time. 

Recall we have Normal priors on $\theta_i$ where $\theta_i \sim \mathcal{N}(\mu_{\text{prior}}, \Sigma_{\text{prior}})$, where $\mu_{\text{prior}} = \boldsymbol{0} \in \mathbb{R}^{3+3+4}$ and $\Sigma_{\text{prior}} = \text{diag}(\sigma_{prior}^2 I_3, \sigma_{prior}^2 I_3, \sigma_{prior}^2 I_4)$. The posterior distribution of the weights given current history $H_{i,t-1}$, $p(\theta_i | H_{i,t-1})$ is conjugate and is also normal. 

$$
\theta_i | H_{i,t-1} \sim \mathcal{N}(\mu_{i,t-1}^{posterior}, \Sigma_{i,t-1}^{posterior}) 
$$
$$\Sigma_{i,t-1}^{posterior} = \bigg(\frac{1}{\eta^2}\Phi_{i,1:t-1}^T \Phi_{i,1:t-1} + \Sigma_{prior}^{-1}\bigg)^{-1}
$$
$$\mu_{i,t-1}^{posterior} = \Sigma_{i,t-1}^{posterior} \bigg(\frac{1}{\eta^2}\Phi_{i,1:t-1}^T \mathbf{R}_{i,1:t-1} + \Sigma_{prior}^{-1}\mu_{prior} \bigg)$$

Note, we fit $\eta^2$ to the ROBAS 2 dataset and fixed it for all of our experiments. In the real study, we want to consider assigning a conjugate prior on $\eta^2$ and updating it at update times.

\subsubsection{Zero-Inflated Poisson Regression Model}
\label{update:zip}
For the 0-Inflated Poisson regression model, the posterior distribution of the weights $\theta_i = \{\alpha_{i,b}, \beta_{i,b}, \alpha_{i,p}, \beta_{i,p} \}$ given data $H_{i,t-1}$, $p(\theta_i | H_{i,t-1})$, does not have a closed form. Therefore we use Metropolis Hastings (MH) with a normal proposal distribution as an approximate posterior sampling method.

\paragraph{\bf{Posterior Density}}
The log-likelihood of the 0-Inflated Poisson Regression Model is:
$$
\log f(R_{i, t} | S_{i, t}, A_{i, t}; \theta_i) = 
\begin{cases} 
\log ((1 - p) + p \exp(-\lambda)) & R = 0 \\
\log p - \lambda + R \log \lambda - \log R! & R = 1, 2, 3,...
\end{cases}
$$
where $p = 1 - \mathrm{sigmoid}\left( m(S_{i, t})^T \alpha_{i,b} + A_{i, t} \cdot f(S_{i, t})^T \beta_{i,b} \right)$, is the probability of the user intending to brush, and $ \lambda = \exp \left( m(S_{i, t})^T \alpha_{i,p} + A_{i, t} \cdot f(S_{i, t})^T \beta_{i,p} \right)$ is the expected Poisson count.

Therefore the log posterior density is:
$$
\log p(\theta_i | H_{i,t-1}) \propto \sum_{n = 1}^N \log f(R_{i, t} | S_{i, t}, A_{i, t}; \theta_i) + \log p(\theta_i)
$$

\paragraph{\bf{Proposal Distribution}}
We choose a Normal distribution for our proposal distribution. At each step of MH, we propose a new sample given the old sample,  $\theta^k_{\text{prop}} \sim \mathcal{N}(\theta^k_{\text{old}}, \gamma^2I)$, where $\theta^k$ denotes the $k$th value of $\theta$.

\paragraph{\bf{Metropolis Hastings Acceptance Ratio}}
The Metropolis Hastings acceptance ratio given a proposed sample $\theta_{\text{prop}}$ and an old sample $\theta_{\text{old}}$ is defined as:

$$
\alpha(\theta_{\text{prop}}, \theta_{\text{old}}) := \min \bigg(1, \frac{p(\theta_{\text{prop}}) / q(\theta_{\text{prop}} | \theta_{\text{old}})}{p(\theta_{\text{old}}) / q(\theta_{\text{old}} | \theta_{\text{prop}})} \bigg)
$$

Since our proposal distribution is symmetric, the log acceptance ratio becomes:

$$
\log \alpha(\theta_{\text{prop}}, \theta_{\text{old}}) := \min(0, \log p(\theta_{\text{prop}}) - \log p(\theta_{\text{old}}))
$$

\subsection{Action Selection at Decision Time}
\label{action_selection}

\subsubsection{\textbf{Posterior Sampling}}
Our action selection scheme at decision time selects action $1$ according to the posterior probability that that arm is optimal. We define these as $\tilde{\pi}_{i,t}$ below. \\

\noindent \underline{Bayesian Linear Regression Model}
Based on the Bayesian Linear Regression model of the reward, specified by Equation~\eqref{eqn:blr}:
$$
\tilde{\pi}_{i,t} = \text{Pr}_{\tilde{\beta} \sim \mathcal{N}(\mu^{post}_{i,t-1}, \Sigma^{post}_{i,t-1}) } \left\{f(S_{i, t})^T \tilde{\beta} > 0 \big| S_{i, t}, H_{i, t - 1} \right\}.
$$
Note that the randomness in the probability above is only over the draw of $\tilde{\beta}$ from the posterior distribution.
\\\\
\noindent \underline{Zero-Inflated Poisson Model}
Based on the Zero-Inflated Poisson model of the reward, specified by Equation~\eqref{eqn:zip_alg}:

$$
\tilde{\pi}_{i,t} = \text{Pr}_{\tilde{\alpha}_{i,b}, \tilde{\alpha}_{i,p}, \tilde{\beta}_{i,b}, \tilde{\beta}_{i,p} } \left\{ \tilde{Z}_{i, t} \tilde{Y}_{i, t} > 0 \big| S_{i, t}, H_{i, t - 1} \right\}
$$
where $\tilde{Z}_{i, t} \sim \text{Bernoulli} \left(1 - \mathrm{sigmoid}(m(S_{i, t})^T \tilde{\alpha}_{i,p} + A_{i, t} \cdot f(S_{i, t})^T \tilde{\beta}_{i,b}) \right)$
and \\
$\tilde{Y}_{i, t} \sim \text{Poisson} \left( \exp \left(m(S_{i, t})^T \tilde{\alpha}_{i,p} + A_{i, t} \cdot f(S_{i, t})^T \tilde{\beta}_{i,p} \right) \right)$.
Note that the randomness in the probability above is only over the draw of $(\tilde{\alpha}_{i,b}, \tilde{\alpha}_{i,p}, \tilde{\beta}_{i,b}, \tilde{\beta}_{i,p})$ from the posterior distribution. \\

\subsubsection{\textbf{Clipping to form Action Selection Probabilities}}

Since we want to facilitate after study analyses, we clip action selection probabilities using action clipping function for some chosen values of $\pi_{\min}, \pi_{\max}$ where $0 < \pi_{\min} \leq \pi_{\max} < 1$ chosen by the scientific team:
\begin{equation}
\label{eqn:clipping_function}
    \mathrm{clip}(\pi) = \min(\pi_{\max}, \max(\pi, \pi_{\min})) \in [\pi_{\min}, \pi_{\max}]
\end{equation}
This means that
$$\pi_{i,t} = \mathrm{clip} \left( \tilde{\pi}_{i,t} \right)$$


\begin{comment}
\subsection{Action Selection For Bayesian Linear Regression Model} 
Based on the Bayesian Linear Regression model of the reward, specified by Equation~\ref{eqn:blr}:

$$
\text{Pr}_{\tilde{\beta} \sim \mathcal{N}(\mu_{post}_i, \Sigma_{post}_i) } \left\{f(S_{i, t})^T \tilde{\beta}_i > 0 \big| S_{i, t}, H_{i, t - 1} \right\}
$$

Notice that $\beta$ is a sample from the posterior for the estimate of the margin, $\mu_{post}, \Sigma_{post}$ are the parameters of the posterior distribution of $\beta$ updated at update time as described in Appendix~\ref{update:blr}.

$$
\pi_t = \phi \bigg(\frac{\sum_{m = 1}^M \mathbb{I}[f(S_{i, t})^T\beta_m > 0]}{M} \bigg)
$$

where $\phi(\pi)$ is the clipping function specified above in Equation~\ref{eqn:clipping_function}.
\\\\
Finally, we sample action $A_t \sim \text{Bern}(\pi_{i, t})$ 

\subsection{Action Selection For Zero-Inflated Poisson Model}

Based on the Zero-Inflated Poisson model of the reward, specified by Equation~\ref{eqn:zip_alg}, we also look at 

$$
\pi_{i,t} = \phi \bigg(\frac{\sum_{m = 1}^M \mathbb{I}[\mathbb{E}[R_{i, t} | S_{i, t}, A_{i, t} = 1; \theta_m] - \mathbb{E}[R_{i, t} | S_t, A_{i, t} = 0; \theta_m] > 0]}{M} \bigg)
$$
where 
$
\mathbb{E}[R_{i, t} | A_{i, t}, S_{i, t}; \theta_m]
= (1 - \sigma(g(S_{i, t})^Tw_b^0 + A_{i, t} * f(S_{i, t})^Tw_b^1)) * \exp(g(S_{i, t})^Tw_p^0 + A_{i, t} * f(S_{i, t})^Tw_p^1)
$, $\theta_m = [w_b^0, w_b^1, w_p^0, w_p^1]$ represents one draw of $\theta$ from the posterior updated at update time as described in Appendix~\ref{update:zip}, and where $\phi(\pi)$ is the clipping function specified above in Equation~\ref{eqn:clipping_function}.

Namely, at update time, we drew $M$ samples of $\theta$, where $\theta_m$ denotes one draw from the posterior. Then at decision time, we sample action $A_{i,t} \sim \text{Bern}(\pi_{i, t})$

\end{comment}

