Estimating Time-to-Death and Determining Risk Predictors for Heart Failure Patients: Bayesian AFT Shared Frailty Models Using INLA Method

1. Introduction

Heart failure (HF) is becoming one of the major health problem across the globe [1]. It is a leading cause of death and poses a rapidly growing public health challenge, currently affecting over 64 million individuals [2,3,4]. With the prevalence of HF increasing with age, the number of patients is rising swiftly in both developed and developing countries with aging populations [5,6]. Additionally, HF is a major contributor to cardiovascular diseases and imposes a significant economic burden on individuals and healthcare systems worldwide [3,7]. A prospective cohort study of 5,823 patients revealed an overall mortality rate of 16.5%, with the highest rates in Africa (34%) and India (23%), intermediate rates in Southeast Asia (15%), and the lowest rates in China (7%), South America (9%), and the Middle East (9%) [8].

Heart failure (HF) poses a significant public health challenge, particularly in developing countries where mortality rates are alarmingly high [9,10,11]. In sub-Saharan Africa (SSA), including Ethiopia, HF is increasingly a leading cause of death, exacerbated by inadequate resources and weak healthcare systems [9,12,13]. Unlike in the US and Europe, HF in SSA often affects younger individuals due to preventable conditions, such as hypertension (43.9%), rheumatic heart disease (15%), and dilated cardiomyopathy (19.5%) [14,15]. In Ethiopia, HF-related hospital deaths account for 24.4% of total fatalities [16,17]. Despite this, there is a notable gap in research regarding the survival rates and duration of HF patients [16,17,18,19]. While advanced survival models like Cox regression [20,21,22] and Bayesian approaches have been employed to analyze HF survival distributions [23], these methods often fall short in addressing heterogeneity and within-cluster dependencies.

Various survival methods, including Cox regression and Accelerated Failure Time (AFT) models such as exponential, Weibull, log-normal, and log-logistic, are used to estimate median survival times and identify prognostic factors [23,24]. While parametric AFT models can offer a more accurate representation of survival data compared to Cox-PH models when the hazard shape is known [25,26], they often lack flexibility in capturing unobserved heterogeneity and associations. Most studies fail to incorporate random effects or frailty terms into these models, which limits their ability to account for cluster-level variability. Frailty models address this limitation by incorporating random effects to handle unobserved heterogeneity and are an extension of the Cox-PH model [23,27,28]. Shared frailty models, where frailties are common within clusters, are particularly suitable for clustered survival time data, offering a more nuanced understanding of survival dynamics within groups [26,29,30,31].

Traditional survival methods, such as Cox regression and parametric AFT models, often struggle with heterogeneity and within-cluster dependence, leading to potential inaccuracies. These methods may also suffer from computational inefficiencies and convergence issues, particularly with Markov-Chain Monte-Carlo (MCMC) approaches [32,33,34].

In contrast, Bayesian analysis with parametric AFT models [35,36], enhanced by the Integrated Nested Laplace Approximations (INLA), overcomes these limitations. The Bayesian approach effectively integrates prior knowledge (prior probability) with data (likelihood), improving model accuracy and handling unobserved heterogeneity [37]. INLA offers faster and more accurate computation than MCMC, making it highly efficient for complex survival data analysis [32]. Thus, Bayesian AFT models with INLA provide a more robust and efficient framework for analyzing survival times [35,38,39], addressing the drawbacks of traditional methods.

In this study, we employed Bayesian Accelerated Failure Time (AFT) shared gamma frailty models using the Integrated Nested Laplace Approximations (INLA) method to estimate the survival time of HF patients. Unlike to the classical survival approaches, our method enjoys the advantage of INLA especially in largely reducing the issue of time complexity. We chose this approach due to the increasing prevalence of HF in hospital settings and the gaps identified in various studies. The proposed method better estimates the parameters and captures the heterogeneity and dependence within clusters, ensuring a more accurate and comprehensive analysis and leading to more effective interventions.

The results from this study can significantly aid medical professionals in identifying high-risk patients and adjusting their treatment plans. Analyzing the survival time of HF patients is crucial for addressing public health issues by identifying factors associated with the survival time in time-to-death based analysis. Additionally, the findings can raise awareness about the factors that contribute to HF patient deaths.

This study can also provide valuable scientific information to the Ministry of Health in Ethiopia, helping policymakers to educate the public about factors that increase the probability of death due to HF. Heart failure is preventable and curable if screened and treated in its early stages with appropriate treatment. By enhancing awareness of these aspects, the study ultimately contributes to reducing mortality rates and improving overall public health outcomes.

Overall, the major contribution of this study are as follows:

1) New Data: This study applies Bayesian AFT shared frailty models using newly compiled data that incorporate recent developments in survival methods for estimating survival time and identifying predictors in heart failure patients. This new data can be public available to supplement the development of other methods and this study would contribute the quality of survival data analysis and improves parameter estimation.

2) New Frailty Model: Unlike the classical survival methods, the proposed frailty model captures unobserved heterogeneity within clusters and further leverages its advantages.

3) Novel INLA Implementation: Addressing the highly complexity of computation, this study makes use of the advanced INLA implementation to accurately estimate the parameters of the survival models.

4) Improved Parameter Estimation: Integrating Bayesian methods into Accelerated Failure Time (AFT) models improves parameter estimation during the heart failure patient data by incorporating prior knowledge and providing a full probability distribution, further enhancing the accuracy and reliability of results. It also offers flexibility in handling random effects and frailty terms, performs better with small samples, and provides a more comprehensive analysis of complex survival data. With this advantage, the proposed method enjoys in handling the issue of heterogeneity and survival analysis based on the heart failure patients.

5) Model Flexibility: The log-normal Bayesian AFT shared gamma frailty model excels in flexibly modeling the new survival data with both individual and group-level random effects, effectively handling unobserved heterogeneity and providing robust estimates. Our results reveal that it outperforms state-of-the-art methods.

The rest of this paper is structured as follows. Section 2 describes the nature of data and methods and Section 3 presents and discusses the data analysis results. Finally, we provide concluding remarks in Section 4.

2. Data and Methods

2.1. Heart Failure Dataset and Study Design

The retrospective study, which was conducted at Jimma University Medical Center (JUMC), Jimma, Ethiopia, focused on collecting data from heart failure (HF) patients over a period of 3 years and 6 months, specifically from July 1, 2015, to January 1, 2019. The primary objective was to gather comprehensive information on these patients by reviewing records from the registration logbook and patient registration cards. Any incomplete data encountered during this process was verified against patient files and excluded if inadequate for analysis. The study aimed to include all HF patients who had received treatment at least once during the study period. Patients were considered eligible if they had complete information on the study variables in the registration book or patient card. Those with insufficient data or who did not commence treatment at the hospital were excluded, as were patients who died due to accidents. Ethical permission was obtained from JUMC prior to data collection. The data extraction was carried out by a team that included a statistician, public health officer, and trained enumerator using a predefined checklist (data extraction form). This rigorous approach ensured that the data collected were reliable and suitable for subsequent analysis. Jimma University Medical Center is a pioneering and great hospital, providing essential services to many patients. This study at JUMC aimed to contribute valuable insights into the management and outcomes of HF patients within the institution during the specified time-frame.

2.2. Study Variables and Definitions

Outcome variable

The response variable was time-to-death of HF patients from the time of diagnosis (measured in months). In this study, we assume that the data follow the right-censored, where patients who die as a result of HF are said to have experienced the event and all others (death not related to HF, dropouts, stopped, and lost to follow-up) are said to be censored (have not experienced the event of interest). The status variable was coded as 0 for censored and 1 for death. Let $Y_{ij}$ be a binary indicator variable such that:

$$Y_{i}j=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{i}^{th}\mathrm{patient}{j}^{th}\mathrm{residence}\left(\mathrm{censored}\right)\hfill \\ 1\hfill & \mathrm{if}{i}^{th}\mathrm{patient}{j}^{th}\mathrm{residence}\left(\mathrm{death}\right)\hfill \end{array}\right.$$

(1)

Predictor variables

The covariates that are considered in this study were as follows: sex, age, smoking cigarettes, diabetes mellitus, hypertension, place of residence, alcohol consumption, history of HF, anemia, treatments taken, etiology of HF, chronic kidney disease, and stages of HF.

The findings of the study will be intended to guide future healthcare interventions and to be well prepared for reduce the influence of heart failure patients.

2.3. Survival Data Analysis

The description of survival data analysis employs non-parametric methods to compare the survival functions of two or more groups. In non-parametric survival analysis, the data are summarized via estimates of the survival function and hazard function. And Kaplan-Meier plots are employed for this purpose [40]. When comparing groups, starting with a graphical display of the data for each group is helpful. A frequency distribution table was employed to summarize the data obtained from the patient registration book based on the study variables at Jimma University Medical Center.

Survival data are censored in the sense that they did not provide complete information since subjects of the study may not have experienced the event of interest [41,67]. Survival analysis is well suited for HF datasets which are very common in medical research since studies in medical areas have a special feature that follow-up studies could start at a certain observation time and could end before all experimental units had experienced an event [26]. The Cox PH model is widely used in survival regression analysis [42]. However, Wei [43] and Cox [44] have highlighted that the AFT model can provide a more realistic interpretation and alternative to the Cox PH model [43,45]. The AFT model is an alternative to the Cox PH and parametric models when the assumptions of the Cox PH model do not hold [46]. Finally, the key differences between the Cox PH and the AFT models are baseline hazard function and ways of estimating coefficients [47].

2.4. Shared Frailty AFT Models

The term `frailty’ was first introduced by [48] for univariate data and later extended to multivariate data by [48,49]. A natural extension of the univariate frailty model would be a multivariate frailty model where allows the individuals in the same cluster to share the same frailty. To account for unobserved heterogeneity, Hougaard [50] introduced the frailty term as an extension of the Cox PH model [50]. Frailty models extend traditional survival models by incorporating random effects (frailties) that act multiplicatively on the baseline hazard function [26,29,51]. The concept of frailty provides a convenient way to introduce random effects in the survival model to account for association and unobserved heterogeneity (frailty) [52]. A shared frailty model is a survival model that includes a random effect term, where the frailties are common or shared among individuals within specific groups [26,53,54]. Shared frailty model is a conditional independence model in which frailty is common to all subjects in a cluster [29,31]. AFT frailty models extend the standard AFT model by incorporating random effects [55,56].

Suppose that there are n clusters and that the $i^{\mathrm{th}}$ cluster has $n_i$ individuals and is associated with an unobserved frailty $u_i$. The log-linear representation of the shared frailty AFT model for the random variable, $T_{ij}$, that is associated with the survival time of the $j^{\mathrm{th}}$ individual in the $i^{\mathrm{th}}$ group or cluster (patient residence). Frailties or random effects $u_i$ are assumed to be identically and independently distributed random variables with a common density function $f(u,\eta )$, where $\eta$ is a parameter of the frailty distribution. Shared frailty models enable individuals within the same group to share a common frailty [31]. This study focused on HF patients in different parts of Ethiopia (urban and rural). The shared frailty model was applied based on their residence. The shared frailty AFT model can be written as follows [55]:

$$log{T}_{ij}=\mu +x_{ij}^{\prime }\beta +u_i+\delta {ϵ}_{ij},i=1,2,\dots ,n,j=1,2,\dots ,n_i$$

(2)

where $T_{ij}$ denotes the time, $\mu$ and $\delta$ are unknown location and scale parameters, $x_{ij}$ are vector of covariates, $\beta$ is a parameter vector, $u_i$ is the frailty term, and ${ϵ}_{ij}$ are independent identically distributed random errors. The conditional hazard function for a given frailty variable, u, at time $t>0$ is,

$$\lambda (t\mid x,u)=u{\lambda }_0\left(t\right)exp\left(x^{\prime }\beta \right)$$

(3)

where $\lambda \left(t\right|x,u)$ is the hazard function at time t for a subject with a set of covariates $x_1,x_2,\dots ,x_p$, ${\lambda }_0\left(t\right)$ is the baseline hazard function, and ${\beta }_1,{\beta }_2,\dots ,{\beta }_p$ are a vector of fixed effect parameters to be estimated. The conditional survival function for a given frailty u at time $t>0$ is,

$$S\left(t\right|x,u)=exp\left(-u{H}_0\left(t\right)exp\left(x^{\prime }\beta \right)\right)$$

(4)

where $S\left(t\right|x,u)$ is the survival function at time t for a subject with a set of covariates $x_1,x_2,\dots ,x_p$, $H_0\left(t\right)$ is the cumulative hazard function.

2.5. Bayesian AFT Shared Gamma Frailty Models

The Bayesian approach is preferred over the frequentist approach in survival analysis because the power of information obtained from the Bayesian approach is much better, as it combines likelihood data with prior information regarding the distribution of the parameters. The Bayesian approach is more useful in clinical data analysis than the frequentist approach and is a more appropriate data analysis technique for clinical researchers [57]. The main reasons to prefer Bayesian statistics include the ability to estimate complex models that traditional statistics cannot, a preference for the Bayesian definition of probability, the incorporation of prior and background knowledge into the analysis, and Bayesian statistics do not rely on large sample sizes (i.e., the central limit theorem). Furthermore, Bayesian statistics allow for incorporating uncertainty about a parameter and updating this knowledge via the prior distribution. The Bayesian approach considers the parameters of the model as random variables and requires specifying prior distributions for them, while the data are considered fixed [58].

Let’s start with the data likelihood function $L(y\mid \theta )$, which is a function that gives the probability of observing the sample data given the current parameters. For a set of unknown parameters in the presence of right censoring, it can be written as:

$$L(y\mid \theta )=\prod _{i=1}^n\prod _{j=1}^{n_i}\left[f{(t_{ij}\mid y_{ij};\theta )}^{{\delta }_{ij}}S{(t_{ij}\mid y_{ij};\theta )}^{1-{\delta }_{ij}}\right]$$

(5)

where ${\delta }_{ij}$ is the censoring indicator (0 = censored and 1 = death) and $f(t_{ij}\mid y_{ij};\theta )$ and $S(t_{ij}\mid y_{ij};\theta )$ are the probability density and survival distributions, respectively [59].

In the Bayesian framework, the prior distribution is a probability distribution that expresses uncertainty about the unknown parameters before incorporating the data. It represents the prior information associated with the parameter of interest. We assumed normal priors with mean and variance for all regression coefficients and a gamma prior with shape and scale parameters for the frailty parameter [35,60]. The posterior probability density function, which summarizes our beliefs about a particular parameter, is obtained using Bayes’ rule. The distribution of a parameter can be updated by combining its prior distribution with the likelihood function, resulting in the posterior density of the parameter. It is obtained by multiplying the prior distribution over all parameters, $\theta$, by the full likelihood function, $L(y\mid \theta )$ [61] which is given by

$$\mathrm{Posterior}={\displaystyle \frac{\mathrm{prior}\times \mathrm{Likelihood}}{\int \mathrm{prior}\times \mathrm{Likelihood}d\theta }}$$

(6)

So, if $L(y\mid \theta )$ is the likelihood function and $p\left(\theta \right)$ is the prior density of the parameter, then the posterior density function, $\pi (\theta \mid y)$, is proportional to the likelihood multiplied by the prior:

$$\pi (\theta \mid y)\propto p\left(\theta \right)L(y\mid \theta )$$

(7)

In our case, the posterior density function, $\pi (\theta \mid y)$, of $\theta$ is given by:

$$\pi (\theta \mid y)={\displaystyle \frac{p\left(\theta \right)\times L(y\mid \theta )}{\int p\left(\theta \right)\times L(y\mid \theta )d\theta }}$$

(8)

The quantity $m\left(y\right)=\int p\left(\theta \right)\times L(y\mid \theta )d\theta$ in the denominator is the normalizing constant of $\pi (\theta \mid y)$ and is often called the marginal distribution of the data or the prior predictive distribution.

Parametric survival models play an important role in Bayesian survival analysis since many Bayesian analyses in practice are carried out using AFT models (Exponential, Weibull, Log-Normal, Log-Logistic) [35,62,63]. We considered shared gamma frailty as the gamma distribution is commonly used for modeling the frailty effect. Gamma frailty is very well suited to failure data in terms of computational and analytical points of view. The most common reason for using the gamma distribution is its mathematical convenience [64,65]. The mean and variance are $E_{\eta }\left(u_i\right)=1$ and ${\mathrm{Var}}_{\eta }\left(u_i\right)=\eta$. Hence, the variance parameter $\eta$ represents the amount of heterogeneity in the risk of the event.

2.6. Integrated Nested Laplace Approximation Method

Survival models can be challenging to fit, especially when dealing with complex censoring schemes. Using the Gibbs sampler and other Markov-Chain Monte-Carlo (MCMC) techniques, fitting complex survival models becomes relatively straightforward, and the availability of software significantly simplifies the implementation [35]. However, MCMC methods have limitations, such as the time burden in approximating the posterior and convergence problems [33,34]. The Integrated Nested Laplace Approximation (INLA) method provides a flexible and fast technique for approximating complex Bayesian models. The Bayesian approach using the INLA method focuses on providing an accurate approximation to the posterior marginal distributions of the parameters in the model [32].

The INLA method was used to estimate the parameters in the Bayesian survival models. According to [32], INLA computes the posterior marginal for each component in the model, from which the posterior expectations and standard deviations can be derived. Additionally, INLA provides both extremely fast and highly accurate approximations to the posterior marginal via the clever use of Laplace approximations and advanced numerical methods, and it can be adapted to fit survival models [66].

2.7. Data Analysis Procedures

The entire data analysis was done using R software to estimate the heart failure patients median survival time while understanding the dataset’s unique characteristics and ensuring methodological precision. First, summary statistics was used to better present the characteristics of the data and Kaplan-Meier curves were fitted to learn the death dynamics, followed by a validation of the basic assumptions of the Cox model and fitting various survival models followed to identify the most suitable one.

For Bayesian models, the Deviance Information Criterion (DIC) was used to compare Bayesian survival models [68]. An alternative is the Watanabe-Akaike Information Criterion (WAIC) [69], which assumes a more fully Bayesian approach to construct a criterion. According to [37], the WAIC is preferable to the DIC. The preferable model is the one with the lowest values of the DIC and WAIC. A lower value indicates a better fit of the model. The most common ways of checking the goodness of fit are the Bayesian Cox-Snell residual plot and predictive distribution. Model-checking and adequacy play an important role in models for survival data. In Bayesian analysis, Chaloner [70] defined the Bayesian version of the residuals.

3. Results and Discussions

In this section, we begin by presenting descriptive results to better understand the nature of our dataset for readers and reviewers of this manuscript. Following this, we present results to clearly demonstrate the advantages of the proposed method. In doing so, we emphasize comparisons with other existing approaches in the literature. Finally, we discuss the analysis results, highlighting the benefits of the method and its broader implications in combating the potential problem of heart failure patients.

3.1. Results

The data for this study was collected from 497 patients who received treatments for HF at least once at Jimma University Medical Center in Jimma, Ethiopia, between July 2015 and January 2019.

3.1.1. Descriptive Results

As summarized in Table 1, the minimum and maximum event time observed from HF patient’s follow-up were 6 and 42 months respectively. Among those HF patients, about 60.56% were censored and the remaining 39.44% have died. Fifty percent of HF patients who received treatment survived 34 months or above.

Almost half, 53.92%, of the HF patients, were female, and the remaining were male during the follow-up study. However, the survival time of male patients seems lower. Most HF patients, about 72.03%, live in rural areas, while the remainder reside in urban areas; however, survival times appear to be shorter for those living in rural areas. The survival time of HF patients seems less as they get older. Among HF patients, 21.93% have ischemic heart disease, 22.54% have rheumatic valvular heart disease, 23.54% have cardiomyopathy, 24.14% have hypertensive heart disease, and 7.85% have other diseases. Regarding smoking status, 79.28% of HF patients were non-smokers and the remaining were smokers. Additionally, 63.38% of HF patients did not use alcohol, while 36.62% were alcohol users.

Moreover, about 18.31% of HF patients are treated in the hospital with a combination of two or more treatments, and 19.72% of HF patients take digoxin. In addition, the remaining 23.14%, 27.56%, and 11.27% of HF patients were treated with spironolactone, atorvastatin, and other treatments respectively. About diabetes, 56.94% of HF patients are non-diabetic, 15.09% have type I diabetes mellitus, and 27.97% have type II diabetes mellitus. Regarding chronic kidney disease among HF patients, approximately 31.19% have chronic kidney disease, while 68.81% do not. HF patients with chronic kidney disease generally have shorter survival times. Most HF patients, 60.16%, do not have hypertension, while the remaining patients have hypertension.

Finally, regarding the stage at which HF patients go to the hospital for treatment, approximately 35.41% are in stage IV, 25.96% in stage III, 21.73% in stage II, and 16.90% in stage I. Most HF patients who pass away, about 54.08%, HF patients go for treatment in the hospital at a later stage, and their survival times lower at this advanced stage.

3.1.2. Kaplan-Meier Estimate for Selective Covariates

The survival distributions of time-to-death of the HF patients were estimated for each group using the Kaplan-Meier (KM) method to compare the survival curves of two or more groups as seen from Figure 1.

As seen from Figure 1, it shows that the overall survival rate at the end of the first year was almost 92.6%, and the overall survival rate at the end of 42 months in this study was 19.2%, 95% confidence interval was (12.8%, 28.8%). Heart failure patients whose age was below 49 years had a higher probability of surviving than patients whose age was 49 to 65 years and those aged 65 years or older. The probability of surviving becomes less for patients aged 65 years or older. The HF patients with stage I had a higher chance of surviving than other stages. The survival curve for patients with stage I was above the survival curve for patients with stage II, stage III and stage IV. The probability of surviving becomes less for HF patients with stage IV. Additionally, HF patients who live in urban areas had a better chance of survival than HF patients who live in rural areas.

3.1.3. Bayesian AFT Shared Gamma Frailty Models

As comparison, we employed Bayesian AFT without frailty and Bayesian AFT with shared gamma frailty models for the HF dataset.

The multivariable survival models of the exponential, Weibull, lognormal, and log-logistic for the baseline hazard function by assuming gamma frailty distributions were fitted for time-to-death of HF patients. In univariable analysis at 5% level of significance, except for sex, treatment, and alcohol consumption, all variables were significant. The covariates significant in univariable analysis entered into multivariable analysis. We assume that all these coefficients have a normal prior with mean 0 and variance of 1000 in the Baysian modeling. In analyzing the frailty parameter (u) via the Bayesian approach, we assume that a gamma prior was applied to the Weibull, log-normal, and log-logistic distributions with shape parameter 1 and inverse scale parameter 0.001 [35].

Table 2 shows the analysis of the HF dataset for model comparison to compare the Bayesian AFT models with frailty with the Bayesian AFT models without frailty. To compare the efficiency of these various models, DIC and WAIC were used, and the one with the smallest value and the best fit was chosen. We noted that adding frailty to the Bayesian AFT models boosted the performance in all AFT models as compared with the Bayesian AFT models without frailty. Then, the Bayesian log-normal shared gamma frailty model (DIC = 1608.812; WAIC = 1607.838) was found to be the best model for survival time of HF patients from the given alternative.

Table 3 shows the final results for the Bayesian log-normal shared gamma frailty model, and as this result shows, the survival time of HF patients is statistically significantly affected by age, chronic kidney disease, diabetes mellitus, history of HF, aetiology of HF, hypertension, anemia, smoking, and stages of HF.

From Table 3, the results of the Bayesian log-normal shared gamma frailty model show that, keeping the effect of other factors constant, the estimated acceleration factor for HF patients aged 49 to 65 years and those aged 65 years or older is 0.7357 [95% CrI: 0.580, 0.931] and 0.6969 [95% CrI: 0.554, 0.873], respectively. Thus, the expected survival time of HF patients decreases by 26.43% and 30.31% for those aged 49 to 65 and 65 or older, respectively, compared to HF patients aged 49 or younger (reference group). The 95% CrI for the acceleration factors in both age groups does not include one, implying that both age groups significantly impact the survival time of HF patients.

For chronic kidney disease, while controlling for other factors, the estimated acceleration factor for HF patients with chronic kidney disease is 0.6630 [95% CrI: 0.569, 0.769], indicating that the expected survival time decreases by 33.70% compared to HF patients without chronic kidney disease. The 95% CrI for the acceleration factor of HF patients with chronic kidney disease does not include one, showing a significant impact on survival time.

Regarding the history of HF, keeping the effect of other factors constant, the estimated acceleration factors for HF patients in medical OPD to be 1.2751 [95% CrI: 1.046, 1.559]. Thus, the expected survival time increases by 27.51% for those in medical OPD compared to new HF patients (reference group). The 95% CrI for this acceleration factor does not include one, showing a significant impact on survival time.

For hypertension, keeping the effect of other factors constant, the estimated acceleration factor for HF patients with hypertension is 0.7320 [95% CrI: 0.626, 0.855], indicating a 26.80% decrease in expected survival time compared to HF patients without hypertension (reference group). The 95% CrI for the acceleration factor of HF patients with hypertension does not include one, implying a significant effect on survival time. Moreover, keeping the effect of other factors constant, the estimated acceleration factor for HF patients who smoked cigarettes is 0.8395 [95% CrI: 0.721, 0.976], indicating a 16.05% decrease in expected survival time compared to non-smokers. The 95% CrI for the acceleration factor of HF patients who smoked cigarettes does not include one, implying that smoking has a significant effect on survival time.

Regarding the aetiologies of heart failure and keeping the effect of other factors constant, the estimated acceleration factors for HF patients with rheumatic valvular heart disease, hypertensive heart disease, cardiomyopathy, and other heart diseases are 0.7211 [95% CrI: 0.568, 0.913], 0.7558 [95% CrI: 0.602, 0.948], 0.7416 [95% CrI: 0.586, 0.936], and 0.6750 [95% CrI: 0.487, 0.933], respectively. Thus, the expected survival time of HF patients decreases by 27.89% for those with rheumatic valvular heart disease, 24.42% for those with hypertensive heart disease, 25.84% for those with cardiomyopathy, and 32.50% for those with other heart diseases, compared to those with ischemic heart disease (reference group). The 95% CrI for these acceleration factors does not include one, implying that HF patients with these aetiologies have a significant impact on survival time.

On the other hand, for diabetes mellitus and keeping the effect of other factors constant, the estimated acceleration factors for HF patients with type I diabetes and type II diabetes are 0.6421 [95% CrI: 0.509, 0.806] and 0.5764 [95% CrI: 0.468, 0.706], respectively. Thus, the expected survival time of HF patients decreases by 35.79% for type I diabetics and 42.36% for type II diabetics compared to non-diabetic HF patients (reference group). The 95% CrI for the acceleration factors of HF patients with both types of diabetes does not include one, implying that diabetes has a significant effect on survival time. Regarding anemia, while controlling for other factors, the estimated acceleration factor for HF patients with anemia is 0.8445 [95% CrI: 0.729, 0.974], indicating that the expected survival time decreases by 15.55% compared to HF patients without anemia. The 95% CrI for the acceleration factor of HF patients with anemia does not include one, implying that anemia significantly affects survival time.

Finally, observing the stages of HF patients while keeping the effect of other factors constant, the estimated acceleration factors for stages II, III, and IV of HF patients are 0.6332 [95% CrI: 0.449, 0.888], 0.7033 [95% CrI: 0.509, 0.968], and 0.6402 [95% CrI: 0.467, 0.875], respectively. This result indicates that the expected survival time decreases by 36.68%, 29.67%, and 35.98% for stages II, III, and IV, respectively, compared to stage I (reference group). The 95% CrI for the acceleration factors of HF patients in stages II, III, and IV does not include one, implying that these stages significantly impact the survival time of HF patients.

From Table 3, the 95% CrI of the acceleration factors for covariates that do not include one are significant at the 5% significance level. The 95% CrI [1.0268, 95%: 1.00598, 1.625] of the frailties in the HF dataset suggests that this variability is significant. The frailty variance ($\eta$) = 0.0264 indicates heterogeneity between residence regions. Additionally, a frailty variance ($\eta$) of 0.0264 less than 1 suggests that HF patients living in rural areas have a higher risk of death (shorter survival time) compared to HF patients living in urban areas. Kendall’s tau ($\tau$ = 0.013) indicates a moderate amount of dependence in the clusters or unobserved heterogeneity between groups (residence region). Therefore, the frailty term had a significant contribution to the model. The Kullback-Leibler divergence (KLD) values for all significant parameters in the Bayesian log-normal shared gamma frailty model were 0, indicating that the posterior distribution was well approximated by a normal distribution. The most efficient algorithm used was a simplified Laplace approximation, which improved efficiency and resulted in faster computation speed.

3.1.4. Bayesian Model Diagnostics

By observing the Bayesian Cox-Snell residual plots in Figure 2, we can see that the Bayesian log-normal shared gamma frailty model best fit HF data-set from the given alternative models, since the plot of Cox-Snell residuals against cumulative hazard function of residuals was approximately a straight line with slope one and Bayesian cox-snell residual plot for Bayesian log-normal shared gamma frailty model was nearest to the line through the origin. In addition, the plot also indicated that the Bayesian log-normal shared gamma frailty model describes the HF data-set well.

Figure 3 includes a 95% credibility interval by observing the posterior density for the parameter’s normal distribution in the HF dataset. In addition, the histograms of the cross-validated Probability Integral Transform (PIT) values show that the posterior predictive p-values are somewhat closer to a uniform distribution, with some outliers observed in the HF data set. The Conditional predictive ordinate (CPO) values are significantly smaller (order of magnitude smaller) than the others, so the observed values would be considered surprising concerning the Bayesian log-normal shared gamma frailty model because the sum of the observations associated with failure flags is equal to zero in the HF dataset. The KLD is a diagnostic that measures the accuracy of the INLA approximation since Table 3 shows that the KLD values for all significant parameters in the Bayesian log-normal shared gamma frailty model were 0.

3.2. Discussion

This study employed the Bayesian AFT shared gamma frailty models on the HF dataset obtained from Jimma University Medical Center. Heart failure is a growing global issue, with its prevalence increasing among the adult population in developing countries with aging populations [5,6]. The study considered 497 HF patients to attain the proposed objectives. The overall survival time observed from HF patient’s follow-up ranged from 6 months to 42 months. The result shows that the estimated median survival time of HF patients who received treatment survived 34 months or above. Moreover, among those HF patients, about 301 (60.56%) were censored and the remaining 196 (39.44%) died. This finding was supported by other studies conducted by [20,23,71].

In the study, among HF patients, 21.93% have ischemic heart disease, 22.54% have rheumatic valvular heart disease, 23.54% have cardiomyopathy, 24.14% have hypertensive heart disease, and 7.85% have other diseases. This finding was consistent with a study conducted by [8,9,16,20,23]. Most HF patients, about 72.03%, live in rural areas, while the remainder reside in urban areas; however, survival times appear to be shorter for those living in rural areas. This finding was consistent with a study conducted by [20,23].

Furthermore, our study results indicate heterogeneity among clusters or groups (residences) at the time of death of heart failure (HF) patients, suggesting residence as a significant clustering factor. In line with our current work, several authors have explored Bayesian frailty models [60,63,64], including parametric models with baseline distributions such as exponential, Weibull, log-normal, and log-logistic, incorporating gamma frailty. Bayesian approaches have been applied to AFT shared gamma frailty models, and efficiency comparisons of different AFT models have utilized DIC and WAIC [68,69]. The Bayesian log-normal AFT shared gamma frailty model emerged as the most effective for describing the HF dataset among the alternatives. Similar findings were reported in previous studies [23]. However, our study, using the INLA method, reveals that age, chronic kidney disease, history of HF, diabetes mellitus, etiology of HF, hypertension, anemia, smoking, and HF stage all significantly impact the survival time of HF patients.

The findings of this study reveal that age significantly impacts the survival time of heart failure (HF) patients, with older individuals (65 years and older) experiencing shorter survival times. This trend is consistent with results from several other studies [19,20,22,23,72,73]. Additionally, HF patients without hypertension tend to have longer survival times compared to those with hypertension, underscoring the significant effect of hypertension on survival. These results are corroborated by studies conducted by [17,19,20,23,73].

On the other hand, the survival time of heart failure (HF) patients who smoke is shorter compared to nonsmokers, consistent with findings from studies by [17,20,23]. Both types of diabetes mellitus significantly impact HF patients’ survival, with those having either type of diabetes experiencing shorter expected survival times compared to HF patients without diabetes. These findings align with studies conducted by [17,19,20,23]. Additionally, chronic kidney disease significantly affects survival time, with HF patients without chronic kidney disease having longer survival times compared to those with the condition. This result is corroborated by studies from [19,20,23,73].

According to findings from [19,20,22,23], anemia significantly impacts the survival time of heart failure (HF) patients, with those having anemia experiencing shorter expected survival times compared to non-anemic HF patients. This study aligns with those results. Additionally, studies [20,23,74] support the finding that HF patients with conditions such as rheumatic valvular heart disease, hypertensive heart disease, cardiomyopathy, and other heart diseases have significantly reduced survival times. The stages of HF also markedly influence survival, with increased stage correlating with decreased survival time. This conclusion is supported by studies [19,20,23], which demonstrate that advancing HF stages significantly impact survival time.

To assess model adequacy, we plotted cumulative hazard plots for the Bayesian Cox-Snell residuals of the Exponential, Weibull, Log-normal, and Log-logistic models, as shown in Figure 2. The Bayesian log-normal shared gamma frailty model approached the line more closely, indicating it was the most suitable for the HF dataset, consistent with previous research [23]. We also used conditional predictive ordinates and probability integral transforms to evaluate the model. Before graphical assessment, we confirmed that no numerical issues arose in the computation of conditional predictive ordinates, as indicated by zero failures. Histograms and scatter plots of the probability integral transform showed that the predictive residuals were reasonably uniformly distributed with some outliers, supporting a good match between the predictive distribution and actual data. These results are in line with studies by [23,66,75], highlighting the robustness of our model.

Diagnostic plots for the Bayesian log-normal shared gamma frailty model, including a 95% credibility interval, show that the posterior density for parameters is normally distributed. The Kullback-Leibler divergence (KLD) values for all significant parameters were zero, demonstrating that the model offers both speed and high accuracy. These findings align with results from [23,66,75]. Despite limitations, such as reliance on potentially incomplete secondary data and missing prognostic factors, this study is pioneering in Ethiopia for its rigorous analysis and critical insights into heart failure patient survival, highlighting the robustness and advantages of the proposed method.

4. Conclusions

The increasing application of survival methods to estimate median survival time and identify associated factors is crucial. However, traditional methods often struggle with overfitting or underfitting, resulting in large standard errors in regression coefficients. They also frequently fail to account for within-cluster dependence and handle heterogeneity effectively. To address these shortcomings, we analyzed the survival times of heart failure patients treated at Jimma University Medical Center. This study employed various advanced survival methods to estimate survival time and identify risk factors for heart failure patients in this region. To achieve the proposed objectives, we employed Bayesian AFT and Bayesian AFT shared frailty models, comparing their performance using the heart failure patient dataset based on DIC and WAIC criteria. Notably, the Bayesian AFT models with shared gamma frailty significantly outperformed the Bayesian AFT models without frailty across all AFT model classes considered in this research.

Among the Bayesian AFT shared frailty models, the log-normal model emerged as the best performer. This finding highlights the importance of incorporating frailty to enhance model accuracy and parameter estimation. Subsequently, we then conducted an extensive data analysis using the model that demonstrated the best performance. This is due to the fact that log-normal Bayesian AFT shared gamma frailty model provides flexible survival data modeling by incorporating both individual-level and group-level random effects, effectively handling unobserved heterogeneity and delivering robust survival estimates. This model consistently outperforms alternatives such as exponential, Weibull and log-logistic AFT models. In this study, 50% of heart failure patients who received treatment survived for 34 months or more. Significant factors affecting survival included age, chronic kidney disease, history of heart failure, diabetes mellitus, etiology, hypertension, anemia, smoking, and heart failure stages. Additionally, survival times varied between rural and urban areas, highlighting the model’s ability to capture important heterogeneity.

More specifically, the study findings suggest that age groups (49 to 65 years and $\ge 65$ years), etiology of heart failure (rheumatic valvular heart disease, cardiomyopathy, hypertensive heart disease, and other conditions), hypertension, anemia presence, history of HF (medical OPD), chronic kidney disease presence, smoking, diabetes mellitus (type I and type II), and heart failure stages (II, III, and IV) reduced the survival time of heart failure patients. Furthermore, residence differences between rural and urban areas significantly influenced HF patient survival times, indicating heterogeneity and highlighting the need for frailty models.

Our results underscore the superior performance of the Bayesian log-normal shared gamma frailty model compared to current state-of-the-art methods.