## Appendix A Simulation Environments
### Imputing Treatment Effect Sizes for Simulation Environments
#### a.5.2 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-\text{sigmoid}(g(S)^{T}w_{i,b}+A\times h(S)^{T}\Delta _{i,B}))\]  
\[D= \begin{cases}ZY^{2},Y\sim\mathcal{N}(g(x)^{T}w_{i,u}+A\times h(x)^{T}\Delta _{i,N},\sigma_{u}^{2})&\text{for hurdle square root}\\ Z\exp(Y),Y\sim\mathcal{N}(g(x)^{T}w_{i,u}+A\times h(x)^{T}\Delta_{i,N},\sigma_ {i,u}^{2})&\text{for hurdle log normal}\\ ZY,Y\sim\text{Poisson}(\exp(g(x)^{T}w_{i,p}+A\times h(x)^{T}\Delta_{i,N}))& \text{for zero-inflated Poisson}\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 A.1, and \(h(S)\) is the feature vector that interacts with the effect size as specified above. Appendix A.5.3. 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). 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). Values of the fitted parameters can be accessed at: [[https://github.com/StatisticalReinforcementLearningLab/pcs-for-rl/tree/main/sim_env_data](https://github.com/StatisticalReinforcementLearningLab/pcs-for-rl/tree/main/sim_env_data)]([https://github.com/StatisticalReinforcementLearningLab/pcs-for-rl/tree/main/sim_env_data](https://github.com/StatisticalReinforcementLearningLab/pcs-for-rl/tree/main/sim_env_data)) (accessed on 1 June 2022). We use these parameters to form the population effect sizes as follows:  
Zero-Inflated Models' Effect Sizes:  
* \(\Delta_{B}=\mu_{B,\text{avg}}\) where \(\mu_{B,\text{avg}}=\frac{1}{4}\sum_{d\in[2:\ 5]}\frac{1}{N}\sum_{i=1}^{N}|w_{i,b}^{(d)}|\). * \(\Delta_{N}=\mu_{N,\text{avg}}\) where \(\mu_{N,\text{avg}}=\frac{1}{4}\sum_{d\in[2:\ 5]}\frac{1}{N}\sum_{i=1}^{N}|w_{i,p}^{(d)}|\). Hurdle Models' Effect Sizes:
* \(\Delta_{B}=\mu_{B,\text{avg}}\) where \(\mu_{B,\text{avg}}=\frac{1}{4}\sum_{d\in[2:\ 5]}\frac{1}{N}\sum_{i=1}^{N}|w_{i,b}^{(d)}|\). * \(\Delta_{N}=\mu_{N,\text{avg}}\) where \(\mu_{N,\text{avg}}=\frac{1}{4}\sum_{d\in[2:\ 5]}\frac{1}{N}\sum_{i=1}^{N}|w_{i,\mu}^{(d)}|\). We use \(w_{i,b}^{(d)},w_{i,p}^{(d)},w_{i,u}^{(d)}\) to denote the \(d^{\text{th}}\) dimension of the vector \(w_{i,b},w_{i,p},w_{i,u}\) respectively; we take the minimum over all dimensions excluding \(d=1\), which represents the weight for the bias/intercept term. 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}(\Delta_{B},\sigma_{B}^{2})\] \[\Delta_{i,N} \sim\mathcal{N}(\Delta_{N},\sigma_{N}^{2})\]  
\(\Delta_{B},\Delta_{N}\) are set to the population-level effect sizes for that base model class as described above. To set \(\sigma_{B}^{2},\sigma_{N}^{2}\), we do the following:  
Zero-Inflated Models:  
* \(\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:\ 5]}|w_{i,b}^{(d)}|\).