Deep Spectral Q-learning with Application to Mobile Health

Precision medicine focuses on providing personalized treatment to patients by taking their personal information into consideration (see e.g., Kosorok and Laber, 2019; Tsiatis et al., 2019). It has found various applications in numerous studies, ranging from the cardiovascular disease study to cancer treatment and gene therapy (Jameson and Longo, 2015). A dynamic treatment regime (DTR) consists of a sequence of treatment decisions rules tailored to each individual patient’s status at each time, mathematically formulating the idea behind precision medicine. One of the major objectives in precision medicine is to identify the optimal dynamic treatment regime that yields the most favorable outcome on average.

With the rapidly development of mobile health (mHealth) technology, it becomes feasible to collect rich longitudinal data through mobile apps in medical studies. A motivating data example is given by the OhioT1DM dataset (Marling and Bunescu, 2020) which contains data from 12 patients suffering from type-I diabetes measured via fitness bands over 8 weeks. Data-driven decision rules estimated from these data have the potential to improve these patients’ health (see e.g., Shi et al., 2020; Zhu et al., 2020; Zhou et al., 2022). However, it remains challenging to estimate the optimal DTR in these mHealth studies. First, the number of treatment stages (e.g., horizon) is no longer fixed whereas the number of patients can be limited. For instance, in the OhioT1DM dataset, only 12 patients are enrolled in the study. Nonetheless, suppose treatment decisions are made on an hourly basis, the horizon is over 1000. Existing proposals in the DTR literature (Murphy, 2003; Zhang et al., 2013; Zhao et al., 2015; Song et al., 2015; Shi et al., 2018; Nie et al., 2021; Fang et al., 2022; Qi and Liu, 2018; Mo et al., 2021; Ertefaie et al., 2021; Guan et al., 2020; Zhang et al., 2018) become inefficient in these long or infinite horizon settings and require a large number of patients to be consistent. Second, patients’ time-varying covariates typically contain mixed frequency data. In the OhioT1DM dataset, some of the variables, such as the continuous glucose monitoring (CGM) blood glucose levels, are recorded every 5 minutes. Meanwhile, other variables, such as the carbohydrate estimate for the meal and the exercise intensity are recorded with a much lower frequency. Concatenating these high-frequency variables over each one-hour interval produces a high-dimensional state vector and directly using these states as input of the treatment policy would yield very noisy decision rules. A naive approach is to use some ad-hoc summaries of the high-frequency data for policy learning. However, this might produce a sub-optimal policy due to the information loss.

Recently, there is a growing line of research in the statistics literature for policy learning and/or evaluation in infinite horizons. Some references include Ertefaie and Strawderman (2018); Luckett et al. (2020); Shi et al. (2022); Liao et al. (2020, 2021); Shi et al. (2021c); Chen et al. (2022); Ramprasad et al. (2022); Li et al. (2022); Xu et al. (2020). In the computer science literature, there is a huge literature on developing reinforcement learning (RL) algorithms in infinite horizons. These algorithms can be casted into as model-free and model-based algorithms. Popular model-free RL methods include value-based methods that model the expected return starting from a given state (or state-action pair) and compute the optimal policy as the greedy policy with respect to the value function (Ernst et al., 2005; Riedmiller, 2005; Mnih et al., 2015; Van Hasselt et al., 2016; Dabney et al., 2018), and policy-based methods that directly search the optimal policy among a parameterized class via policy gradient or action critic methods (Schulman et al., 2015a, b; Koutník et al., 2013; Mnih et al., 2016; Wang et al., 2017b). Model-based algorithms are different from the model-free algorithms in the sense that they model the transition dynamics of the environment and use the model of environment to derive or improve policies. Popular model-based RL methods include Guestrin et al. (2002); Janner et al. (2019); Lai et al. (2020); Li et al. (2020), to name a few. See also Arulkumaran et al. (2017); Sutton and Barto (2018); Luo et al. (2022) for more details. These methods cannot be directly applied to datasets such as OhioT1DM as they haven’t considered the mixed frequency data.

