A non-mixture cure rate model for analyzing the survival lifetimes of cardiovascular heart failure patients

Abstract

The present study considers non-mixture models based on the discrete Burr XIII distribution to model recurrent event data in the presence of a cure fraction. In this context, we provide an alternative to the standard Cox proportional hazards model using a discretized distribution to analyze lifetime data assuming a non-mixture structure for cure rates. In a Bayesian setting, the proposed methodology was considered for analyzing a real dataset from a retrospective cohort study that aimed to evaluate specific clinical conditions that affect the lifetimes of 299 heart failure patients admitted to the Institute of Cardiology and Allied Hospital – Faisalabad, Pakistan (April-December, 2015). The model validation process was addressed using the Cox-Snell residuals, which allowed us to identify the suitability of the proposed non-mixture cure rate model.

Keywords

Bayesian approach Burr XIII distribution heart failure non-mixture model survival regression analysis

1. Introduction

In medical research, many studies are developed to evaluate how some subject-specific factors impact the probability of an exposed patient experiencing the event of interest (such as death, for example). Conversely, researchers may also be interested in understanding which conditions explain the existence of an underlying cluster of right-censored patients. In this case, it would be essential to access the mechanism responsible for characterizing “cured” patients. Standard survival analysis techniques, as the Cox proportional hazards model (Cox, 1972), provide no straightforward estimation for cure rates, which we acknowledge as the primary motivation for adopting mixture or non-mixture cure fraction model. In this way, several approaches have been proposed for modeling cure rates on univariate lifetime data (Farewell, 1982; De Angelis et al., 1999; Cancho & Bolfarine, 2001; Price & Manatunga, 2001; Yu et al., 2004; Yin & Ibrahim, 2005; Lambert et al., 2006; Lu, 2010; Othus et al., 2012; Achcar et al., 2012). The bivariate setting was treated by Wienke et al. (2003) and Wienke et al. (2006); in each case, the authors have introduced model-based structures for the analysis of bivariate time-to-event data cure rates. Besides, a general multivariate cure rate model was derived, and well-discussed by Cancho et al. (2016).

Lambert (2007) describes a cure rate model that incorporates the expected (or background) mortality for each subject, and so one can estimate cure fractions when some patients die of anything other than the investigated cause. Such an approach may be advantageous in practice since the information of cure at the individual level will rarely be available. Lambert also points out that, usually, the relative survival curve stays at a plateau level after several years in oncology studies. This plateau effect occurs when the mortality rate of diseased subjects is the same as the expected mortality rate in the general population, or equivalently, the excess mortality rate is equal to zero; that is, there is a cure fraction in the target population.

According to Vahidpour (2016), there are two approaches to modeling cure fractions in lifetime data. The first one is based on the mixture cure rate model, also known as the standard cure rate model (De Angelis et al., 1999; Tsodikov et al., 2003; Lambert et al., 2006), which has been widely applied in the analysis of survival datasets with some subjects that do not experience the event of interest before the end of the study. The second approach is based on non-mixture cure rate models, which are not very popular in survival analysis (Achcar et al., 2012; Vahidpour, 2016). In this way, let us denote $T$ as a random variable representing the time until observing the event of interest in a survival study. Following Tsodikov et al. (2003), the non-mixture cure rate model defines an asymptote of the cumulative hazard function of the baseline distribution, and hence for the cure rate. In this case, the survival function (sf) of the non-mixture cure rate model is given by

$\displaystyle S(t)=\rho^{F_{0}(t)}=\exp\{\log(\rho)F_{0}(t)\},t\in\mathbb{R}_{% +},$ (1)

where $\rho\in(0,1)$ is the probability of cure and $F_{0}(t)=1-S_{0}(t)$ is the (proper) cumulative distribution function (cdf) of the baseline distribution that naturally accounts for the non-cured or susceptible subjects in the target population. From Eq. (1), one can derive the overall probability density function (pdf) as

$\displaystyle f(t)=-\log(\rho)f_{0}(t)\exp\{\log(\rho)F_{0}(t)\},t\in\mathbb{R% }_{+},$

and so the hazard function (hf) of $T$ can be expressed as

$\displaystyle h(t)=-\log(\rho)f_{0}(t),t\in\mathbb{R}_{+},$

where $f_{0}(t)$ is the pdf of the baseline distribution.

The primary goal of this paper is to explore the structure of non-mixture cure fraction models for analyzing long-term lifetimes depending on potential risk factors. Specifically, we aim to introduce a non-mixture cure fraction model based on the Discrete Burr XIII (DB-XIII) distribution (Krishna & Pundir, 2009; Akdoğan et al., 2015, Kamari et al., 2015), whose continuous counterpart is widely used in medical applications nowadays, mainly due to the great flexibility of its hazard rate curves. As pointed out by Mazucheli et al. (2018), the DB-XIII distribution preserves some properties from its origin, and so we acknowledge that it may be interesting to consider such a model to analyze discrete lifetime datasets.

Several statistical techniques are commonly used in medical research to analyze univariate survival datasets with censored observations. Beyond the Cox proportional hazard model, we mention, in this context, the Kaplan-Meier estimator (Kaplan & Meier, 1958) for empirical survival functions and the non-parametric hypothesis tests for treatments comparison, such as Wilcoxon, Log-Rank, or Cox-Mantel tests. However, these methods may not provide accurate results when analyzing discrete lifetimes whose underlying population has long-term survivors. We acknowledge this fact as another excellent motivation for considering the DB-XIII model as the baseline distribution for characterizing lifetimes in survival studies.

