Some applications of phase-type distributions in recurrent events

Multi-state models are the most common models used for the description of longitudinal survival data. A multi-state model is a model for a stochastic process, which is characterized by a set of states and the possible transitions among them [11]. The states represent different situations of the individual (healthy, diseased, etc) along a follow-up. Special multi-state models that have been widely used in biomedical applications are the three-state progressive model, the illness-death model, or the bivariate model [9]. One of the important features of multiple state models is reversible illness-death model that is the model we, in this paper, drive some formulas in a Markovian context.
In many biomedical studies, the event of interest can occur more than once over the follow-up time for patients. Such events are called recurrent events [2]. The recurrent events are repeated events, which are of the same type such as acute exacerbations in asthmatic children, seizures in epileptics, cancer recurrences, myocardial infarctions, migraine pain, and ear infections. Several statistical models have been proposed in the literature to analyze recurrent events including [14],[15],[16] and [17]. [18] describes the R package TPmsm which aims at implementing nonparametric and semiparametric estimators for the transition probabilities in progressive illness-death models. Also another R package, TP.idm is introduced in [19] that implements a novel non-Markovian estimator for the transition probability matrix in the progressive illness-death model under right censoring. [12] reviews several modeling approaches following the methodology of the multi-state models, and focusing on the estimation of several quantities such as the transition probabilities and survival probabilities.
In this paper, we use a finite-state continuous time Markov process with a single absorbing state. The state space in the Markov process is partitioned to some subspaces. Each subspace will represent a state of a human beings health state. For example, in cancer case, each subspace is interpreted as either a stage of cancer or recovery. Each subspace will include states such that the semi-Markov property can be reflected in the models. See [23] for more details.

As in our model the time to death follows a phase-type (PH) distribution, we can take advantages of the PH properties. It is good to remind that the set of PH is dense in the class of all distributions defined on the non-negative real numbers, see [7]. There are closed-form expressions for the distribution and density functions, and this also applies to the Laplace–Stieltjes transform and we are able to obtain the expected value and the all non-central moments of the PH random variable by successive differentiation of the Laplace transform.
Many authors recently have used PH distribution in their researches. [1] uses PH for modeling disability insurance. The model represents the aging process as the passage through a number of phases of decreasing vitality. When disabled, individuals additionally pass through several stages that represent duration of disability. [23] uses PH for modeling of skin cancer patients in the United States and estimate parameters related to the aging process that can be useful for comparing the physiological age processes of patients with cancer and healthy people. [8] used phase-type distribution for mortality analysis, and obtained conditional survival probabilities of the time of death and the actuarial present values of the whole life insurance and annuity. [26] used phase-type distribution for analyzing data on lengths of stay of hospital patients. [4] uses PH to model the interval time between first and second birth. A special case of PH is Coxian distribution that can be used to represent survival times in terms of phases through which an individual may progress until leaves a system, such as hospital stay or time till death. [5] used the Coxian to model the patients stay in hospital. The contribution of this paper are presented in Theorem 3.1 and by developing formulas, in the PH context, presented in examples. It can be seen that how some formulas in recurrent events context are calculated tractably in a Markov chain, by using the Markov property.
This paper is organized as follows. Section 2 provides an introduction to PH as well as the notational convention that will be used in describing our model. In Section 3, we describe the structure of our model and present the main theorem which is the main contribution of this paper. In this section two examples will be presented. Section 4 closes the paper with some concluding remarks and a brief discussion of possible future studies.

Consider a continuous-time Markov process $\{J_{t},t\geq 0\}$ with the space state, $\Gamma=\{0,1,\ldots,m\}$, $m\in\mathbb{N}$, where the state $0$ is absorbing and the rest are transient. Assume the intensity matrix and initial probability vector associated with $J$, are denoted by Q and $\boldsymbol{\beta}$, respectively. Clearly, we have that