In the RL literature, a few works have considered dimension reduction to handle the high dimensional state system. In particular, Murao and Kitamura (1997) proposed to segment the state space and learn a cluster representation of the states. Whiteson et al. (2007) proposed to divide the state space into tilings to represent each state. Both papers proposed to discretize the state space for dimension reduction. However, it can lead to considerable information loss (Wang et al., 2017a). Sprague (2007) proposed a iterative dimension reduction method using neighborhood components analysis. Their method uses a linear basis function to model the Q-function and cannot allow more general nonlinear function approximation. Recently, there are a few works that employ principal components analysis (PCA) for dimension reduction in RL (Curran et al., 2015, 2016; Parisi et al., 2017). However, none of the aforementioned papers formally established the theoretical guarantees for their proposals. Moreover, these methods are motivated by applications in games or robotics and their generalization to mHealth applications with mixed frequency data remains unknown.

In the DTR literature, a few works considered mixed frequency data which include both scalar and functional covariates. Specifically, McKeague and Qian (2014) proposed a functional regression model for optimal decision making with one functional covariate. Ciarleglio et al. (2015) and Ciarleglio et al. (2016) extended their proposal to a more general setting with multiple scalar and functional covariates. Ciarleglio et al. (2018) considered variable selection to handle the mixed frequency data. Laber and Staicu (2018) applied functional PCA to the functional covariates for dimension reduction. All these works considered single stage decision making. Their methods are not directly applicable to the infinite horizon settings.

Our contributions are as follows. Scientifically, mixed frequency data frequently arise in mHealth applications. Nonetheless, it has been less explored in the infinite horizon settings. Our proposal thus fills a crucial gap, and greatly extends the scope of existing approaches to learning DTRs. Methodologically, we propose a deep spectral Q-learning algorithm for dimension reduction. The proposed algorithm achieves a better empirical performance than those that either directly use the original mixed frequency data or its ad-hoc summaries as input of the treatment policy. Theoretically, we derive an upper error bound for the regret of the proposed policy and decompose this bound into the sum of approximation error and estimation error. Our theories offer some general guidelines for practitioners to select the number of principal components in the proposed algorithm.

The rest of this paper is organized as follows. We introduce some background about DTR and the mixed frequency data in Section 2. We introduce the proposed method to estimate the optimal DTR in Section 3 and study its theoretical properties in Section 4. Empirical studies are conducted in Section 5. In Section 5.3, we apply the proposed method to the OhioT1DM dataset. Finally, we conclude our paper in Section 6.

Neural network with ReLU activation functions in 2.3 is commonly used in value-based reinforcement learning algorithms. However, in medical studies, the training dataset is often of limited size, with a few thousands or 10 thousands of observations in total (see e.g., Marling and Bunescu, 2020; Liao et al., 2021). Meanwhile, the data contains high frequency state variables, which yields a high-dimensional state system. Directly using these states as input will procedure a very noisy policy. This motivates us to consider dimension reduction in RL.

A naive approach for dimension reduction is to use some summary statistics of the high frequency state as input for policy learning. For instance, on the OhioT1DM dataset, the average of CGM blood glucose levels between two treatment decision points can be used as the summary statistic, as in Zhu et al. (2020); Shi et al. (2022); Zhou et al. (2022). In this paper, we propose to use principal component analysis to reduce the dimensionality of $\{Z_{i,t,j}\}_{j=1}^{J}$. We expect that using PCA can preserve more information than some ad-hoc summaries (e.g., average).

To apply PCA in the infinite horizon setting, we need to impose some stationarity assumptions on the concatenated high dimensional variables $Z_{i,t}^{T}\in\mathbb{R}^{m}$: $E[Z_{i,t}]=\bf\mathbf{\mu}$ and $Cov[Z_{i,t}]=\mathbf{G}$ for some mean vector $\mathbf{\mu}\in\mathbb{R}^{m}$ and covariance matrix $\mathbf{G}\in\mathbb{R}^{m\times m}$ that are independent of $t$. In real data application, we can test whether the concatenated high-frequency variable $Z_{i,t}$ is weak stationary (see e.g., Dickey and Fuller, 1979; Said and Dickey, 1984; Kwiatkowski et al., 1992). If it is weak stationary, the concatenated high frequency covariate $Z_{i,t}$ will automatically satisfy the two assumptions above. Similar assumptions have been widely imposed in the literature (see e.g., Shi et al., 2021b; Kallus and Uehara, 2022). Without loss of generality, we assume $\mathbf{\mu}=0_{m}$ for simplicity of notations. For the covariance matrix $\mathbf{G}$, it is generally unknown. In practice, we recommend to use the sample covariance estimator $\hat{\mathbf{G}}$.

