Bayesian Estimation of Multivariate Autoregressive Hidden Markov Model with Application to Breast Cancer Biomarker Modeling

Abstract

In this work, a first-order autoregressive hidden Markov model (AR(1)HMM) is proposed. It is one of the suitable models to characterize a marker of breast cancer disease progression essentially the progression that follows from a reaction to a treatment or caused by natural developments. The model supposes we have observations that increase or decrease with relation to a hidden phenomenon. We would like to discover if the information about those observations can let us learn about the progression of the phenomenon and permit us to evaluate the transition between its states (supposed discrete here). The hidden states governed by the Markovian process would be the disease stages, and the marker observations would be the depending observations. The parameters of the autoregressive model are selected at the first level according to a Markov process, and at the second level, the next observation is generated from a standard autoregressive model of first order (unlike other models considering the successive observations are independents). A Markov Chain Monte Carlo (MCMC) method is used for the parameter estimation, where we develop the posterior density for each parameter and we use a joint estimation of the hidden states or block update of the states.

1. Introduction

The main motivation behind this work is to characterize progression in breast cancer. In fact, disease progression cannot be assessed correctly without the use of biomarkers, which would effectively monitor the evolution of the patient health state, and this is the case for breast cancer. The major challenge in this matter for researchers and clinicians is to unravel the stage of the disease, so as to tailor the treatment for each patient and to monitor the response of a patient to a treatment.

Currently, studies have shown that there is a correlation between the levels of certain markers such as cancer antigen CA15-3, carcinoembryonic antigen (CEA), and serum HER2 Neu with the stage of the disease [1]. This gives an opportunity of using a hidden Markov model (HMM) to predict the stage of the disease based on biomarker data and to address the effectiveness of the treatments in their influence on the transition of the cancer from one state to another. In HMM, we have two constituents: the Markovian hidden process suitable to represent the breast cancer stage and the observation process given by the biomarker data. By the way, we can learn about the disease transition rates and how it progresses from primary breast cancer to advanced cancer stage, for example.

Indeed, HMM is a useful tool for tackling numerous concrete problems in many fields but some possible applications of HMM are in speech processing [2], in biology [3], in disease progression [4], in economics [5, 6], and in gene expression [7]. For a complete review of HMM, the reader is referred to Zucchini and MacDonald [8], in which properties and definitions of HMM are presented in a plausible way with both classical estimation by maximum likelihood method and expectation maximization (EM) algorithm and the new Bayesian inference is addressed.

The model we consider here is a variation of the regular hidden Markov model, since we will use extensions to incorporate dependence among successive observations, suggesting autoregressive dependence among continuous observations. Consequently, we have relaxed the conditional independence assumption from a standard HMM, because we would like to add some dynamics to the patient disease progression and because in reality the current patient biomarker observation is dependent on the past one. In fact, the autoregressive assumption in HMM has shown its advantage over regular HMM that cannot catch the strong dependence between successive observations (e.g., Ref. [9]). A similar model to ours can be found in Ref. [10]. This kind of models, which were first proposed in Ref. [11] to describe econometrics time series, is generalization of both HMM and autoregressive models, will be effective in representing multiple heterogeneous dynamics such as the disease progression dynamics, and can be even generalized to a regime-switching ARMA models such as in Ref. [12].

Moreover, Our model can also be viewed as an extension of the multivariate double-chain Markov model (DCMM) developed by Ref. [13], where there are two discrete Markov chains of first order: the first Markov chain is observed and the second one is hidden. In contrast to this DCMM, our multivariate first-order autoregressive hidden Markov model (MAR(1)HMM) will lead to continuous observations, where each observation conditional on the hidden process will depend on the previous observation according to an autoregressive process of first order. This dynamic is promising for continuous observed disease biomarkers.

Parameter estimation is very challenging for HMM family models since the likelihood is not available in a closed form most of the time. Thus, we call for a Markov Chain Monte Carlo (MCMC) procedure instead of a maximum likelihood-based approach. This choice rises from the fact that the Bayesian analysis uses prior knowledge about the process being measured, and it allows direct probability statements and an approximation of posterior distributions for the parameters. Instead in the maximum likelihood approach, we cannot have declared prior or have exact distribution for the parameters when the likelihood is untractable or when we have missing data (e.g., Refs. [14–16]).

Since the realization of HMM includes two separate entities: the parameters and the hidden states, the Bayesian computation is carried out after augmenting the likelihood by the missing hidden states [17]. The hidden states are sampled using a Gibbs sampler adopting a joint estimation of the hidden states or block update of the states (instead of a single update of each state separately) by means of a forward filtering/backward smoothing algorithm. Given the hidden states, we can compute the autoregressive parameters and the transition probabilities of the Markov chain by Gibbs sampler from their posterior densities after specifying conjugate priors for the parameters. Hence, the MCMC algorithm will alternate between simulating the hidden states and the parameters. Finally, we can obtain posteriors statistics such as the means, standard deviations and confidence intervals after assessing the convergence of the MCMC algorithm.

This chapter is organized as follows: after a preliminary on HMM, a description of the model is given in Section 3. In Section 4, we will give the Bayesian estimation of the parameters and the hidden states and provide the details of the MCMC algorithm, before presenting the results of a simulation studies in Section 5 and we will finish by a conclusion.

2. Preliminary

Since the model suggested is of the HMM type, we will describe HMM in more detail: an HMM is a stochastic process ${\{X_t,Y_t\}}_{t=0}^T$ , where ${\left\{X_t\right\}}_{t=0}^T$ is a hidden Markov chain (unobservable) and ${\left\{Y_t\right\}}_{t=0}^T$ is a sequence of observable independent random variables such that Y _t depends only on X _t for the time t = 0,1,…T. Here the process ${\left\{X_t\right\}}_{t=0}^T$ evolves independently of ${\left\{Y_t\right\}}_{t=0}^T$ and is supposed to be a homogeneous finite Markov chain with probability transition matrix Π of dimension a × a, where a indicates the number of the hidden states and Π₀ = (Π₀₁,…,Π_0a ) is the initial state distribution.

