Abstract
Survival analysis is a collection of statistical techniques which examine the time it takes for an event to occur, and it is one of the most important fields in biomedical sciences and other variety of scientific disciplines. Furthermore, the computational rapid advancements in recent decades have advocated the application of Bayesian techniques in this field, giving a powerful and flexible alternative to the classical inference. The aim of this study is to consider the Bayesian inference for the generalized log-logistic proportional hazard model with applications to right-censored healthcare data sets. We assume an independent gamma prior for the baseline hazard parameters and a normal prior is placed on the regression coefficients. We then obtain the exact form of the joint posterior distribution of the regression coefficients and distributional parameters. The Bayesian estimates of the parameters of the proposed model are obtained using the Markov chain Monte Carlo (McMC) simulation technique. All computations are performed in Bayesian analysis using Gibbs sampling (BUGS) syntax that can be run with Just Another Gibbs Sampling (JAGS) from the R software. A detailed simulation study was used to assess the performance of the proposed parametric proportional hazard model. Two real-survival data problems in the healthcare are analyzed for illustration of the proposed model and for model comparison. Furthermore, the convergence diagnostic tests are presented and analyzed. Finally, our research found that the proposed parametric proportional hazard model performs well and could be beneficial in analyzing various types of survival data.
1. Introduction
The healthcare domain has evolved significantly in recent years as a result of computational developments. The use of Bayesian statistics in healthcare has encouraged the application of computational developments, providing a powerful and versatile alternative to traditional methodologies used in healthcare [1]. The progress of Bayesian approaches in healthcare aims to make an individual’s life more affordable and comfortable, similar to how smartphones have made life easier [2]. Despite the fact that the idea of applying computational Bayesian statistics to survival analysis dates back to the 19th century, McMC techniques are now garnering more attention in the literature because of abundant and cheap computation [3]. The application of deep learning to the context of parametric survival models was discussed by [4]. Through an efficient training process, the quality of the developed system improves. Data portioning is done three times to confirm the trained algorithms (training-testing-validation) [5]. The main goal of this article is to present the Bayesian parametric proportional hazard model using BUGs syntax.
The statistical analysis of survival data is an essential topic in many fields, including medicine, biology, environmental science, healthcare, economics, engineering, social science, and epidemiology, among others. Probability distributions serve as the foundation for survival models. The family of distributions can be parametric, semiparametric, or nonparametric. The parametric survival models lead to more efficient and smaller standard errors of the estimates than semiparametric and nonparametric models [6] if the distributional assumption is correct, to be more specific.
In analyzing survival data, parametric survival models are crucial. The benefits of using parametric survival models include the following: (1) handling all types of censored data (left, right, interval, double, and middle); (2) application of survival analysis in a healthcare care problem and (3) producing better estimation when you have a theoretical expectation of the baseline hazard; also, (4) they can apply random effects—frailty models—and can also be used to estimate expected lives, not only hazard ratios like the accelerated failure time models [7].
The proportional hazards (PH) model, in which covariates affect the hazard rate function, and the accelerated failure time (AFT) model, in which covariates affect both the hazard rate and time scale, are the two most common methods for developing parametric regression models for survival data [8]. However, other class of models have also been proposed such as the accelerated hazard (AH) model [9] and the proportional odds (PO) model [10].
One of the first steps in using a parametric approach to model survival data is to choose a suitable baseline distribution that can capture significant features of the observations of interest. Certain probability distributions are widely used in the modelling of survival data. Only a few are closed under the proportional hazard model, and none are flexible enough to describe a wide range of survival data [11]. Most of the distributions closed under the PH assumption fails to model a nonmonotone (i.e., bathtub and unimodal) survival data sets.
The log-logistic (LL) distribution has a wide range of applications in survival data analysis and can accommodate unimodal survival data sets. The distribution is closed under both proportionality odds (PO) and multiplication of failure time (AFT) frameworks [7]. It is not a PH model, but an AFT model. However, when the log-logistic distribution is generalized, it has the appealing feature of being a member of all classes of parametric hazard-based regression models of the survival analysis because its failure rate function is quite versatile and its cumulative hazard function (chf) has a tractable form.
Extensive efforts have been made over the last decades to extend classical distributions to use as a baseline distribution for parametric hazard-based regression models. Many modifications to the LL distribution have been introduced to make it more adaptable to a wide range of hazard shapes [12–16]. The generalized log-logistic distribution (GLL) is one such model, which modifies the log-logistic distribution by inducing an additional shape parameter [17]. The model is tractable and closed under the PH assumption and can account for both nonmonotone and monotone hazard rates [11]. On the other hand, recent computational advances have advocated for the use of Bayesian techniques in the field of survival and reliability analysis.
The motivating ideas behind our work on Bayesian parametric proportional hazard (PH) model with GLL baseline hazard are as follows: (i) despite the fact that there are some classical distributions closed under the PH framework, none of which is flexible enough to incorporate both monotone and nonmonotone hazard rate; (ii) Bayesian inference does not rely on asymptotic approximation for statistical inference; (iii) the availability of software makes Bayesian implementation for hazard-based complicated models relatively more straightforward and simple than classical inference [18]; (iv) parametric PH model may lead to more precise estimates than the semiparametric PH model; and, last but not the least, (v) the use of generalized distributions that can capture both monotone and nonmonotone hazard rate functions is what makes our work unique and more appealing to biostatisticians, epidemiologists, healthcare workers, and other applied researchers in multiple disciplines.
To the best of author’s knowledge, no Bayesian inferences study has been conducted on the PH model with generalized log-logistic baseline hazard. As a result, in this paper, we consider the Bayesian inference for the generalized log-logistic proportional hazard model, beginning with the PH model formulation and assumptions, revising the generalized log-logistic distribution, and verifying that the GLL distribution is closed under the PH framework. In addition, we discuss the inferential procedures and how to obtain the classical and Bayesian estimators for the model’s parameters. We also compare the proposed model to other existing distributions closed under the PH framework, and one interesting feature of this model is that it can incorporate different hazard rate shapes. Hence, the formulation of the parametric PH model and its lifetime function, the inferential procedures using both classical and Bayesian approaches, and the development of the computational algorithms to fit the proposed PH model and its competing models using RJAGS in R software are the novelty of this study.
The article is structured as follows: the PH model formulation, assumptions, and its probabilistic functions are discussed in Section 2. Section 3 revises the most common probability distributions closed under the PH model. Section 4 presents the proposed baseline hazard function which is a generalized log-logistic (GLL) distribution. The GLL distribution under the PH model is presented in Section 5. Section 6 discusses the inferential procedures of the proposed model. In Section 7, we present an McMC simulation study to assess the performance of the proposed model. Section 8 presents the application of the proposed model to two right-censored cancer data sets with monotone and nonmontone hazard rates. In addition, the convergence diagnostics of the McMC techniques were discussed. The Bayesian model selection criterion is presented in Section 9. Finally, in the final portion, the article’s concluding notes are offered, and some future works are mentioned.
2. PH Model Formulation and Assumptions
In many real-life applications, survival times are affected by explanatory variables. The explanatory variable vector is related to response variable through a regression model. An important aspect of survival modelling is the inclusion of explanatory variables. The hazard-based regression models can be formulated in a number of ways. One of the most frequently used method is the proportional hazard (PH) model formulation.
PH models play an essential role in analyzing time-to-event data and are broadly used in survival and reliability analysis as well as in joint modelling of survival and longitudinal data [7]. It is the most popular parametric model in medical studies and clinical trials because of the existence of a semiparametric PH hazard model which is robust against the distributional assumption of the survival time. The parametric PH model is given with the similar form to the Cox PH model. It is the parametric form of the Cox PH models [6].
2.1. PH Model Formulation
The parametric proportional hazard (PH) models are formulated using a defined baseline hazard and a link function for the covariates which is defined as follows:
The most commonly found option for the link function is the exponential (or log-linear) function. In this work, we define the PH model with the assumption that
2.2. PH Assumptions
The PH model assumption is that the effect of covariates is to increase or decrease the hazard rate function by a proportionate amount which does not depend on t. The assumption of the PH model can be defined as follows:where is called the baseline hazard.
Simplifying, we get,
The main difference between the Cox PH model and the parametric PH model is that the baseline hazard function is assumed to follow a specific distribution when it is fitted to the data. Using equation (2), we can see that the hazard ratio (HR) comparing any two specifications of the covariates, for example, is
The above equation shows us that the baseline hazards cancel each other from this ratio, so the hazard rate for one individual is proportional to the hazard rate for any other individual. On the other hand, the proportionality constant is independent of time which makes the main assumption of this model [6]. As a result, the model is known as the proportional hazard (PH) model in the literature.
Unlike most parametric regression models including accelerated failure time (AFT) models, PH models does not include an intercept [19]. More properly, the vector in the PH model is not assumed to have . An intercept would get confounded with the baseline hazard function
2.3. Lifetime Functions Describing the PH Model
The five frequent representatives of a lifetime distribution function that are used to characterize the PH model are addressed in this section.
2.3.1. Hazard Rate Function of the PH Model
In the PH analysis, one of the most important lifetime functions is the concept of the hazard rate function (hrf). The hazard rate function , abbreviated by hrf, also called the instantaneous failure rate or as force of mortality of a PH model is of the form:
2.3.2. Cumulative Hazard Function of the PH Model
The hazard or survival functions, rather than the cumulative distribution or probability density function, are typically used in the PH analysis of survival data. The hazard rate function is used to interpret the most common survival regression models; however, the cumulative hazard function (chf), also known as the integrated hazard rate function, can be easily written down. Hence, the chf of a PH model takes the following form:
2.3.3. Survival Function of the PH Model
The survivor function (sf) for a PH model can be derived using the following relationship between survival function and the hazard rate function. Hazard function is given by
Cumulative hazard function:
Using the above expressions, we can easily find
2.3.4. Cumulative Distribution Function of the PH Model
The cdf of the PH model, also known as the lifetime distribution function, is given by
2.3.5. Probability Density Function of the PH Model
The pdf or the failure density function of the PH model is defined as
The five representatives used here were chosen for their special meaning for lifetime data, their intuitive appeal, their utility in survival data analysis, and, last but not the least, their popularity in probability theory and statistics.
The PH model can be formulated without assuming a probability distribution for survival times, and this leads to the well-known Cox PH model [20]. On the other hand, assuming a probability distribution for survival times leads to the fully parametric PH model. The most common parametric survival models used are as follows: exponential, Weibull, Gompertz, log-logistic, log-normal, gamma, and the generalized gamma distributions. Only the exponential, Weibull, and Gompertz distributions are used for the PH model. The log-logistic and the log-normal distributions are not closed under the PH framework. Weibull distribution is the only one that is closed under both parametric AFT and PH models.
3. Distributions Closed under PH Framework
In this section, we present most common parametric distributions that are closed under the PH framework and are used to analyze survival data. These distributions have been studied and used in various contexts in the literature.
3.1. Exponential PH Model
Exponential distribution is a continuous probability distribution with only one unknown parameter It is the simplest distribution for lifetime distribution models. The distribution is not flexible enough to describe commonly encountered hazard rate shapes for survival data. The pdf, cdf, sf, hrf, and chf of the exponential random variable are, respectively, as follows.
Let ,where is the scale parameter . A short value of shows low risk and long survival, where a large value shows high risk and short survival. For the PH model, the exponential baseline hazard is
So, according to the formulation of the PH framework, the hazard rate for an individual with covariate vector and link function is
Applying the log-linear function we can simplify into
In this equation, the hrf has the exponential distribution with scale parameter which indicates that the PH assumption is satisfied with the exponential distribution. It is worth mentioning that the exponential distribution is often found to be inadequate to describe survival data. This makes the applicability of this distribution fairly limited.
The other lifetime distributions of the exponential PH model are as follows.
The survival function of the exponential PH model is
The pdf of the exponential PH model is
The cdf of the exponential PH model is
The chf of the exponential PH model is
3.2. Gompertz PH Model
Gompertz distribution is named after Benjamin Gompertz, a British mathematician and actuary, who developed it in 1825. It is a continuous probability distribution used for modelling adult life spans and other application under different disciplines such as actuarial science, demography, survival, and reliability analysis. This distribution is flexible and can be skewed both in right and in left. The pdf, cdf, sf, hrf, and chf of the exponential random variable are, respectively, as follows.
Let ,where is the rate parameter, . When the survival time then has an exponential distribution; therefore, Gompertz distribution is a generalization of exponential distribution. For the PH model, the Gompertz baseline hazard rate function is given by
So, according to the formulation of the PH framework, the hazard rate for an individual with covariate vector and link function is
Applying the log-linear function we can simplify into
In the above equation, it is straightforward that the PH property is satisfied. However, the Gompertz PH model is rarely used in the real-life applications.
The other lifetime distributions of the Gompertz PH model are as follows: the survival function of the Gompertz PH model is
The pdf of the Gompertz PH model is
The cdf of the Gompertz PH model is
The chf of the Gompertz PH model is
3.3. Weibull PH Model
Weibull distribution is a generalization of the exponential distribution. It is a versatile distribution that can take on the characteristics of other types of continuous distributions. It has an additional parameter compared to the exponential. The additional parameter describes the shape of the hazard functions, based on the value of the shape parameter [21]. The pdf, cdf, sf, hrf, and chf of the Weibull random variable are, respectively, as follows.
Let ,where is the shape parameter and is the rate parameter. The hazard rate function increases when decreases for and constant for When , the Weibull distribution reduces to exponential. It is worth mentioning that the Weibull distribution does not accommodate nonmonotone (i.e., unimodal or bathtub) hazard rates.
For the PH model the Weibull baseline hazard is
So, according to the formulation of the PH framework, the hazard rate for an individual with covariate vector and link function is
Applying the log-linear function we can simplify into
In this equation, the model has the Weibull distribution with rate parameter and shape parameter which indicates that the PH assumption is satisfied with the Weibull distribution with constant
The other lifetime distributions of the PH Weibull model are as follows: the survival function of the Weibull PH model is
The pdf of the Weibull PH model is
The cdf of the Weibull PH model is
The chf of the Weibull PH model is
4. Parametric Baseline Hazard
The parametric baseline hazard function is essential because it determines which hazard shapes can be captured by the proportional hazard (PH) model. Most classical distributions that are closed under the PH framework, such as the exponential, Weibull, and Gompertz distributions, are incapable of accommodating unimodal hazard shapes. As a result, it is worth looking into some modifications to the classical distributions that can account for both monotone and nonmonotone hazard rates.
In this paper, we consider the Bayesian inference for the parametric PH models with generalized log-logistic (GLL) baseline. The GLL is a flexible survival distribution proposed by [11]. This distribution has a characteristic similar to those of the log-logistic distribution. Also, the advantage of the GLL distribution is that it approaches to Weibull in the limit. These properties allowed the GLL to handle both monotone and nonmonotone hazard functions, and also it makes to be a baseline distribution that is closed under both AFT and PH model [22] like the Weibull distribution. The distribution is adaptable, and the two shape parameters enable a wide range of hazard shapes. It also includes a variety of important distributions such as the exponential, Weibull, Burr XII, and log-logistic distributions. In addition, when compared to competitors, it is relatively tractable. We refer to, for more information on the distribution and its properties, [17].
For a positive-valued random variable T, the hrf of the GLL distribution with three unknown parameters is given by
The chf of the GLL distribution is given by
The distribution function of the GLL model is of the form:
The survival function (sf) of the GLL model is given bywhere are parameters and .
The quantile function of the GLL model is given by
The reverse cumulative hazard rate function is expressed as follows:
Figure 1 illustrates shapes that the failure rate function can accept such as constant, increasing, decreasing, V-shape, and unimodal among others.
5. The Proposed PH Model
For the PH model, the generalized log-logistic baseline hazard is
So, according to (2), the hazard rate for an individual with covariate vector and link function is
Applying the log-linear function we can simplify into
In this equation, the hrf can be recognized as a generalized log-logistic distribution as well, but contrary to (36), the rate parameter is and shape parameters are which indicates that the PH assumption is satisfied with the GLL distribution and the proposed model is closed under the PH framework.
The other lifetime distribution functions for the GLL PH model are as follows: the survivor function of the GLL PH model is
The pdf of the GLL PH model is
The cdf of the GLL PH model is
The chf of the GLL PH model is
6. Model Inference
We discuss the classical approach (using maximum likelihood (MLE)) and Bayesian approach (assuming noninformative priors) estimation techniques for the proposed parametric PH model parameters in this section.
6.1. MLE for Right Censored Survival Data
We examine the challenge of estimating the proposed model’s distributional parameters and regression coefficients for right-censored survival data in this section. Because of its appealing qualities, such as consistency, asymptotic efficiency, asymptotic unbiasedness, and asymptotic normality, MLE is one of the most common strategies for estimating the parameters of hazard-based regression models. Let there be individuals with lifetimes represented by Assuming that the data are subject to right censoring, we observe ), where corresponds to a potential censoring time for individual Allow ) that equals 1 if and 0 otherwise.
Suppose that a right-censored random sample with data , is available, where is the censoring time or a survival time according to whether respectively, and is an column vector of external covariates for the ith individual, , and is the vector of regression coefficients. When the parametric PH model is considered, the censored likelihood function can be expressed as
An iterative optimization procedure (e.g., Newton–Raphson algorithm) can be used to obtain the maximum likelihood estimation . Hypothesis testing and interval estimations of model parameters are possible due to the MLEs’ approaching normality [7]. The natural logarithm of the likelihood function, so-called log-likelihood function can be written as follows:where is a vector of the regression coefficients and is the vector of the baseline distributional parameters.
In our case, if we assume that and . Use (36) for and note that is the baseline cumulative hazard rate function as given by (37). The full log-likelihood function of the GLL PH model can be expressed as follows:
To obtain the MLE’s of and , we can maximize (51) directly with respect to or we can solve the nonlinear equations below or the 1st derivative of the log-likelihood function. The 1st derivatives of the log-likelihood function are
To maximize log-likelihood functions, many software packages are available including proven optimization algorithms.
6.2. Bayesian Inference
In this section, Bayesian inference was used to estimate distributional parameters and regression coefficients using objective (or noninformative) priors to obtain proper posterior distributions.
6.2.1. Priors for the Model Parameters
The specification of a prior distribution is a crucial aspect of any Bayesian inference. In parametric survival regression models, this is especially true. As a result, the prior scenario is built in this study using a noninformative independent prior for the parameters. The marginal prior distribution for every regression coefficient is prompted as a normal distribution centred at zero and with a small precision, ; on the other hand, a gamma distribution, is chosen as the marginal prior distribution for the parameters of the GLL PH model due to the versatility of gamma distribution that include the noninformative priors (uniform) on the shape parameters. Many research publications in the literature, such as Danish and Aslam [23, 24], considered the assumption of the gamma priors for the baseline hazard parameters of PH models. Alvares et al. [1] took the assumption of independent gamma priors for the baseline hazard parameters of eight different parametric survival models. Muse et al. [22] used the assumption of independent gamma priors for the baseline hazard parameters of the of the generalized log-logistic AFT model, and other researchers take these priors into account.
For the baseline parameters of the GLL-PH model, we assume independent gamma priors.
Prior to that, we had the regression coefficients (assuming a normal distribution).
The density function of the combined prior distribution of all unknown parameters and the regression coefficients are given as
6.2.2. The Likelihood Function
Unfortunately, the likelihood function of this generalized model is not implemented in BUGS and JAGS syntax. To generate the likelihood function, we use the “zero’s trick” method that become popular in survival analysis and relies on Poisson modelling of expanded or reconstructed data [1]. The zero's trick approach works on the assumption that perhaps the contribution of a Poisson observable of zero is ; if we set with observable data as a vector of , we receive the right contributions of the proposed model [18].
6.2.3. The Posterior Distribution
The joint posterior density function is equal to the multiplication of the prior distribution and the likelihood function the joint posterior density function of the parameters of GLL PH model given the data can be expressed using Bayes’ theorem aswhere the first four terms on the equation represent the prior specification for the unknown parameters and are assumed to be independent and is the likelihood function expressed as follows:
The marginal distributions of the model parameters and the normalising joint posterior density function are difficult to calculate analytically, requiring high-dimensional integration and no close form inferences. To obtain estimates, we use McMC simulation methods, which involve sampling from the posterior distribution through using the Metropolis–Hastings Algorithm.
7. Simulation Study
In this section, we undertake an extensive simulation investigation to demonstrate the proposed parametric proportional hazard model’s good Bayesian features. The parameter values are chosen to construct situations that mimic cancer population studies using a cancer that is severe (with a lower five-year survival rate), such as lung cancer [9, 25]. We demonstrate parameter estimation, the effect of censoring proportions, and sample sizes on inference in more detail.
7.1. Generating Survival Data from the PH Model
To simulate survival data for the GLL PH model, we use the inversion technique [40, 41] to generate survival data. This strategy is based on the link between a survival random variable’s cumulative hazard rate function and a standard uniform random variable. When the cumulative hazard rate function has a closed form expression, it may be immediately applied, inverted, and readily implemented with R [26]. The censoring rates were estimated using administrative censorship at (1) Tc = 5 years, which resulted in around 20% censoring in all sets, and (2) Tc = 3 years, which resulted in about 30% censoring in all sets.
For the purposes of this simulation, we assume that survival times are distributed using the generalized log-logistic distribution (, , ). Using the reverse chf given in equation (41), lifetime data can be simulated as follows:where , , > 0.
7.2. Simulation Design
The simulation analysis was carried out by conducting a series of simulations with different sample sizes (n = 100, and 300) sets and censoring proportions (Tc = 20 and 30 percentages), all based on the PH model in equation (1). The GLL PH model’s true parameter vector is set as follows: (1) set I: distributional parameter values (, , and ) and covariates (2) set II: distributional parameter values , , and and covariates .
The values of the covariates were simulated as follows: (1) combination of uniform distributions with 0.25 probability on (30, 65), 0.35 probability on (65, 75), and 0.40 probability on (75, 85) years old was used to simulate the continuous covariate “age,” and (2) the binary covariates “treatment” and “gender” were both simulated using a 0.5 binomial distribution. We recommend that the reader can refer [9] for further details.
7.3. Posterior Analysis of the Simulated Data
We fitted the proposed PH model with GLL baseline hazard to assess its Bayesian properties in the simulation sets. With all censoring rates and different sample sizes, each simulation set was used to estimate the proposed PH model. JAGS software [27] was used to approximate posterior distributions using three parallel chains with 40,000 iterations each plus another 3,000 for the burn-in period. To minimize autocorrelation in the sequences, the chains were thinned further by storing every 10th draw.
7.4. Measures of Performance
The actual mean, standard deviation (SD), Naive standard error, bias, percentage of bias, coverage probability (CP), potential scale reduction factor and the effective number of separate simulation draws were used to test the posterior distribution stability for the suggested PH model.
7.4.1. Evaluating the Performance of the Estimators
We calculate the bias of the estimators using:
An underestimation is indicated by a negative bias, whereas an overestimation is shown by a positive bias.
7.4.2. Accuracy of the Estimators
The mean square error (MSE) is a good indicator of overall accuracy and is calculated as follows:
This metric determines how accurate the estimates are as follows. The lower the MSE, the more accurate the estimations of impacts.
The Naive standard error, which is calculated by dividing the posterior standard deviation by the square root of the sample size, is another accuracy metric. As a result, the smaller the standard error, the larger the sample size. The Naïve SE incorporates simulation error rather than posterior uncertainty.
7.4.3. Coverage
The 95 percent coverage probability (CP) is the percentage of N simulated data sets in which the true estimates were included in the 95 percent confidence interval. The more precise the estimations are, the closer the outcome is to a 95 percent coverage probability. The following is how CP is expressed:
7.4.4. Convergence Diagnostics
Quantitatively, Gelman et al. [28] recommended that the acceptable limit of multivariate potential scale reduction factor (MPSRF) and potential scale reduction factor (PSRF) be near 1 , and the effective number of sample size simulation draws be greater than or equal to 100 for checking the convergence of McMC simulations. It is clear from the summary characteristics (Tables 1–4) that the PSRF is less than 1.1, that number of sample size simulation draws is larger than 100, and that Naive SE is smaller than the standard deviations (SD) for all of the distributional parameters and regression coefficients, as expected, indicating that the McMC algorithm has converged to the posterior distribution. Trace plots, autocorrelation plots, and Gelman plot diagnostics are the most common ways to judge the convergence of a McMC simulation graphically [28]. The McMC simulation has been achieved as evidenced by the trace plot, density plot, autocorrelation plot, and Gelman diagnostic plots for each distributional parameter and regression coefficients. That is, the McMC simulation for the GLL PH model explores the target posterior distribution appropriately.
7.5. Simulation Results
Tables 1–4 shows the simulation results for the posterior mean, bias, Naive standard error (SE), mean square error (MSE), coverage probability (CP), Gelman–Rubin diagnostic (, and the number of sample size simulation draws (no. of Eff) of the proposed PH model, and Figures 2–5 shows the visual summary for the convergence diagnostics.
Based on these findings, we may deduce that, as the sample size grows, the biases and MSE of the estimators decrease; also, the censoring proportion impacts the bias and MSE of the estimators, with larger censoring rates increasing the bias and MSE. The Gelman–Rubin diagnostic, on the other hand, as well as the number of efficiency sample size draws show that convergence has been attained. However, the estimators’ coverage probability was close to 95%.
8. Practical Illustrations
In this section, two real-life survival data sets dealing with right-censored cancer data sets were considered to demonstrate the flexibility and applicability of the proposed GLL PH in modelling different survival data sets with different hazard rate shapes.
8.1. Data Set I: Lung Cancer Survival Data
8.1.1. Data Description
In this section, we reanalyse the data set reported in [29] which is available in the R package survival. The Veterans Administration Lung Cancer Study Group followed up on n = 137 patients in this data set. For this clinical investigation, the censorship rate is around 6.5 percent (9 observations out of 137 were censored). The response and exploratory factors in this clinical trial are the time until death (in days), the number of months from diagnosis to study enrolment (diagt), age (in years), a history of previous lung cancer therapy (prior), and the trt = (treatment = conventional chemotherapy).
8.1.2. Hazard Rate Shape
The hazard rate function appears to be unimodal or decreasing in Figure 6 based on the TTT plot (careful inspection reveals a slight indication of unimodality). The data could be evaluated with a model like the log-logistic distribution, which can accommodate decreasing or unimodal hazard rate forms. However, because the classical LL distribution is not closed under the PH framework, we employ the GLL distribution, which is closed and can encompass various hazard rate shapes. The box plot, histogram, and TTT plots are shown in Figure 6.
8.1.3. Proportionality Assumption
There are two widely used methods for assessing the PH assumption: (1) graphical diagnostics based on (a) time-dependent variables [7] and (b) standardized Schoenfeld residuals [30] and (2) statistical tests. The standardized Schoenfeld residuals are used in this section to evaluate the PH assumption of the Cox model for each covariate included in the model. Based on Figure 7 and the significance threshold of 5%, there is no evidence to reject the proportional hazards assumption. As a result, we anticipate that the GLL PH model will provide a good fit when compared to the other existing parametric PH model employed in this study.
8.1.4. Posterior Analysis
In this paper, we assume the noninformative independent framework with a normal prior for (regression coefficients) and an independent gamma prior for the distributional parameters with hyperparameter values (.
(1) Numerical Summary. We looked at various quantities of interest and their numerical values using the McMC sample of posterior properties for the generalized log-logistic proportional hazard model using the lung cancer data in this section.
The posterior summaries for the generalized log-logistic PH model parameters using Veterans lung cancer data sets are illustrated in Table 5. The probability that the corresponding parameter is +ve is given in the last row of Table 5. A posterior probability of 0.5 indicates that a positive parameter value is as likely as a negative one. Once we’ve saved the posterior sample for each model parameter, we can compute the posterior probability, for example, for , using mean
(2) Visual Summary. We looked at density strip plots, trace plots, Gelman–Rubin diagnostic plots, Ergodic mean plots, and autocorrelation diagnostic plots in this section to get a visual description of the posterior properties. These plots and graphs provide a nearly comprehensive representation of the parameters’ posterior uncertainty for the application of the lung cancer data sets.
(3) Density Plots. Density can be compared to the fundamental shapes associated with typical analytic distributions, and density plots can reveal behaviour in the tails, skewness, existence of multimodal behaviour, and data outliers. Figure 8 illustrates the density plots for the GLL PH model parameters.
(4) Time-Series Plots. One of the most common diagnostics of an McMC simulation is a time series plot (or trace plot) [28]. Figure 9 shows that the McMC sampling process converges to the joint posterior distribution with no periodicity. As a result, we can say that the chains have converged.
(5) Brooks–Gelman–Rubin (BGR) Convergence Diagnostic. Gelman and Rubin [31] propose a convergence diagnostic technique to check the McMC algorithms simulation and is based on within chain variance and between chain variance. Gelman et al. [28] suggested that the limit of acceptance of potential scale reduction factor (PSRF) to be less than 1.1. Figure 10 shows us that both PSRF and MPSRF are less than 1.1.
(6) Running Mean Plots. The running mean (also referred to Ergodic mean) is a well-known convergence diagnostic for McMC algorithms. The Ergodic mean is defined as the mean of all simulated sample values of up to a specific iteration [32]. Ergodic mean is used to observe the convergence pattern of the McMC chains. Figure 11 shows us the Ergodic mean plots for the regression coefficients and the baseline hazard parameters. It is quite clear from the running mean time-series plots that the chains converge after N iterations to their mean values. However, these plots display only at the mean of the baseline hazard parameters and the regression coefficients and therefore are inadequate.
(7) Autocorrelation Plots. Although the autocorrelation plot is not strictly a convergence diagnostic tool, it does aid in indirectly assessing the convergence of the McMC simulation process [33]. Figure 12 shows the autocorrelation plots for all parameters and regression coefficients.
8.1.5. Convergence of McMC Algorithm for the Veterans Lung Cancer Data Set
Computational developments in the previous few decades have recently emerged as a very useful instrument for employing McMC approaches [34] and fitting Bayesian survival regression models in time-to-event analysis. The complicated posterior distribution is sampled using the McMC algorithm. As a result, when an algorithm converges to the target posterior distribution, the Markov chain is stationary, and adding more samples will not change the shape and position of the posterior distribution’s density in a meaningful way and hence will not change the estimations or other relevant outcomes.
(1) Common Statistical Tests for Convergence Diagnostics. The convergence of the McMC algorithm was checked quantitatively using conventional statistical tests for convergence diagnostics: (1) Brooks–Gelman–Rubin diagnostics [28]; (2) Raftery and Lewi diagnostics [35]; (3) Heidelberger and Welch’s diagnostic tests [36]; and (4) Geweke diagnostics [37]. For more information about these tests, we can refer to [34]. Table 6 indicates the Geweke, Raftery–Lewis, and Heidelberger–Welch diagnostics for the GLL PH model parameters.
(2) Graphical Techniques for Convergence Diagnostics. Convergence diagnostics of an McMC algorithm can be examined graphically, including: (1) time series plot; (2) autocorrelation plot; (3) running mean plot; and (4) Gelman–Rubin plots. See Figures 9–12.
8.2. Data Set II: Larynx Cancer Data Sets
8.2.1. Data Description
Lifetimes for 90 patients with larynx-cancer, according to the stage of cancer tumour (stages I–IV) are given in Table 7. The study time or time to death are recorded in months (where, shows us the censored time). Alvares et al. [1]; Wang et al. [8]; and Christensen et al. [19] discussed the data from different aspects under different hazard-based regression models, and the data were first reported by [38]. The survival times (in months) of patients is illustrated in Table 7.
The other covariates of the data are as follows: (1) age (in years) at diagnosis and (2) the year of diagnosis. One goal of this study was to see if the age, year of diagnosis, and stage of cancer were associated with the death of patients with laryngeal cancer.
8.2.2. Hazard Rate Shape
Based on the TTT plot, the hazard rate function is an increasing hazard in Figure 13. The data could be analyzed using a model such as the Weibull distribution, which can handle monotone hazard rate forms. We adopt the GLL distribution, which would be represented by the PH framework and can accommodate a variety of hazard rate shapes to see its applicability of the monotone (increasing) hazard rates. Figure 13 shows the box plot, histogram, and TTT plots.
8.2.3. Proportionality Assumption
We investigated if the proportional hazards model could be used with this data set. The underlying assumption of the Cox model for each explanatory variable utilized in the model is depicted in Figure 14. With a significance level of 5%, there is no evidence to reject the PH assumption. As a result, we anticipate that the parametric PH model will provide a strong fit.
8.2.4. Posterior Analysis
In this paper, we assume the noninformative independent framework with for (regression coefficients) and an independent gamma prior for the distributional parameters with hyperparameter values (
(1) Numerical Summary. We looked at various quantities of importance as well as their numerical values using the McMC sample of posterior properties for the generalized log-logistic proportional hazard model considering the larynx data in this section.
The posterior summaries for the GLL-PH model parameters using larynx cancer data are illustrated in Table 8. The probability that the corresponding parameter is +ve is given in the last row of Table 8.
(2) Visual Summary. We looked at density strip plots (Figure 15), trace plots (Figure 16), Ergodic mean plots (Figure 17), autocorrelation plots (Figure 18), and Gelman–Rubin diagnostic plots (Figure 19), in this section, to get a visual description of the posterior properties. These plots and graphs provide a nearly comprehensive representation of the parameters’ posterior uncertainty.
8.2.5. Convergence Diagnostic Tests for the Larynx Cancer Data Using GLL PH Model
(1) Statistical Tests. Table 9 indicates the Geweke, Raftery–Lewis, and Heidelberger–Welch diagnostics for the GLL PH model parameters.
(2) Graphical Techniques. Convergence diagnostics of an McMC algorithm for the larynx cancer data set are presented in Figures 16–19.
8.2.6. Hazard Ratio (HR)
One of the most intriguing aspects of PH models is that the regression coefficients can be interpreted using the hazard ratio, which is preferred by many clinicians.
A key feature for PH models is the hazard ratio (HR), also known as the relative risk, between two individuals with covariate vectors . The HR is defined aswhich does not depend on time t. the hazard function in the numerator is equal to this constant HR times the hazard in the denominator, i.e.,
Hence, the name “proportional hazards model” [19]. For example, the posterior distributions of the HR between two individuals of the same age and diagyr (year of diagnosis) but in different stages can be easily summarized.
Table 10 depicts the posterior characteristics of the hazard ratio between two men of the same age and diagnosis year (diagyr) but in different stages.
9. Bayesian Model Selection
In this study, we will use the deviance information criterion (DIC) to distinguish between the proposed models. DIC is a popular Bayesian model selection criterion. This criterion is available in most McMC packages [39]. The DIC is computed as follows:where denotes the deviance’s posterior mean and is a goodness of fit test for parametric survival models and calculates as the difference between , and it is denoted the effective number of proposed model parameters.
9.1. Data Set I
Table 11 displays some posterior characteristics for the three PH models (generalized log-logistic, Gompertz, and Weibull). Even though the estimates of the regression coefficient are significant compared, the flexibility provided by the GLL distribution’s additional shape parameter contributes to its ultimate superiority over the Gompertz and Weibull models and the DIC shows us its goodness-of-fit and versatility comparing to the competing parametric PH models.
9.2. Data Set II
Table 12 displays some posterior characteristics for the three PH models (generalized log-logistic, Gompertz, and Weibull). Even though the estimates of the regression coefficient are significant compared, the flexibility provided by the GLL distribution’s additional shape parameter contributes to its ultimate superiority over the Gompertz and Weibull models and the DIC demonstrates us its goodness-of-fit and versatility comparing to the competing parametric PH models.
10. Conclusion and Future Work
In this paper, we explored how to derive Bayesian estimates of the baseline hazard parameters and the regression coefficients of the parametric proportional hazard model with generalized log-logistic baseline hazard using right-censored survival data utilizing McMC approaches. The McMC techniques offer an alternative technique for estimating the parameters of the proposed model that is more flexible than frequentist techniques such as maximum likelihood estimation. Bayesian inference was performed with a variety of priors, and the convergence pattern was investigated using various diagnostic procedures.
To test the performance of the proposed parametric PH model, a comprehensive McMC simulation study was conducted. According to the simulation results, the PH model produces better results, with fewer absolute biases and MSEs for most regression coefficients and baseline distributional parameters. The behavior of the PH model in a generic PH regression situation comprising numerous covariates was also examined using synthetic right-censored data sets. Our findings indicate that the PH model performs well when handling with multiple factors. The paper’s final analysis focused on a real-world application involving two well-known right-censored survival data sets for lung cancer and laryngeal cancer patients. In conclusion, the findings of the proposed parametric PH model show that it performs better and is superior to the other competing PH model, as well as indicating significant distributional parameters and regression coefficients.
Furthermore, for both simulation and real-data analysis, the regression coefficients were assumed to have a normal prior, and the baseline distribution parameters were assumed to have an independent gamma prior to compute the quantities of importance derived from the proposed model’s posterior distribution. It has been attempted to create a visual summary and other essential graphs to aid in the interpretation of results and decision making. Finally, we hope that this paper will be an extension of the work of Khan and Khosa [11] and will encourage researchers who employ parametric hazard-based regression models to conduct their analyses using the Bayesian approach from the BUGs codes with the help of the R software’s RJAGS package.
In terms of future work, we intend to produce an R package to fit the most prevalent parametric hazard-based regression models, including the PH model. The method given in this study can also be applied to multiple event scenarios, such as the competing risk model, and to survival data with a cure fraction rate. It can also be applied to joint model frameworks. Other types of censored and truncated observations, such as left censoring, interval censoring, and double censoring, could be used in future research. This is outside the scope of this study and will be addressed in future ones.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
All authors contribute significantly to this manuscript.
Acknowledgments
The first author thanks the Pan African University, Institute for Basic Sciences, Technology and Innovation (PAUSTI), Nairobi, Kenya, for supporting his work. This work was supported by Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2022R299), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.