We describe our procedure as follows. By the spectral decomposition, we have $\mathbf{\hat{G}}=\sum_{k=1}^{m}{\hat{\lambda}_{k}\hat{U}_{k}\hat{U}_{k}^{T}}$, where $\lambda_{1}\geq\lambda_{2}\geq...\geq 0$ are the eigenvalues and $\hat{U}_{k}$’s are the corresponding eigenvectors. This allows us to represent ${Z}_{i,t}$ as $\sum_{k=1}^{m}{\hat{\lambda}_{k}^{1/2}\hat{V}_{k}^{(i,t)}\hat{U}}_{k}$, where $\hat{V}_{k}^{(i,t)}$s are the estimated principal component scores, given by $\hat{\lambda}_{k}^{-1/2}{Z}_{i,t}^{T}\hat{U}_{k}$. For any $\kappa$, the estimated principal component scores $\hat{\mathbf{V}}_{i,t,\kappa}=(\hat{V}_{1}^{(i,t)},\hat{V}_{2}^{(i,t)},...,\hat{V}_{\kappa}^{(i,t)})$ correspond to the $\kappa\leq m$ largest eigenvalues of the concatenated high frequency variable $Z_{i,t}$. When $\kappa=m$, using these principal component scores is equivalent to using the original high-frequency variable ${Z}_{i,t}$. We will approximate $Q^{\pi}({X}_{i,t},\mathbf{V}_{i,t,m},{A}_{i,t})$ by $Q^{\pi}({X}_{i,t},\hat{\mathbf{V}}_{i,t,\kappa},{A}_{i,t})$ and propose to use neural fitted Q-iteration algorithm by Riedmiller (2005) to learn the estimated optimal policy. We detail our procedure in Algorithm 1.

In Algorithm 1, we fit $|\mathcal{A}|$ neural networks corresponding to each $a$ in $Q(s,a)$. This is reasonable in settings where the action space is small. When the dataset is small (such as the OhioT1DM dataset), we recommend to set $n$ to $N(T-1)$ such that all the data transactions (instead of a random subsample) will be used in each iteration.

Similar to the original neural fitted Q-iteration algorithm in Riedmiller (2005), the intuition of this algorithm is also based on the Bellman optimality equation (1). In each step $k$ of Algorithm 1, $\tilde{Q}_{k}$ estimates $Q^{\pi^{*}}$ and response $R_{i,t}(\tilde{Q}_{k})=Y_{i,t}+\gamma\operatorname*{max}_{a\in\mathcal{A}}\tilde{Q}_{k}(X_{i,t+1},\hat{\mathbf{V}}_{i,t+1,\kappa},a)$ corresponds to the right hand side of equation (1). Therefore, fitting the regression of $R_{i,t}$ with $\tilde{Q}_{k+1}$ is to solve the Bellman optimality equation. The key difference between Algorithm 1 and original neural fitted Q-iteration algorithm is that the high dimensional inputs $Z_{i,t}^{T}=[Z_{i,t,1}^{T},Z_{i,t,2}^{T},\ldots,Z_{i,t,J}^{T}]$ is involved in the state space, and is mapped to a lower dimensional vector $\hat{\mathbf{V}}_{i,t,\kappa}$ during the learning process, so the neural networks $\tilde{Q}_{0}(x,v_{\kappa},a)$ takes principle components $\hat{\mathbf{V}}_{i,t,\kappa}$ rather than original high dimensional $Z_{i,t}$ as input.

Before discussing the asymptotic properties of our proposed Q-learning method, we introduce some notations.

