A Flexible Bayesian Parametric Proportional Hazard Model: Simulation and Applications to Right-Censored Healthcare Data

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.


Introduction
e healthcare domain has evolved significantly in recent years as a result of computational developments. e 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]. e 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]. e application of deep learning to the context of parametric survival models was discussed by [4]. rough an efficient training process, the quality of the developed system improves. Data portioning is done three times to confirm the trained algorithms (training-testingvalidation) [5]. e main goal of this article is to present the Bayesian parametric proportional hazard model using BUGs syntax. e 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.
e family of distributions can be parametric, semiparametric, or nonparametric. e 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. e 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]. e 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. e log-logistic (LL) distribution has a wide range of applications in survival data analysis and can accommodate unimodal survival data sets. e 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][13][14][15][16]. e generalized log-logistic distribution (GLL) is one such model, which modifies the log-logistic distribution by inducing an additional shape parameter [17]. e 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. e 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. e 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. e 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. e 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.

PH Model Formulation and Assumptions
In many real-life applications, survival times are affected by explanatory variables. e explanatory variable vector is related to response variable through a regression model. An important aspect of survival modelling is the inclusion of explanatory variables. e 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-toevent 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. e 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].

PH Model Formulation.
e parametric proportional hazard (PH) models are formulated using a defined baseline hazard and a link function ψ(x ′ β) for the covariates which is defined as follows: is a monotone function that has a one -to-one correspondence, (1) e most commonly found option for the link function ψ(x ′ β) is the exponential exp(x ' β) (or log-linear) function. In this work, we define the PH model with the assumption that ψ(x ′ β) � exp(x ′ β).

PH Assumptions.
e 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. e assumption of the PH model can be defined as follows: (2) where h 0 (t) is called the baseline hazard.
Simplifying, we get, e 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, (x and x * ) is e 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 X in the PH model is not assumed to have x ≡ 1. An intercept would get confounded with the baseline hazard function h 0 .

Lifetime Functions Describing the PH Model.
e five frequent representatives of a lifetime distribution function that are used to characterize the PH model are addressed in this section.

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 ). e hazard rate function h(t|x), abbreviated by hrf, also called the instantaneous failure rate or as force of mortality of a PH model is of the form:

Cumulative Hazard Function of the PH Model.
e hazard or survival functions, rather than the cumulative distribution or probability density function, are typically used in the PH analysis of survival data. e 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: Cumulative hazard function: Using the above expressions, we can easily find

Cumulative Distribution Function of the PH Model.
e cdf of the PH model, also known as the lifetime distribution function, is given by

Probability Density Function of the PH Model.
e pdf or the failure density function of the PH model is defined as e 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. e 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. e 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. e 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.

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. ese distributions have been studied and used in various contexts in the literature.