By using the model to be proposed, the second goal of this paper is to model data from a retrospective cohort study conducted by Ahmad et al. (2017), which consists of 299 patients with cardiovascular heart failure, who were admitted to the Institute of Cardiology and Allied Hospital – Faisalabad, Pakistan (April-December, 2015). The records include patients aged 40 years or above, having left ventricular systolic dysfunction, belonging to the New York Heart Association (NYHA) classes III and IV. We have considered data of such nature since heart failure includes coronary heart disease, diabetes, high blood pressure, and other conditions that arouse great interest in the medical field.

According to Denolin et al. (1983), heart failure is the pathophysiologic state in which an abnormality of cardiac function is responsible for the failure of the heart to pump blood at a rate commensurate with the requirements of the metabolizing tissues. It is essential to point out that not all conditions that lead to cardiovascular heart failure can be reversed, but treatments can improve the signs and the related symptoms. One of these conditions, for example, is the myocardial ischemia (Buja, 2005) that results from severe impairment of coronary blood supply and produces a spectrum of clinical syndromes. For more details of those irreversible heart conditions and the epidemiology of heart failure, the interested reader should consult (Kannel & Belanger, 1991).

This paper is organized as follows. In Section 2, we present the non-mixture DB-XIII cure rate model and its regression structure. We also describe the adopted Bayesian approach, the associated numerical procedures considered for inferential purposes, and the model validation approach based on the Cox-Snell residuals. In Section 3, we assess the empirical properties of the Bayesian estimators through an intensive Monte Carlo simulation study. In Section 4, we present an exploratory data evaluation of the cardiovascular heart failure dataset, and, in the following, the proposed cure rate model is applied for identifying (possible) risk factors for heart failure. General comments and concluding remarks are addressed in Section 5.

2. Materials and methods

2.1 Non-mixture DB-XIII cure rate model

A discrete random variable $T$ is said to follow a DB-XIII distribution if its probability mass function (pmf) can be written as

$\displaystyle P(T=t;\phi,\alpha)=\phi^{\log(1+t^{\alpha})}-\phi^{\log(1+(t+1)^% {\alpha})},t\in\mathbb{N}_{0},$ (2)

for $\alpha\in\mathbb{R}_{+}$ and $\phi\in(0,1)$ . Noticeably, Eq. (2) refers to a simplified version of the DB-XIII distribution, which is achieved by assuming that the scale parameter of the Burr XIII distribution is equal to 1 (Mazucheli et al., 2018).

A non-mixture cure rate model based on the DB-XIII distribution has baseline cdf given by

$\displaystyle F_{0}(T\leqslant t;\phi,\alpha)=P(T\leqslant t;\phi,\alpha)=1-% \phi^{\log(1+t^{\alpha})},t\in\mathbb{N}_{0},$

and so we have

$\displaystyle S(T>t;\phi,\alpha,\rho)=\exp\{\log(\rho)[1-\phi^{\log(1+t^{% \alpha})}]\},t\in\mathbb{N}_{0},$ (3)

where $\rho$ represents the probability of cure.

In order to evaluate the effect of covariates (and risk factors) on the probability of a subject being a long-term survivor, we have to specify a monotonic, invertible, and twice differentiable link function, say $g$ , in which $\rho_{i}=g^{-1}(\bm{x}_{i}^{\intercal}\bm{\beta})(i=1,\ldots,n)$ is well defined into $(0,1)$ , where $\bm{\beta}=(\beta_{0},\beta_{1},\ldots,\beta_{q})$ is the vector of $q+1$ regression coefficients and $\bm{x}_{i}^{\intercal}=(1,x_{1i},\ldots,x_{qi})$ is a vector of covariates that may include, for example, dummy variables, cross-level interactions, and polynomials. For the requested purpose, one may choose any suitable mapping $g\hskip-2.845276pt:\mathbb{R}\rightarrow(0,1)$ . The popular choice in this context is the $\log\text{it}$ link function

$\displaystyle\log\text{it}(\rho_{i})=\log\left(\frac{\rho_{i}}{1-\rho_{i}}% \right)=\bm{x}_{i}^{\intercal}\bm{\beta},$

but the users can adopt other specifications (e.g., probit, log-log, cauchit, among others).

2.2 Inference

Let $\bm{T}=(T_{1},\ldots,T_{n})$ be a random sample of size $n$ taken from Eq. (2). In practical situations, we have that each $T_{i}$ may not be completely observed as it might be subjected to right censoring. In this context, let $C_{i}$ be the censoring time of the $i$ th subject. Then, we could actually observe $T_{i}=\min\{T_{i},C_{i}\}$ and $\delta_{i}=\mathcal{I}(T_{i}<C_{i})$ , that is, $\delta_{i}=1$ if $T_{i}$ was an observed lifetime and $\delta_{i}=0$ if it was right-censored. Thus, given the vector $\bm{t}=(t_{1},\ldots,t_{n})$ of observed/censored lifetimes, the log-likelihood function of $\bm{\zeta}=(\phi,\alpha,\bm{\beta})$ can be written as

$\displaystyle\ell(\bm{\zeta};\bm{t},\bm{\delta})=\sum_{i=1}^{n}\delta_{i}\log% \{-\log[g^{-1}(\bm{x}_{i}^{\intercal}\bm{\beta})]\}{}+\sum_{i=1}^{n}\delta_{i}% \log[\phi^{\log(1+t_{i}^{\alpha})}-\phi^{\log(1+(t_{i}+1)^{\alpha})}]{}+\sum_{% i=1}^{n}\log[g^{-1}(\bm{x}_{i}^{\intercal}\bm{\beta})]-\sum_{i=1}^{n}\log[g^{-% 1}(\bm{x}_{i}^{\intercal}\bm{\beta})]\phi^{\log(1+t_{i}^{\alpha})},$