We define $L_{\mathcal{F}_{0}}=\sup_{f\in\mathcal{F}_{0}}{\sup_{x\neq y}{\frac{{|f(y)-f(x)|}^{2}}{\left\|y-x\right\|^{2}}}}$ as the Lipschitz constant for the sparse ReLU class $\mathcal{F}_{SReLU}(L,\{d_{j}\}_{i=0}^{L+1},s,V_{max})$ used in $\mathcal{F}_{0}$. Note that this class $\mathcal{F}_{0}$ is used in the original neural fitted Q-iteration algorithm to model the Q-function, where the dimension of high frequency part $Z$ in state is not reduced through PCA. We further define a function class $\mathcal{F}_{2}$ such that it models the Q-function by first converting high frequency part $Z$ into its principal component scores and then use a sparse ReLU neural network to obtain the resulting Q-function. More specifically, $\mathcal{F}_{2}$ is a set of functions $\{f_{2}\}$, such that $f_{2}(x,z,a)=f_{0}(x,\hat{v}_{\kappa},a)$, where $f_{0}\in\mathcal{F}_{0}$ and $\hat{v}_{\kappa}$ is the vector containing first $\kappa$ principle components of $Z$. Note that $\mathcal{F}_{2}$ is the function class that we use in Algorithm 1 to model the Q-function. That is, $\tilde{Q}_{k}\in\mathcal{F}_{2},k\in\{1,2,...,K\}$ (formal definition of $\mathcal{F}_{2}$ and another function class $\mathcal{F}_{1}$ not mentioned here are in Section A of Appendix).

In addition, we introduce the Hölder smooth function class by $\mathcal{C}_{r}(\mathcal{D},\beta,H)$ with $\mathcal{D}\in\mathbb{R}^{r}$ to be the set of function input. The definition is:

We next construct a network structure $\mathcal{G}(\{p_{j},t_{j},\beta_{j},H_{j}\}_{j\in[q]})$ with the component function on each layer of this network belonging to Holder smooth function class $\mathcal{C(\mathcal{D},\beta,H)}$, which is called composition of Holder Smooth functions. This composition network contains $q$ layers, with each layer being $g_{j}:{[a_{j},b_{j}]}^{p_{j}}\rightarrow{[a_{j+1},b_{j+1}]}^{p_{j+1}}$, such that $g_{jk}$ the $k$th component ($k\in[p_{j+1}]$) in layer $j$ satisfies that $g_{jk}\in\mathcal{C}_{t_{j}}({[a_{j},b_{j}]}^{t_{j}},\beta_{j},H_{j})$ with $t_{j}\leq p_{j}$ inputs. Similar to the definition of $\mathcal{F}_{0}$, we can define the function class $\mathcal{G}_{0}(\{p_{j},t_{j},\beta_{j},H_{j}\}_{j\in[q]})$(simply denoted as $\mathcal{G}_{0}$) on $\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ such that each function $g\in\mathcal{G}_{0}$ satisfies that $g(.,a)\in\mathcal{G}(\{p_{j},t_{j},\beta_{j},H_{j}\}_{j\in[q]})$ for $\forall a\in\mathcal{A}$. The relation between function class $\mathcal{G}_{0}$ and the network structure $\mathcal{G}$ is similar to the relation between function class $\mathcal{F}_{0}$ and the neural network $\mathcal{F}_{SReLU}$ in Definition 4.1. See Definition 4.1 of Fan et al. (2020) for more details on $\mathcal{G}_{0}(\{p_{j},t_{j},\beta_{j},H_{j}\}_{j\in[q]})$.

Next we will introduce the three major assumptions for our theorems:

Among the three assumptions, the exponential decaying structure of eigenvalues in Assumption 4.3 can be commonly found in the literature of high-dimensional and functional data analysis (see e.g., Reiß and Wahl, 2020; Crambes and Mas, 2013; Jirak, 2016). This assumption is to control the information loss caused by using the first $\kappa$ principal component scores of $Z_{i,t}$ only. Assumption 4.4 is about the consistency of the estimators $\hat{\mathbf{G}}$ and similar assumptions are imposed in the literature of functional data analysis (see e.g., Laber and Staicu, 2018; Staicu et al., 2014). Using similar arguments in proving Theorem 5.2 of Zhang and Wang (2016), we can show that such an assumption holds in our setting as well. It is to bound the error caused by the estimation of the covariance matrix. Assumption 4.5 is referred to as the completeness assumption in the literature (see e.g., Chen and Jiang, 2019; Uehara et al., 2022, 2021). This assumption is automatically satisfied when the transition kernel and the reward function satisfy certain smoothness conditions.

Our first theorem is concerend with the nonasymptotic convergence rate of $\tilde{Q}_{K}$ in Algorithm 1.

