Off-Policy Evaluation with Irregularly-Spaced, Outcome-Dependent Observation Times

Consider a study where treatment decisions are made at irregularly-spaced time points $\{T_{0}=0\,,T_{1}\,,T_{2},\ldots,T_{k},\ldots\}$. For each $k\geq 1$, let $X_{k}=T_{k}-T_{k-1}$ denote the gap times, and let $\bm{S}_{k}\in\mathbb{S}$ be a vector of state variables that summarizes the subject’s information collected up to and including time $T_{k}$. At each $T_{k}$, investigators can observe the status of subject $\bm{S}_{k}=S(T_{k})$ and the gap time $X_{k}=T_{k}-T_{k-1}$, and then take an action $A_{k}=A(T_{k})\in\mathbb{A}$ based on $\bm{S}_{k}$ and $X_{k}$. The reward of action $A_{k}$ is observed at $T_{k+1}$, and thus denoted as $R(T_{k+1})$. For simplicity, we assume that the state space $\mathbb{S}$ is a subspace of $\mathbb{R}^{d}$ where $d$ is the number of state vectors, and the action space $\mathbb{A}$ is a discrete space $\{0,1\,,\ldots\,,m-1\}$ with $m$ denoting the number of actions. We also assume that $\bm{S}_{0}$, $X_{0}$, and $A_{0}$ exist and are observable at $T_{0}$. An illustration of the data structure is given in Figure 1.

Let $\pi(\cdot|\cdot,\cdot)$ denote a policy that maps the state and gap time $(\bm{s},x)\in\mathbb{S}\times[0,\infty]$ to a probability mass function $\pi(\cdot|\bm{s},x)$ on $\mathbb{A}$. Under such a policy $\pi$, a decision maker will take action $A_{k}=a$ with probability $\pi(a|\bm{S}_{k},X_{k})$ for $a\in\mathbb{A}$ and $k\geq 0$. Aiming at off-policy evaluation for any given dynamic policy $\pi$, we first define the corresponding value functions. Unlike the existing literature (Shi et al., 2021), we treat $\{T_{k}\}_{k\geq 0}$ as jump points of a counting process $N(t)$ and propose the following two value functions

	$\displaystyle V_{k}^{\pi}(\bm{s},x,t)=$	$\displaystyle\ \operatorname{E^{\pi}}\left[\int_{t}^{\infty}\gamma^{u-t}R(u)\ \mathrm{d}N(u)\nonscript\;\middle|\nonscript\;\bm{S}_{t}=\bm{s}\,,X_{t}=x\,,T_{k}=t\right]$		(1)
	$\displaystyle V_{\mathcal{I},k}^{\pi}(\bm{s},x,t)=$	$\displaystyle\ \operatorname{E^{\pi}}\left[\int_{t}^{\infty}\gamma^{u-t}R(u)\ \mathrm{d}u\ \nonscript\;\middle|\nonscript\;\bm{S}_{t}=\bm{s}\,,X_{t}=x\,,T_{k}=t\right]\,.$		(2)

Here, $0<\gamma<1$ is a discounted factor that balances the immediate and long-term effects of treatments; $R(t)$, an extension of the aforementioned notation $R(T_{k})$, is redefined as a reward process only observed at $\{T_{k}\}_{k\geq 1}$, but could be potentially available at time $t\geq 0$; and $E^{\pi}$ denotes the expectation when all the actions follow policy $\pi$. By the definition of $N(t)$, $V_{k}^{\pi}(\bm{s},x,t)$ defined in (1) can be rewritten as an expectation of discounted cumulative reward, given as

$\displaystyle V^{\pi}_{k}(\bm{s},x,t)$

$\displaystyle\equiv\operatorname{E^{\pi}}\left[\sum_{j\geq k+1}\gamma^{T_{j}-T_{k}}R(T_{j})\nonscript\;\middle|\nonscript\;\bm{S}_{k}=s,X_{k}=x,T_{k}=t\right]$

(3)

Therefore, we refer to $V^{\pi}_{k}(\bm{s},x,t)$ as the value function based on cumulative reward and $V^{\pi}_{\mathcal{I},k}(\bm{s},x,t)$ as the value function based on integrated reward. Similarly, we define $Q_{k}^{\pi}(\bm{s},x,a,t)=\operatorname{E^{\pi}}\left[\sum_{j\geq k+1}\gamma^{T_{j}-T_{k}}R(T_{j})\nonscript\;\middle|\nonscript\;\bm{S}_{k}=\bm{s}\,,X_{k}=x\,,A_{k}=a\,,T_{k}=t\right]$ as the action-value function (or Q-function) based on cumulative reward, while the one based on integrated reward is defined as $Q_{\mathcal{I},k}^{\pi}(\bm{s},x,a,t)=\operatorname{E^{\pi}}\left[\int_{t}^{\infty}\gamma^{u-t}R(u)\ \mathrm{d}u\ \nonscript\;\middle|\nonscript\;\bm{S}_{k}=\bm{s}\,,X_{k}=x\,,A_{k}=a\,,T_{k}=t\right]\,.$

Remark 1. Both value functions, whether based on cumulative or integrated reward, may depend not only on the current value of state $\bm{S}_{k}=\bm{s}$, but also on the time $T_{k}=t$ at which the policy is evaluated. This differs from the value functions developed under regularly-spaced observations and introduces challenges to the policy evaluation process. To address this issue, we incorporate $X_{k}$ into the value functions and will demonstrate in subsequent sections that, under modified Markov and time-homogeneous assumptions, the value function at $T_{k}=t$ can be independent of $k$ and $t$, given the current state $\bm{S}_{k}$ and gap time $X_{k}$.

Remark 2. The two value functions differ in terms of quantity, interpretation, application scenarios, and estimation procedures. First, the value function based on cumulative reward, $V^{\pi}_{k}(\bm{s},x,t)$, can be highly correlated with the rate of observation process $N(t)$. For example, assuming $R(t)$ is a non-negative function, a higher rate of $N(t)$ would yield more decision stages $T_{j}$, leading to a higher policy value. In contrast, $V^{\pi}_{\mathcal{I},k}(\bm{s},x,t)$, which integrates the reward function over the follow-up period, provides a more principled evaluation of the policy as it remains independent of the number of observed decision stages. As a result, the two value functions are suited to different application scenarios, and the choice of value function should depend on how the reward function is defined. Taking our motivating PD treatment study as an example, when the outcome of interest is whether a patient’s probing pocket depths (PPD) are reduced after treatment, the reward $R(T_{k})$ is defined as an indicator of PPD reduction, making the value function based on cumulative reward preferable. Conversely, when the reward function $R(t)$ is defined as the PPD measurement at time $t$ for $t\in[0,\infty]$, policy evaluation should be independent of the measurement times, and the value function based on integrated reward provides a comprehensive assessment of the patient’s oral health status. Finally, since $R(t)$ is unobservable when $t\neq T_{j}$, additional assumptions about $R(t)$ and $N(t)$ are required to ensure that the expectation of $\int_{t}^{\infty}\gamma^{u-t}R(u)\mathrm{d}u$ is estimable. This makes the evaluation process for $V^{\pi}_{\mathcal{I},k}(\bm{s},x,t)$ different from that for $V^{\pi}_{k}(\bm{s},x,t)$; further details are provided in Section 4.

Now, we posit assumptions on the states, actions, rewards, and gap times.

• (A1)

  (Markov and time-homogeneous) For all $k\geq 0$, and for any measurable set $\mathcal{B}\subset\mathbb{S}\times\mathbb{R}^{+}$ with $\bm{S}_{k}$ and $X_{k}$ as the current state and gap times, respectively, the next state variable $\bm{S}_{k+1}$ and gap time $X_{k+1}$ satisfy

  $\displaystyle P\bigg{[}\begin{pmatrix}\bm{S}_{k+1}\\ X_{k+1}\end{pmatrix}\in\mathcal{B}\bigg{|}\{\bm{S}_{j}\}_{0\leq j\leq k},\{X_{j}\}_{0\leq j\leq k},\{A_{j}\}_{0\leq j\leq k},\{T_{j}\}_{0\leq j\leq k}\bigg{]}=\mathcal{P}(\mathcal{B};\bm{S}_{k},X_{k},A_{k})\,.$

  Here, $\mathcal{P}$ denotes the joint distribution of the next state and next gap time conditional on the current state, action, and gap time.

• (A2)

  (Conditional mean independence) For all $k\geq 0$ and for some bounded function $r$, the reward $R(T_{k+1})$ satisfies

  $\displaystyle E\left[R(T_{k+1})|\{R(T_{j}),\bm{S}_{j},X_{j},A_{j}\}_{1\leq j\leq k},\bm{S}_{k+1}\,,X_{k+1}\right]=r(\bm{S}_{k},X_{k},A_{k},\bm{S}_{k+1},X_{k+1})\,.$

Assumption A1 implies that, given the current state, action, and gap time, the joint distribution of $(\bm{S}_{k+1},X_{k+1})$ does not depend on the historical observations $\{\bm{S}_{j},X_{j},A_{j},T_{j}\}_{j<k}$, or the current observation time $T_{k}$. Instead of assuming the marginal distributions of the states and the gap times separately, we present assumption A1 based on the joint distribution of the next state and the next gap time. This indicates that, not only the marginal distributions of $\bm{S}_{k+1}$ and $X_{k+1}$, but also the dependence between $\bm{S}_{k+1}$ and $X_{k+1}$, should be independent of the historical observations and the current observation time $T_{k}$, conditionally on $\{\bm{S}_{k},X_{k},A_{k}\}$. For instance, in the PD treatment study, the recall interval $X_{k+1}$ often correlates strongly with the PD status at the subsequent visit $\bm{S}_{k+1}$. This correlation, along with the values of $\bm{S}_{k+1}$ and $X_{k+1}$, may be influenced by the current PD status of patients and the treatments provided by the dentists. However, it’s reasonable to assume conditional independence of these factors from historical states, treatments, and recall intervals, given the latest observations $\{\bm{S}_{k},X_{k},A_{k}\}$.