where the matrix T is $m\times m$, sub–intensity matrix ( matrix $\mathbf{B}=(b_{ij})_{i,j=1,\dots,m}$ is called a sub–intensity matrix if $b_{ii}\leq 0$, $b_{ij}\geq 0,\textrm{and}\,\sum_{k=1}^{m}b_{ik}\leq 0$, with strick inequality for at least on $i$, for  $i,j=1,2,\dots,m;~{}i\neq j$). The vector $\textbf{t}_{0}$ is the transition rates from the transient states to the absorbing state, 0. Therefore we have $\textbf{T}\textbf{1}+\textbf{t}_{0}=\textbf{0}$, (where 0, 1 are column vectors of zeros and ones with proper dimensions, respectively). It can be easily proved that

( matrix exponential of the squared matrix $\mathbf{B}$ is defined as $e^{\mathbf{B}}=\sum_{l=0}^{\infty}\frac{\textbf{B}^{l}}{l!}$).
We also assume that $P(J_{0}\in\{0\})=\beta_{0}=0$. As a result, $\boldsymbol{\beta}=(\beta_{0},\boldsymbol{\alpha})$ where we have that $\boldsymbol{\alpha}\textbf{1}=1$.
If we define, $Y=\inf\{t;J_{t}=0\}$, then it is said that $Y$ is a PH random variable with representation $(\boldsymbol{\alpha},\textbf{T})$. For a continuous random variable $Y$, the event $Y>y$ implies that the process started from any transient state has not reached the absorbing state by time $y$, thus the distribution function of $Y$ can be obtained as follows:

Thus the survival function of $Y$ is

By taking derivatives of $F(y)$, one can obtain the probability density function of $Y$ as follows ,

The Laplace-transform and $k^{th}$ moments are given in the following.

For a complete review of PH distributions refer to [3].

Based on the structure of a representation of the PH, several classes of phase-type distributions can be distinguished. The structure of a PH representation often has an impact on its application, as some structures allow more efficient solutions. The most important distinction is the one into Acyclic and General Phase-type distributions: Every acyclic phase-type distribution has at least one Markovian representation without cycles in the sub–generator, while for general phase-type distributions cycles are allowed. An application of the general PH is in the reversible illness-death model. In this paper we consider the reversible illness-death model in a Markovian environment.

Since the Markov process has only a single absorbing state, the time of death (the time of absorption) follows a PH distribution. Note also that the time staying in recovery or disease statuses can be obtained by using the Markov property.
Furthermore, the sub-intensity matrix, T, of the Markov process in (1) is given by

where $\textbf{T}_{i,j},i\neq j,~{}i,j=1,2,...,k$ are $n\times n$ matrices containing the transition rates from the stage $i$ to $j$. The matrices $\textbf{T}_{i}$’s, $i=1,2,...,k$ are the sub-intensity matrices (with non-zero main and upper diagonal elements) for the Markov chain describing a sojourn in stage $i$. The mortality rates, $\mathbf{t}_{0}=-\textbf{T}\boldsymbol{1}$, is an $n\times k$-dimensional column vector containing the rate of death from each of the $n\times k$ alive states.
In this paper, we assume that the person at diagnose has age of $x$ and we will remove $x$ from the notations.
From now on, it as also assumed that $\boldsymbol{\alpha}=(\boldsymbol{\alpha_{1}},\boldsymbol{\alpha_{2}},\dots,\boldsymbol{\alpha_{k}})$, where $\boldsymbol{\alpha_{i}}$’s are all $1\times n$ vectors corresponding to initial probabilities to start from stage $i$ ($\boldsymbol{\alpha_{1}}$ corresponds to the recovery which can be assumed to be $\boldsymbol{0}$). The normalized and the transposed vectors of a vector $\boldsymbol{\beta}$ are denoted by $\hat{\boldsymbol{\beta}}$ and $\boldsymbol{\beta}^{\prime}$, respectively.
Now with the notations above and assumptions, by using the Markov property we are able to calculate the probability of various number of transitions until time $t$.

The descriptive statistics of the data by Age (intervals of 5-year), Year (from 1967 to 1974, intervals of 1) and Surgery (yes=1 or no=0) are presented in Tables 2, 3 and 4, respectively.