where $\bm{\delta}=(\delta_{1},\ldots,\delta_{n})$ is a vector whose components indicates whether the corresponding value on $\bm{t}$ is observed or right-censored.

In a Bayesian framework, one may assume that no specialized information is available to justify the choice of informative priors for the model parameters. In this context, we specify prior distributions such that, even for moderate sample sizes, the information provided by the data should dominate the prior information. The non-informative prior distributions adopted in this work are given by

$\displaystyle\phi\sim\text{Beta}(1,1),\alpha\sim\text{Gamma}(0.001,0.001),% \text{ and }\bm{\beta}\sim\text{Normal}_{q+1}(\bm{0},10^{2}\bm{\mathcal{I}}_{q% }),$

where $\bm{\mathcal{I}}_{q+1}$ is a identity matrix of size $q+1$ . Thus, a Bayesian approach for the non-mixture cure rate model based on the DB-XIII distribution can be considered by writing the unnormalized joint posterior distribution of the vector $\bm{\zeta}$ as

$\displaystyle\pi(\bm{\zeta};\bm{t},\bm{\delta})\propto\exp\{\ell(\bm{\zeta};% \bm{t},\bm{\delta})\}\pi(\alpha)\pi(\phi)\pi(\bm{\beta}).$ (4)

From this point of view, inferences for vector $\bm{\zeta}$ are entirely based on the marginal posterior densities, which can be obtained by integrating Eq. (4), apart from the normalizing constant, which we assume to be finite. However, deriving analytical expressions for these densities is infeasible, mainly due to the complexity of the associated log-likelihood function. In this case, we may resort to a suitable Markov Chain Monte Carlo (MCMC) methods (Gelfand & Smith, 1990; Chib & Greenberg, 1995; Achcar & Leandro, 1998) to draw pseudo-random samples from the marginal posterior densities. In this work, we have adopted the Metropolis-within-Gibbs (MwG) sampler, whose steps are described in the following.

•

Step 1: Choose an initial value $\bm{\zeta}^{(0)}$ for $\bm{\zeta}$ . Denote $\bm{\zeta}$ at the $k$ th step as $\bm{\zeta}^{(k)}$ ;

•

Step 2: Generate $\phi^{(k)}$ from $\pi(\phi;\bm{t},\bm{\delta},\alpha^{(k-1)},\bm{\beta}^{(k-1)})$ ;

•

Step 3: Generate $\alpha^{(k)}$ from $\pi(\alpha;\bm{t},\bm{\delta},\phi^{(k)},\bm{\beta}^{(k-1)})$ ;

•

Step 4: Generate $\bm{\beta}^{(k)}$ from $\pi(\bm{\beta};\bm{t},\bm{\delta},\alpha^{(k)},\phi^{(k)})$ ;

•

Step 5: Repeat Steps 2–4 $N$ times;

•

Step 6: Compute the component-wise Monte Carlo posterior estimate of $\bm{\zeta}$ as

$\displaystyle\hat{\bm{\zeta}}={\displaystyle\frac{1}{N-b}}\,\sum_{k=b}^{N-1}% \bm{\zeta}^{(k+1)},$

where $b$ is the burn-in period.

To apply the proposed Bayesian approach to fit the non-mixture DB-XIII cure rate model, we have adopted the MwG algorithm for MCMC sampling. For each generated sample, a chain with $N=$ 200,000 values was generated for each component of $\bm{\zeta}$ , considering a burn-in period of 5% of the chain’s size. To obtain pseudo-independent samples from the posterior distribution (4), one out every 100 generated values was kept, resulting in chains of size 2,000 for each parameter. All computations were performed using the R2jags package, which is available in the R environment.

2.3 Model validation

Model validation procedures play an essential role when evaluating the suitability of any fitted model. In general, residual metrics are widely used in such a context, being those measures responsible for indicating departures from the underlying model assumptions by quantifying the data variability that the fitted model did not accommodate. Assessing model suitability using residual metrics is common nowadays; however, deriving appropriate residuals is not always easy for non-normal models typically adopted for describing cure rates in survival studies. In this way, we will consider here a popular residual metric proposed by Cox and Snell (Cox & Snell, 1968), which can be straightforwardly used in our context to assess the appropriateness of the proposed model when used in the analysis of real datasets. The Cox-Snell residuals are defined by

$\displaystyle r_{i}=-\log[S(t_{i})].$