We denote the probability density function of Y _t = y _t given X _t = k for k ∈ {1, …, a} with $P_{x_t}(y_t,{\theta }_k)$ , where θ _k refers to the parameters of P when X _t = k. We suppose further that the processes Y _t |X _t and Y _t′ |X _t′ are independent for t ≠ t′. Let Θ = (θ ₁,…, θ _a ) and θ = (Π₀, Π, Θ), and then, the HMM can be described as follows: First, the likelihood of the observations and the hidden states can be decomposed to $P(y_0,\dots ,y_T,x_0,\dots ,x_T,\theta )=P(y_0,\dots ,y_T|x_0,\dots ,x_T,\theta )P(x_0,\dots ,x_T,\theta ).$ Since ${\left\{X_t\right\}}_{t=0}^T$ is a Markov chain, $P(x_0,\dots ,x_T,\theta )={\text{Π}}_0(x_0){\displaystyle {\displaystyle \prod _{t=1}^T}}\text{Π}(x_t|x_{t-1})$ . Under the conditional independence of the observations given the hidden states, $P(y_0,\dots ,y_T|x_0,\dots ,x_T,\theta )=P_{x_0}(y_0|{\theta }_{x_0}){\displaystyle {\displaystyle \prod _{t=1}^T}}{P}_{x_t}(y_t|{\theta }_{x_t})$ . Consequently, the likelihood function for the hidden states and the observations is given by $P(y_0,y_1,\dots ,y_T,x_0,x_1,\dots ,x_T,\theta )=\text{Π}(x_0)P_{x_0}(y_0|{\theta }_{x_0}){\displaystyle {\displaystyle \prod _{t=0}^T}}\text{Π}(x_t|x_{t-1})P_{x_t}(y_t|{\theta }_{x_t})$ .

3. Model description and specification

The MAR(1)HMM model we consider in this work is a hidden Markov model, where conditionally on the latent states, the observations are not independent like it is the case for a regular hidden Markov model. Instead, the current observation is allowed to depend on the previous observation according to an autoregressive model of first order. As in an HMM model, the latent states evolve according to a discrete first-order time homogeneous Markov model. We consider data of n continuous random variables observed over time, each of potentially different lengths, i.e., for each individual i = 1,2, …, n, we observe a vector $y_{i\mathrm{,.}}={(y_{i,u_i},\dots ,y_{i,m_i})}^T$ , with u _i < m _i .

Define $u_0={\displaystyle \underset{1\le i\le n}{\min }}\left\{u_i\right\}$ and $M={\displaystyle \underset{1\le i\le n}{\max }}\left\{m_i\right\}$ and note that the times u _i and m _i may vary over the entire observation period from u ₀ to M with the restriction that u _i – m _i ≥ 1, for i = 1,2,…,n.

We assume, for i = 1,2,…,n for integer time t = u _i ,…,m _i , that the random variable Y _i,t taking nonnegative values depends only on the states X _t and the previous observation Y _i,t–1, and based on the model developed by Farcomeni and Arima [10], we get the following model:

$Y_{i,t}{|}_{X_t=x_t}={\beta }^{(x_t)}{Y}_{i,t-1}+{\mu }^{(x_t)}+{\epsilon }_{i,t}\mathrm{.}$ E1

The choice of the autoregressive part of the model is motivated by the fact that successive biomarker observations are most of the time correlated from many diseases unlike the hypothesis of independence between observations in HMMs.

We interpret x as the vector of the hidden health states of the patients; in the case of breast cancer, those states would be localized or advanced metastatic breast cancer for example, while y is the vector of the biomarkers observed and measured for the patients. The ε _i , _t are normal variables with mean 0 and variance σ ² such that ε _i , _t and ε _i′ , _t′ are uncorrelated, (i, t) ≠ (i′, t′).

The parameters ${\beta }^{\left(x_t\right)}$ and ${\mu }^{\left(x_t\right)}$ are parameters taking values in $\mathbb{R}$ for each hidden state and σ ² ∈ ${\mathbb{R}}^{+}$ .

Similar to Ref. [13], the transition matrix of the Markov chain Π is time homogeneous with dimension a × a where a is the number of hidden states, and Π = (Π _gh , g = 1,…,a; h = 1,…,a) where Π _gh = P(X _t = h|X _t–1 = g), for g = 1,2,…,a; h = 1,2,…, a; and t = u ₀₊₁,…, M. We let the first state $X_{u_0}$ to be selected from a discrete distribution with vector of probabilities r = (r ₁,…,r _a ). Also we consider the time of initial observation u _i , the initial observed state $y_{i,u_i}$ , and the number of consecutive time points that were observed m _i – u _i + 1. Let μ = (μ ⁽¹⁾,…, μ ^(a)), β = (β ⁽¹⁾,…,β ^(a)), and θ = (μ, β, σ ², r, Π) be the set of all parameters in the model. We suppose that the individuals, i.e., Y _i , _t , behave independently conditionally on X. Therefore, for i = 1, …, n, $P(y_{i\mathrm{,.}}|y_{i,u_i},x,\theta )={\displaystyle {\displaystyle \prod _{t=u_i+1}^{m_i}}}P(y_{i,t}|y_{i,t-1},x_t,\text{Θ})$ and $P(x|\theta )=P(x_{u_0}){\displaystyle {\displaystyle \prod _{t=u_0+1}^M}}P(x_t|x_{t-1},\text{Π})$ , where $P(x_t|x_{t-1},\text{Π})=P(X_t=x_t|X_{t-1}=x_{t-1},\text{Π})={\text{Π}}_{x_{t-1},x_t}\mathrm{.}$ Then, the likelihood density for the observations of all individuals y = (y ₁,…,y _n ) given first time vector of observations $y_0=(y_{\mathrm{1,}{u}_1},\dots ,y_{n,u_n})$ , x, and θ is $P(y|y_0,x,\theta )={\displaystyle {\displaystyle \prod _{i=1}^n}}P(y_{i\mathrm{,.}}|y_{i,u_i},x,\theta ),$

This is due to the conditional independence of the y _i , given x and θ. The joint mass of each y _i ,. and x given $y_{u_i}$ and θ can be written as follows: $P(y_{i\mathrm{,.}},x|y_{i,u_i},\theta )=P(y_{i\mathrm{,.}}|y_{i,u_i},x,\theta )\times P(x|y_{u_i},\theta ).$ Using the Markov property of the hidden process, we have after simplification $P(x|y_{i,u_i},\theta )\propto P(y_{i,u_i}|x,\theta )P(x|\theta )=P(y_{i,u_i}|x_{u_i},\theta )r_{x_{u_0}}{\text{Π}}_{x_{u_0},x_{u_0+1}}\times \cdots \times {\text{Π}}_{x_{M-1},x_M}\mathrm{.}$ In addition, $P(y_{i\mathrm{,.}}|y_{i,u_i},x,\theta )={\displaystyle {\displaystyle \prod _{t=u_i+1}^{m_i}}}P(y_{i,t}|y_{i,t-1},x,\theta )$ , and consequently,