In Table 2, the average age of heart patients is 45.2 years. The percent of transplanted patients with 45-50 years old is maximum with 78% transplanting. The rate of death for transplanted patients($\frac{\#\text{death of transplanted patients}}{\#\text{ transplanted patients}}$) increases with age.

In Table 3, the number of heart patients’ enrollment is about the same for every year but the percent of transplanted patients for the fifth year of study is maximum with 84% transplanting. The mortality rate for transplanted patients has decreased over the years so it is expected that the influence of year of acceptance on hazard is negative.

The heart patients are in the illness status when the disease is diagnosed. At any time thereafter, they may become transplanted or deceased. As explained in Section 1, rates of transition from the illness state depend on the duration of disease. Therefore the model must be a semi-Markov model.
We propose a model involving $2$ stages (Illness and Transplant) as shown in Figure 2. In each stage, the patients may move from left to right in the figure as they age. After a heart transplant, patients move right from a state in the illness stage to a state in the transplant stage. Of course, patients may die while in any state, as indicated by the arrows to death state in Figure 2.
In order to describe how probabilities values are calculated, we define some notations for our model. For $t\geq 0$, let $J_{t}$ represent the status of a patient at time $t$. We use the finite-state continuous-time Markov process with $2n$ alive states: $n$ transient states for disease stage, $n$ transient states for transplant stage, and one absorbing state which is the death state. Suppose that the space state is partitioned into the set $E_{1}$ of $n$ illness states and the set $E_{2}$ of $n$ transplant states and $D$, for death.
Hence the sub-intensity matrix (4) will be reduced to

where the $\textbf{T}_{1}$ and $\textbf{T}_{2}$ are the sub-intensity (with non-zero main and upper diagonal elements) matrices for the Markov chain describing a sojourn in $E_{1}$, and $E_{2}$, respectively. The matrix $\mathbf{T}_{12}$ contains the transition rates from the illness stage to the transplant stage. Let $\boldsymbol{\alpha}$ be an 2n-dimensional column vector providing the initial state probabilities.

The elements of the proposed sub-intensity matrix T are:

where $\lambda_{0}$ is the rate of transition from one state to the next, $\lambda_{1}$ is the rate of transplant from disease. The rate of mortality for the stages decease and transplant are shown by $\mathbf{q}^{1}_{i}$ and $\mathbf{q}^{2}_{i}$, respectively.

for $i=1,2,\dots,n$. Where the constant $a$ is interpreted as a background rate that a general reflection of the living environment, $q$ is a scale parameter and $p$ is a measure of the relative impact of the occurrence of aging. By descriptive statistics from heart data, rate of death for heart patients, before transplant is higher than after. So we add the parameter $b$ to the elements of $\mathbf{q}^{1}_{i}$. The three regression coefficients, $\gamma_{1},\gamma_{2}$, and $\gamma_{3}$ are included in the model to represent Age, Year, Surgery effects.
Also it was assumed that patients at acceptance will be in the first phase of the model. That is, $\boldsymbol{\alpha}=(1,0,..,0)^{\prime}$.

We aim to estimate the parameters by the maximum likelihood method. In the heart transplant data, there are four situations for a patient: He stays in disease stage (stage 1) until the end of the study period, dies before the transplant, dies after the transplant or he is alive after transplant (stage 2) until end of the follow-up. Since the Markov process has only a single absorbing state, the time of death (the time of absorbing) follows a PH distribution. In order to estimate the parameters of the model, by using the phase-type density function (3) with representation ($\boldsymbol{\alpha},\mathbf{T}$), we can construct expressions for the contribution of the likelihood function.

1. The probability of staying in disease state for a period of time $t$:

2. The probability of staying in disease state for $s$ time units and then stay in transplant state until time $t$, will be:

3. The probability of death in disease state after $t$ time units is:

where $\textbf{t}_{10}=-\textbf{T}_{1}\textbf{1}$.

4. The probability of death at time $t$ from the transplant after staying in disease state for $s$ time units is:

where $\mathbf{t}_{20}=-\textbf{T}_{2}\textbf{1}$. We estimate the parameters in (3) by maximizing the likelihood function,

where $n_{k}$s, $k=1,2,3,4,$ represent the number of patients in each scenarios explained in (12)-(15) and $\theta=(a,~{}b,~{}q,~{}p,~{}\lambda_{0},~{}\lambda_{1},~{}\gamma_{1},~{}\gamma_{2},~{}\gamma_{3})$.
The FMINCON function in MATLAB Software has been used for finding the maximum likelihood estimates. We estimated the parameters for different numbers of $n$. Based on the results shown in Table 5, $n=3$ was chosen. Note that the our criterion is maximizing $l(\theta)=log(L(\theta))$.

The parameters values, the standard deviations and the confidence intervals at level 95% are shown in Table 6.
In the fitted model, the influence of Age at acceptance on hazard is positive, while effects of Year and Surgery are both negative as they were expected.
The standard deviation of each parameter is obtained by the bootstrap technique. The main benefit of the bootstrap is that it allows statisticians to set the confidence intervals of the parameters without having to make unreasonable assumptions. It creates multiple resamples (with replacement) from a single set of observations, and computes the effect size of interest on each of these resamples. The bootstrap resamples of the effect size can then be used to determine the $95\%$ confidence interval. See [21] for more about bootstrap techniques.

It was our interest to know about the significance of the covariates and the parameter $b$ in the model. To this end, we have conducted hypothesis testing for some null tests. In Table 7 the results are presented.

The likelihood ratio tests ($LRT$) has been used for the testing. For example, for the null hypothesis $H_{0}:\forall i,\gamma_{i}=0$ against the alternative $H_{1}:\exists i,\gamma_{i}\neq 0$ can be expressed:

which is greater than $\chi^{2}_{0.95,3}=7.8$ and therefore the null hypothesis is disproved. As seen in Table 7, we can obtain the values of $LRT$ for the model without every variable separately. $\Theta$ is the parameter space and $\Theta_{0}$ is the parameter under the null hypothesis.
For testing the null hypothesis $H_{0}:b=0$ against the alternative $H_{1}:b\neq 0$, $LRT$ is $39.06$, by comparing with $\chi^{2}_{0.95,1}=3.8$, we reject $H_{0}$. Hence, it can be said that mortality before heart transplant is greater than after transplant.

Now using the value of parameters, the probability of the number of transition until time $t$ can be calculated. These probabilities are shown in Table 8 for different Age ( $30$ and $50$), Year (years $3$ and $5$) and Surgery ($0$ and $1$). As we can seen in this table, the probability without transition decreases as $t$ increases because the patient may die or transplant after more time of staying in disease stage. The event of no transition decreases with age, increases with year and surgery.
As we said before, mortality increases with age, so probability of staying in the disease stage for older patients are reduced, and he/she is more likely to die. Therefore, it is expected that the probability of one transition increases with age. Mortality decreases with year and surgery, so the probability for no transition also increases as the year and surgery increase. The probability of one transition increases until 6 months but decreases thereafter because the probability of transplant decreases in the later months and the probability of death increases. The probability of two transitions constantly increases. It is obvious that this probability increases with age because of its positive influence on death, decreases with year and surgery because of their negative influence on death.

Given $J_{u}\in E_{i}$, the time of staying continuously in stage $i$ of disease has a PH distribution with representation $(\boldsymbol{\hat{\alpha}}^{u}_{i},\textbf{T}_{i})$, where

$\textbf{I}_{E\,E_{c}}$, $c=1,\dots,k$, is $nk\times n$ matrix with $(n(c-1)+l,l)$ entry equal to 1, $l=1,\dots,n$ and all other entries are 0, and $\textbf{1}_{c}$ is a column vector with $nk$ entries such that elements $n(c-1)+l$, $l=1,\dots,n$ equal to $1$ and all other entries $0$.

The integrand in (17) can be interpreted as the survival function of a PH distribution with representation of $(\boldsymbol{\hat{\alpha}}^{u}_{i},\textbf{T}_{i})$. Given starting from stage illness (stage 1), expected times without transition for $u=0$, are shown in Table 9 for different ages. It is expected that the 30-year-old patient spends 24 days in disease stage during one month. This expected time decrease as age increases because of transition to death state.

Finally, the probability of transition from stage $i$ to stage $j$ in a period of time $t$ can be calculated easily as follows

For Age$=30$, Year$=3$ and Surgery$=0$, and $u=0$ (18) can be seen in Table 10. The formula for the probability of transition from disease stage to transplant is greater than death. Mortality from disease and transplant stages increase over time but mortality before transplant is greater than after.

The rates of transition from one state to the next are given by

where $\lambda$ is assumed to be $0.2$.
The rates of recovery from disease stage, defined based on the degree of malignancy of the disease, are

where $\gamma$ is the coefficient of recovery; assume $\gamma=0.1$
and the rates of transition from one stage of disease to the stage before, defined based on the degree of malignancy of the disease, are given by

where $\beta$ is assumed to be $0.2$.
The rates of transition from one stage of disease to the next are

with $a=10^{-3}$, $q=10^{-6}$ and $p=4.5$.
The rates of transition from recovery stage to each of the disease stages are given by

whereas the mortality rates are given by

where

and for $l=0,1,...,4$

So other entries of T are zero except the diagonal entries which are:

Finally, the initial probability is $\boldsymbol{\alpha}^{\prime}=(\boldsymbol{\alpha}^{\prime}_{1},\boldsymbol{\alpha}^{\prime}_{2},...,\boldsymbol{\alpha}^{\prime}_{6})$ where $\boldsymbol{\alpha}_{i}$’s are $5$-dimensional column vector. After the diagnosis of cancer, first stage, $i^{\textrm{th}}$, is also known so $\boldsymbol{\alpha}^{\prime}_{i}$ is $(1,0,0,0,0)$ and other elements of $\boldsymbol{\alpha}$ will be zero.

Also we can obtain the probability of one transition from stage $i$ to stage $j$, $P[N_{ij}(t)=1]$, by using Equation (18).

As seen in Table 13, probability of one transition from the lower stages of cancer to recovery stage is higher than the higher stages, this means that the lower the number, the less the cancer has spread and the higher number of cancer has spread more. In higher stages probability of transition to recovery is near to zero. Also probability of transition from the lower stages of cancer to death is lower than the higher stages and probability of transition from higher stages to death is near to one. Other possible probabilities from one stage to another is shown in this table.

In this article, the multi-state models, that are widely used in the medical studies for many diseases, has been developed via taking advantages of phase–type distribution, to model the recurrent events. It is important for decision makers and patients to be aware of information about their illness, such as the probability of relapse, the time stay in each stage of recovery or disease, the probability of recovery, and so on. Using Markov’s properties and phase-type distributions, we present a formula for calculating the probability of the number of transitions between the stage os diseases. To calculate these probabilities, which are interdependent and defined as recursive relations, we need to use differential equations to calculate them. In this article, we used two examples to calibrate our model. Example $1$ was for 103 Stanford heart patients who have two stages: heart disease and heart transplant. The mortality rate was defined for these stages, which have several parameters. We estimate these parameters by the maximum likelihood method to estimate the matrix T. The standard deviations and the confidence intervals for the parameters are obtained by the bootstrap technique. Finally, the probability of the number of transition until time $t$, the expected time of continuously staying in each stage and the probability of transition from stage $i$ to stage $j$ in a period of time $t$ are obtained. Another example is a simulated cancer that has one recovery stage and five cancer stages. In the simulated example, which had a more complex model, similar calculations have been done.

Data availability statement: NA

In this section we prove Theorem 3.1. The Equation (5) is obvious. To prove the rest, we divide $[0,t]$ into $m$ equal subintervals, each length will be $ds=\frac{t}{m}$. First, we assume that any transition between the stages will occur in a multiplier of $ds$. Assume also that $s_{\kappa}=\kappa\,ds$, $\kappa=0,1,\dots$, then we have the followings: