A general piecewise multi-state survival model: application to breast cancer

Abstract

Multi-state models are considered in the field of survival analysis for modelling illnesses that evolve through several stages over time. Multi-state models can be developed by applying several techniques, such as non-parametric, semi-parametric and stochastic processes, particularly Markov processes. When the development of an illness is being analysed, its progression is tracked periodically. Medical reviews take place at discrete times, and a panel data analysis can be formed. In this paper, a discrete-time piecewise non-homogeneous Markov process is constructed for modelling and analysing a multi-state illness with a general number of states. The model is built, and relevant measures, such as survival function, transition probabilities, mean total times spent in a group of states and the conditional probability of state change, are determined. A likelihood function is built to estimate the parameters and the general number of cut-points included in the model. Time-dependent covariates are introduced, the results are obtained in a matrix algebraic form and the algorithms are shown. The model is applied to analyse the behaviour of breast cancer. A study of the relapse and survival times of 300 breast cancer patients who have undergone mastectomy is developed. The results of this paper are implemented computationally with MATLAB and R.

Appendices

Appendix A

The parameters of the model are estimated by a maximum likelihood function. These parameters are the matrices T_u (or parameters inside these matrices), the regression covariate vectors β^u, for u = 1,…, k and the cut-points, all of them estimated jointly. We assume that n items are observed, all beginning in state 1, and item i is observed at m_i change times, the last time being death or censorship. Given that the item is observed at change times, then for any item, the value of the covariate vector and the corresponding state is observed. Therefore, a sequence of times, states and values of the covariate vector is achieved for each item i: $0 = t_{i,1} < t_{i,2} < \cdots < t_{{i,m_{i} }}$, $1 = x_{1}^{i} , \ldots , \, x_{{m_{i} }}^{i}$ and ${\mathbf{z}}_{{l_{1} }}^{i} , \ldots ,{\mathbf{z}}_{{l_{{m_{i} }} }}^{i}$, respectively. ${\mathbf{z}}_{{l_{s} }}^{i}$ corresponds to the covariate vector for the interval that contains the time $t_{i,s}$ for item i and for $s = 1, \ldots ,m_{i}$.

We assume k − 1 unknown positive integer cut-points, c₀ = 0 < c₁ < ··· < c_k−1 < c_k = ∞. The likelihood function for estimating the parameters is given by

$$L\left( {c_{1} , \ldots ,c_{k - 1} ,{\mathbf{T}}_{u} ,{\varvec{\upbeta}}^{u} ,u = 1, \ldots ,k} \right) = \prod\limits_{i = 1}^{n} {\prod\limits_{s = 2}^{{m_{i} }} {h_{{x_{s - 1}^{i} ,x_{s}^{i} }} \left( {\left. {{\mathbf{T}}_{u} ,{\varvec{\upbeta}}^{u} ,u = 1, \ldots ,k} \right|t_{i,s - 1} ,t_{i,s} ,{\mathbf{z}}_{{l_{s - 1} }}^{i} , \ldots ,{\mathbf{z}}_{{l_{s} }}^{i} } \right)} } .$$

For the calculations, we define the intervals $I_{q} = \left[ {c_{q - 1} ,c_{q} } \right[;J_{q} = \left] {c_{q - 1} ,c_{q} } \right] ,\, \, j = 1, \ldots ,k$. Let $f_{x}^{q} \left( {t,{\mathbf{z}}_{q}^{i} ;{\mathbf{T}}_{q} ,{\varvec{\upbeta}}^{q} } \right)$ be the sojourn time probability in state x at time t calculated by using the matrix ${\mathbf{P}}_{q} \left( {{\mathbf{z}}_{q}^{i} } \right)$. Given that the state at any cut-point is known, then the factors in the likelihood function have the following expressions,

• 1. If t_i,s−1 and t_i,s belong to intervals I_j and J_j, respectively,

  $$h_{{x_{s - 1}^{i} ,x_{s}^{i} }} \left( {\left. {{\mathbf{T}}_{j} ,{\varvec{\upbeta}}^{j} } \right|t_{i,s - 1} ,t_{i,s} ,{\mathbf{z}}_{{l_{s - 1} }}^{i} , \ldots ,{\mathbf{z}}_{{l_{s} }}^{i} } \right) = f_{{x_{s - 1}^{i} }}^{j} \left( {t_{i,s} - t_{i,s - 1} - 1,{\mathbf{z}}_{j}^{i} ;{\mathbf{T}}_{j} ,{\varvec{\upbeta}}^{j} } \right)T_{{x_{s - 1}^{i} x_{s}^{i} }}^{j} \left( {{\mathbf{z}}_{j}^{i} } \right) .$$

• 2. If t_i,s−1 and t_i,s belong to interval I_j−1, J_j, respectively,

  $$\begin{aligned} h_{{x_{s - 1}^{i} ,x_{s}^{i} }} \left( {\left. {{\mathbf{T}}_{u} ,{\varvec{\upbeta}}^{u} ,u = j - 1,j} \right|t_{i,s - 1} ,t_{i,s} ,{\mathbf{z}}_{{l_{s - 1} }}^{i} , \ldots ,{\mathbf{z}}_{{l_{s} }}^{i} } \right) = & f_{{x_{s - 1}^{i} }}^{j - 1} \left( {c_{j - 1} - t_{i,s - 1} ,{\mathbf{z}}_{j - 1}^{i} ;{\mathbf{T}}_{j - 1} ,{\varvec{\upbeta}}^{j - 1} } \right) \\ & \quad \times f_{{x_{s - 1}^{i} }}^{j} \left( {t_{i,s} - c_{j - 1} - 1,{\mathbf{z}}_{j}^{i} ;{\mathbf{T}}_{j} ,{\varvec{\upbeta}}^{j} } \right)T_{{x_{s - 1}^{i} ,x_{s}^{i} }}^{j} \left( {{\mathbf{z}}_{j}^{i} } \right). \\ \end{aligned}$$

• 3. If $t_{i,s - 1} \in I_{j} \;{\text{and}}\;t_{i,s} \in J_{q} \;{\text{with}}\;q - j \ge 2$,

  $$\begin{aligned} h_{{x_{s - 1}^{i} ,x_{s}^{i} }} \left( {\left. {{\mathbf{T}}_{u} ,{\varvec{\upbeta}}^{u} ,u = j, \ldots ,q} \right|t_{i,s - 1} ,t_{i,s} ,{\mathbf{z}}_{{l_{s - 1} }}^{i} , \ldots ,{\mathbf{z}}_{{l_{s} }}^{i} } \right) = & f_{{x_{s - 1}^{i} }}^{j} \left( {c_{j} - t_{i,s - 1} ,{\mathbf{z}}_{j}^{i} ;{\mathbf{T}}_{j} ,{\varvec{\upbeta}}^{j} } \right) \\ & \quad \times \prod\limits_{u = j + 1}^{q - 1} {f_{{x_{s - 1}^{i} }}^{u} \left( {c_{u} - c_{u - 1} ,{\mathbf{z}}_{u}^{i} ;{\mathbf{T}}_{u} ,{\varvec{\upbeta}}^{u} } \right)} f_{{x_{s - 1}^{i} }}^{q} \left( {t_{i,s} - c_{q} - 1,{\mathbf{z}}_{q}^{i} ;{\mathbf{T}}_{q} ,{\varvec{\upbeta}}^{q} } \right)T_{{x_{s - 1}^{i} ,x_{s}^{i} }}^{q} \left( {{\mathbf{z}}_{q}^{i} } \right). \\ \end{aligned}$$

The likelihood function is maximized by considering several restrictions. The matrices ${\mathbf{P}}_{q}$ and ${\mathbf{P}}_{q} \left( {{\mathbf{z}}_{q}^{i} } \right)$ associated with the model should be stochastic matrices for any covariate vector ${\mathbf{z}}_{q}^{i}$. This restriction will not allow probabilities less than zero or greater than one for any values of the parameters.

Then, the cut-points are estimated, and the optimum values $c_{1} , \ldots ,c_{k - 1}$ are the values that verify

$$c_{1} , \ldots ,c_{k - 1} \in {\rm N}\,{\text{such}}\,{\text{that}}\,L\left( {c_{1} , \ldots ,c_{k - 1} ,{\hat{\mathbf{T}}}_{u}^{{c_{1} , \ldots ,c_{k - 1} }} ,{\hat{\mathbf{\beta }}}_{u}^{{c_{1} , \ldots ,c_{k - 1} }} ,u = 1, \ldots ,k} \right) = \mathop {\hbox{max} }\limits_{{v_{j} }} \left\{ {L\left( {v_{1} , \ldots ,v_{k - 1} ,{\hat{\mathbf{T}}}_{u}^{{v_{1} , \ldots ,v_{k - 1} }} ,{\hat{\mathbf{\beta }}}_{u}^{{v_{1} , \ldots ,v_{k - 1} }} ,u = 1, \ldots ,k} \right)} \right\} ,$$

subject to $0 < v_{j} < v_{j + 1} \,{\text{for}}\, \, j = 1, \ldots ,k - 2$ and $v_{k - 1} < \mathop {\hbox{max} }\limits_{i} \left\{ {t_{{i,m_{i} }} } \right\}$, where v_j belongs to the set of natural numbers for any j with the corresponding restrictions. $\left( {{\hat{\mathbf{T}}}_{u}^{{v_{1} , \ldots ,v_{k - 1} }} ,{\hat{\mathbf{\beta }}}_{u}^{{v_{1} , \ldots ,v_{k - 1} }} ,u = 1, \ldots ,k} \right)$ are the maximum likelihood estimates of $\left( {{\mathbf{T}}^{u} ,{\varvec{\upbeta}}^{u} ,u = 1, \ldots ,k} \right)$ for $\nu_{1} , \ldots ,\nu_{k - 1}$.

The likelihood function has been implemented computationally with Matlab and it is maximized by using the function fmincon of this programme. This function is used to find the minimum of a constrained nonlinear multivariable function by using the interior-point algorithm.

Appendix B

See Tables 11 and 12.

Table 11 Contingency table of observed and expected counts for the homogeneous model

Table 12 Contingency table of observed and expected counts for the piecewise model