According to Cox and Snell (1968), if the obtained model fit is adequate, then the Cox-Snell residuals should follow an Exponential distribution with mean equals 1 (if $T$ has survival distribution $S(t)$ , then $-\log[S(T)]\sim\text{Exp}(1)$ . However, if a survival time is right-censored, then the corresponding Cox-Snell residual, say $r_{i}^{+}$ , is lower than $r_{i}$ , which was obtained from an uncensored observation with the same lifetime. In this way, Crowley and Hu (1977) proposed the modified Cox-Snell residuals, based on the mean and the median of an $\text{Exp}(1)$ distribution. These modified residuals were derived by assuming that the difference between the cumulative hazard functions, $H(t_{i})$ and $H(t_{i}^{+})$ , also follow $\text{Exp}(1)$ distributions. Thus, the modified Cox-Snell residuals for censored observations are defined by

$\displaystyle r_{i}^{+}=1-\log[S(t_{i})]\text{ or }r_{i}^{+}=\log(2)-\log[S(t_% {i})].$

In this work, we have adopted

$\displaystyle\kappa_{i}=\left\{\begin{array}[]{ll}-\log[S(t_{i})],&\text{ if $% t_{i}$ was observed},\\ 1-\log[S(t_{i}),&\text{ if $t_{i}$ was right-censored}.\\ \end{array}\right.$

The overall procedure for obtaining the Cox-Snell residuals is summarized in the following.

•
Step 1: Fit the adopted model;
•
Step 2: Estimate the Cox-Snell residuals, $\hat{\kappa}_{i}$ , for each observation;
•
Step 3: Estimate the empirical survival function of the Cox-Snell residuals, $\widehat{S}_{e}(\hat{\kappa}_{i})$ , via Kaplan-Meier;
•
Step 4: Compute $-\log[\widehat{S}_{e}(\hat{\kappa}_{i})]$ ;
•
Step 5: Plot $\hat{\kappa}_{i}$ vs. $-\log[\widehat{S}_{e}(\hat{\kappa}_{i})]$ .

After following those steps, one should evaluate the results carefully. In general, if the displayed dots are close to a straight line with unit slope and intercept equals 0, then the fitted model can be considered adequate.
3. Simulation study

The empirical properties of an estimator can be accessed through Monte Carlo simulations. In this way, we have conducted an intensive simulation study aiming to evaluate the performance of the proposed Bayesian approach in some specific situations. The simulation process was carry out by generating 500 pseudo-random samples of sizes $n=$ 25, 50, 75 and 100 from the proposed non-mixture DB-XIII cure rate model defined from Eq. (3), with $\phi=\alpha=0.5$ and varying the right-censored observations and cure rate proportions.

The inverse-transform method for discrete distributions (Rubinstein & Kroese, 2008) was considered in the process of generating the pseudo-random samples. As previously described, the adopted Bayesian point estimator is the posterior expected value, and here we have evaluated the performance of such an estimator by assessing its bias (B) and its mean squared error (MSE). Thus, using samples of the model parameters, one may obtain approximate Monte Carlo estimators for these measures as

$\displaystyle\widehat{\text{B}}(\hat{\theta})\approx{\displaystyle\frac{1}{500% }}\sum_{i=1}^{500}(\hat{\theta}_{i}-\theta)\text{ and }\widehat{\text{MSE}}(% \hat{\theta})\approx{\displaystyle\frac{1}{500}}\sum_{i=1}^{500}(\hat{\theta}_% {i}-\theta)^{2},$

where $\theta=\phi$ , $\alpha$ , or $\rho$ , and $\{\hat{\theta}\}_{i=1}^{500}$ is the sequence of posterior means based on chains of size 2,000. The variance of $\hat{\theta}$ can be estimated as the difference between MSE and the square of the bias. Besides, the coverage probability (CP) of the 95% highest posterior density (HPD) intervals was estimated by the approximation

$\displaystyle\widehat{\text{CP}}(\hat{\theta})\approx{\displaystyle\frac{1}{50% 0}}\sum_{i=1}^{500}\delta_{i}(\theta),$

where $\delta_{i}(\theta)$ assumes 1 if the $i$ th HPD interval contains the actual value $\theta$ and zero otherwise.

The obtained results are presented in Table 1. We have noticed in our study that the biases of $\hat{\phi}$ and $\hat{\alpha}$ are positive, and $\hat{\rho}$ has a negative bias in most cases. However, the overall accuracy of the Bayesian estimators improved with increasing sample sizes since the estimated MSEs have decreased as $n$ increased. Additionally, one can notice that the estimated CPs of the HPD intervals for $\phi$ and $\alpha$ were converging to the nominal level (95%), and the CPs of the HPD interval for $\rho$ were always greater than 95%.

4. Cardiovascular heart failure

4.1 Exploratory data analysis

As previously described, one of the goals of this work is to analyze data from a retrospective cohort study (Ahmad et al., 2017) related to 299 patients with heart failure comprising of 105 women (35%) and 194 men (65%) with more than 40 years old and had left ventricular systolic dysfunction, lying in the NYHA classes III and IV. The follow-up time was 4–285 days, with an average of 130 days. As potential risk factors to explain mortality caused by cardiovascular heart disease, the following variables were considered: age, serum sodium, serum creatinine, gender, smoking, blood pressure (BP), ejection fraction (EF), anemia, platelets, creatinine phosphokinase (CPK) and diabetes. The baseline characteristics for dead and censored patients diagnosed with cardiovascular heart disease are summarized in Table 2.

Table 1
Results of the simulation study for the non-mixture DB-XIII cure rate model

$n$	Par.	Bias	MSE	CP (%)	$n$	Par.	Bias	MSE	CP (%)
Cure rate: 70% – Censoring: 40%					Cure rate: 50% – Censoring: 50%
25	$\phi$	$\phantom{-}0.3821$	$0.1461$	$96.89$	25	$\phi$	$\phantom{-}0.3868$	$0.1497$	$97.00$
	$\alpha$	$\phantom{-}0.4710$	$0.2219$	$98.75$		$\alpha$	$\phantom{-}0.4696$	$0.2205$	$98.45$
	$\rho$	$-0.6460$	$0.4175$	$99.74$		$\rho$	$-0.4475$	$0.2004$	$98.23$
50	$\phi$	$\phantom{-}0.3666$	$0.1345$	$96.01$	50	$\phi$	$\phantom{-}0.3711$	$0.1378$	$96.75$
	$\alpha$	$\phantom{-}0.4628$	$0.2142$	$98.15$		$\alpha$	$\phantom{-}0.4604$	$0.2119$	$98.02$
	$\rho$	$-0.6277$	$0.3944$	$98.69$		$\rho$	$-0.4278$	$0.1833$	$97.69$
75	$\phi$	$\phantom{-}0.3414$	$0.1168$	$95.92$	75	$\phi$	$\phantom{-}0.3450$	$0.1192$	$96.02$
	$\alpha$	$\phantom{-}0.4471$	$0.1999$	$97.02$		$\alpha$	$\phantom{-}0.4438$	$0.1970$	$97.54$
	$\rho$	$-0.5974$	$0.3578$	$98.21$		$\rho$	$-0.3935$	$0.1557$	$96.21$
100	$\phi$	$\phantom{-}0.2914$	$0.0854$	$95.15$	100	$\phi$	$\phantom{-}0.2909$	$0.0850$	$95.70$
	$\alpha$	$\phantom{-}0.4085$	$0.1671$	$96.20$		$\alpha$	$\phantom{-}0.4010$	$0.1611$	$96.87$
	$\rho$	$-0.5389$	$0.2932$	$97.35$		$\rho$	$-0.3213$	$0.1063$	$95.35$
Cure rate: 30% – Censoring: 60%					Cure rate: 10% – Censoring: 70%
25	$\phi$	$\phantom{-}0.3764$	$0.1418$	$96.29$	25	$\phi$	$0.2633$	$0.1321$	$96.10$
	$\alpha$	$\phantom{-}0.4668$	$0.2179$	$98.75$		$\alpha$	$0.3619$	$0.1911$	$96.45$
	$\rho$	$-0.2241$	$0.0545$	$98.20$		$\rho$	$0.1831$	$0.0396$	$95.94$
50	$\phi$	$\phantom{-}0.3592$	$0.1292$	$95.79$	50	$\phi$	$0.2440$	$0.1118$	$95.75$
	$\alpha$	$\phantom{-}0.4569$	$0.2088$	$98.02$		$\alpha$	$0.3500$	$0.1826$	$96.07$
	$\rho$	$-0.1985$	$0.0399$	$97.69$		$\rho$	$0.0972$	$0.0118$	$95.69$
75	$\phi$	$\phantom{-}0.3309$	$0.1097$	$95.37$	75	$\phi$	$0.2132$	$0.0983$	$95.35$
	$\alpha$	$\phantom{-}0.4378$	$0.1918$	$97.89$		$\alpha$	$0.3287$	$0.1639$	$96.00$
	$\rho$	$-0.1539$	$0.0252$	$96.11$		$\rho$	$0.0485$	$0.0033$	$95.26$
100	$\phi$	$\phantom{-}0.2749$	$0.0760$	$95.12$	100	$\phi$	$0.1521$	$0.0640$	$95.12$
	$\alpha$	$\phantom{-}0.3902$	$0.1528$	$97.10$		$\alpha$	$0.1757$	$0.0420$	$95.08$
	$\rho$	$-0.0762$	$0.010$	$95.05$		$\rho$	$0.0161$	$0.0006$	$95.01$

Table 2

Baseline characteristics for dead and censored patients diagnosed with cardiovascular heart disease

Continuous variables			Categorical variables
Risk factor	Non-censored (96)	Censored (203)	Risk factor	Non-censored (96)	Censored (203)
	Mean $\pm$ Std. Dev.	Mean $\pm$ Std. Dev.		Count (%)	Count (%)
Creatinine	1.8358 $\pm$ 1.4686	1.1849 $\pm$ 0.6541	Male	62 (64%)	132 (65%)
Sodium	135.38 $\pm$ 5.0016	137.22 $\pm$ 3.9829	Smoking	30 (31%)	66 (32%)
CPK	670.20 $\pm$ 1316.6	540.05 $\pm$ 753.80	Diabetes	40 (42%)	85 (42%)
Age	65.215 $\pm$ 13.215	58.762 $\pm$ 10.638	BP	40 (42%)	66 (32%)
Platelets	256381 $\pm$ 98525	266657 $\pm$ 97531	Anemia	54 (56%)	83 (40%)
EF	33.469 $\pm$ 12.525	40.267 $\pm$ 10.860

Figure 1.

Kaplan-Meier plots for the OS function (left-panel) and the EF-based survival estimates (right-panel).

Figure 1 depicts the Kaplan-Meier plots of the overall survival (OS) function and the survival of patients with EF $<$ 40 and EF $\geqslant$ 40. Clinically speaking, lower EF values may indicate cardiovascular heart failure (Florea et al., 2016). In both settings, one can notice that the survival function forms a plateau, which insightfully characterize long-term survivors in the available sample. Noticeably, after a long period (approximately 250 days), the survival probability does not reach zero, which means that for a specific portion of the patients, cardiovascular heart disease does not evolve. The Kaplan-Meier estimates for the survival plateau is 57.6% for the OS and 49.6% (EF $\geqslant$ 40) against 74.4% (EF $<$ 40) for the EF-based survival functions.

The estimated empirical survival functions illustrate a typical case where one should consider cure rate models in the data analysis. Thus, in the present study, we aimed to identify which risk factors are more likely to affect the chance of cure, describing the nature of the probabilities using the non-mixture DB-XIII cure rate model, which accommodates long-term survivors. In this context, accounting for cure fraction is crucial as it allows us to evaluate the different influences of risk factors on the behavior of cured patients (evidenced by a significant odds ratio) and which risk levels may lead patients to death (evidenced by a non-significant hazard ratio). Therefore, for the $i$ th $(i=1,\ldots,299)$ subject, the adopted regression structure for analyzing the covariates’ effect in the probability of cure is given by

$\displaystyle\log\left({\displaystyle\frac{\rho_{i}}{1-\rho_{i}}}\right)=\beta% _{0}+\sum_{j=1}^{11}\beta_{j}x_{ji},$ (5)

where

•

$x_{1i}$ : Gender (0: female; 1: male);

•

$x_{2i}$ : Smoking Habit (0: no; 1: yes);

•

$x_{3i}$ : Diabetes (0: no; 1: yes);

•

$x_{4i}$ : BP (0: low; 1: high);

•

$x_{5i}$ : Anemia (0: no; 1: yes);

•

$x_{6i}$ : Creatinine (0: $\leqslant$ 1.5; 1: $>$ 1.5);

•

$x_{7i}$ : Age (continuous);

•

$x_{8i}$ : EF (continuous);

•

$x_{9i}$ : Sodium (continuous);

•

$x_{10i}$ : Pletelets (continuous – logarithmic scale);

•

$x_{11i}$ : CPK (continuous – logarithmic scale).

4.2 Results and discussion

This section is dedicated to presenting and discussing the main results obtained from the fit of the non-mixture DB-XIII cure rate model adopted for analyzing the cardiovascular heart disease dataset, using the regression structure provided by Eq. (5). Before fitting the proposed model, we were concerned about the multicollinearity issue, so we have performed a preliminary analysis to evaluate whether the covariates are highly correlated to each other. Thus, in Fig. 2 we present a graphical display of the covariates’ correlation matrix based on the Spearman rank coefficient. One can notice that Gender and Smoking Habit are the most correlated covariates, but the coefficient’ magnitude ( $\approx$ 0.446) does not suggest such a high association degree. Besides, we verified that none of the eigenvalues of the correlation matrix are close to zero (the lower is $\approx$ 0.523), which would indicate highly correlated variables.

Table 3
Parameter estimates and 95% HPD intervals from the fitted non-mixture DB-XIII cure rate regression model

Parameter	Mean	Std. Dev.	95% HPD interval	Hazard ratio
$\beta_{0}$ (Intercept)	$-$ 0.0821	3.1532	( $-$ 6.1793, 6.1813)	0.9212
$\beta_{1}$ (Gender – Male)	1.6262	2.8236	( $-$ 3.9306, 7.2808)	5.0847
$\beta_{2}$ (Smoking – Yes)	1.1013	2.8552	( $-$ 4.9610, 6.4169)	3.0081
$\beta_{3}$ (Diabetes – Yes)	0.0931	2.6313	( $-$ 5.4754, 4.9251)	0.9111
$\beta_{4}$ (BP – High)	1.2406	2.7023	( $-$ 3.9963, 6.8367)	3.4576
$\beta_{5}$ (Anemia – Yes)	1.6302	2.7442	( $-$ 3.7180, 7.1356)	5.6302
$\beta_{6}$ (Creatinine $>$ 1.5)	$-$ 4.5238	3.1297	( $-$ 10.4291, 1.6493)	0.0108
$\beta_{7}$ (Age) ${}^{*}$	1.7459	1.4566	(0.1713, 4.9503)	5.7311
$\beta_{8}$ (EF) ${}^{*}$	$-$ 0.9912	0.8712	( $-$ 2.8232, $-$ 0.0543)	0.3711
$\beta_{9}$ (Sodium)	$-$ 0.3115	0.3966	( $-$ 1.1408, 0.4059)	0.7323
$\beta_{10}$ (log (Pletelets))	$-$ 0.6893	2.7167	( $-$ 6.0090, 4.9000)	0.5019
$\beta_{11}$ (log (CPK))	$-$ 2.8531	2.3598	( $-$ 7.8892, 1.6788)	0.0577
$\phi$	2.1609	0.8652	(0.9414, 3.8532)	–
$\alpha$	0.9987	0.0011	(0.9965, 0.9999)	–

${}^{*}$ Significant at 5%.

Figure 2.

Graphical display of the covariates’ Spearman Rank correlation matrix.

After running the MwG algorithm and performing tests for diagnosing convergence, the stationarity of the generated chains (after burn-in) was revealed. The obtained Bayesian results are summarized in Table 3, where the posterior parameter estimates and the 95% highest posterior density (HPD) intervals are presented. From the results displayed in Table 3, one can notice that age and EF are statistically significant risk factors for developing cardiovascular heart disease. In this setting, the obtained results suggest that increasing age represents a higher risk for the disease since the estimated hazard ratio is 5.731 $\gg$ 1. On the other hand, higher EF values indicate lower risk since $\hat{\beta}_{8}$ is negative, and the hazard ratio is below 1. Figure 3 (Half-Normal probability plot) suggests that the model fit is adequate since all the estimated Cox-Snell residuals are lying within the simulated envelope, which also indicates that there is no severe violation of model assumptions. Additionally, one can notice that the proposed non-mixture cure rate model fits considerably well since the estimated residuals are close to the dashed median line.

Figure 3.

Half-Normal plot with simulated envelope for the Cox-Snell residuals.

Our findings are corroborated by Solomon et al. (2005) that investigated the influence of ejection fraction on cardiovascular upshots in heart failure patients. In their study, the authors showed that left ventricular ejection fraction (LVEF) is a powerful predictor for heart failure across a broad spectrum of ventricular functions; however, such a measure does not further assess cardiovascular risks of heart failure patients with EF $>$ 45. Based on a study of the prevalence of right ventricular (RV) dysfunction in patients with preserved LVEF comparing RV functions of patients with preserved EF and reduced EF, Puwanant et al. (2009) concluded that the prevalence of RV systolic and diastolic dysfunctions were not uncommon in patients with heart failure with preserved EF. However, patients with preserved EF presented a milder degree of RV systolic and diastolic dysfunctions than those with reduced EF.

Watson et al. (2015) used MicroRNA signatures to differentiate preserved from reduced EF, independent of an echocardiography exam, which constitutes a significant challenge for the medical community heart studies. On the other hand, Reddy et al. (2018) developed noninvasive diagnostic criteria for preserved EF on euvolemic patients with dyspnea since diagnosis of cardiovascular heart failure is admittedly challenging in these cases.

Finally, compared to the results obtained by Ahmad et al. (2017) using the standard Cox proportional hazards model, our results are notably different when it comes to the effect of the risk factors. In the Cox model, we have the following covariates as statistically significant: age, EF, serum creatinine, serum sodium, and anemia. Ahmad et al. (2017) interpreted the results regarding the increasing age similarly to our’s previous conclusion. As for EF, they have pointed out that the hazard rate between patients with EF $\leqslant$ 30 was 67% and 59% higher than patients with 30 $<$ EF $\leqslant$ 45 and EF $>$ 45 respectively. For serum creatinine, the death risk was more than double for a unit increase of such chemical compound; for serum sodium, they noticed that one unit (meq/L) decreases the hazard by 6%; and, for anemia, they concluded that an anemic patient had 76% more chances of death when compared to a non-anemic one. Nonetheless, we may consider that, as the standard Cox model has no direct form to incorporate the cure rate parameters, the obtained results might not adequately indicate the statistically significant risk factors related to the occurrence of cardiovascular heart disease.

5. Concluding remarks

This work introduced a non-mixture cure fraction model based on the DB-XIII distribution as an alternative to the standard Cox proportional hazard model. The proposed model was considered to analyze a medical survival dataset related to long-term heart failure patients with left ventricular systolic dysfunction. We found promising results in discovering the significant risk factors affecting the survival times and the subjects’ cure rate.

We have considered a fully Bayesian approach for model estimation. The adopted method was based on the MwG algorithm for sampling pseudo-random values from the posterior distribution of model parameters. We acknowledge the remarkable computational simplicity of estimating cure rate models baseline by the DB-XIII distribution. Furthermore, Bayesian techniques allow incorporating prior information from experts, leading to much more insightful inferential results. Besides, such methods also provide straightforward interpretations based on the estimated model parameters, which is a prominent concern in medical applications.

Using the standard Cox model for the same dataset, we found the results obtained by Ahmad et al. (2017) quite different from ours. This fact is directly related to adopting a regression structure for cure fractions based on a flexible discrete distribution. Our study identified only two significant factors for cardiovascular heart failure: age and EF. Noticeably, the observed differences highlight the importance of choosing adequate statistical models, especially in clinical studies, as the widely used Cox proportional hazards model does not accommodate cure fractions. In general, accounting for cure rates may provide much more accurate model-based results when the empirical survival function forms a plateau, so indicating the existence of long-term survivors.

Footnotes

Acknowledgments

We thank the Co-Editor-in-Chief and the anonymous referee for the careful reading and thoughtful suggestions for improving this work’s content.

References

Achcar

J.A.

Coelho-Barros

E.A.

, & Mazucheli

(2012). Cure fraction models using mixture and non-mixture models. Tatra Mountains Mathematical Publications, 51(1), 1-9.

Achcar

J.A.

, & Leandro

R.A.

(1998). Use of Markov Chain Monte Carlo methods in a Bayesian analysis of the Block and Basu bivariate exponential distribution. Annals of the Institute of Statistical Mathematics, 50(3), 403-416.

Ahmad

Munir

Bhatti

S.H.

Aftab

, & Raza

M.A.

(2017). Survival analysis of heart failure patients: A case study. PloS One, 12(7), e0181001.

Akdoğan

Kuş

, & Kınacı

(2015). Point estimation of parameters in discrete Burr distribution based on type I censored sample. İstatistik Türk İstatistik Derneği Dergisi, 7(3), 80-86.

Buja

L.M.

(2005). Myocardial ischemia and reperfusion injury. Cardiovascular Pathology, 14(4), 170-175.

Cancho

V.G.

, & Bolfarine

(2001). Modeling the presence of immunes by using the exponentiated-Weibull model. Journal of Applied Statistics, 28(6), 659-671.

Cancho

V.G.

Dey

D.K.

, & Louzada

(2016). Unified multivariate survival model with a surviving fraction: An application to a Brazilian customer churn data. Journal of Applied Statistics, 43(3), 572-584.

Chib

, & Greenberg

(1995). Understanding the Metropolis-Hastings algorithm. The American Statistician, 49(4), 327-335.

Cox

D.R.

(1972). Regression models and life tables. Journal of the Royal Statistical Society, Series B, 34, 187-220.

10.

Cox

D.R.

, & Snell

E.J.

(1968). A general definition of residuals. Journal of the Royal Statistical Society. Series B (Methodological), 30(2), 248-275.

11.

Crowley

, & Hu

(1977). Covariance analysis of heart transplant survival data. Journal of the American Statistical Association, 72, 27-36.

12.

De Angelis

Capocaccia

Hakulinen

Soderman

, & Verdecchia

(1999). Mixture models for cancer survival analysis: Application to population-based data with covariates. Statistics in Medicine, 18(4), 441-454.

13.

Denolin

Kuhn

Krayenbuehl

Loogen

, & Reale

(1983). The defintion of heart failure. European Heart Journal, 4(7), 445-448.

14.

Farewell

V.T.

(1982). The use of mixture models for the analysis of survival data with long-term survivors. Biometrics, 1041-1046.

15.

Florea

V.G.

Rector

T.S.

Anand

I.S.

, & Cohn

J.N.

(2016). Heart failure with improved ejection fraction: Clinical characteristics, correlates of recovery, and survival: Results from the Valsartan heart failure trial. Circulation: Heart Failure, 9(7), e003123.

16.

Gelfand

A.E.

, & Smith

A.F.

(1990). Sampling-based approaches to calculating marginal densities. Journal of the American Statistical Association, 85(410), 398-409.

17.

Kamari

Bevrani

, & Kamary

(2015). Bayesian estimate of discrete Burr distribution with two parameters. Research & Reviews: Journal of Statistics and Mathematical Sciences, 1(2), 62-68.

18.

Kannel

W.B.

, & Belanger

A.J.

(1991). Epidemiology of heart failure. American Heart Journal, 121(3), 951-957.

19.

Kaplan

E.L.

, & Meier

(1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457-481.

20.

Krishna

, & Pundir

P.S.

(2009). Discrete Burr and discrete Pareto distributions. Statistical Methodology, 6(2), 177-188.

21.

Lambert

P.C.

(2007). Modeling of the cure fraction in survival studies. Stata Journal, 7(3), 351.

22.

Lambert

P.C.

Thompson

J.R.

Weston

C.L.

, & Dickman

P.W.

(2006). Estimating and modeling the cure fraction in population-based cancer survival analysis. Biostatistics, 8(3), 576-594.

23.

(2010). Efficient estimation for an accelerated failure time model with a cure fraction. Statistica Sinica, 20, 661.

24.

Mazucheli

Oliveira

R.P.

Peralta

, & Emanuelli

I.P.

(2018). Application of discrete Burr XII distribution in the analysis of animal production data. Ciência e Natura, 40, 25.

25.

Othus

Barlogie

LeBlanc

M.L.

, & Crowley

J.J.

(2012). Cure models as a useful statistical tool for analyzing survival. Clinical Cancer Research, 18(14), 3731-3736.

26.

Price

D.L.

, & Manatunga

A.K.

(2001). Modelling survival data with a cured fraction using frailty models. Statistics in Medicine, 20(9-10), 1515-1527.

27.

Puwanant

Priester

T.C.

Mookadam

Bruce

C.J.

Redfield

M.M.

, & Chandrasekaran

(2009). Right ventricular function in patients with preserved and reduced ejection fraction heart failure. European Journal of Echocardiography, 10(6), 733-737.

28.

Reddy

Y.N.

Carter

R.E.

Obokata

Redfield

M.M.

, & Borlaug

B.A.

(2018). A simple, evidence-based approach to help guide diagnosis of heart failure with preserved ejection fraction. Circulation, 138(9), 861-870.

29.

Rubinstein

R.Y.

, & Kroese

D.P.

(2008). Simulation and the Monte Carlo method. Wiley Series in Probability and Statistics. John Wiley & Sons, Hoboken, NJ, 2nd edition.

30.

Solomon

S.D.

Anavekar

Skali

McMurray

J.J.

Swedberg

Yusuf

Granger

C.B.

Michelson

E.L.

Wang

Pocock

, et al. (2005). Influence of ejection fraction on cardiovascular outcomes in a broad spectrum of heart failure patients. Circulation, 112(24), 3738-3744.

31.

Tsodikov

Ibrahim

, & Yakovlev

(2003). Estimating cure rates from survival data: An alternative to two-component mixture models. Journal of the American Statistical Association, 98(464), 1063-1078.

32.

Vahidpour

(2016). Cure Rate Models. PhD thesis, École Polytechnique de Montréal.

33.

Watson

C.J.

Gupta

S.K.

O’Connell

Thum

Glezeva

Fendrich

Gallagher

Ledwidge

Grote-Levi

McDonald

, et al. (2015). MicroRNA signatures differentiate preserved from reduced ejection fraction heart failure. European Journal of Heart Failure, 17(4), 405-415.

34.

Wienke

Lichtenstein

, & Yashin

A.I.

(2003). A bivariate frailty model with a cure fraction for modeling familial correlations in diseases. Biometrics, 59(4), 1178-1183.

35.

Wienke

Locatelli

, & Yashin

A.I.

(2006). The modelling of a cure fraction in bivariate time-to-event data. Austrian Journal of Statistics, 35(1), 67-76.

36.

Yin

, & Ibrahim

J.G.

(2005). Cure rate models: A unified approach. Canadian Journal of Statistics, 33(4), 559-570.

37.

Tiwari

R.C.

Cronin

K.A.

, & Feuer

E.J.

(2004). Cure fraction estimation from the mixture cure models for grouped survival data. Statistics in Medicine, 23(11), 1733-1747.

A non-mixture cure rate model for analyzing the survival lifetimes of cardiovascular heart failure patients

Abstract

Keywords

1. Introduction

2.1 Non-mixture DB-XIII cure rate model

4. Cardiovascular heart failure

4.1 Exploratory data analysis

Table 1 Results of the simulation study for the non-mixture DB-XIII cure rate model

Table 3 Parameter estimates and 95% HPD intervals from the fitted non-mixture DB-XIII cure rate regression model

Footnotes

Acknowledgments

References

Table 1
Results of the simulation study for the non-mixture DB-XIII cure rate model

Table 3
Parameter estimates and 95% HPD intervals from the fitted non-mixture DB-XIII cure rate regression model