By denoting $\widetilde{\bm{S}}_{k+1}=(\bm{S}_{k+1}^{\top},X_{k+1})^{\top}$ as an extended vector of state variables, it is not hard to see that A1 takes a similar form as the Markov and time-homogeneous assumption commonly adopted in the RL literature. However, in the following sections, we will demonstrate that the gap times $\{X_{k}\}_{k\geq 1}$ and the state variables $\{\bm{S}_{k}\}_{k\geq 1}$ play different roles in the estimation procedure and theoretical properties. This distinction explains why we do not simply include $X_{k+1}$ in the state vector.

For illustration, we present two data-generating schemes that satisfy assumption A1, while A1 itself does not specify any structure on the relationship between $\bm{S}_{k+1}$ and $X_{k+1}$. The two schemes bridge the proposed framework with modulated renewal process models and Markov decision process models, as well as providing guidelines for constructing specific models on the states, actions, and gap times.

Scheme 1. Consider a stochastic process where, given $\bm{S}_{k}=\bm{s}$, $X_{k}=x$ and $A_{k}=a$, the next state $\bm{S}_{k+1}$ is generated as

$\displaystyle\text{(S1-s)}\quad P\left(\bm{S}_{k+1}\in\mathcal{B}_{s}|\bm{S}_{k}=\bm{s},X_{k}=x,A_{k}=a,T_{k},\{\bm{S}_{j},X_{j},A_{j},T_{j}\}_{j<k}\right)=\mathcal{P}_{S}(\mathcal{B}_{s};\bm{s},x,a)$

for any measurable set $\mathcal{B}_{s}\subset\mathbb{S}$, while given the value of next state $\bm{S}_{k+1}=\bm{s^{\prime}}$, the next gap time $X_{k+1}$ is generated as

$\displaystyle\text{(S1-x)}\quad P\left(X_{k+1}\leq x^{\prime}|\bm{S}_{k+1}=\bm{s^{\prime}},\bm{S}_{k}=s,X_{k}=x,A_{k}=a,T_{k},\{\bm{S}_{j},X_{j},A_{j},T_{j}\}_{j<k}\right)=\mathcal{P}_{X}(x^{\prime};\bm{s^{\prime}},\bm{s},x,a)$

for any $x^{\prime}\in\mathbb{R}^{+}$. Here, $\mathcal{P}_{S}$ denotes the transition function of the next state conditional on the current state, action, and gap time, and $\mathcal{P}_{X}$ denotes the cumulative distribution function of next gap time conditional on the current state, gap time, action, as well as the next state.

Remark 3. As a special case, when $\mathbb{S}$ consists of finite or countable number of states and $\mathcal{P}_{S}(\mathcal{B}_{s};\bm{s},x,a)$ is independent of the value of current gap time $x$, (S1-s) and (S1-x) together form a time-homogenous semi-Markov decision process (Bradtke and Duff, 1994; Parr, 1998; Du et al., 2020), where $\{\bm{S}_{k}\}_{k\geq 0}$ is the Markov chain and $X_{k+1}$ is the sojourn time in the $k$th state $\bm{S}_{k}$.

Remark 4. Considering that $\{X_{k}\}_{k\geq 1}$ are the gap times of a counting process, we can restate the generation procedure in the model (S1-x) using the notations of the counting process. Let $\{N(t)\}_{t\geq 0}$ denote a right-continuous counting process with jumps at $\{T_{k}\}_{k\geq 1}$. Then, we have $T_{N(t)}=T_{k}$, $X_{N(t)}=X_{k}$, $A_{N(t)}=A_{k}$, $\bm{S}_{N(t)}=\bm{S}_{k}$, and $\bm{S}_{N(t)+1}=\bm{S}_{k+1}$ for all $T_{k}\leq t<T_{k+1}$. According to (S1-s), the distribution/generation of next state $\bm{S}_{N(t)+1}$ is fully determined by $\{\bm{S}_{N(t)},X_{N(t)},A_{N(t)}\}$, and therefore could be assumed available at time $t$. Let $\mathcal{F}(t)=\left\{\{N(u)\}_{u\leq t},\{\bm{S}_{N(u)}\}_{u\leq t},\{X_{N(u)}\}_{u\leq t},\{A_{N(u)}\}_{u\leq t},\bm{S}_{N(t)+1}\right\}$ denote all information available at $t$. It is not hard to observe that if $N(t)$ satisfies

	$\displaystyle\text{(S1-x${}^{\prime}$)}\qquad E\left[\mathrm{d}N(t)|\mathcal{F}(t-)\right]=$	$\displaystyle\lambda\{t-T_{N(t-)};\bm{S}_{N(t-)+1},\bm{S}_{N(t-)},X_{N(t-)},A_{N(t-)}\}\mathrm{d}t$
	$\displaystyle=$	$\displaystyle\lambda\{t-T_{k};\bm{S}_{k+1},\bm{S}_{k},X_{k},A_{k}\}\mathrm{d}t\,,\text{for}\ k=\max\{j;T_{j}<t\}\quad$