Exponential PH Model.
Exponential distribution is a continuous probability distribution with only one unknown parameter k. It is the simplest distribution for lifetime distribution models. e distribution is not flexible enough to describe commonly encountered hazard rate shapes for survival data. e pdf, cdf, sf, hrf, and chf of the exponential random variable are, respectively, as follows.
Let X ∼ exponential(k), where k > 0 is the scale parameter and t ≥ 0. A short value of k 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 x and link function ψ(x) is Applying the log-linear function ψ(x ′ β) � exp(x ′ β) , we can simplify into (15) In this equation, the hrf has the exponential distribution with scale parameter k. exp(x ′ β) 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. is makes the applicability of this distribution fairly limited. e other lifetime distributions of the exponential PH model are as follows.
e survival function of the exponential PH model is e pdf of the exponential PH model is e cdf of the exponential PH model is e chf of the exponential PH model is 4 Journal of Healthcare Engineering

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. is distribution is flexible and can be skewed both in right and in left. e pdf, cdf, sf, hrf, and chf of the exponential random variable are, respectively, as follows.
Let X ∼ Gompertz(k, α), where k > 0 is the rate parameter, α > 0 is the shape parameter, and t ≥ 0. When k ≥ 0, 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 x and link function ψ(x) is Applying the log-linear function ψ(x ′ β) � exp(x ′ β) , we can simplify into h GoPH (t) � αk.e tk . exp x ′β � αk.e tk . exp β 1 x 1 + β 2 x 2 + . . . + β p x p . (23) 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. e other lifetime distributions of the Gompertz PH model are as follows: the survival function of the Gompertz PH model is e pdf of the Gompertz PH model is e cdf of the Gompertz PH model is e chf of the Gompertz PH model is

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. e additional parameter describes the shape of the hazard functions, based on the value of the shape parameter [21]. e pdf, cdf, sf, hrf, and chf of the Weibull random variable are, respectively, as follows. Let X ∼ Weibull(k, α), where α > 0 is the shape parameter and k > 0 is the rate parameter. e hazard rate function increases when α > 1, decreases for α < 1, and constant for α � 1. When α � 1, 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 x and link function ψ(x) is Applying the log-linear function ψ( In this equation, the model has the Weibull distribution with rate parameter k. exp(x ' β) and shape parameter α which indicates that the PH assumption is satisfied with the Weibull distribution with constant α. e other lifetime distributions of the PH Weibull model are as follows: the survival function of the Weibull PH model is e pdf of the Weibull PH model is Journal of Healthcare Engineering e cdf of the Weibull PH model is e chf of the Weibull PH model is

Parametric Baseline Hazard
e 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. e GLL is a flexible survival distribution proposed by [11]. is 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. ese 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. e 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 k > 0, α > 0, η > 0 is given by e chf of the GLL distribution is given by e distribution function of the GLL model is of the form: e survival function (sf ) of the GLL model is given by where k > 0, α > 0, and η > 0 are parameters and θ � (k, α, η) ' . e quantile function of the GLL model is given by e 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.

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 x and link function ψ(x) is Applying the log-linear function In this equation, the hrf can be recognized as a generalized log-logistic distribution as well, but contrary to (36), the rate parameter is k * � k. exp (x ′ β) 1/α and shape parameters are α and η which indicates that the PH assumption is satisfied with the GLL distribution and the proposed model is closed under the PH framework. e other lifetime distribution functions for the GLL PH model are as follows: the survivor function of the GLL PH model is e pdf of the GLL PH model is e cdf of the GLL PH model is e chf of the GLL PH model is

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.

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 n individuals with lifetimes represented by T 1 , T 2 , . . . , T n .
Assuming that the data are subject to right censoring, we Suppose that a right-censored random sample with data D � (t i , δ i , x i ), i � 1, 2, . . . , n, is available, where t i is the censoring time or a survival time according to whether δ i � 0 or 1, respectively, and x i � x 1 , x 2 , . . . , x n is an n × 1 column vector of external covariates for the i th individual, ϑ is the vector of parameters associated with the baseline distribution, 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 ϑ ofϑ. Hypothesis testing and interval estimations of model parameters are possible due to the MLEs' approaching normality [7]. e 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 ϑ ′ � (k, α, η) is the vector of the baseline distributional parameters.
In our case, if we assume that a � n i�1 δ i ,p i � exp(x i ′ β) and q i � (ηt i ) α . Use (36) for h 0 (.) and note that H 0 (t; θ) � t 0 h(u)du is the baseline cumulative hazard rate function as given by (37). e full log-likelihood function of the GLL PH model can be expressed as follows:  Journal of Healthcare Engineering To obtain the MLE's of θ ′ � (k, α, η) and β ′ , we can maximize (51) directly with respect to (k, α, η) and β ′ or we can solve the nonlinear equations below or the 1 st derivative of the log-likelihood function. e 1 st derivatives of the loglikelihood function are To maximize log-likelihood functions, many software packages are available including proven optimization algorithms.

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.

Priors for the Model Parameters.
e 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. e marginal prior distribution for every regression coefficient β m, m � 1, . . . , 5, is prompted as a normal distribution centred at zero and with a small precision, N(0, 0.001); on the other hand, a gamma distribution, gamma (10,10), 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 loglogistic 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).
e density function of the combined prior distribution of all unknown parameters and the regression coefficients are given as Journal of Healthcare Engineering 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]. e zero's trick approach works on the assumption that perhaps the contribution of a Poisson (λ) observable of zero is exp(− λ); if we set λ � − log(f(t i |ϑ, β, x)) with observable data as a vector of 0 ′ s, we receive the right contributions of the proposed model [18].

e Posterior Distribution.
e joint posterior density function is equal to the multiplication of the prior distribution p(α, k, η, β ′ ) and the likelihood function the joint posterior density function of the parameters α, k, η, an d β ' of GLL PH model given the data can be expressed using Bayes' theorem as where the first four terms on the equation represent the prior specification for the unknown parameters and are assumed to be independent and L(α, k, η, β ' ) is the likelihood function expressed as follows: e 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.

Simulation Study
In this section, we undertake an extensive simulation investigation to demonstrate the proposed parametric proportional hazard model's good Bayesian features. e 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.

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. is 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]. e 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 (α, η, k). Using the reverse chf given in equation (41), lifetime data can be simulated as follows: where a, η, and k > 0.

Simulation
Design. e 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). e GLL PH model's true parameter vector is set as follows: (1) 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.

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 Journal of Healthcare Engineering 9 [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.

Measures of Performance.
e actual mean, standard deviation (SD), Naive standard error, bias, percentage of bias, coverage probability (CP), potential scale reduction factor (R), and the effective number of separate simulation draws were used to test the posterior distribution stability for the suggested PH model.

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.

Accuracy of the Estimators.
e mean square error (MSE) is a good indicator of overall accuracy and is calculated as follows: is metric determines how accurate the estimates are as follows.
e lower the MSE, the more accurate the estimations of impacts.
e 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. e Naïve SE incorporates simulation error rather than posterior uncertainty.

Coverage.
e 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. e more precise the estimations are, the closer the outcome is to a 95 percent coverage probability. e following is how CP is expressed:

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 R < 1.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]. e 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. at is, the McMC simulation for the GLL PH model explores the target posterior distribution appropriately.

Simulation Results
. Tables 1-4 shows the simulation results for the posterior mean, bias, Naive standard error    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. e 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%.

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.

Data Set I: Lung Cancer Survival Data
In this section, we reanalyse the data set reported in [29] which is available in the R package survival. e 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). e 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).

Hazard Rate Shape.
e 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). e 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. e box plot, histogram, and TTT plots are shown in Figure 6.

Proportionality Assumption.
ere 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. e 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.

Posterior Analysis.
In this paper, we assume the noninformative independent framework with a normal prior N(0, 0.001) for β ′ s (regression coefficients) and an independent gamma prior for the distributional parameters α ∼ G(a 1 , b 1 ), η ∼ G(a 2 , b 2 ), and k ∼ G(a 3 , b 3 (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. e posterior summaries for the generalized log-logistic PH model parameters using Veterans lung cancer data sets are illustrated in Table 5. e 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 β 1 , using mean (β 1 > 0).
(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. ese 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. e running mean (also referred to Ergodic mean) is a well-known convergence diagnostic for McMC algorithms. e 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.     Chain   22000  32000  42000  52000 2000  12000  22000  32000  42000  52000 2000  12000  22000  32000  42000  52000   Iteration   1  2  3 alpha beta [1] beta [2] beta [3] beta [4] eta kappa 1 2 3 Figure 11: e running mean plots for the baseline distributional parameters and regression coefficients for the Veterans lung cancer data set.
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. e 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.  Table 7. e 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]. e survival times (in months) of patients is illustrated in Table 7.
e 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. beta [2] beta [4] beta [3] eta kappa 1 2 3 Figure 12: Autocorrelation plots for all the baseline distributional parameters and regression coefficients for the Veterans lung cancer data set.

Hazard Rate Shape.
Based on the TTT plot, the hazard rate function is an increasing hazard in Figure 13. e 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.

Proportionality Assumption.
We investigated if the proportional hazards model could be used with this data set. e 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.

Posterior Analysis.
In this paper, we assume the noninformative independent framework with N(0, 0.001) for β ′ s (regression coefficients) and an independent gamma prior for the distributional parameters α ∼ G(a 1 , b 1 (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. e posterior summaries for the GLL-PH model parameters using larynx cancer data are illustrated in Table 8. e 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. ese plots and graphs provide a nearly comprehensive representation of the parameters' posterior uncertainty.

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.   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 x 1 and x 2 . e HR is defined as 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.

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. is criterion is available in most McMC packages [39]. e DIC is com-  beta [4] beta [5] eta kappa Chain 1 2 3 Figure 18: Autocorrelation plots for all the baseline hazard parameters and regression coefficients for the larynx cancer data. eta kappa beta [5] beta [4] beta [3] beta [2] Parameter beta [ Figure 19: PSRF of the baseline hazard parameters and the regression coefficients for the larynx cancer data.

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. e 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. e 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. e 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. e 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. is is outside the scope of this study and will be addressed in future ones.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.