More details of neural network structure and sample size assumptions for Theorem 4.6 can be found in Section B of Appendix. Theorem 4.6 provides an error bound on the estimated Q-function $\tilde{Q}_{K}$. Based on this theorem, we further establish the regret bound of the estimated policy $\pi_{K}$ obtained via Algorithm 1. Toward that end, we need another assumption:

The margin type condition Assumption 4.7 is commonly used in the literature. Specifically, in classification, the margin conditions are imposed to bound the excess risk (Tsybakov, 2004; Audibert and Tsybakov, 2007). In dynamic treatment regime, a similar assumption is introduced for proving the convergence of State-Value function in a finite horizon setting (Qian and Murphy, 2011; Luedtke and van der Laan, 2016). In RL, these assumptions were introduced by (Shi et al., 2022) to obtain sharper regret bound for the estimated optimal policy.

The proofs of the two theorems are included in Section C of Appendix. We summarize our theoretical findings below. First, we notice that the convergence rate of regret in Theorem 4.8 is faster than the convergence of estimated Q-functions in Theorem 4.6. This is due to the margin type Assumption 4.7 which enables us to obtain a sharper error bound. Similar results have been established in the literature. See e.g., Theorem 3.3 in Audibert and Tsybakov (2007), Theorem 3.1 in Qian and Murphy (2011), Theorems 3 and 4 in Shi et al. (2022).

From Theorem 4.8, it can be seen that the regret bound is mainly determined by three parameters: the sample size $n$, the number of principal components $\kappa$ and the number of iterations $K$. Here, the first term $|\mathcal{A}|(n^{\alpha^{*}}{(\log^{\xi^{*}}{n})}+d_{1}^{*}\kappa)n^{-1}\log^{\xi^{*}+1}{n}$ on the right-hand-side corresponds to the estimation error, which decreases with $n$ and increases with $\kappa$. The second term $L_{\mathcal{F}_{0}}(e^{-\zeta\kappa}-e^{-\zeta m}+n^{-2\Delta})$ corresponds to the approximation error, which decreases with both $n$ and $\kappa$. The remaining term $\frac{\gamma^{K+1}}{(1-\gamma)^{2}}R_{max}$ is the optimization error which will decrease as the iteration number $K$ in Algorithm 1 grows.

Compared with the existing results on the convergence rate of deep fitted Q-iteration algorithm (see Theorem 4.4 of Fan et al., 2020), our theorems additionally characterize the dependence upon the number of principal components. Specifically, selecting the first $\kappa$ principal components induces the information loss (e.g., bias) that is of the order $e^{-\zeta\kappa}-e^{-\zeta m}$, but reduces the model complexity caused by high frequency variables from $d_{1}^{*}m$ to $d_{1}^{*}\kappa$ and hence the variance of the policy estimator. This represents a bias-variance trade-off. Notice that the bias decays at an exponential order, when the training data is small, reducing the model complexity can be more beneficial. Thus, our algorithm is expected to perform better than the original fitted Q-iteration algorithm in small samples, as shown in our numerical studies.

Finally, The number of principal components shall diverge with $n$ to ensure the consistency of the proposed algorithm. Based on the two theorems, the optimal $\kappa^{*}$ that balances the bias and variance trade-off shall satisfy $\kappa^{*}\asymp\log(n)$ (details are given in Section D of Appendix). Thus, when $n$ goes to infinity, we will eventually take $\kappa^{*}=m$ and our Algorithm 1 will be equivalent to the original neural fitted Q-iteration. This is just an asymptotic guideline for selecting the number of principal components. We provide some practical guidelines in the next section.

In summary, we propose a deep spectral fitted Q-iteration algorithm to handle mixed frequency data in infinite horizon settings. The algorithm relies on the use of PCA for dimension reduction and the use of deep neural networks to capture the non-linearity in the high-dimensional system. In theory, we establish a regret bound for the estimated optimal policy. Our theorem provides an asymptotic guideline for selecting the number of principal components. In empirical studies, we demonstrate the superiority of the proposed algorithm over baseline methods without dimension reduction or use ad-hoc summaries of the state. We further offer practical guidelines to select the number of principal components. The proposed paper is built upon on PCA and deep Q-learning. It is worthwhile to investigate the performance of other dimension reduction and RL methods. We leave it for future research.