or some intensity function $\lambda$, then for each $k\geq 0$, the gap time $X_{k+1}$ with the hazard function $h(x)=\lambda\{x;\bm{S}_{k+1},\bm{S}_{k},X_{k},A_{k}\}$ satisfies model (S1-x) too. On the other hand, if the distribution function $\mathcal{P}_{X}$ in the model (S1-x) is absolutely continuous wrt. the Lebesgue measure, then the counting process defined by $N(t)=\sum_{j\geq 1}I(T_{j}\leq t)$ also satisfies (S1-x^′). As a matter of fact, by viewing $\{\bm{S}_{k},A_{k},\bm{S}_{k+1}\}_{k\geq 0}$ as the covariates available at $\{T_{k}\}_{k\geq 0}$, (S1-x^′) takes a similar form to the modulated renewal process model. The assumption is natural for user-initiated visits and has been popularly adopted in the statistical literature (Cox, 1972; Lin et al., 2013; Chen et al., 2018).

Scheme 2. As an alternative, scheme 2 generates the gap time first and then determines the value of the next state. Specifically, assume that the next gap time $X_{k+1}$ satisfies

$$\text{(S2-x)}\quad P\left(X_{k+1}\leq x^{\prime}|\bm{S}_{k}=s,X_{k}=x,A_{k}=a,T_{k},\{\bm{S}_{j},X_{j},A_{j},T_{j}\}_{j<k}\right)=\widetilde{\mathcal{P}}_{X}(x^{\prime};\bm{s},x,a)\qquad$$

for $x^{\prime}\in\mathbb{R}^{+}$. Then, given the value $X_{k+1}=x^{\prime}$, the next state $\bm{S}_{k+1}$ satisfies

$$\text{(S2-s)}\quad P\left(\bm{S}_{k+1}\in\mathcal{B}_{s}|X_{k+1}=x^{\prime},\bm{S}_{k}=\bm{s},X_{k}=x,A_{k}=a,T_{k},\{\bm{S}_{j},X_{j},A_{j},T_{j}\}_{j<k}\right)=\widetilde{\mathcal{P}}_{S}(\mathcal{B}_{s};x^{\prime},\bm{s},x,a)\,,$$

for any measurable set $\mathcal{B}_{s}\subset\mathbb{S}$. Here, $\widetilde{\mathcal{P}}_{X}$ denotes the cumulative distribution function of the next gap time conditional on the current state, gap time, and action, while $\widetilde{\mathcal{P}}_{S}$ denotes the transition function of the next state conditional on the current state, action, gap time, as well as the next gap time. Let $\mathcal{\widetilde{F}}(t)=\left\{\{N(u)\}_{u\leq t},\{\bm{S}_{N(u)}\}_{u\leq t},\{X_{N(u)}\}_{u\leq t},\{A_{N(u)}\}_{u\leq t}\right\}$. With arguments similar to Remark 4, it can be shown that (S2-x) is equivalent to

$$\text{(S2-x${}^{\prime}$)}\quad E\left[\mathrm{d}N(t)\big{|}\mathcal{\widetilde{F}}(t-)\right]=\lambda(t-T_{N(t-)};\bm{S}_{N(t-)},X_{N(t-)},A_{N(t-)})=\lambda(t-T_{k};\bm{S}_{k},X_{k},A_{k})\,,$$

for $k=\max\{j;T_{j}<t\}$, and takes the form of a modulated renewal process.

To estimate $V^{\pi}_{k}(\bm{s},x,t)$ and $Q^{\pi}_{k}(\bm{s},x,a,t)$, we first present the following two lemmas, whose proof is available in Section A.1 of the Supplementary Material.

Lemma 1. Under assumptions A1 and A2, the action-value function based on cumulative reward $Q^{\pi}_{k}(\bm{s},x,a,t)$ exists and satisfies $Q^{\pi}_{k}(\bm{s},x,a,t)=Q^{\pi}_{0}(\bm{s},x,a,0)$ for all $k\geq 0$, $t\geq 0$, $s\in\mathbb{S}$, $x\in\mathbb{R}^{+}$ and $a\in\mathbb{A}$.

By Lemma 1, the value of $Q^{\pi}_{k}(\bm{s},x,a,t)$ is time-homogeneous, meaning the expected value of a policy $\pi$ does not depend on the number of past observations $k$ or the current observation time $t$. Therefore, we can omit both $k$ and $t$, and denote the action-value function simply as $Q^{\pi}(\bm{s},x,a)$ for brevity. Similarly, the state-value function of policy $\pi$ can also be presented as $V^{\pi}(\bm{s},x)=\sum_{a\in\mathbb{A}}Q^{\pi}(\bm{s},x,a)\pi(a|\bm{s},x)\,,$ and the value function of policy $\pi$ under some reference distribution $\mathcal{G}$ can be presented as $V^{\pi}(\mathcal{G})=\int\sum_{a\in\mathbb{A}}Q^{\pi}(\bm{s},x,a)\pi(a|\bm{s},x)d\mathcal{G}(\bm{s},x)\,,$ both of which are independent of $k$ and $t$. Furthermore, we derive a Bellman equation as formulated in Lemma 2.

Lemma 2. Under assumptions A1 and A2, the action-value function based on cumulative reward $Q^{\pi}(\bm{s},x,a)$ satisfies a Bellman equation

$\displaystyle Q^{\pi}(\bm{S}_{k},X_{k},A_{k})=$

$\displaystyle E\left[\gamma^{X_{k+1}}R(T_{k+1})+\gamma^{X_{k+1}}\sum_{a\in\mathbb{A}}Q^{\pi}(\bm{S}_{k+1},X_{k+1},a)\pi(a|\bm{S}_{k+1},X_{k+1})\bigg{|}\bm{S}_{k},X_{k},A_{k}\right]$

for all $k\geq 0$.

To solve the Bellman equation, let $Q(\bm{s},x,a;\bm{\theta})$ denote a model for $Q^{\pi}(\bm{s},x,a)$ indexed by $\bm{\theta}\in\Theta$, and satisfies $Q(\bm{s},x,a;\bm{\theta}_{\pi})=Q^{\pi}(\bm{s},x,a)$ for some $\bm{\theta}_{\pi}\in\Theta$. Define

$\displaystyle\delta_{k+1}(\bm{\theta})=$

$\displaystyle\gamma^{X_{k+1}}R(T_{k+1})+\gamma^{X_{k+1}}\sum_{a\in\mathbb{A}}Q(\bm{S}_{k+1},X_{k+1},a;\bm{\theta})\pi(a|\bm{S}_{k+1},X_{k+1})-Q(\bm{S}_{k},X_{k},A_{k};\bm{\theta})\,.$

According to Lemma 2, we have $E[\delta_{k+1}(\bm{\theta}_{\pi})|\bm{S}_{k},X_{k},A_{k}]=0$ for all $k\geq 0$. Assume the map $\bm{\theta}\rightarrow Q(\bm{s},x,a;\bm{\theta})$ is differentiable everywhere for each fixed $s$, $x$ and $a$, and let $\bm{\dot{Q}}_{\bm{\theta}}$ denote the gradient of $Q$ with respect to $\bm{\theta}$. Then, we construct an estimating function for $\bm{\theta}_{\pi}$ as

$$\bm{U}(\bm{\theta})=\frac{1}{\sum_{i}K_{i}}\sum_{i=1}^{n}\sum_{k=0}^{K_{i}-1}\bm{\dot{Q}}_{\bm{\theta}}(\bm{S}_{i,k},X_{i,k},A_{i,k};\bm{\theta})\widehat{\delta}_{i,k+1}(\bm{\theta})\,,$$

(4)

where $\{X_{i,k}\,,\bm{S}_{i,k}\,,A_{i,k}\,,R_{i}(T_{i,k+1});1\leq k\leq K_{i},1\leq i\leq n\}$ are $n$ independent and identically distributed trajectories of $\{\bm{S}_{k},X_{k},A_{k},R(T_{k+1})\}_{k\geq 0}$. $K_{i}\in\mathbb{Z}_{+}$ denotes the total number of decisions observed on the $i$th sequence, and $\delta_{i,k+1}(\bm{\theta})=\gamma^{X_{i,k+1}}R_{i}(T_{i,k+1})+\gamma^{X_{i,k+1}}\sum_{a\in\mathbb{A}}Q(\bm{S}_{i,k+1},X_{i,k+1},a;\bm{\theta})\pi(a|\bm{S}_{i,k+1},X_{i,k+1})-Q(\bm{S}_{i,k},X_{i,k},A_{i,k};\bm{\theta})\,.$

