Bayesian inference for the Markov-modulated Poisson process with an outcome process

Abstract

In medical research, understanding changes in outcome measurements is crucial for inferring shifts in health conditions. However, traditional methods often struggle with large, irregularly longitudinal data and fail to account for the tendency of individuals in poorer health to interact more frequently with the healthcare system. Additionally, clinical data can lack information on terminating events like death. To address these challenges, we start from the continuous-time hidden Markov model which models observed data as outcomes influenced by latent health states. Our extension incorporates a point process to account for the impact of health states on observation timings and includes a ‘death’ state to model unobserved terminating events through a Poisson process, where transition rates depend on the latent health state. This approach captures both the severity of the disease and the timing of healthcare interactions. We present an exact Gibbs sampler procedure that alternates between sampling the latent health state paths and the model parameters. By including the ‘death’ state, we mitigate biases in parameter estimation that would arise from solely modelling ‘live’ health states. Simulation studies demonstrate that the proposed Gibbs sampler performs effectively. We apply our method to Canadian healthcare data, offering valuable insights for healthcare management.

claim data, Gibbs sampler, hidden Markov models, longitudinal analysis, Markov-modulated Poisson process

1 Introduction

In medical research, observed changes in outcome measurement are used to infer changes in the underlying condition of the patient. Data must be collected over time (longitudinally) on a patient to characterize the effect of risk factors and to evaluate the effectiveness of treatments. For instance, understanding the progression of chronic diseases, such as chronic obstructive pulmonary disease (COPD), is important to inform early diagnosis, personalized care, and health system management. With no existing cure, the severity of COPD increases over time, producing more frequent exacerbations, and increasing use of healthcare (Burge and Wedzicha, 2003). This worsening disease trajectory defines COPD and many other chronic illnesses. Despite this natural history of worsening disease, health system indicators currently describe COPD from a cross-sectional perspective. These traditional measures of service volume, such as hospitalization rates, neglect the temporal disease progression within individuals that is central to managing COPD. Data from clinical and administrative systems have the potential to advance this understanding, but traditional methods for modelling disease progression are not well-suited to analysing high-dimensional data collected at irregular intervals. Many cases of severe illness can be traced back to a lack of appropriate and ongoing follow-up health care. To help understand the long-term causes of severe and preventable illness, we need to understand not only how many patients have severe illness but also what their ‘care paths’ or ‘trajectories’ were before they ended up with severe illness. Private and public healthcare systems now have the ability to maintain longitudinal digital health records. These data from sources, such as electronic health records, health systems and data from mobile health applications, allow the study of subject-specific trajectories, which can inform clinical and public health decision-making. However, data are typically recorded only when a subject interacts with a provider, resulting in irregularly spaced longitudinal observations, and the patterns of clinical interactions vary from subject to subject. Although the data allow dynamic monitoring of the underlying disease progression that governs the data recorded in the health system, this disease progression cannot be observed directly, and so inferential methods are needed to ascertain the latent progression. The statistical analysis of such data, especially when the cohort is large and observations are frequent, is complex.

A multi-state model captures the status of an individual longitudinally as a time-discretization of a continuous-time Markov process, and in many settings the common assumption that the health status of a patient falls into one of a finite number of categories (usually ordered by severity) is a plausible one. Inference for continuous-time Markov processes observed discretely and either precisely or partially has been well studied (see, for example, Bladt and Sørensen, 2005; Pfeuffer, 2017). However, health system records typically do not include direct information about underlying health status; instead, it is inferred indirectly through a proxy measure. Such a process can be modelled as a continuous-time hidden Markov model (CTHMM, e.g. Cappé et al., 2005), which has been studied recently in medical applications such as disease progression (Luo et al., 2021) and cancer screening (Klausch et al., 2023; Lange et al., 2015). While a CTHMM can address the issue of non-equidistant longitudinal observations, in medical research, it is possible that the timing of healthcare utilization itself is informative in the assessment of health status progression. For example, patients who are in more severe health states typically interact with the healthcare system more frequently. A common approach to modelling informative observation times is to use a Poisson process (e.g. Cook and Lawless, 2018). However, if the intensity of the Poisson process depends on the state of a hidden Markov process, it is often referred to as a Markov-modulated Poisson process (MMPP) (e.g. Fearnhead and Sherlock, 2006; Rydén, 1996). There have been studies modelling this observation time process using frequentist approaches, such as Lange et al. (2015) and Mews et al. (2023).

In addition, data from the claim and healthcare system provide no information for terminating events, such as death. However, appropriate parameter inference will only follow if we distinguish a lack of interaction with the healthcare system due to, for example, mild illness from a lack of interaction due to the subject being dead. Nevala et al. (2024) included the death data from a different data system to model colorectal cancer progression. However, to the best of our knowledge, there is no existing method to account for the terminating event, such as death, when it is not available in a dataset and event times are informative. In this article, we develop a Bayesian inference procedure that takes into account the observation process by modelling the observed data as an outcome process and the time of observations as a point process jointly alongside the CTHMM, given the underlying latent health state. Moreover, we also consider the death time as missing data, and model the death as an informative censoring. This extension allows us to model disease severity and death not only from the types of care received but also from the times between different observed events and their respective frequency. It also provides insights on how long patients tend to survive in each state.

The article is organized as follows. Section 2 demonstrate the model structure, followed by the likelihood construction in Section 2.1 and the forward–backward method to calculate the probability of each state at the observation time in Section 3. Section 4 layouts the Gibbs sampler to draw inference for the model, and some empirical studies, including simulation and real data analysis, are in Section 5. The article concludes in Section 6 with a summary of findings and discusses areas for future research.

2 A MMPP with an outcome process

In this section, we describe the continuous-time hidden Markov model with the observation time following a Poisson process. In the sequel, for any collection of scalars, $a_1,\dots ,a_T$ or $a^1,\dots ,a^T$⁠, we denote the concatenations by $a_{1:T}$ and $a^{1:T}$⁠, respectively. Thus, for a particular individual we observe the data $o_{1:T}:=\{o_1,\dots ,o_T\}$ at time points ${\tau }_{1:T}:=\{{\tau }_1,\dots ,{\tau }_T\}$⁠. We assume that the data arise as a consequence of a latent continuous-time Markov chain (CTMC) $\{X_s{\}}_{s\in [0,\tau ]}$ taking values on the finite integer set $\mathcal{K}:=\{1,2,\dots ,K\}$⁠. We denote the generator of the latent process by $\mathbf{Q}$ and the initial distribution by ν. At observation time ${\tau }_t$ (⁠$t=1,\dots ,T$⁠), the latent process is in state $X_{{\tau }_t}$⁠, and the observation is $o_t$⁠.

Conditional on the latent process being in state $k\in \mathcal{K}$⁠, at time ${\tau }_t$⁠, the outcome observation $o_t$ is drawn from a distribution with a probability mass or density function $f_t(o_t|X_{{\tau }_t}=k)\equiv$  $f_k(o_t|z_t,{\beta }_k)$⁠, where $z_t\in {\mathbb{R}}^D$ is a vector of covariates (including an intercept) and β is a D-vector of covariate effects. In all of our simulations, including the real-data example, the prior for ${\beta }_k$⁠, is chosen to be conditionally conjugate to $f_k$ so that it is straightforward to sample from the conditional posterior for each ${\beta }_k$ given $X_{{\tau }_1},\dots ,X_{{\tau }_T}$⁠. It is possible to extend inference to the non-conjugate case using the Metropolis–Hastings algorithm; see the online supplementary material for a brief discussion. Define the matrix $\mathbf{B}=({\beta }_{d,k})$ for $d=1,\dots ,D$ and $k=1,\dots ,K$ as the regression coefficient matrix.

In addition, we allow the observation process itself (that is, the times that the observations are made) to be informative about the latent trajectory. Specifically, let $N_s$ be the number of observations (events of the Poisson process) on $[0,s]$⁠, $N_s=\sum _{t=1}^{T}1({\tau }_t\le s)$⁠, so that $N_0=0$ and $N_{{\tau }_t}=t$⁠. We assume that $N_s$ follows a Poisson process whose intensity is ${\lambda }_k$ when $X_s=k$⁠. The joint process $(N_s,X_s)$ is a MMPP, and inference approaches have been extensively studied in, for example, Rydén (1996) and Fearnhead and Sherlock (2006). To simplify notation, we write $X_{[0,\tau ]}$ for $\{X_s{\}}_{s\in [0,\tau ]}$ and $N_{[0,\tau ]}$ for $\{N_s{\}}_{s\in [0,\tau ]}$⁠.

Conditional on $X_s$⁠, the observations $O_{1:T}$ are assumed to be independent of each other and of $N_s$⁠; also $N_s$ is assumed to evolve independently of $X_s$ except that its intensity is ${\lambda }_{X_s}$⁠. Therefore, the observation time ${\tau }_t$ is generated from the MMPP $(N_s,X_s)$ and the outcome $O_t$ is then generated from $f_{X_{{\tau }_t}}(o_t|z_t,{\beta }_{X_{{\tau }_t}})$⁠. Figure 1 provides a schematic of the presumed data generating structure for one subject. Conditional on ν, $\mathbf{Q}$⁠, $\mathit{\lambda }:=({\lambda }_1,\dots ,{\lambda }_k)$ and $\mathbf{B}$⁠, data for each subject are assumed to be generated independently of the data for all other subjects. Ideally we observe the process over some time interval $[0,\tau ]$ for some $\tau >{\tau }_T$ and with ${\tau }_1>0$⁠. However, for each patient, our data start with a diagnosis of COPD, so ${\tau }_1=0$⁠.

Figure 1.

Schematic of the presumed data generating mechanism for the Markov-modulated Poisson process with an outcome process $O_{{\tau }_t}$⁠. $X_s$ represents the process underlying the observed data, with observation (event) time points denoted ${\tau }_t,t=1,\dots ,T$ for a state-dependent Poisson process $N_s$⁠; $X_s$ represents the hidden Markov process.

2.1 Likelihood for the MMPP with an outcome process

We are interested in inference for the data generated from the MMPP with an outcome process. We now construct the complete-data likelihood for one subject over the time interval $[0,\tau ]$⁠. Since subjects are independent, the likelihood for all subjects is the product of the individual likelihoods.

The likelihood function for $\mathbf{Q}$ and ν is (Bladt and Sørensen, 2005)

$$\mathcal{L}(\mathbf{Q},\nu |x_{[0,\tau ]})=f(x_{[0,\tau ]}|\mathbf{Q},\nu )={\nu }_{x_0}{\displaystyle \prod _{l=1}^K}{\displaystyle \prod _{m\ne l}}{q_{l,m}}^{N_{l,m}}\exp \left(-q_{l,m}{R}_l\right),$$

where $N_{l,m}$ is the number of transitions from state l to state m in the time interval $[0,\tau ]$⁠, and $R_l$ is the total time that the process has spent in state l in $[0,\tau ]$⁠. The values of $q_{l,m}$ are the elements in $\mathbf{Q}$⁠. In addition, the likelihood function for $\mathit{\lambda }$ is

$$\mathcal{L}(\mathit{\lambda };x_{[0,\tau ]},{\tau }_{1:T})=f({\tau }_{1:T}|x_{[0,\tau ]},\mathit{\lambda })={\displaystyle \prod _{i=1}^K}{{\lambda }_i}^{N_i}\exp (-{\lambda }_i{R}_i),$$

where $N_i$ is the number of events at state i. Note that the quantities $N_{l,m}$⁠, $R_l$⁠, and $N_i$ all are unobserved, but they can be computed when the underlying process $X_{[0,\tau ]}$ is known.

The complete likelihood, derived from $o_{1:T}$⁠, ${\tau }_{1:T}$ and $X_{[0,\tau ]}$ can be factorized:

$$\mathcal{L}(\mathbf{B},\mathit{\lambda },\mathbf{Q},\nu )\equiv f(o_{1:T},{\tau }_{1:T},x_{[0,\tau ]}|\mathbf{B},\mathit{\lambda },\mathbf{Q},\nu )=f(o_{1:T}|\mathbf{B},x_{[0,\tau ]})f({\tau }_{1:T}|x_{[0,\tau ]},\mathit{\lambda })f(x_{[0,\tau ]}|\nu ,\mathbf{Q}).$$

Thus,

$$\mathcal{L}(\mathbf{B},\mathit{\lambda },Q,\nu )={\displaystyle \prod _{t=1}^T}{f}_{X_{{\tau }_t}}\left(o_t|z_t,{\beta }_{x_{{\tau }_t}}\right)\times {\nu }_{x_0}{\displaystyle \prod _{l=1}^K}{\displaystyle \prod _{m\ne l}}{q_{l,m}}^{N_{l,m}}\exp \left(-q_{l,m}{R}_l\right){\displaystyle \prod _{i=1}^K}{{\lambda }_i}^{N_i}\exp (-{\lambda }_i{R}_i).$$

(1)

If independent Gamma-distributed priors are assigned to all of the components ${\lambda }_i$ and $q_{l,m}$ (see Section 4.1) then, conditional on the complete chains for all individuals, the posteriors for these parameters also have independent Gamma distributions. This permits the use of the Gibbs steps described in Section 4.2. An alternative partial imputation, which necessitates Metropolis–Hastings updates for these parameters, is discussed in the online supplementary material.

2.2 Data-induced window restriction

Our data are of interactions between patients with COPD and the health service during a 17-year window between 1st January 1998 and 31st December 2014. If a patient was diagnosed with COPD within the window then the first observation for that patient is the time of diagnosis. The patient may have interacted with the health service within the 17-year window and before they were diagnosed with COPD, but we do not know about these interactions. Therefore, for any given patient, we set ${\tau }_1=0$⁠; we start the clock at the first interaction. Since we know there will always be an interaction at time 0, this observation does not contribute to the count $N_k(\tau )$ in (1). Of necessity, therefore, the initial distribution ν represents an average of the distribution of the state when an individual is first diagnosed with COPD and, for those already had with COPD at the start of 1998, the distribution of the state at their first interaction within the window.

For each patient, we do not have the time of death or, even, an indicator of whether or not death occurred before the end of 2014. If the last interaction with the health service was just before the end of the window then it seems likely the patient was still alive at the end of the window; however, if the last interaction was well before the end of the window and there had been a high frequency of interactions before this then it is likely that the patient died or moved away from the province shortly after the final interaction. In Section 5, for the data analysis and one of the simulation studies, we circumvent this issue by introducing an additional state for the hidden chain, corresponding to death (or moving away). This is an absorbing state and there are no observations whilst in this state.

Figure 2 shows three example patients’ health trajectories. The observation window is between 0 and τ, and assume that all patients have an observation at time 0. Patient 1 died right before τ, and we have the record for the last observation but do not have information for the death. Patient 2, on the other hand, has one additional interaction after the observation window and died afterwards. For these patients, the death happens close to the observation window. Modelling death does not have much impact on estimates of $\mathbf{Q}$ or $\mathit{\lambda }$⁠. However, Patient 3 has an early death, and the trajectory for observed data is only half of the observation window. For these patients who die early, in particular, we need to model death instead of, for example, assuming the patient is alive until the end of the window.

Figure 2.

Example of patients’ health trajectories. Circle represents the observation, and cross represents the death (missing). The end of the observation window is at time τ.

3 The forward–backward algorithm

Before we discuss inference for the model, we first derive the specifics of the general forward–backward algorithm (Baum and Petrie, 1966) during the window $[0,\tau ]$ which enable us to simulate the unobserved Markov process $X_{[0,\tau ]}$⁠.

In this section, it will be clearer to denote the end of the observation window by ${\tau }_e$ rather than τ and to define $X_{{\tau }_{T+1}}\equiv X_{{\tau }_e}$⁠. Note that ${\tau }_{1:T}$ are random observation times but ${\tau }_e$ is the fixed end of the window. Just as the event ${\tau }_{1:T}$ is equivalent to having observation-process events at times ${\tau }_{1:T}$ and at no other times in $[{\tau }_1,{\tau }_T]$⁠, the event ${\tau }_{1:T},{\tau }_e$ implies, additionally, that there are no observation process events in the interval $({\tau }_T,{\tau }_e]$⁠. Recall, also, that our data format forces us to set ${\tau }_1=0$⁠, so that the fact that an event happens at time 0 does not form part of the data; ${\tau }_1$ now has the same purpose as ${\tau }_e$⁠, defining the boundary of the observation window. However, since there was an event at time ${\tau }_1$ there is an observation, $o_1$⁠, which does need to be accounted for.

We use the notation $j:k$ to indicate all indices between and including j and k; for example $o_{j:k}\equiv (o_j,o_{j+1},\dots ,o_k)$⁠. Thus, we first define the forward variable,

$${\alpha }_{t,k}=\mathbb{P}(o_{1:t},{\tau }_{1:t},X_{{\tau }_t}=k),t=1,\dots ,T,\text{and}{\alpha }_{e,k}=\mathbb{P}(o_{1:T},{\tau }_{1:T},{\tau }_e,X_{{\tau }_e}=k).$$

In our case, ${\tau }_1=0$ is guaranteed, so ${\alpha }_{1,k}={\nu }_{k}f(o_1|X_{{\tau }_1}=k)$⁠. Then for $t=2,\dots ,T$⁠,

$$\begin{array}{rl}{\alpha }_{t,k}& =\mathbb{P}\left(o_t|X_{{\tau }_t}=k,o_{1:t-1},{\tau }_{1:t}\right)\mathbb{P}(o_{1:t-1},{\tau }_{1:t},X_{{\tau }_t}=k)\\ & =f_t\left(o_t|X_{{\tau }_t}=k\right)\times {\displaystyle \sum _{i=1}^K}\mathbb{P}(o_{1:t-1},{\tau }_{1:t-1},{\tau }_t,X_{{\tau }_t}=k,X_{{\tau }_{t-1}}=i)\\ & =f_t\left(o_t|X_{{\tau }_t}=k\right)\times {\displaystyle \sum _{i=1}^K}\mathbb{P}(o_{1:t-1},{\tau }_{1:t-1},X_{{\tau }_{t-1}}=i)\mathbb{P}({\tau }_t,X_{{\tau }_t}=k|X_{{\tau }_{t-1}}=i)\\ & =f_t\left(o_t|X_{{\tau }_t}=k\right)\times {\displaystyle \sum _{i=1}^K}{\alpha }_{t-1,i}\mathbb{P}({\tau }_t,X_{{\tau }_t}=k|X_{{\tau }_{t-1}}=i).\end{array}$$

Further, ${\alpha }_{e,k}=\sum _{i=1}^K{\alpha }_{T,i}\mathbb{P}({\tau }_e,X_{{\tau }_e}=k|X_{{\tau }_T}=i)$⁠.

To enact the above recursion we must calculate the quantities, $\mathbb{P}({\tau }_t,X_{{\tau }_t}=k|X_{{\tau }_{t-1}}=i)$⁠, $t=2,\dots ,T$⁠, and $\mathbb{P}({\tau }_e,X_{{\tau }_e}=k|X_{{\tau }_T}=i)$⁠. For $t=2,\dots ,T$⁠, this is the probability of the next event after time ${\tau }_{t-1}$ being at time ${\tau }_t$ and the chain being in state k at this time, all given that the chain was in state i at time ${\tau }_{t-1}$⁠. The final quantity differs from this in that there has been no event from just after ${\tau }_T$⁠, up to and including ${\tau }_e$⁠.

Consider a new process $\{Y_t{\}}_{t\ge {\tau }_{t-1}}$ on $\{1,\dots ,K,✓\}$⁠. The absorbing state $✓$ occurs as soon as there has been at a new event in the Poisson process $N_s$⁠. Specifically, at time ${\tau }_{t-1}$ the process is in state i; i.e. $Y_{{\tau }_{t-1}}=X_{{\tau }_{t-1}}$⁠. Up until the next event in the Poisson process, $N_s$⁠, we continue to have $Y_s=X_s$⁠, but from the time of the next Poisson event onwards, $Y_s=✓$⁠. The process $\{Y_s{\}}_{s\ge {\tau }_{t-1}}$ is a CTMC (Mark and Ephraim, 2013), and its generator is

$$\left(\begin{array}{cc}Q-\mathbf{\Lambda }& {\mathit{\lambda }}^{\mathrm{\top }}\\ {\mathbf{0}}^{\mathrm{\top }}& 0\end{array}\right),$$

where $\mathbf{\Lambda }=\text{diag}({\lambda }_1,\dots ,{\lambda }_K)$⁠. Writing $Y_{{\tau }^{-}}$ for $\underset{s\uparrow \tau }{lim}{Y}_s$⁠, and defining ${\eta }_{t-1,t,i,k}:=$  $[\exp \{(Q-\mathbf{\Lambda })({\tau }_t-{\tau }_{t-1})\}{]}_{i,k}$⁠, we see that

$$\mathbb{P}({\tau }_t,X_{{\tau }_t}=k|X_{{\tau }_{t-1}}=i)=\mathbb{P}(Y_{{\tau }_t^{-}}=k|Y_{{\tau }_{t-1}}=i){\lambda }_k={\left[\exp \{(Q-\mathbf{\Lambda })({\tau }_t-{\tau }_{t-1})\}\right]}_{i,k}{\lambda }_k={\eta }_{t-1,t,i,k}{\lambda }_k.$$

Similarly,

$$\mathbb{P}(X_{{\tau }_e}=k,{\tau }_e|X_{{\tau }_T}=i)={\left[\exp \{(Q-\mathbf{\Lambda })({\tau }_e-{\tau }_T)\}\right]}_{i,k}=:{\eta }_{T,e,i,k}.$$

Therefore, the recursions for the forward variable are

$${\alpha }_{t,k}=f_t\left(o_t|X_{{\tau }_t}=k\right)\times {\displaystyle \sum _{i=1}^K}{\alpha }_{t-1,i}{\eta }_{t-1,t,i,k}{\lambda }_k(t=2,\dots ,T),$$

with ${\alpha }_{e,k}=\sum _{i=1}^K{\alpha }_{T,i}{\eta }_{T,e,i,k}$⁠.

From its definition, $a_{e,k}\propto \mathbb{P}(X_{{\tau }_e}=k|o_{1:T},{\tau }_{1:T},{\tau }_e)$ so we may sample a value for $x_{{\tau }_e}$ with the probabilities of the individual states proportional to the elements of $a_e$⁠. Suppose that we have simulated $x_{{\tau }_e},x_{{\tau }_T},\dots ,x_{{\tau }_{t+1}}$ from their joint distribution given $o_{1:T}$⁠, ${\tau }_{1:T}$ and ${\tau }_e$⁠. We now describe how to simulate $x_{{\tau }_t}$ backward from its conditional distribution given $x_{{\tau }_e},x_{{\tau }_T},\dots ,x_{{\tau }_{t+1}}$⁠, $o_{1:T}$⁠, ${\tau }_{1:T}$ and ${\tau }_e$⁠. From the hidden Markov structure, conditional on $x_{{\tau }_{t+1}}$⁠, all random variables with time indices before $t+1$ are independent of those with indices after $t+1$⁠, so

$$p_{t,j,k}:=\mathbb{P}(X_{{\tau }_t}=j|x_{{\tau }_e},x_{{\tau }_T},\dots ,x_{{\tau }_{t+2}},X_{{\tau }_{t+1}}=k,o_{1:T},{\tau }_{1:T},{\tau }_e)=\mathbb{P}(X_{{\tau }_t}=j|X_{{\tau }_{t+1}}=k,o_{1:t},{\tau }_{1:t+1}),$$

where ${\tau }_{T+1}$ is understood to be ${\tau }_e$⁠. Hence

$$p_{t,j,k}\propto b_{t,j,k}:=\mathbb{P}(o_{1:t},{\tau }_{1:t+1},X_{{\tau }_t}=j,X_{{\tau }_{t+1}}=k)=\mathbb{P}(o_{1:t},{\tau }_{1:t},X_{{\tau }_t}=j)\mathbb{P}({\tau }_{t+1},X_{{\tau }_{t+1}}=k|X_{{\tau }_t}=j),$$

again using the conditional independence structure. Thus, for $t=1,\dots ,T-1$⁠,

$$b_{t,j,k}={\alpha }_{t,j}{\eta }_{t,t+1,j,k}{\lambda }_k\propto {\alpha }_{t,j}{\eta }_{t,t+1,j,k},$$

since k is known. For $t=T$⁠, the constant ${\lambda }_k$ term does not appear in the first place, leading to the same result. Finally,

$$\mathbb{P}(X_{{\tau }_t}=j|x_{{\tau }_e},x_{{\tau }_T},\dots ,x_{{\tau }_{t+2}},X_{{\tau }_{t+1}}=k,o_{1:T},{\tau }_{1:T},{\tau }_e)=\frac{b_{t,j,k}}{{\displaystyle \sum _{i=1}^K}{b}_{t,i,k}}.$$

The above recursions enable us to simulate $X_{{\tau }_t}$⁠, $t=e,T,\dots ,1$⁠, conditional on $\mathbf{B},\nu ,Q,\mathit{\lambda }$⁠, $o_{1:T}$ and ${\tau }_{1:T}$ and (no further observations up to) ${\tau }_e$⁠. To simulate from $X_{[0,\tau ]}$⁠, we also need to fill in the gaps between the observation times. The conditional dependence structure of the model means that given $X_{{\tau }_{t-1}}$ and $X_{{\tau }_t}$⁠, $X_{[{\tau }_{t-1},{\tau }_t]}$ is independent of $O_{1:T}$ and of $X_{{\tau }_1},\dots ,X_{{\tau }_{t-2}}$ and $X_{{\tau }_{t+1}},\dots ,X_{{\tau }_T}$⁠. Furthermore, given ${\tau }_{t-1}$ and ${\tau }_t$⁠, $X_{[{\tau }_{t-1},{\tau }_t]}$ is independent of ${\tau }_{1:T}$⁠. However, the fact that two neighbouring event times, ${\tau }_{t-1}$ and ${\tau }_t$ bracket the time interval of interest means that there are no $N_s$ events in $({\tau }_{t-1},{\tau }_t)$⁠.

Combining the above information, to simulate from $X_{({\tau }_{t-1},{\tau }_t)}$ given $o_{1:T},{\tau }_{1:T}$ and $X_{{\tau }_1},\dots ,X_{{\tau }_T}$ with $X_{{\tau }_{t-1}}=j$ and $X_{{\tau }_t}=k$⁠, we must simulate from $Y_{({\tau }_{t-1},{\tau }_t)}$ given $Y_{{\tau }_{t-1}}=j$ and $Y_{{\tau }_t}=k$⁠. We employ the uniformization method of Fearnhead and Sherlock (2006) and Rao and Teh (2013).

4 Bayesian inference for the MMPP with an outcome process

We are interested in Bayesian analysis of the MMPP with an outcome process. Bayesian inference for CTHMMs has been explored under various likelihood formulation, such as simulation-based (Luo et al., 2021) and approximation (Williams et al., 2020) methods. To the best of our knowledge, it has never been implemented for an MMPP model with an outcome process under a full Bayesian paradigm. We outline the steps for a Metropolis-within-Gibbs Markov chain Monte Carlo (MCMC) algorithm to draw posterior samples given suitable prior distributions for the parameters. Our scheme uses the complete data log-likelihood in (1) in the spirit of Luo et al. (2021) and Fearnhead and Sherlock (2006). Throughout, we assume that the number of states, K, is fixed and known.

4.1 Prior distributions

As stated in Section 2, the model for the outcome process when $X_{{\tau }_t}=k$⁠, $f_k(o_t|z_t,{\beta }_k)$ is such that a family of conditionally conjugate priors exists; we choose the prior for each ${\beta }_k$ from this family, and set ${\pi }_0(\mathbf{B})=\prod _{k=1}^K{\pi }_{0,k}({\beta }_k)$⁠. The initial distribution, ν, is given a Dirichlet prior, so that, conditional on $X_0=x_0$⁠, its posterior is also Dirichlet.

For $\mathit{\lambda }$ and $\mathbf{Q}$⁠, we choose

$${\pi }_0(\mathit{\lambda })={\displaystyle \prod _{k=1}^K}\mathsf{G}\mathsf{a}\mathsf{m}({\lambda }_k;a_k^{\lambda },b_k^{\lambda })\text{and}{\pi }_0(Q)={\displaystyle \prod _{j=1}^K}{\displaystyle \prod _{k\ne j}}\mathsf{G}\mathsf{a}\mathsf{m}(q_{j,k};a_{j,k}^q,b_{j,k}^q),$$

where $\mathsf{G}\mathsf{a}\mathsf{m}(x;a,b)$ is the density of a $\mathsf{G}\mathsf{a}\mathsf{m}\mathsf{m}\mathsf{a}(a,b)$ random variable, evaluated at x. From (1), conditional on $X_{[0,\tau ]}$ and ${\tau }_{1:T}$⁠, the posterior distributions for each ${\lambda }_k$ and $q_{j,k}$ (⁠$j\ne k$⁠) are mutually independent and all have Gamma distributions.

For $\mathit{\lambda }$ and $\mathbf{Q}$⁠, it is tempting to represent the absence of knowledge about each parameter through a prior with a low shape parameter, $a<1$⁠. However, any model with K states includes a model with $K-1$ states, with a particular state, i, never visited. In this case, the data, $o_{1:T}$ and ${\tau }_{1:T}$⁠, provide no further information on ${\lambda }_i$ and $q_{i,k}$ (⁠$k\ne i$⁠). Thus, the posteriors for all $\mathit{\lambda }$ and $\mathbf{Q}$ parameters are a mixture including at least one component that is simply the prior. If the prior shape parameter for that parameter is below 1, the prior density rises to $\mathrm{\infty }$ as the parameter tends to 0, and, because the posterior is a mixture density that includes the prior, the posterior density also rises to $\mathrm{\infty }$ as each parameter $\downarrow 0$⁠, which usually does not represent our prior belief and, moreover, can lead to poor mixing of the MCMC chain.

4.2 Exact Gibbs sampler

Each step of our Gibbs sampler simulates a parameter or parameter vector, or an aspect of the latent process, $X_{[0,\tau ]}$⁠, from its conditional distribution given all of the other information, including ${\tau }_{1:T}$ and $o_{1:T}$⁠.

The presentation in Sections 2 and 3 focused on a single individual. To make the full process clear, let there be N subjects, and let the nth have $T^n$ observations $o_{1:T^n}^n$ at times ${\tau }_{1:T^n}^n$⁠. For subject n, we shift time so that ${\tau }_1^n=0$⁠, set ${\tau }^n=\tau -{\tau }_1^n$⁠, where τ corresponds to the end of 2014 and define the latent path to be $X_{[0,{\tau }^n]}^n$⁠. Each latent path $X_{[0,{\tau }^n]}^n$ gives rise to summary statistics $N_{l,m}^n$⁠, $R_l^n$⁠, and $N_i^n$ as defined for a specific subject in Section 2.1. We write $\mathbf{\Theta }=(\mathbf{B},\mathit{\lambda },\mathbf{Q},\nu )$ and ${\mathcal{D}}^n=\{T^n,{\tau }_{1:T^n}^n,O_{1:T^n}^n\}$⁠.

Following Fearnhead and Sherlock (2006); Luo et al. (2021), the Gibbs sampler has the following distinct steps.

1. For each n in $1,\dots ,N$⁠, simulate $(x_{{\tau }_1^n}^n,\dots ,x_{{\tau }_{T^n}^n}^n,x_{{\tau }^n}^n)$ from its conditional distribution given $\mathbf{\Theta }$⁠, ${\mathcal{D}}^n$ via the steps in Section 3.

2. For each n in $1,\dots ,N$ and each t in $1,\dots ,T^n$⁠, simulate $x_{({\tau }_t^n,{\tau }_{t+1}^n)}$ (identifying ${\tau }_{T+1}^n$ with ${\tau }^n$⁠) via the direct-sampling method of Hobolth and Stone (2009) or uniformization in Fearnhead and Sherlock (2006) and Rao and Teh (2013). This provides $x_{[0,{\tau }^n]}^n$ and, hence, the summaries $R_l^n$⁠, $N_{l,m}^n$ and $N_l^n$⁠, $n=1,\dots ,N$⁠.

3. Update ν from its conditional Dirichlet posterior given $x_0^{1:N}$⁠.

4. For $l=1,\dots ,K$ and $m\ne l$⁠, update $q_{l,m}$ from its conditional Gamma posterior given $R_l^{1:N}$ and $N_{l,m}^{1:N}$⁠.

5. For $k=1,\dots ,K$⁠, update each ${\lambda }_k$ from its conditional Gamma posterior given $R_k^{1:N}$ and $N_k^{1:N}$⁠.

6. For $k=1,\dots ,K$⁠, update ${\beta }_k$ from its conditional posterior given $x_{{\tau }_1^1}^1,\dots ,x_{{\tau }_{T^1}^1}^1,\dots ,$  $x_{{\tau }_1^N}^N,\dots ,x_{{\tau }_{T^N}^N}^N$⁠.

Our first simulation study involves a two-state chain and allows us to demonstrate the utility of both forms of information: the event times and values observed at these times. Our second simulated example, and our real-data example consider a three-state chain, with the third state being death and where the state can only increase. The non-zero transition rates are $q_{1,2}$⁠, $q_{1,3}$⁠, and $q_{2,3}$⁠, and the non-zero observation-process rates are ${\lambda }_1$ and ${\lambda }_2$⁠.

5 Empirical studies

5.1 Example 1

In this example, we demonstrate the performance of the proposed Gibbs sampler, and investigate the convergence, mixing and inferences on parameters under different scenarios.

For each of $N=50$ subjects, we simulate a realization of a two-state latent process without the death, $X_{[0,5]}$⁠, using a generator of

$$\mathbf{Q}=\left(\begin{array}{cc}-1& 1\\ 3& -3\end{array}\right),$$

and an initial distribution of $\nu =(0.8,0.2)$⁠. The outcome process is Gaussian, $O_t\sim \mathsf{N}({\beta }_k,1)$⁠, where ${\beta }_k$ is a scalar, so that $\mathbf{B}=({\beta }_1,{\beta }_2)$⁠. For each individual, the time window is $[0,5]$⁠; individuals are started from the initial distribution, ν, and the trajectory, events and observations simulated for the entire window. Thus, here, in contrast with our real-data scenario, there is no observation at time 0, and the first forward variable is ${\alpha }_{0,k}={\nu }_k$⁠.

We consider four scenarios in this example.

• Scenario I: $\mathbf{B}=(-1,1)$ and the point process rate $\lambda =(4,12)$⁠.