$P(y_{i\mathrm{,.}},x|y_{i,u_i},\theta )\propto r_{x_{u_0}}P(y_{i,u_i}|x_{u_i},\theta ){\displaystyle {\displaystyle \prod _{t=u_0+1}^M}}{\text{Π}}_{x_{t-1},x_t}{\displaystyle {\displaystyle \prod _{t=u_i+1}^{m_i}}}P(y_{i,t}|y_{i,t-1},x,\theta ).$ Finally, under the hypothesis of normal error distribution for the autoregressive parameters of the model (Eq. (1)) and the Chapman-Kolmogorov property, the joint distribution of y _i ,. and x given $y_{i,u_i}$ and θ can be simplified to:

$\begin{array}{c}P(y_{i\mathrm{,.}},x|y_{i,u_i}\theta )\propto P(y_{i,u_i}|x_{u_i},\theta ){\displaystyle {\displaystyle \prod _{h=1}^a}}{r}_h^{{\chi }_{\left\{x_{u_0}\right\}}(h)}{\displaystyle {\displaystyle \prod _{t=u_0+1}^M}}{\displaystyle {\displaystyle \prod _{g=1}^a}}{\displaystyle {\displaystyle \prod _{h=1}^a}}{\text{Π}}_{g,h}{}^{{\chi }_{\{x_t,x_{t-1}\}}(g,h)}\\ \times {\displaystyle {\displaystyle \prod _{t=u_i+1}^{m_i}}}{\displaystyle {\displaystyle \prod _{h=1}^a}}{\left[\frac{1}{\sigma }\varphi \left(\frac{y_{i,t}-{\mu }^{(h)}-{\beta }^{(h)}{y}_{i,t-1}}{\sigma }\right)\right]}^{{\chi }_{\left\{x_t\right\}}(h)},\end{array}$ E2032

where ϕ denotes the density of a standard normal distribution $\mathcal{N}(\mathrm{0,1})$ and ${\chi }_{\left\{A\right\}}(x)$ is the usual indicator function of a set A. Finally, the joint distribution of y and x has the following form:

$\begin{array}{l}P(y,x|y_0,\theta )\propto {\displaystyle {\displaystyle \prod _{h=1}^a}}{r}_h^{{\chi }_{\left\{x_{u_0}\right\}}(h)}{\displaystyle {\displaystyle \prod _{t=u_0+1}^M}}{\displaystyle {\displaystyle \prod _{g=1}^a}}{\displaystyle {\displaystyle \prod _{h=1}^a}}{\text{Π}}_{g,h}{}^{{\chi }_{\{x_t,x_{t-1}\}}(g,h)}\\ \times {\displaystyle {\displaystyle \prod _{i=1}^n}}{\displaystyle {\displaystyle \prod _{l=1}^a}}{\left[\frac{1}{\sigma }\varphi \left(\frac{y_{i,u_i}-{\mu }^{(l)}}{\sigma }\right)\right]}^{{\chi }_{\left\{x_t\right\}}(l)}\times {\displaystyle {\displaystyle \prod _{i=1}^n}}{\displaystyle {\displaystyle \prod _{t=u_i+1}^{m_i}}}{\displaystyle {\displaystyle \prod _{h=1}^a}}{\left[\frac{1}{\sigma }\varphi \left(\frac{y_{i,t}-{\mu }^{(h)}-{\beta }^{(h)}{y}_{i,t-1}}{\sigma }\right)\right]}^{{\chi }_{\left\{x_t\right\}}(h)}\mathrm{.}\end{array}$ E2

4. Bayesian estimation of the model parameters

We will use a Bayesian approach to estimate the model parameters. Inference in the Bayesian framework is obtained through the posterior density, which is proportional to the prior multiplied by the likelihood. The posterior distribution for our model, as in most cases, cannot be derived analytically, and we will approximate it through MCMC methods specifically designed for working with the augmented likelihood with the hidden states. In fact, MCMC methods start by specifying the prior density Π(θ) for the parameters. Since the data Y are available, the general sampling methods work recursively by alternating between simulating the full conditional distributions X given y and θ given x and y.

4.1. Prior distributions

Under the assumption of independence between the parameters $\theta =(\mu ,\beta ,{\sigma }^2,r,\text{Π})$ , the prior density could be written as $P(\theta )=P(r)P(\text{Π})P(\mu )P(\beta )P({\sigma }^2)$ . r is the parameters of a multinomial distribution; hence, the natural choice for the prior would be a Dirichlet distribution $r\sim \mathbb{D}({\alpha }_{01},\dots ,{\alpha }_{0a})$ . Later on, ${\displaystyle {\displaystyle \sum _{j=1}^a}}{\text{Π}}_{ij}=1$ , and we assume that ${\text{Π}}_i\sim \mathbb{D}({\delta }_{i1},\dots ,{\delta }_{ia})$ for each row i of the transition matrix. This choice of the Dirichlet prior can be even the default $\mathbb{D}(\mathrm{1,}\dots \mathrm{,1})$ as recently discussed in Ref. [18]. In fact, a Dirichlet prior is justified because the posterior density of each row of the transition matrix is proportional to the density of a Dirichlet distribution, and hence, choosing a Dirichlet prior would give a posterior Dirichlet. This can be justified as follows for a given set of parameters $\lambda =({\lambda }_1,\dots ,{\lambda }_a)$ from a discrete or from a multinomial density:

$\pi (x_1,\dots ,x_a,{\lambda }_1,\dots ,{\lambda }_a)=\frac{n!}{x_1!\dots x_a!}{\lambda }_1^{x_1}\dots {\lambda }_a^{x_a}$ for the nonnegative integers $x_1,\dots ,x_a$ , with ${\displaystyle {\displaystyle \sum _{i=1}^a}}{x}_i=n\mathrm{.}$ This probability mass function can be expressed, using the gamma function Γ, as $\pi (x_1,\dots ,x_a,{\lambda }_1,\dots ,{\lambda }_a)=\frac{\Gamma \left({\displaystyle {\displaystyle \sum _{i=1}^a}}{x}_i+1\right)}{{\displaystyle {\displaystyle \prod _{i=1}^a}}\Gamma (x_i+1)}{\displaystyle {\displaystyle \prod _{i=1}^a}}{\lambda }_i^{x_i}\mathrm{.}$ This form shows its resemblance to the Dirichlet distribution, and by starting from supposing the prior $\lambda \propto \mathbb{D}({\alpha }_0,\dots ,{\alpha }_a)$ , the posterior is $P(\lambda |x)\propto P(\lambda )P(x|\lambda )\propto {\displaystyle {\displaystyle \prod _i}}{\lambda }_i^{x_i}{\displaystyle {\displaystyle \prod _i}}{\lambda }_i^{{\alpha }_i-1}\propto {\displaystyle {\displaystyle \prod _i}}{\lambda }_i^{x_i+{\alpha }_i-1}\propto \mathbb{D}(x_1+{\alpha }_1,\dots ,x_a+{\alpha }_a).$

Furthermore, concerning the priors for parameters of the autoregressive model, we suppose for $h=\mathrm{1,}\dots ,a$ : ${\mu }^{(h)}\sim \mathcal{N}({\alpha }_h,{\tau }_h)$ , ${\beta }^{(h)}\sim \mathcal{N}(b_h,c_h)$ , and inverse gamma ( $\mathbb{I}\mathbb{G}$ ) prior for ${\sigma }^2\sim \mathbb{I}\mathbb{G}(\epsilon ,\zeta )$ . ${\alpha }_h,{\tau }_h,b_h,c_h,ϵ,\zeta$ are hyperparameters to be specified. For more details on Bayesian inference and prior selection in HMM, the reader is referred to Ref. [19]. In our case, prior distributions for the autoregressive parameters were proposed by Ref. [20] for a mixture autoregressive model, who points out that they are conventional prior choices for mixture models.

4.2. Sampling the posterior distribution for the hidden states

Chib [21] developed a method for the simulation of the hidden states from the full joint distribution for the univariate hidden Markov model case. We will describe his full Bayesian algorithm for the univariate hidden Markov model before a generalization to our MAR(1)HMM.

4.2.1. Chib’s algorithm for the univariate hidden Markov model for estimation of the states

Suppose we have an observed process $Y_n=(y_1,\dots ,y_n)$ and the hidden states $X_n=(x_1,\dots ,x_n)$ , θ are the parameters of the model. We adopt for simplicity $X_t=(x_1,\dots ,x_t)$ the history of the states up to time t and $X^{t+1}=(x_{t+1},\dots ,x_n)$ the future from t + 1 to n. We use the same notation for Y _t and Y ^t+1.

For each state $x_t\in \{\mathrm{1,2,}\dots ,a\}$ for $t=\mathrm{1,2,}\dots ,n$ , the hidden model can be described by a conditional density given the hidden states $\pi (y_t|Y_{t-1},x_t=k)=\pi (y_t|Y_{t-1},{\theta }_k)$ , $k=\mathrm{1,}\dots ,a$ , with x _t depending only on x _t–1 and having transition matrix Π and initial distribution Π₀, and the parameters for π(.) are $\theta =({\theta }_1,\dots ,{\theta }_a)$ .

Chib [21] shows that it is preferable to simulate the full latent data $X_n=(x_1,\dots ,x_n)$ from the joint distribution of $x_1,\dots ,x_n|Y_n,\theta$ , in order to improve the convergence property of the MCMC algorithm because instead of n additional blocks if each state is simulated separately, only one additional block is required. First, we write the joint conditional density as

$P(X_n|Y_n,\theta ,\text{Π})=P(x_n|Y_n,\theta )P(x_{n-1}|Y_n,x_n,\theta ,\text{Π})\times \cdots \times P(x_1|Y_n,X^2,\theta ,\text{Π}).$ E4

For sampling, it is sufficient to consider the sampling of x _t from $P(x_t|Y_n,X^{t+1},\theta ,\text{Π})$ . Moreover, $P(x_t|Y_n,X^{t+1},\theta ,\text{Π})\propto P(x_t|Y_t,\theta ,\text{Π})P(x_{t+1}|x_t,\text{Π})$ . This expression has two ingredients: the first is $P(x_{t+1}|x_t,\text{Π})$ , which is the transition matrix from the Markov chain. The second is $P(x_t|Y_t,\theta ,\text{Π})$ that would be obtained by recursively starting at t = 1.

The mass function $P(x_{t-1}|Y_{t-1},\theta ,\text{Π})$ is transformed into $P(x_t|Y_t,\theta ,\text{Π})$ , which is in turn transformed into $P(x_{t+1}|Y_{t+1},\theta ,\text{Π})$ and so on. The update is as follows: for $k=\mathrm{1,}\dots ,a$ , we could write

$P(x_t=k|Y_t,\theta ,\text{Π})=\frac{P(x_t=k|Y_{t-1},\theta ,\text{Π})\pi (y_t|y_{t-1},{\theta }_k)}{{\displaystyle {\sum }_{l=1}^{a}P(x_t=l|Y_{t-1},\theta ,\text{Π})\pi (y_t|y_{t-1},{\theta }_l)}}\mathrm{.}$ E5

These calculations are initialized at t = 0, by setting $P(x_1|Y_0,\theta )$ to be the stationary distribution of the Markov chain. Precisely, the simulation proceeds for $k=\mathrm{1,}\dots ,a,$ recursively by first simulating $P(x_1=k|Y_0,\theta ),$ from the initial distribution Π₀(k) and $P(x_1=k|Y_1,\theta ,\text{Π})\propto P(x_1=k|Y_0,\theta ,\text{Π})\pi (y_1|Y_0,{\theta }_k).$ Then, we get by forward calculation $P(x_t=k|Y_{t-1},\theta )={\displaystyle {\sum }_{l=1}^a}{\text{Π}}_{lk}P(x_{t-1}=l|Y_{t-1},\theta )$ , for each $t=\mathrm{2,}\dots ,n$ , where Π _lk is the transition probability and $P(x_t=k|Y_t,\theta )\propto P(x_t=k|Y_{t-1},\theta ,P)\text{Π}(y_t|Y_{t-1},{\theta }_k).$ The last term in the forward computation $P(x_n=k|Y_n,\theta )$ would serve as a start for the backward pass, and we get recursively for each $t=n-\mathrm{1,}\dots \mathrm{,1.}$ ; $P(x_t=k|Y_n,X^{t+1},\theta )\propto P(x_t|Y_{t-1},\theta ,\text{Π})P(x_{t+1}|x_t=k,\text{Π}),$ which permits the obtention of $X_n=(x_1,\dots ,x_n)$ .

4.2.2. Simulating the hidden states for the MAR(1)HMM

Returning to our model, and adopting notations and algorithm developed by Fitzpatrick and Marchev, f will denote the observation density for the MAR(1)HMM, and for $u_0<t<M$ ; $x_{-t}=(x_{u_0},\dots ,x_t),$ $x^t=(x_t,\dots ,x_M),$ $y(t)=(y_{i,t},i=\mathrm{1,2,}\dots ,n)$ , $y_{,t}={\displaystyle \bigcup _{i:u_i<t}}\{y_{i,u_i},\dots ,y_{i,\min \{t,m_i\}}\}$ , and $y^t={\displaystyle {\displaystyle \bigcup }_{i:t<m_i}}\{y_{i,\max \{t+\mathrm{1,}{u}_i\}},\dots ,y_{i,m_i}\}$ . The posterior distribution of the hidden state could be written as: $P(x_{-M}|y_{,M},\theta )=P(x_M|y_{,M},\theta )\times \cdots \times P(x_{u_0}|y_{,M},x^{u_0+1},\theta )$ . So we could sample the whole sequence of states by sampling from $P(x_t|y_{,M},x^{t+1},\theta ).$ Hence, the estimation of the hidden states is performed recursively by first initializing

$\begin{array}{l}P(x_{u_0}|y_{,u_0},\theta )\propto P(y_{,u_0}|x_{u_0})P(x_{u_0}|r);y_{,u_0}=\{y_{i,u_i,},u_i=u_0,i=\mathrm{1,}\dots ,n\}.\\ P(x_{u_0+1}=k|y_{,u_0},\theta )={\displaystyle {\displaystyle \sum _{l=1}^a}}{\text{Π}}_{lk}P(x_{u_0}=l|Y_{,u_0});k=\mathrm{1,}\dots ,a\mathrm{.}\\ P(x_{u_0+1}=k|y_{,u_0+1},\theta )\propto P(x_{u_0+1}=k|y_{,u_0},\theta )f(y(u_0)|y_{,u_0},{\theta }_k).\end{array}$ E6

We perform a similar calculation for every state at time t, and we conclude by calculating $P(x_M=k|y_{,M-1},\theta )={\displaystyle {\sum }_{l=1}^a}{\text{Π}}_{lk}P(x_{M-1}=l|Y_{,M-1},\theta ),$ and $P(x_M=k|y_{,M},\theta )\propto P(x_M=k|y_{,M-1},\theta )f(y(M)|y_{,M-1},{\theta }_k)$ . Later on, we get $P(x_M=k|y_{,M},\theta )$ , which permits the simulation of $P(x_M|y_{,M},\theta ).$ Finally, by backward calculation, we simulate from the probabilities $P(x_t|y_{,M},x^{t+1},\theta )\propto P(x_{t+1}|x_t,\text{Π})P(x_t|y_{,t},\theta )$ for each time $t=M-\mathrm{1,}\dots ,u_0$ . Those backward probabilities would permit the simulation of the latent states.

4.3. Sampling from P(θ|x, y)

4.3.1. Sampling Π

Under the prior assumption of Dirichlet prior for each row of the transition matrix $P({\text{Π}}_i)\propto \mathbb{D}({\delta }_{i1},\dots ,{\delta }_{ia})$ , and the independence assumption between those rows, the posterior distribution for Π _i can be developed using Eq. (2) as follows: Let n _ij denote the number of single transitions from state i to state j, so

$\begin{array}{l}P({\text{Π}}_i|y,x)\propto P({\text{Π}}_i){\displaystyle {\displaystyle \prod _{t=u_0+1}^M}}{\displaystyle {\displaystyle \prod _{j=1}^a}}{\text{Π}}^{{\chi }_{\{x_{t-1},x_t\}}(i,j)}\propto P({\text{Π}}_i){\displaystyle {\displaystyle \prod _{j=1}^a}}{\text{Π}}_{ij}^{n_{ij}}\propto {\displaystyle {\displaystyle \prod _{j=1}^a}}{\text{Π}}_{ij}^{{\delta }_{ij}+n_{ij}-1}\\ \propto \mathbb{D}({\delta }_{i1}+n_{i1},\dots ,{\delta }_{ia}+n_{ia}).\end{array}$ E7

4.3.2. Sampling posterior distribution for initial distribution

Let $n_{0l}={\chi }_{x_{u_0}}(l)$ , for $l=\mathrm{1,}\dots ,a$ . Using (2), under Dirichlet prior $\mathbb{D}({\delta }_{01},\dots ,{\delta }_{0a})$ for the parameter r, we obtain $P(r|x,y)\propto P(r){\displaystyle {\displaystyle \prod _{l=1}^a}}{r}_l^{{\chi }_{\left\{x_{u_0}\right\}}(l)}\propto {\displaystyle {\displaystyle \prod _{l=1}^a}}{r}_l^{{\delta }_{0l}+n_{0l}-1}\propto \mathbb{D}({\delta }_{01}+n_{01},\dots ,{\delta }_{0a}+n_{0a}).$

4.3.3. Sampling posterior distribution for the autoregressive parameters $\mu ,\beta ,{\sigma }^2$

When a complete conditional distribution is known such as the normal distribution or beta distribution, we use the Gibbs sampler to draw the random variable. This is the case for our model. Let us define $n_{u_i}(l)={\displaystyle {\displaystyle \sum _{i=1}^n}}{\chi }_{\{x_{u_i}=l\}},n_l={\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+1}^{m_i}}}{\chi }_{\{x_t=l\}},N={\displaystyle {\displaystyle \sum _{l=1}^a}}{n}_l,n_{0l}={\chi }_{\left\{x_{u_0}\right\}}(l).$ So for $l=\mathrm{1,2,}\dots ,a;$ by supposing $\mathcal{N}({\alpha }_l,{\tau }_l)$ as prior distribution and using Eq. (2), the conditional posterior distribution of μ ^(l) is:

$\begin{array}{c}P({\mu }^{(l)}|y,x)\propto P({\mu }^{(l)}){\displaystyle {\displaystyle \prod _{i=1}^n}}{\left[\frac{1}{\sigma }\varphi \left(\frac{y_{i,u_i}-{\mu }^{(l)}}{\sigma }\right)\right]}^{{\chi }_{\left\{x_{u_i}\right\}}(l)}\times {\displaystyle {\displaystyle \prod _{i=1}^n}}{\displaystyle {\displaystyle \prod _{t=u_i+1}^{m_i}}}{\left(\frac{1}{\sigma }\varphi (\frac{y_{i,t}-{\mu }^{(l)}-{\beta }^{(l)}{y}_{i,t-1}}{\sigma })\right)}^{{\chi }_{\left\{x_t\right\}}(l)}\mathrm{.}\\ \propto \exp \frac{-1}{2}\{\frac{{({\mu }^{(l)}-{\alpha }_l)}^2}{{\tau }_l}+{\displaystyle {\displaystyle \sum _{i=\mathrm{1,}{x}_{u_i}=l}^n}}{\left(\frac{y_{i,u_i}-{\mu }^{(l)}}{\sigma }\right)}^2+\\ {\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+\mathrm{1,}{x}_t=l}^{m_i}}}{\left(\frac{y_{i,t}-{\mu }^{(l)}-{\beta }^{(l)}{y}_{i,t-1}}{\sigma }\right)}^2\}.\end{array}$ E8

then ${\mu }^{(l)}/y,x\sim \mathcal{N}({{\displaystyle \tilde{\tau }}}_l,{{\displaystyle \tilde{\alpha }}}_l)$ with inverse mean ${{\displaystyle \tilde{\tau }}}_l{}^{-1}=\frac{n_{u_i}(l)+n_l}{{\sigma }^2}+\frac{1}{{\tau }_l}$ and variance

${{\displaystyle \tilde{\alpha }}}_l={{\displaystyle \tilde{\tau }}}_l\left(\frac{{\displaystyle {\displaystyle \sum _{i=\mathrm{1,}{x}_{u_i}=l}^n}}{y}_{i,u_i}+{\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+\mathrm{1,}{x}_t=l}^{m_i}}}(y_{i,t}-{\beta }^{(l)}{y}_{i,t-1})}{{\sigma }^2}+\frac{{\alpha }_l}{{\tau }_l}\right)\mathrm{.}$ E9

For ${\beta }^{(l)},l=\mathrm{1,}\dots ,a$ , and similar to μ ^(l), $\mathcal{N}(b_l,c_l)$ was proposed as prior choice to obtain:

$P({\beta }^{(l)}|y,x)\propto P({\beta }^{(l)}){\displaystyle {\displaystyle \prod _{i=1}^n}}{\displaystyle {\displaystyle \prod _{t=u_i+1}^{m_i}}}{\left[\frac{1}{\sigma }\varphi \left(\frac{y_{i,t}-{\mu }^{(l)}-{\beta }^{(l)}{y}_{i,t-1}}{\sigma }\right)\right]}^{{\chi }_{\left\{x_t\right\}}(l)},$ E10

and therefore, ${\beta }^{(l)}/y,x\sim \mathcal{N}({{\displaystyle \tilde{c}}}_l,{{\displaystyle \tilde{b}}}_l)$ with inverse mean ${{\displaystyle \tilde{c}}}_l{}^{-1}=\frac{1}{c_l}+\frac{{\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+\mathrm{1,}{x}_t=l}^{m_i}}}{y}_{i,t-1}^2}{{\sigma }^2}$ and variance

${{\displaystyle \tilde{b}}}_l={{\displaystyle \tilde{c}}}_l\left(\frac{b_l}{c_l}+\frac{{\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+\mathrm{1,}{x}_t=l}^{m_i}}}(y_{i,t}-{\mu }^{(l)})y_{i,t-1}}{{\sigma }^2}\right)\mathrm{.}$ E11

For the posterior distribution of σ ², by supposing $\mathbb{I}\mathbb{G}\left(\epsilon ,\zeta \right)$ as prior, we deduce from Eq. (2)

$\begin{array}{c}P({\sigma }^2|y,x)\propto ({\sigma }^2{)}^{-(\epsilon +1)}\exp (-\frac{\zeta }{{\sigma }^2}){\displaystyle {\displaystyle \prod _{i=1}^n}}\left[\frac{1}{\sigma }\varphi \left(\frac{y_{i,u_i}-{\mu }^{(x_{u_i})}}{\sigma }\right)\right]\\ \times {\displaystyle {\displaystyle \prod _{i=1}^n}}{\displaystyle {\displaystyle \prod _{t=u_i+1}^{m_i}}}\left[\frac{1}{\sigma }\varphi \left(\frac{y_{i,t}-{\mu }^{(x_t)}-{\beta }^{(x_t)}{y}_{i,t-1}}{\sigma }\right)\right],\end{array}$ E12

consequently ${\sigma }^2/y,x\sim \mathbb{I}\mathbb{G}({\displaystyle \tilde{\epsilon }},{\displaystyle \tilde{\zeta }})withparameters{\displaystyle \tilde{\epsilon }}=\frac{n_{u_i}+N}{2}+\epsilon$ and

${\displaystyle \tilde{\zeta }}=\frac{{\displaystyle {\displaystyle \sum _{i=1}^n}}{(y_{i,u_i}-{\mu }^{(x_{u_i})})}^2+{\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+1}^{m_i}}}{(y_{i,t}-{\mu }^{(x_t)}-{\beta }^{(x_t)}{y}_{i,t-1})}^2}{2}+\zeta \mathrm{.}$ E13

Finally, the algorithm is ran for d = 1,…,D iterations by alternating between the following steps, where in each step we compute a conditional posterior for the given parameter:

The MCMC algorithm:

1. For $h=\mathrm{1,2,....,}a$ , give reference values for the hyperparameters α _h , τ _h , a _h , b _h , δ _0h , and δ _ih for $i=\mathrm{1,2,....,}a\mathrm{.}$

2. Initialization (Step d = 1 of the MCMC iterations): Initialize Π⁽¹⁾, r ⁽¹⁾, μ ⁽¹⁾, β ⁽¹⁾, and σ ²⁽¹⁾.

3. Simulation of the hidden states:

   1. Initialization of forward simulation: $P(x_{u_0}^{(d)}|y_{,u_0},\theta )\propto P(y_{,u_0}|x_{u_0}^{(d)})P(x_{u_0}^{(d)}|r^{(d)})$ , with $y_{,u_0}=\{y_{i,u_i,},u_i=u_0,i=\mathrm{1,\dots ,}n\}.$

   2. Forward simulation: For $k=\mathrm{1,\dots .,}a$ and $t=u_0+\mathrm{1,\dots .,}M:$ $P(x_t^{(d)}=k|y_{,t-1},\theta )={\displaystyle {\sum }_{l=1}^a}{\text{Π}}_{lk}^{(d)}P(x_{t-1}^{(d)}=l|Y_{,t-1},\theta )$ and

      $P(x_t^{(d)}=k|y_{,t},\theta )=\frac{P(x_t^{(d)}=k|y_{,t-1},{\theta }_k)f(y(t)|y_{,t-1},{\theta }_k)}{{\displaystyle {\sum }_{l=1}^a}P(x_t^d=l|y_{,t-1},\theta )f(y(t)|y_{,t-1},{\theta }_l)}$ .

   3. Initialization of backward simulation: For $k=\mathrm{1,\dots .,}a$ , given

      $P(x_M^{(d)}=k|y_{,M},\theta )$ from forward simulation, we get $P(x_M^{(d)}|y_{,M},\theta ).$

   4. Backward simulation: For $k=\mathrm{1,\dots .,}a$ and $t=M-\mathrm{1,\dots ,}{u}_0:$

      $P(x_t^{(d)}|y_{,M},x^{t+1(d)},\theta )\propto P(x_{t+1}^{(d)}|x_t^{(d)},\pi )P(x_t^{(d)}|y_{,t},\theta ).$

4. Estimation of the initial distribution and the transition distribution

   1. for $l=\mathrm{1,\dots .,}a$ , $k=\mathrm{1,\dots .,}a$ . Calculate $n_{0l}={\chi }_{\left\{x_{u_0}^{(d)}\right\}}(l)$ and $n_{kl}={\Sigma }_{t=u_0+1}^M{\chi }_{\{x_{t-1}^{(d)},x_t^{(d)}\}}(k,l).$

   2. Sample $(r_1^{(d+1)}\mathrm{,\dots .,}{r}_a^{(d+1)})\propto \mathbb{D}({\delta }_{01}+n_{01}\mathrm{,\dots .,}{\delta }_{0a}+n_{0a}).$

   3. For $i=\mathrm{1,\dots .,}a$ ; sample $({\text{Π}}_{i1}^{(d+1)}\mathrm{,\dots .,}{\text{Π}}_{ia}^{(d+1)})\propto \mathbb{D}({\delta }_{i}1+n_{i1}\mathrm{,\dots .,}{\delta }_{i}a+n_{ia}).$

5. Simulation of μ: For $l=\mathrm{1,\dots .,}a$ ,

   1. ${{\displaystyle \tilde{\tau }}}_l{}^{-1}=\frac{n_{u_i}(l)+n_l}{{\sigma }_{(d)}^2}+\frac{1}{{\tau }_l}\mathrm{.}$

   2. ${{\displaystyle \tilde{\alpha }}}_l={{\displaystyle \tilde{\tau }}}_l\left(\frac{{\displaystyle {\displaystyle \sum _{i=\mathrm{1,}{x}_{u_i}=l}^n}}{y}_{i,u_i}+{\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+\mathrm{1,}{x}_t=l}^{m_i}}}(y_{i,t}-{\beta }_{(d)}^{(l)}{y}_{i,t-1})}{{\sigma }_{(d)}^2}+\frac{{\alpha }_l}{{\tau }_l}\right)\mathrm{.}$

   3. Simulate ${\mu }_{(d+1)}^{(l)}/y,x\sim \mathcal{N}({{\displaystyle \tilde{\alpha }}}_l,{{\displaystyle \tilde{\tau }}}_l).$

6. Simulation of β: For $l=\mathrm{1,\dots .,}a$ ,

   1. ${{\displaystyle \tilde{c}}}_l{}^{-1}=\frac{1}{c_l}+\frac{{\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+\mathrm{1,}{x}_t=l}^{m_i}}}{y}_{i,t-1}^2}{{\sigma }_{(d)}^2}\mathrm{.}$

   2. ${{\displaystyle \tilde{b}}}_l={{\displaystyle \tilde{c}}}_l\left(\frac{b_l}{c_l}+\frac{{\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+\mathrm{1,}{x}_t=l}^{m_i}}}(y_{i,t}-{\mu }_{d+1}^{(l)})y_{i,t-1}}{{\sigma }_{(d)}^2}\right)$ .

   3. Simulate ${\beta }_{(d+1)}^{(l)}/y,x\sim \mathcal{N}({{\displaystyle \tilde{b}}}_l,{{\displaystyle \tilde{c}}}_l).$

7. - Simulation of σ ²:

   1. ${\displaystyle \widetilde{ϵ}}=\frac{n_{u_i}+N}{2}+ϵ$ .

   2. ${\displaystyle \tilde{\zeta }}=\frac{{\displaystyle {\displaystyle \sum _{i=1}^n}}{(y_{i,u_i}-{\mu }_{(d+1)}^{(x_{u_i})})}^2+{\displaystyle {\displaystyle \sum _{i=1}^n}}{\displaystyle {\displaystyle \sum _{t=u_i+1}^{m_i}}}{(y_{i,t}-{\mu }_{(d+1)}^{(x_t)}-{\beta }_{(d+1)}^{(x_t)}{y}_{i,t-1})}^2}{2}+\zeta$ .

   3. Simulate ${\sigma }_{(d+1)}^2/y,x\sim \mathbb{I}\mathbb{G}({\displaystyle ϵ},{\displaystyle \tilde{\zeta }})$ .

5. Simulation study

In this section, we apply our results to the breast cancer model discussed earlier. The main reason behind our work is that the progression of breast cancer cannot be seen directly unless we use observations related to the disease that could characterize its progression; those observations here are quantities which could be measured; they are called biomarkers, where the word biomarker is used to designate any objective indication of a biological process or disease condition including during treatment and should be measurable. Furthermore, biomarkers are increasingly used in the management of breast cancer patients. One example is reported in Ref. [22], stating that there is correlation between elevation of CEA and/or CA15-3 and disease progression, in breast cancer patients. Also we use the autoregressive dependence among the observations to add more dynamics to the model unlike conventional HMMs where the successive observations given the Markov process are independent. We used the classification of breast cancer in three states: local where the disease is confined within the breast, the regional phase when the lymph nodes are involved, and the distant stage where the cancer is found in other parts of the body. We restrict ourselves to these three stages unlike other stage classifications that divide the progression in more than three stages such as the TNM (tumor, node, and metastasis) system. By lack of finding data about breast cancer biomarkers, we will confine ourselves to simulate an MAR(1)HMM model for observation time M = 24, and a number of individuals n = 210, a = 3 for Markov states number, with the length observation time for each individual selected uniformly between 2 and M. The simulation process supposes we have for the autoregressive means $\mu =({\mu }^{(1)},{\mu }^{(2)},{\mu }^{(3)})=(\mathrm{12,24,36})$ , since markers such as CA15–3 increase as the disease advances toward metastatic breast cancer. In addition, CA15–3 increase rapidly between successive observations, and thus, we take in the simulation the parameters $\beta =({\beta }^{(1)},{\beta }^{(2)},{\beta }^{(3)})=(\mathrm{0.2,0.4,0.8})$ .

The algorithm of simulation works as follows:

1. For each individual $i=\mathrm{1,}\dots ,n$ , choose m _i the length of observation for that individual i.

2. Generate each discrete disease state x _t using transition matrix $\text{Π}=(\mathrm{0.7,0.2,}0.1;\mathrm{0.1,0.6,0.3};\mathrm{0.2,0.3,0.5})$ for $t=u_{0+1},\dots ,M$ .

3. Generate the observations y _i,t for all individuals using our model 1.

We choose a prior ${\sigma }^2\sim \mathbb{I}\mathbb{G}(\mathrm{0.001,0.001})$ , a $\mathbb{D}(\mathrm{1,}\dots \mathrm{,1})$ prior for each row of Π, and Gaussian noninformative priors for the μs and the βs. Having the hidden states and the observations, we ran our algorithm for 8000 MCMC iterations. MCMC algorithm convergence was assessed by analyzing MCMC iterations mixing plots that are shown in Figure 1 , autocorrelation sample graphs checking as illustrated in Figure 2 , and inspecting histograms of posterior densities for the parameters of the models in Figure 3 . All parameters show good mixing of chains, autocorrelations that decay immediately after a few lags, and perfect posterior densities fitting. Also the Gelman [23] potential scale reduction factor (PSRF) was plot. The PSRF is measured for more than two MCMC chains (three chains in this works are considered), and it is measured for each parameter of the model; it should show how the chains have forgotten their initial values and that the output from all chains is indistinguishable. It is based on a comparison of within-chain and between-chains variances and is similar to a classical analysis of variance; when the PSRF is high (perhaps greater than 1.1 or 1.2), then we should run our chains out longer to improve convergence to the stationary distribution. Each PSRF declines to 1 as the number of iterations approaches infinity to confirm convergence. All the parameters have shown a PSRF less than 1.1 as the number of iteration increases and by the way a good sign of convergence ( Figure 4 ). Moreover, we should point out that in the family of Markov switching model there is the so-called label switching problem (e.g., Ref. [24]) which arises identifiability problem, and hence, we would not estimate perfectly the parameters. In addition, the posterior densities could show evidence of multimodality. some authors postprocess the output of the MCMC to deal with the issue (e.g., [25]), while other uses a random permutation of the parameters in each iteration of the MCMC algorithm (e.g., [26]) or one can call for an invariant loss function method (e.g., [27]). In our case, no identifiability issue is noticed since we used well-separated prior hyperparameters. Even when we start from different initial values for the parameters, our algorithm converges immediately after a few iterations.

Figure 1.

Markov chain mixing for each parameter through MCMC algorithm simulation.

Figure 2.

Autocorrelation sample plots for parameters of the model.

Figure 3.

Posterior densities for the parameters of the model (after 8000 iterations).

Figure 4.

Potential scale reduction factor convergence to less than 1.02 with more iterations.

Finally and before giving our results, we should report that the simulation of the Dirichlet posterior was carried out following ([28, p. 22], [29, p. 155]) who reported that the posterior Dirichlet parameters should be simulated using the beta distribution approach. Table 1 shows how the posterior values estimated from algorithm are very close to the true ones.

Parameter	True value	Posterior statistics
		Mean	Standard deviation	Confidence interval (5%)
μ ^1	12	11.929	0.047	(11.851–12.005)
μ ^2	24	23.923	0.047	(23.847–24.000)
μ ^3	36	35.843	0.070	(35.729–35.959)
β ^1	0.2	0.2016	0.0012	(0.1997–0.2035)
β ^2	0.4	0.4018	0.0009	(0.4004–0.4032)
β ^3	0.8	0.8022	0.0010	(0.8005–0.8038)
π 11	0.7	0.688	0.068	(0.5715–0.797)
π 12	0.2	0.223	0.062	(0.129–0.332)
π 13	0.1	0.090	0.042	(0.032–0.17)
π 21	0.1	0.091	0.035	(0.041–0.154)
π 22	0.6	0.607	0.059	(0.507–0.701)
π 23	0.3	0.302	0.055	(0.214–0.397)
π 31	0.2	0.153	0.053	(0.075–0.250)
π 32	0.3	0.368	0.071	(0.257–0.488)
π 33	0.5	0.479	0.073	(0.358–0.599)
σ ^2	2	2.023	0.032	(1.970–2.077)

Table 1.

Posterior inference for the parameters of the MAR(1)HMM model.

6. Conclusion

We have extended the method of Chib [21] for block update estimation of the states to a MAR(1)HMM model. Furthermore, we would like to point out that our model can easily be extended to include missing observations, as we should only add an extra step in each MCMC iteration to estimate the missing observations. Also, we can estimate the autoregressive model for different values of the autoregressive order, p ≥ 1, by evaluating the Bayesian information criterion to select the best order that fits the observations of the model. Our model would capture the complexity and the dynamics of the evolution of breast cancer by introducing the latent states; the probabilities of transition between the latent states allow to compare among the effects of treatments on slowing or accelerating the transition of the disease from one health stage to another the autoregressive parameter mean values corresponding to different stages of the disease would guide medical doctors and scientists to monitor patients in different phases of the disease. The model incorporates individual observations with different lengths.

Last but not least, we like to mention the utilities of switching diffusion processes in addressing and analyzing many complicated applications such as in finance and risk management. Our future work would be to apply these processes to explore disease progression, because they are characterized by the coexistence of continuous dynamics and discrete events as well as their interactions.

Acknowledgments

We would like to thank the editorial staff for the comments that helped in improving this work. Also, we would like to thank the supporters of this work: The Lalla Salma Foundation Prevention and Treatment of Cancer, Rabat, Morocco; and the Germano-Morrocan Program for Scientific Research PMARS 2015-060.