Generally, (4) can be solved by numerical methods, but in certain cases, it may have a closed form. For example, consider a B-spline model $Q(\bm{s},x,a;\bm{\theta})=\bm{\xi}(\bm{s},x,a)^{\top}\bm{\theta}\,,$ where $\bm{\xi}(\bm{s},x,a)=\{\Phi_{L}^{\top}(\bm{s},x)\mathbb{I}(a=0),\Phi_{L}^{\top}(\bm{s},x)\mathbb{I}(a=1),\ldots,\Phi_{L}^{\top}(\bm{s},x)\mathbb{I}(a=m-1)\}^{\top}\,,$ $\bm{\theta}=(\bm{\theta}_{0}^{\top},\bm{\theta}_{1}^{\top},\cdots,\bm{\theta}_{m-1}^{\top})^{\top}\,,$ $\Phi_{L}(\cdot)=\{\Phi_{L,1}(\cdot),\cdots,\Phi_{L,L}(\cdot)\}^{\top}$ is a vector consisting of $L$ B-spline basis functions, and $\bm{\theta}_{a}$ is a $L$-dimensional parameter vector for each $a\in\mathbb{A}$. Denote $\bm{\xi}_{i,k}=\bm{\xi}(\bm{S}_{i,k},X_{i,k},A_{i,k})$ and $\bm{\zeta}_{\pi,i,k+1}=\sum_{a\in\mathbb{A}}\bm{\xi}(\bm{S}_{i,k+1},X_{i,k+1},a)\pi(a|\bm{S}_{i,k+1},X_{i,k+1})\,.$ Then, the estimating equation takes the form $\bm{U}(\bm{\theta})=\frac{1}{\sum_{i}K_{i}}\sum_{i=1}^{n}\sum_{k=0}^{K_{i}-1}\bm{\xi}_{i,k}\Big{\{}\gamma^{X_{i,k+1}}R_{i}(T_{i,k+1})+(\gamma^{X_{i,k+1}}\bm{\zeta}_{\pi,i,k+1}-\bm{\xi}_{i,k})^{\top}\bm{\theta}\Big{\}}=0\,,$ and can be solved as

$\displaystyle\widehat{\bm{\theta}}_{\pi}=\biggl{\{}\frac{1}{\sum_{i}K_{i}}\sum_{i=1}^{n}\sum_{k=0}^{K_{i}-1}\bm{\xi}_{i,k}(\bm{\xi}_{i,k}-\gamma^{X_{i,k+1}}\bm{\zeta}_{\pi,i,k+1})^{\top}\biggr{\}}^{-1}\bigg{\{}\frac{1}{\sum_{i}K_{i}}\sum_{i=1}^{n}\sum_{k=0}^{K_{i}-1}\bm{\xi}_{i,k}\gamma^{X_{i,k+1}}R_{i}(T_{i,k+1})\bigg{\}}\,.$

Therefore, we propose to estimate $Q^{\pi}(\bm{s},x,a)$ by $\widehat{Q}^{\pi}(\bm{s},x,a)=Q(\bm{s},x,a;\widehat{\bm{\theta}}_{\pi})=\bm{\xi}(\bm{s},x,a)^{\top}\widehat{\bm{\theta}}_{\pi}\,,$ estimate state-value function $V^{\pi}(\bm{s},x)$ by $\widehat{V}^{\pi}(\bm{s},x)=\sum_{a\in\mathbb{A}}\pi(a|\bm{s},x)\bm{\xi}(\bm{s},x,a)^{\top}\widehat{\bm{\theta}}_{\pi}$, and estimate $V^{\pi}(\mathcal{G})$ by $\widehat{V}^{\pi}(\mathcal{G})=\bm{\zeta}_{\pi,\mathcal{G}}^{\top}\widehat{\bm{\theta}}_{\pi}\,,$ where $\bm{\zeta}_{\pi,\mathcal{G}}=\sum_{a\in\mathbb{A}}\int\pi(a|\bm{s},x)\bm{\xi}(\bm{s},x,a)\mathrm{d}\mathcal{G}(\bm{s},x)$.

In this subsection, we further propose an alternative policy evaluation procedure that relies on a specified model for the observation process. Here, the motivation for modeling the observation process arises from two aspects. First, as demonstrated in Subsection 2.3, under the revised time homogeneity and Markov assumptions for states and gap times, the irregularly spaced and outcome-dependent observation times can be represented as the occurrence times of a renewal process modulated by the state and action sequences (Cox, 1972). Therefore, positing a modulated renewal process model on the observation process is consistent with the structure of the data and helps us make more accurate inferences. Second, by employing an intensity model for the observation process, the integration of rewards can be transformed into a cumulative form. This enables us to leverage policy evaluation techniques based on cumulative reward for estimating the value function based on integrated reward. More details are given in Section 4.

For brevity, we focus on Scheme 1 given in Section 2.3; similar derivations can also be applied to Scheme 2. The modification begins with the Bellman equation outlined in Lemma 2. Under assumption (S1-x), we can express the equation as

		$\displaystyle Q^{\pi}(\bm{S}_{k},X_{k},A_{k})$
	$\displaystyle=$	$\displaystyle E\left[\gamma^{X_{k+1}}R(T_{k+1})+E\big{[}\gamma^{X_{k+1}}\sum_{a\in\mathbb{A}}\pi(a|\bm{S}_{k+1},X_{k+1})Q^{\pi}(\bm{S}_{k+1},X_{k+1},a)\big{|}\bm{S}_{k+1},\bm{S}_{k},X_{k},A_{k}\big{]}\bigg{|}\bm{S}_{k},X_{k},A_{k}\right]$
	$\displaystyle=$	$\displaystyle E\left[\gamma^{X_{k+1}}R(T_{k+1})+\int\gamma^{x^{\prime}}\sum_{a\in\mathbb{A}}\pi(a|\bm{S}_{k+1},x^{\prime})Q^{\pi}(\bm{S}_{k+1},x^{\prime},a)d\mathcal{P}_{X}(x^{\prime};\bm{S}_{k+1},\bm{S}_{k},X_{k},A_{k})\big{|}\bm{S}_{k},X_{k},A_{k}\right]\,.$