• Scenario II: $\mathbf{B}=(-1,1)$ and the point process rate $\lambda =(8,8)$⁠.

• Scenario III: $\mathbf{B}=(0.8,1)$ and the point process rate $\lambda =(4,12)$⁠.

• Scenario IV: We use the same $\mathbf{B}$ and λ as Scenario III but perform Bayesian inference using a CTHMM (Luo et al., 2021), ignoring the point process.

In Scenario I, both the Poisson intensity and the Gaussian expectation are distinct between states, whereas in Scenario II, while the Gaussian expectations are distinct, the Poisson intensities are the same. In Scenario III, the Gaussian expectations are similar but the intensity rates distinct, while Scenario IV is extremely challenging, with no information from the point process and very little information from the outcome process.

We place non-informative priors on all the parameters in this study. Each of ${\lambda }_1$⁠, ${\lambda }_2$⁠, $q_{12}$⁠, and $q_{21}$ have independent $\mathsf{G}\mathsf{a}\mathsf{m}\mathsf{m}\mathsf{a}(1,1/8)$ priors. We place a $\mathcal{N}(0,\mathrm{10,000})$ prior on each element of $\mathbf{B}$⁠, and a $\mathsf{B}\mathsf{e}\mathsf{t}\mathsf{a}(1,1)$ prior for ν.

Trace plots, auto-correlations plots, and histograms for ${\lambda }_1$ and $q_{12}$ are shown in Figures 3 and 4; Integrated auto-correlation times (IACTs) for ${\lambda }_1$⁠, ${\lambda }_2$⁠, $q_{1,2}$⁠, $q_{2,1}$⁠, ${\beta }_1$⁠, and ${\beta }_2$ are provided in Table 1. In Scenario I, where the data are generated from a model characterized by clear separation between the states, the Gibbs sampler performs well, with fast convergence and mixing. In Scenario II, where the timing of the events provides no information about the state, the Gibbs sampler’s performance declines slightly, though it still remains promising. The posterior distribution of ${\lambda }_1$ is slightly off the true value. Challenges arise when the data are generated from a less well-separated outcome process, leading to notable auto-correlations in the posterior samples, particularly those associated with $\mathbf{Q}$⁠. The reduction in mixing efficiency is considerable; however, the chain does still mix adequately over the 20,000 iterations shown. By contrast, in Scenario IV, when the point-process information is removed from Scenario III, the exploration of the posterior over the 20,000 iterations is so poor that we would judge that the chain has perhaps not yet even converged. Moreover, the posterior for $q_{1,2}$ includes values that are not strongly supported by the data. This illustrates the most important point: especially in cases where the outcome process is not particularly informative, the information from the timings of the events provides an important contribution to the posterior and, possibly because of the tighter posterior and more clearly defined paths for each $X_{[0,\tau ]}$⁠, improves the mixing of the chain.

Figure 3.

Example 1: Trace, ACF plots, and histogram for λ using the exact Gibbs sampler for different scenarios.

Figure 4.

Example 1: Trace, ACF plots, and histogram for $q_{12}$ using the exact Gibbs sampler for different scenarios.

Table 1.

Example 1: Integrated auto-correlation times (IACTs) of the posterior sample for each parameter in each scenario

Scenario	${\lambda }_1$	${\lambda }_2$	$q_{1,2}$	$q_{2,1}$	${\beta }_{1,1}$	${\beta }_{1,2}$
I	3.4	5.0	7.8	8.4	5.4	4.6
II	3.9	4.8	12.4	10.7	7.1	7.1
III	27.4	22.9	52.7	47.7	13.3	14.5
IV	n/a	n/a	(257.0)	(261.4)	(276.2)	(384.0)

Note. Bracketed terms indicate that the effective sample size was $<100$⁠, so there is reduced confidence in the estimate of the IACT.

5.2 Example 2

In this example, we demonstrate valid inference for simulated data with a similar pattern to our real data. We generate the data from a three-state model which is similar to that used for our real data, with $N=500$ subjects. For each subject, we simulate a realization of a three-state latent process with the death as the absorbing state, $X_{[0,10]}$⁠,

$$\mathbf{Q}=\left(\begin{array}{ccc}-0.21& 0.20& 0.01\\ 0.00& -0.05& 0.05\\ 0.00& 0.00& 0.00\end{array}\right),$$

so trajectories can only progress forward with no chances to transit back to previous states. The outcome process is $O_t\sim \mathsf{P}\mathsf{o}\mathsf{i}\mathsf{s}\mathsf{s}\mathsf{o}\mathsf{n}({\mu }_t)$ with $\log {\mu }_t=B_{1,k}+z_t{B}_{2,k}$ when in state k, and where the time-varying covariate is $Z_t\sim$Bernoulli$(0.65)$⁠; we set

$$\mathbf{B}=\left(\begin{array}{cc}-0.69& 0.77\\ -0.13& -0.39\end{array}\right).$$

This is equivalent to $O_t|(Z_t=1)\sim \mathsf{P}\mathsf{o}\mathsf{i}\mathsf{s}\mathsf{s}\mathsf{o}\mathsf{n}({\mu }_{1,k})$ and $O_t|(Z_t=0)\sim \mathsf{P}\mathsf{o}\mathsf{i}\mathsf{s}\mathsf{s}\mathsf{o}\mathsf{n}({\mu }_{0,k})$⁠, with $\log {\mu }_{1,k}=B_{1,k}+B_{2,k}$ and $\log {\mu }_{0,k}=B_{1,k}$⁠. In this case, the outcome expectations (marginalized over the covariate) are $0.35{\mu }_{0,k}+0.65{\mu }_{1,k}$⁠, which are $0.46$ and $1.71$ for the two states, respectively. The point process rate is set as $\lambda =(6,10,0{)}^{\mathrm{\top }}$⁠. State 3 is the death state without any observed outcome. If the patient progresses to the absorbing state before the observation window, i.e. death, then the last observation is the final one before death. Otherwise, the last observation is the one right before the observation window, $\tau =10$⁠. We assumed that the first observation is made at time 0 according to the initial distribution. Similar to Example 1, we place non-informative priors on all the parameters in this study. Each of ${\lambda }_1$⁠, ${\lambda }_2$⁠, $q_{12}$⁠, and $q_{21}$ have independent $\mathsf{G}\mathsf{a}\mathsf{m}\mathsf{m}\mathsf{a}(1,1/8)$ priors. We place independent $\mathsf{G}\mathsf{a}\mathsf{m}\mathsf{m}\mathsf{a}(0.1,0.1)$ priors on each of the four Poisson expectations, ${\mu }_{i,k}$⁠, $i=0,1$⁠, $k=1,2$⁠, and a $\mathsf{B}\mathsf{e}\mathsf{t}\mathsf{a}(1,1)$ prior for ${\nu }_1$⁠.

Table 2 presents the integrated auto-correlation times (IACTs) for $\mathit{\lambda }$⁠, $\mathbf{Q}$⁠, and $\mathbf{B}$ over 20,000 iterations and shows that all of these parameters mix well. Figure 5 presents a histogram of the marginal posterior for each of these parameters, with each true parameter value falling within the main posterior mass. This example demonstrate that by simulating the data-induced window similar to the real data, we can achieve valid inference.

Figure 5.

Example 2: marginal histograms of $q_{12}$⁠, $q_{1,3}$⁠, $q_{2,3}$⁠, ${\lambda }_1$⁠, and ${\lambda }_2$⁠; the associated vertical red lines depict the true parameter value.

Table 2.

Example 2: Integrated auto-correlation times (IACTs) of the posterior sample for each parameter

${\lambda }_1$	${\lambda }_2$	$q_{1,2}$	$q_{1,3}$	$q_{2,3}$	${\beta }_{1,1}$	${\beta }_{1,2}$	${\beta }_{2,1}$	${\beta }_{2,2}$
3.0	3.0	3.2	3.5	3.1	3.4	3.7	3.2	3.0

5.3 Application: Canadian healthcare system data

In this section, we delve into the proposed methodology designed to enhance our understanding of the progression of chronic diseases, particularly COPD. The significance of this understanding lies in its potential to help with effective management of healthcare systems and pave the way for methods of early diagnosis and personalized care, as the model can help the management team understand the ‘care path’ or ‘trajectory’ of the cohort before they become severely ill. While clinical and administrative data harbour valuable insights, traditional methods for modelling disease progression face limitations when dealing with data collected at irregular intervals and interaction frequency, such as when patients choose to engage with the healthcare system.

To address these challenges, our proposed methodology is applied to a COPD cohort derived from healthcare administrative data. An open, dynamic cohort was established, comprising a 25% random sample of Montreal residents in Canada. In 1998, a similar random sample was drawn from the Régie de l’assurance maladie du Québec (RAMQ)-registered population within the census metropolitan area of Montreal. Subsequently, at the onset of each subsequent year, 25% of individuals born in Montreal or relocating to Montreal within the past year were sampled to maintain cohort representativeness. Follow-up concluded upon death or a change in residential address outside Montreal; however, neither death nor the time of address change is recorded in our data. This administrative database encompasses outpatient diagnoses and procedures submitted through RAMQ billing claims, as well as inpatient claims for procedures and diagnoses. Drug dispensing data are available for individuals covered by drug insurance through RAMQ, which includes roughly half the population, especially all those aged 65 and above. All data are interconnected through an anonymized version of the RAMQ number, and administered by the Surveillance Lab at McGill Clinical & Health Informatics research group.

Utilizing established case definitions based on diagnostic codes (Lix et al., 2018), a total of 76,888 COPD patients were enrolled, with an incident event (ICD-9 491x, 492x, 496x; ICD-10 J41–J44) occurring after a minimum of 2 years at risk with no events. The observation period spanned from January 1998, commencing with the patients’ first diagnosis, until December 2014. Physicians observed these patients solely during medical visits, which transpired when patients chose to interact with the healthcare system. In these medical visits, pertinent medical information was gathered, and from this observed information, we aim to infer the patients’ unobserved disease status, which is represented using a discrete disease state model. In this study, the disease status was indirectly gauged through a proxy: the number of prescribed medications at the time of observation. Notably, this information was available exclusively for patients with drug insurance, leading us to focus on patients aged over 65 years, constituting the study population for our analysis. This full dataset has been previously analysed in Luo et al. (2021). Due to the data restriction, we are allowed to randomly sample 1,000 patients to demonstrate our proposed method. Figure 6 shows one COPD patient’s follow-up from this dataset. As shown in the figure, we observe the number of drugs prescribed as the observation process. In addition, we have also observe variations in the type of encounters, i.e. general practitioner visit (GP), hospitalization (HOSP), specialist visit (SPEC), and emergency department visit (ED), which will be incorporated as time-varying covariates in our model. These changes in the observation point process, transitioning from sparse to dense, reflect the evolving severity of the chronic disease under study. In the dataset, the average number of drugs prescribed from a visit for GP, ED, HOSP, and SPEC is 5.06, 5.57, 5.28, and 5.17 respectively, suggesting little deviation among the types of healthcare utilization. Over the 1,000 patients, the average time from the first recorded interaction until the last observation is about 13.71 years.

Figure 6.

COPD patient’s trajectory from the Canadian healthcare system data.

We apply the model from Section 2 to this sample dataset, with the inference methodology from Section 4. We focus on a three-state model with State 3 (the unobserved death state) as an absorbing state. We place independent $\mathsf{G}\mathsf{a}\mathsf{m}\mathsf{m}\mathsf{a}(1,1/8)$ priors on each ${\lambda }_i$ and all elements in $\mathbf{Q}$⁠, and a $\mathsf{G}\mathsf{a}\mathsf{m}\mathsf{m}\mathsf{a}(0.1,0.1)$ prior on each expected number of drugs dispensed for each healthcare utilization, and a $\mathsf{B}\mathsf{e}\mathsf{t}\mathsf{a}(1,1)$ prior for $({\nu }_1,{\nu }_2)$⁠, with ${\nu }_3=0$⁠. The outcome process is modelled with a Poisson distribution, with the types of healthcare utilization (HOSP, SPEC, GP, and ED) as a time-dependent covariate; equivalently, for each of States 1 and 2 there are four parameters, indicating the expected number of drugs dispensed for each of the four values of the covariate.

We ran the algorithm for 20,000 iterations; for all parameters, the estimated IACTs were <5, indicating excellent mixing of the MCMC algorithm. Table 3 displays the estimates for all of the model parameters. The narrow 95% credible intervals associated with these estimated parameters indicate low uncertainty within the estimated densities. The posterior median for the number of prescribed drugs rises by about a factor of four from State 1 to State 2, reflecting the worsening severity of patients’ health conditions. In State 2, expected drug dispensations for emergency department visits and hospitalizations are slightly higher than other types of visits. Point process rates are roughly 35% higher in State 2 than in State 1, suggesting increased healthcare system interactions during poorer health conditions. This indicates a correlation between disease severity and healthcare utilization, with State 2 indicating the higher level of healthcare engagement among patients. In the state transition dynamics, the median stay in State 1 is roughly $2.3$ years and the probability that a patient is still in State 1 after 4 years is approximately $0.30$⁠. When a transition from State 1 occurs, the probability that it will be directly to the death state is $0.23$⁠. The estimates for State 2 suggest a more persistent state (median stay being over 4 years) and potentially more severe state of illness. Patients in State 1 are less likely to die (or move from the province) in any given year, compared to those in State 2. The median time to death in State 2 is 4.3 years. We also notice that ${\nu }_1\approx {\nu }_2$ meaning that conditional on a healthcare visit, patients are almost equally likely to start in each of the two states in the healthcare system.

Table 3.

Posterior median of MMPP parameters associated with 95% credible intervals

Parameter	State 1	State 2
λ (Marked point process rate)	5.07 (4.97, 5.16)	6.82 (6.72, 6.93)
ν (initial state distribution)	0.54 (0.51, 0.58)	0.46 (0.42, 0.49)
Expected number of drugs dispensed (GP)	2.12 (2.10, 2.14)	7.28 (7.26, 7.31)
Expected number of drugs dispensed (ED)	1.72 (1.67, 1.78)	7.47 (7.43, 7.53)
Expected number of drugs dispensed (HOSP)	1.77 (1.72, 1.81)	7.36 (7.31, 7.41)
Expected number of drugs dispensed (SPEC)	1.79 (1.75, 1.85)	7.16 (7.10, 7.21)
$q_{1.}$ (Transition rate from State 1)	$-$	0.23 (0.21, 0.25)
$q_{.3}$ (Transition rate to death state)	0.07 (0.06, 0.08)	0.16 (0.14, 0.17)

We also ran the Gibbs sampler without the point process as with Scenario IV in Example 1. The parameter estimates are reproduced in the online supplementary material. They show that ignoring the point process information makes patients seem to remain in State 1 for longer. The other substantial difference is that the transition rate to the death state from State 1 is slightly greater than that from State 2. We also fitted a four-state model, again with only transitions to higher states allowed, and now with State 4 (the unobserved death state) as an absorbing state. We compared a Laplace-based approximation of its marginal likelihood with a similar approximation to the marginal likelihood of the three-state model; see the online supplementary material for details. The estimated marginal likelihoods for three-state and four-state models are $-2,702.93$ and $-2,687.05$⁠, respectively, leading to corresponding BICs of $5,502.57,5,612.19$⁠. This shows that the three-state model is preferable; results for the four-state model are also presented in the online supplementary material.

This inferred state dynamic provides valuable insights for healthcare management, offering guidance on the performance of current COPD interventions. The tendency for patients in State 1 to remain in that state initially before transitioning to State 2 underscores the need for interventions aimed at preventing or delaying disease progression, particularly during the early stages. For patients already in State 2, the observed persistence suggests a potentially more challenging treatment landscape, emphasizing the importance of comprehensive and long-term management strategies. Understanding these state transition dynamics informs healthcare providers and policymakers about the effectiveness of existing interventions and guides the development of targeted strategies to improve patient outcomes. By identifying critical transition points and patient populations at higher risk of disease progression, healthcare resources can be allocated more effectively to provide timely and appropriate interventions, ultimately improving the management and care of individuals with COPD. However, it is important to note that there is no guarantee that the latent state and trajectory groupings will directly correspond to the underlying biologic processes or disease stages as our initial model is built without clinical measurements. Therefore, the identified states correspond to the general health conditions of COPD patients, implying that a higher disease level or intensity of healthcare encounters does not strictly indicate a deterioration in COPD but could signify any temporary or prolonged period of poor health condition. By using this Canadian cohort dataset, the initial latent groupings are expected to generate hypotheses about how the patients are being managed.

6 Discussion

In this article, we have considered a generalization of both the hidden Markov model and the MMPP to an MMPP with an outcome process, and we have proposed an exact Gibbs sampler to facilitate Bayesian inference. Our data model also incorporates an unobserved terminating event (death or moving away from the area).

This advance is driven by the need to analyse longitudinal data and model an unobserved death so as to comprehend the progression of a disease processes. A critical decision revolves around whether observation times are considered non-informative or informative, on the state of the Markov chain, which relates to the disease severity. Neglecting to account for an informative observation process in the disease dynamics or to incorporate an unobserved terminating event, such as death, can potentially lead to biased parameter estimates. With the large amount of claim data and healthcare data available, there is a pressing need to work on appropriate statistical modelling framework for longitudinal healthcare studies. The proposed technique for inferring death prevents the bias that would arise if a model without death were employed with a shortened observation window only up to the final event time. Our approach offers a framework for analysing complex healthcare data, providing insights into disease progression and patient health/death states that can inform clinical decision-making and healthcare interventions.

Since our algorithm is a Gibbs sampler, it requires no tuning parameters and could be used ‘off the shelf’ by practitioners. Simulation studies indicate the effectiveness of the Gibbs sampler and show how the event times contribute vital information, especially when the outcome process lacks a clear separation between the states. By integrating observation times and death into the model, we enhance our understanding of disease dynamics and healthcare utilization patterns, thus advancing the methodology for analysing longitudinal data in healthcare research. In our application, we imposed a constraint that the process can only progress forward since COPD gets progressively worse over time (Cleveland Clinic, 2024). However, one could make no restrictions on the allowed transitions in the state space. In the absence of other information, the hidden state would then be non-identifiable with all permutations of the state labels equivalent. To avoid the ensuing identifiability issue, meaningful priors could be placed on the outcome models for each state or the outcome process event rates—that is, we would distinguish the states via their impact on the observable quantities or imposing certain restrictions on $\mathbf{Q}$ (see Section 4 in Luo et al., 2021 for a detailed discussion).

Although the proposed model possesses a flexible structure, there remains some scope for future research to explore several possibilities in analysing such irregularly spaced longitudinal data. The current setup assumes time-homogeneous processes $X_s$ and $N_s$⁠. Both processes may plausibly depend on time. For instance, within any given health state, older patients might engage more frequently with the healthcare system. However, introducing complex structures to model the transition in $X_s$ poses computational challenges in both frequentist and Bayesian inference (see for example, Kendall et al., 2024, and a discussion in the online supplementary material of this article). In addition, the issue of identifying the number of hidden states in MMPPs has been a long-standing challenge. In the context of HMMs, this topic has been explored using approaches such as marginal likelihood, cross-validated likelihood, and complete-data likelihood (see, Celeux and Durand, 2008; Luo and Stephens, 2021; Pohle et al., 2017); however, model selection criteria based on these quantities become unreliable as the dimensionality increases beyond 50 (Chen and Chen, 2008). Moreover, while the current model categorizes patients into discrete states, which is appropriate for some chronic diseases thought to have discrete stages of progression, relaxing the assumption of discrete states could yield insights. Relaxing the assumption that patients evolve through discrete states of underlying disease progression, would allow estimation of a continuum of severity that may better represent the trajectory of COPD and other chronic diseases. Substantial advances in statistical methodology would be required to extend this approach to a continuous state space and to allow the latent state distribution to have an interpretation as a continuous index or score. For example, one could use the features included in comorbidity indices to measure multimorbidity in terms of the ability to predict future mortality and health services use. The inference problem for this type of assumption involves multivariate diffusion processes, using discrete-time data that may be incomplete and affected by measurement errors (e.g. Golightly and Sherlock, 2022).

Acknowledgments

The authors would like to thank Prof. David Buckeridge from the Surveillance Lab at McGill Clinical & Health Informatics research group for granting us access to the real data that enabled us to demonstrate the empirical study in this article. The authors would also like to thank two referees for their constructive comments, which have improved the content and clarity of this article.

Data availability

Real data used in this article to illustrate our findings are not shared as restrictions apply to the availability of these data, which were used under license for this study. Simulated data related to Example 1 are generated by the code supplied in the online supplementary material.

Supplementary material

Supplementary material is available online at Journal of the Royal Statistical Society: Series C.

Author notes

Conflicts of interest The authors have no conflict of interest to disclose.