Let $\widehat{\mathcal{P}}_{X}$ denote an estimate of the distribution function $\mathcal{P}_{X}$. Similar to (4), we construct $\widetilde{\bm{U}}(\bm{\theta})=\frac{1}{\sum_{i}K_{i}}\sum_{i=1}^{n}\sum_{k=0}^{K_{i}-1}\bm{\dot{Q}}_{\bm{\theta}}(\bm{S}_{i,k},X_{i,k},A_{i,k};\bm{\theta})\widetilde{\delta}_{i,k+1}(\bm{\theta})\,,$ with $\widetilde{\delta}_{i,k+1}(\bm{\theta})=\gamma^{X_{i,k+1}}R_{i}(T_{i,k+1})+\int\gamma^{x^{\prime}}\sum_{a\in\mathbb{A}}\pi(a|\bm{S}_{i,k+1},x^{\prime})Q(\bm{S}_{i,k+1},x^{\prime},a;\bm{\theta})\mathrm{d}\widehat{\mathcal{P}}_{X}(x^{\prime};\bm{S}_{i,k+1},\bm{S}_{i,k},X_{i,k},A_{i,k})-Q(\bm{S}_{i,k},X_{i,k},A_{i,k};\bm{\theta})\,.$ Further approximating action-value functions with B-splines $Q(\bm{s},x,a;\bm{\theta})=\bm{\xi}(\bm{s},x,a)^{\top}\bm{\theta}$, then the solution to $\widetilde{\bm{U}}(\bm{\theta})=0$ takes the form of

$\displaystyle\widetilde{\bm{\theta}}_{\pi}=\biggl{\{}\frac{1}{\sum_{i}K_{i}}\sum_{i=1}^{n}\sum_{k=0}^{K_{i}-1}\bm{\xi}_{i,k}(\bm{\xi}_{i,k}-\widetilde{\mathcal{U}}_{\pi,i,k+1})^{\top}\biggr{\}}^{-1}\bigg{\{}\frac{1}{\sum_{i}K_{i}}\sum_{i=1}^{n}\sum_{k=0}^{K_{i}-1}\bm{\xi}_{i,k}\gamma^{X_{i,k+1}}R_{i}(T_{i,k+1})\bigg{\}}$

(5)

with $\widetilde{\mathcal{U}}_{\pi,i,k+1}=\int\gamma^{x^{\prime}}\sum_{a\in\mathbb{A}}\pi(a|\bm{S}_{i,k+1},x^{\prime})\bm{\xi}(\bm{S}_{i,k+1},x^{\prime},a)d\widehat{\mathcal{P}}_{X}(x^{\prime};\bm{S}_{i,k+1},\bm{S}_{i,k},X_{i,k},A_{i,k})\,.$

Then, the estimator for $Q^{\pi}(\bm{s},x,a)$ is $\widetilde{Q}^{\pi}(\bm{s},x,a)=\bm{\xi}(\bm{s},x,a)^{\top}\widetilde{\bm{\theta}}_{\pi}$, the estimator for state-value function $V^{\pi}(\bm{s},x)$ is $\widetilde{V}^{\pi}(\bm{s},x)=\sum_{a\in\mathbb{A}}\pi(a|\bm{s},x)\bm{\xi}(\bm{s},x,a)^{\top}\widetilde{\bm{\theta}}_{\pi}$ and the estimator for the value function $V^{\pi}(\mathcal{G})$ under given reference distribution $\mathcal{G}$ is $\widetilde{V}^{\pi}(\mathcal{G})=\bm{\zeta}_{\pi,\mathcal{G}}^{\top}\widetilde{\bm{\theta}}_{\pi}$. Here, extra effort is required in estimating $\mathcal{P}_{X}$, the distribution of the next gap time $X_{k+1}$ conditional on $\bm{S}_{k+1},\bm{S}_{k},X_{k}$ and $A_{k}$. As an illustration, we propose a Cox-type counting process model that satisfies (S1-$x^{\prime}$), and takes the form

$\displaystyle E[\mathrm{d}N(t)|\mathcal{F}(t-)]=\lambda(t-T_{k};\bm{S}_{k+1},\bm{S}_{k},X_{k},A_{k})\ \mathrm{d}t=\lambda_{0}(t-T_{k})\exp(\bm{\beta}_{0}^{\top}\bm{Z}_{k})\ \mathrm{d}t\,,$

(6)

Here, $k=\max\{j;T_{j}<t\}$, $\bm{Z}_{k}=\bm{\phi}(\bm{S}_{k+1},\bm{S}_{k},X_{k},A_{k})\in\mathbb{R}^{q}$ for some pre-specified $q$-dimensional function $\bm{\phi}$, $\bm{\beta}_{0}$ is a vector of unknown parameters and $\lambda_{0}$ is an unspecified baseline intensity function. Let $\{N_{i}(t)\}_{1\leq i\leq n}$ denote $n$ independent and identically distributed trajectories of $N(t)$. For $i=1,\ldots,n$ and $k=1,\ldots,K_{i}$, define $N_{i,k}(x)=N_{i}(T_{i,k-1}+x)-N_{i}(T_{i,k-1})$, $\bm{Z}_{i,k}=\bm{\phi}(\bm{S}_{i,k+1},\bm{S}_{i,k},X_{i,k},A_{i,k})$. Then, we propose to estimate $\bm{\beta}_{0}$ by solving the equation

$\displaystyle\bm{U}_{z}(\beta)=\frac{1}{\sum_{i}K_{i}}\sum_{i=1}^{n}\sum_{k=0}^{K_{i}-1}\int_{0}^{\tau}\{\bm{Z}_{i,k}-\overline{\bm{Z}}(x;\bm{\beta})\}\ \mathrm{d}N_{i,k+1}(x)=0\,,$

(7)

where $\overline{\bm{Z}}(x;\bm{\beta})={\{\sum_{i,k}I(X_{i,k+1}\geq x)\exp(\bm{\beta}^{\top}\bm{Z}_{i,k})\bm{Z}_{i,k}\}^{-1}}{\{\sum_{i,k}I(X_{i,k+1}\geq x)\exp(\bm{\beta}^{\top}\bm{Z}_{i,k})\}}\,,$ and $\tau$ is a positive constant satisfying $P(X_{k}\geq\tau)>0$ for each $k$. Let $\widehat{\bm{\beta}}$ denote the obtained estimator. Then, we can further estimate $\Lambda_{0}(x)=\int_{0}^{x}\lambda_{0}(u)du$ by $\widehat{\Lambda}_{0}(x)=\int_{0}^{x}\{\sum_{i}\sum_{k}I(X_{i,k+1}\geq u)\exp(\widehat{\bm{\beta}}^{\top}\bm{Z}_{i,k})\}^{-1}\{\sum_{i}\sum_{k}\mathrm{d}N_{i,k+1}(u)\}$, and can estimate $\mathcal{P}_{X}$ by $\widehat{\mathcal{P}}_{X}(x^{\prime};\bm{s^{\prime}},\bm{s},x,a)=1-\exp\{-\widehat{\Lambda}_{0}(x^{\prime})\exp(\widehat{\bm{\beta}}^{\top}\bm{z})\}$ for $\bm{z}=\bm{\phi}(\bm{s^{\prime}},\bm{s},x,a)$.

We first introduce the notations and impose necessary assumptions. Denote $r(\bm{s},x,a)=E[\gamma^{X_{k+1}}R(T_{k+1})|\bm{S}_{k}=s,X_{k}=x,A_{k}=a]$. Let $p(\bm{s^{\prime}},x^{\prime};\bm{s},x,a)$ be the conditional density function satisfying $\mathcal{P}(\mathcal{B}|\bm{s},x,a)=\int_{\mathcal{B}}p(\bm{s^{\prime}},x^{\prime};\bm{s},x,a)\mathrm{d}\bm{s^{\prime}}\mathrm{d}x^{\prime}$ for the joint distribution $\mathcal{P}(\cdot|\bm{s},x,a)$ defined in (A1). Then, (A3) is a regularity condition ensuring the action function $Q^{\pi}(\bm{s},x,a)$ to be continuous (Shi et al., 2021, Lemma 1).

• (A3)

  There exists some $p_{0}>0$, such that $r(\cdot,\cdot,a)$ and $p(\cdot,\cdot;\bm{s},x,a)$ are Hölder continuous functions with exponent $p_{0}$ and with respect to $(\bm{s},x)\in\mathbb{S}\times\mathbb{R}^{+}$ .

As discussed in Section 2.2, the transition trajectories of states and gap times $\{\widetilde{S}_{k}=(\bm{S}_{k}^{\top},X_{k})^{\top}\}_{k\geq 0}$ form a time-homogeneous Markov chain under assumption (A1). Suppose $\{(\bm{S}_{k}^{\top},X_{k})^{\top}\}_{k\geq 0}$ has a unique invariant joint distribution with density function $\mu(\cdot,\cdot)$ under some behavior policy $b$. Similarly, when the reference distribution $\mathcal{G}$ is absolutely continuous with respect to the Lebesgue measure, there exists a density function $v_{0}(\bm{s},x)$ satisfying $\mathrm{d}\mathcal{G}(\bm{s},x)=v_{0}(\bm{s},x)\mathrm{d}s\mathrm{d}x$, which is also a joint distribution for the states and gap time. We posit the following assumptions on $\mu$ and $v_{0}$.

• (A4)

  $\mu(\cdot,\cdot)$ and $v_{0}(\cdot,\cdot)$ are uniformly bounded away from 0 and $\infty$.