Abstract
This article investigates a survival analysis under randomly censored mortality distribution. From the perspective of frequentist, we derive the point estimations through the method of maximum likelihood estimation. Furthermore, approximate confidence intervals for the parameters are constructed based on the asymptotic distribution of the maximum likelihood estimators. Besides, two parametric bootstraps are implemented to construct the approximate confidence intervals for the unknown parameters. In Bayesian framework, the Bayes estimates of the unknown parameters are evaluated by applying the Markov chain Monte Carlo technique, and highest posterior density credible intervals are also carried out. In addition, the Bayes inference based on symmetric and asymmetric loss functions is obtained. Finally, Monte Carlo simulation is performed to observe the behavior of the proposed methods, and a real data set of COVID-19 mortality rate is analyzed for illustration.
1. Introduction
An item is said to be randomly censored when the information on time occurrence of an event is not available due to either loss to follow-up, withdrawal, or nonoccurrence of the outcome event before the end of the study. In general, random censoring can be defined as a situation when the items in the experiment are lost or removed from the experiment randomly before its failure. For example, in reliability engineering, an item may have to be removed from the test before its complete failure because of its breakage or for saving time and money. Practically, in a medical study, some patients may leave the course of treatment or die before its completion. Gilbert [1] first introduced the definition of random censoring when items are removed from the test at different random time points.
Thereafter, several authors have studied the statistical inference for different models using random censoring, including Ghitany and Al-Awadhi’s study [2] which discussed maximum likelihood estimation of Burr XII distribution parameters using random censoring. Friesl and Hurt [3] proposed some basic ideas of both the construction and investigation of the properties of the Bayesian estimates of certain parametric functions of the exponential distribution under the model of random censoring. Danish and Aslam [4] discussed the Bayesian inference of the unknown parameters of the randomly censored Weibull distribution. Krishna et al. [5] presented the maximum likelihood and the Bayes estimators for the unknown parameters of Maxwell distribution under random censoring. Krishna and Goel [6] dealt with the geometric distribution under randomly censored data. Danish et al. [7] presented the Bayesian inference for the parameters of randomly censored Burr-type XII distribution. Krishna and Goel [8] dealt with the classical and Bayesian estimation for two-parameter exponential distribution having scale and location parameters with randomly censored data. Kumar [9] discussed classical and Bayesian estimation in log-logistic distribution under random censoring. Kumar and Kumar [10] dealt with the estimation of the parameters and reliability characteristics in inverse Weibull distribution based on the random censoring model. EL-Sagheer et al. [11] obtained the point and interval estimations for the parameters of the three parameter Burr XII distribution based on randomly censored data.
In 1825, the British Benjamin Gompertz presented the law of geometrical progression pervades large portions of different tables of mortality for humans, and statisticians later called this law the Gompertz “mortality” distribution [12]. This distribution occupies an important position in modeling human mortality, income, actuarial chart, and different scientific disciplines fields such as biological and marketing science in addition to reliability analysis and life testing. Many authors have contributed to the statistical methodology and characterization of this distribution by several studies including Wu et al. [13] discussed estimation of the parameters of the Gompertz distribution (GD) under the first-failure-censored sampling plan. Chang and Tsai [14] obtained estimation results concerning a progressively type-II censored sample from the Gompertz distribution. Wu et al. [15] discussed the expected test time for the two-parameter GD under progressive censoring with binomial removals. Soliman et al. [16] investigated Bayes and frequentist estimators for the two-parameter GD as well as the reliability and hazard rate functions using progressive first-failure censoring plan. Mohie El-Din et al. [17] studied the estimation of the parameters for GD using maximum likelihood and Bayesian methods based on type-II progressively hybrid censored data.
Recently, Abu-Moussa et al. [18] combined the adaptive progressive type-II censoring model with the general progressive model, to obtain the estimates for the parameters of GD. Wang and Gui [19] discussed the estimation of the parameters for GD and prediction using general progressive type-II censoring. The probability density function (PDF), cumulative distribution function (CDF), reliability function , and hazard rate function of the GD are given, respectively, bywhere and are the shape and scale parameters, respectively. The GD has an increasing hazard rate function (see Figure 1) and positive skewness; in addition to this, it has a unimodal PDF. It is worth noting that Willekens [20] discussed the relations between other distributions and GD thoroughly, for instance, exponential, Weibull, and extreme value distributions. Schematically, the random censoring can be described as follows: assume items are put on experiment with their lifetimes as which are independent and identically distributed (IID) random variables with PDF and CDF . Also, assume that be their random censoring times with respective PDF and CDF and . Let ’s and ’s be mutually independent so that one observes IID random pairs , where , and is the indicator variable can be defined as if and otherwise . Thus, the joint PDF of and iswhere and . The random variables and satisfy the proportional hazards model with proportionality constant if
From equation (5) and (6), the joint PDF of and can be written as
Hence, from (1), (2) and (7), we obtain
In this article, based on the randomly censored data, we discuss the classical and Bayesian estimation procedures for the unknown parameters and some lifetime indices such as reliability and hazard rate functions of GD. This paper proceeds as follows: Section 2 deals with the maximum likelihood estimate and asymptotic confidence intervals. Two parametric bootstrap methods are presented in Section 3. Bayesian estimates using MCMC technique are provided in Section 4. In Section 5, a simulation study is conducted to evaluate the quality of these approaches. A real data set of COVID-19 is presented to illustrate the application of the proposed inference procedures in Section 6. Finally, the work is concluded in Section 7.
2. Maximum Likelihood Estimation
Let be a randomly censored sample from GD as discussed in Section 1; from (8), the likelihood function can be written as
Thus, the log-likelihood function becomes
Taking partial derivatives of with respect to , , and , respectively, and equating each result to zero, the following likelihood equations can be obtained:
The maximum likelihood estimates (MLEs) of , , and say , , and are the solutions of equations (11)–(13). We plot the profile log-likelihood function of , , and as in Figure 2, to prove the uniqueness property as a simple proof. It is a unimodal function for all parameters. For the solution of the system of these equations, Newton–Raphson iteration method is used to obtain the desired MLEs in such situations. The algorithm is described as follows:(1)Use the method of moments or any other methods to estimate the parameters , , and as starting point of iteration, denote the estimates as , and set .(2)Calculate and the observed Fisher information matrix , given in Section 2.1.(3)Update as(4)Set and then go back to Step 1.(5)Continue the iterative steps until is smaller than a threshold value. The final estimates of , , and are the MLE of the parameters, denoted as , , and .
Once we get desired MLEs of , , and using the invariant property of MLEs, see Zehna, the MLEs of and for given mission time can be obtained as
2.1. Asymptotic Confidence Intervals
In this section, we derive the Fisher information matrix (FIM) for the construction of asymptotic confidence intervals (ACIs) of the parameters , see Zheng and Gastwirth [21]. The FIM is defined by the negative expectation of the partial second derivative of the log-likelihood functionwith
Since MLE has asymptotic normality property under certain regularity conditions, the estimator has asymptotic distribution , see Lawless [22], and the inverse matrix of is
Therefore, two-sided equal tail ACIs for the parameters are given bywhere is the percentile of the standard normal distribution with right-tail probability .
2.2. Delta Method
In order to construct the ACIs of and , we need to find the variances of them. So, Greene [23] implemented the delta method for this purpose as follows.
Assume that be the first partial derivative of , where
Then, the asymptotic variances of are given bywhere is the transpose matrix of . Thus, ACIs for can be written as
3. Bootstrap Technique
When the sample size is small, it is known that the normal approximation is inappropriate. Thus, we propose a resampling technique on the basis of the bootstrap algorithm, including the percentile bootstrap (BP) algorithm and the bootstrap-t (BT) algorithm using the same steps mentioned in DiCiccio and Efron [24].
3.1. Percentile Bootstrap Algorithm
(1)Based on the original data , calculate the MLEs of the parameters , , and expressed as , , and by solving equations (11)–(13).(2)Use , , and to generate a bootstrap sample from GD.(3)Similar to Step (1), based on the bootstrap sample, compute the MLEs of , , and expressed as , and .(4)Due to the invariance property, obtain the MLEs of the reliability and hazard functions using equations (3) and (4), denoted as and , .(5)Repeat Steps (2)–(4) times to obtain based on bootstrap samples, where .(6)Arrange all components in the ascending order, the bootstrap estimates are .(7)Let be the CDF of . Define for given . Then, a two-side confidence intervals of are given by3.2. Bootstrap-t Algorithm
The bootstrap-t algorithm is slightly more complicated compared to the bootstrap-p algorithm, but more accurate when the sample size is small. It also gives better results if we look at the resulting mean squared error. The following steps can be used to obtain the bootstrap-t confidence intervals. (1)–(4): the same as the percentile bootstrap algorithm.(5)Define the statistics , where and obtained using the FIM for , , , and delta method for and .(6)Repeat Steps (2)–(5) times; then, we have .(7)Arrange in the ascending order to obtain .(8)Let be the CDF of . Define for given . Thus, two-side bootstrap-t confidence intervals of or are given by
4. Bayesian Estimation
Bayes estimation characterizes the problems more rationally and reasonably quite different from maximum likelihood and bootstrap methods because it takes into consideration both the information from observed sample data and the prior information. In this section, we consider the Bayesian estimators of , , , , , and the corresponding credible intervals (CRIs) using the Markov chain Monte Carlo (MCMC) technique under both symmetric (squared error (SE)) and asymmetric (linear-exponential (LINEX)) loss functions. In case where underestimation is more serious than overestimation or in other words, when positive and negative estimation errors have different results, in such situations, a strong fine should be put on underestimation. This fine can be obtained by using an asymmetric loss function called LINEX loss function. The LINEX loss function has the following form:where is constant. Then, we havewhere , , , , and . The sign and magnitude of represent the direction and degree of symmetry, respectively. When , the larger loss is incurred for overestimation and the overestimation is more serious than underestimation, and vice-versa, see Calabria and Pulcini [25] for more details. For small values of , the LINEX loss function is approximated as a squared loss function. Thus, selecting appropriate values of , one can assign desired unequal weights to the positive and negative estimation errors.
Regarding variation in values of , it can be seen that, as the loss function departs from symmetry, the risk under the LINEX loss function increases and it may be noted that overestimation results in slightly higher risk than underestimation. For problems on how to choose the LINEX loss parameter values, see [25]. A real problem for the analyst is how to choose the value of the shape parameter of the selected loss function, so that the asymmetry of the loss function reflects the situation. A solution to this problem is when the selected function is the LINEX loss, search the value of that satisfieswhere is the value of the ratio of the loss for overestimation and the loss for underestimation, the error being the same absolute value but with different sign.
The Bayesian deduction requires appropriate choice of priors for the parameters. Arnold and Press [26] pointed out that, from a strict Bayesian viewpoint, there is clearly no way in which one can say that one prior is better than any other. Presumably one has one’s own subjective prior and must live with all of its lumps and bumps. But if we have enough information about the parameter(s), then it is better to make use of the informative prior(s) which may certainly be preferred over all other choices. Otherwise, it may be suitable to resort to use noninformative or vague priors, see [27].
It is known that the family of gamma distributions is flexible enough to cover a large variety of prior beliefs of the experimenter, see [28]. Therefore, we assume that the unknown parameters , , and are stochastically independently distributed with conjugate gamma prior, and the conjugacy is with respect to equation (9), as
Here, the hyperparameters and , , are chosen to reflect the prior knowledge about the unknown parameters. The elicitation of hyperparameter values when the informative priors are taken into account is very important. Notice that the hyperparameter values can be selected based on the past available data. In this regard, suppose that we have number of samples available from G distribution; furthermore, let the associated MLEs of are , . Now, on equating the mean and variance of , with the mean and variance of considered priors, the hyperparameter values can be obtained. In this paper, we have considered the gamma prior , for which mean and variance , ; therefore, for , we have , ; for , we have ; and for , we have . Therefore, on equating the mean and variance of , , with the mean and variance of gamma priors, we get
Now, on solving the above equations, the estimated hyperparameters turn out to have in the following form:
One may also refer to a similar type of discussion presented by Singh and Tripathi [29]. Therefore, the joint prior of the parameters , , and is obtained by
Via Bayes’ theorem across combining the likelihood function given by equation (9) with the joint prior given by equation (33), the joint posterior density function of , , and can be written as follows:
Therefore, the Bayes estimator for any function of , , and , say under SE loss function is the posterior expectation of and is given by
The full conditional posterior densities of , , and can be obtained from equation (34) as follows:
Obviously, it is easy to see that samples of and can be implemented using any gamma generating routine. By observing equation (37), we know that the conditional posterior distribution of cannot be reduced analytically to a common distribution. In addition, we cannot sample directly from it in the usual way. So, Metropolis–Hastings (M-H) algorithm within Gibbs sampling is implemented to generate random samples from the posterior distribution using a proposal density, for more details, see Tierney [30]. The M-H algorithm within Gibbs sampling is set up in the following algorithm steps:(1)Set the initial values and set (2)Generate from gamma distribution (3)Generate from gamma distribution (4)Using M-H algorithm, generate from with the proposal distribution, where represents the variance of in a variance-covariance matrix(i)Generate from .(ii)Evaluate the acceptance probability(iii)Generate a from a Uniform distribution.(iv)If , accept the proposal and set , else set .(5)Calculate and as(6)Set .(7)Repeat Steps times.(8)Under SE and LINEX Loss functions the Bayes estimate of , where , , , and can be obtained by where is the burn-in period, i.e., a number of iterations in MCMC before the stationary distribution is achieved and , , , and .
To compute the highest posterior density CRIs of , order as . Then, symmetric CRIs of are
5. Numerical Explorations
In this section, we conduct a simulation study to compare the performance of the estimates and confidence intervals developed in the preceding sections using different methods. Here, we report the simulation results in the case of , , and . Then, the true values of and at time are evaluated to be 0.952007 and 0.149182, respectively. The performance of estimators are evaluated in terms of mean square error (MSE) which is computed as , where , , , , , , and for the point estimates, also average lengths (ALs) and coverage probability (CPs), which computed as the number of CIs that covered the true values divided by 1000, for interval estimates.
Bayes estimates and the highest posterior density CRIs are computed based on MCMC samples and discard the first values as burn-in. In addition, we assume the informative gamma priors for , , and that is, when the hyperparameters are and , . Moreover, CRIs were computed for each simulated sample. In our study, we consider different sample sizes , and 100. The results of means, MSEs, ALs, and CPs of estimates are displayed in Tables 1–10.
From the results, we observe the following:(i)It is clear that, from Tables 1–10, as sample size increases, the MSEs and ALs decrease, as expected.(ii)Bayes estimates are performed better than the MLEs and bootstrap methods because Bayes estimates have the smallest MSEs and ALs for all parameters and reliability characteristics.(iii)Bayes estimate under LINEX loss function with provides better estimates because of having the smallest MSEs. Hence, the analytical ease with which results can be obtained using asymmetric loss functions makes them attractive for use in applied problems.(iv)Bayes estimates under LINEX loss function for the choice are performed better than their estimates for the choice in the sense of having smaller MSEs.(v)Bayes estimates under SE loss function are equal or slightly larger Bayes estimates under LINEX loss function for the choice in terms of MSEs and ALs.(vi)Bootstrap methods are performed better than the ML method in terms of MSEs and ALs. Also, BT performs better than BP in terms of MSEs and ALs.(vii)The coverage probabilities of the asymptotic confidence intervals and highest posterior density credible intervals are all close to the desired level of 0.95. In addition, the Bayesian credible intervals under LINEX loss function for the choice have the highest coverage probabilities.
6. Real Data Analysis
In this section, for illustrative purposes, a real-life data for COVID-19 mortality rate is applied to inspect the estimation procedures discussed in Section 5. We consider the set of a real-life data, which is reported in Almongy et al. [31]. This data represents a COVID-19 mortality rate belonging to Italy of 59 days, which is recorded from 27 February to April 2020. The observations with sign indicate censored times. Data are shown in Table 11.
Before going further, we fit randomly censored GD for this data set. For the goodness of fit test, we compute the Kolmogorov–Smirnov (KS) distances between the empirical distribution and the fitted distribution functions. The KS is 0.1077, and the associated value is 0.5683. Therefore, according to the value, we can say that the GD fits quite well to the above data. Empirical, and plots are shown in Figure 3, which clears that the GD fits the data very well. Based on the previous sample of randomly censored, the maximum likelihood and two parametric bootstrap estimates for , , , , and are determined to be as in Table 12. Moreover, Table 13 shows the approximate confidence intervals for , , , , and under MLEs and bootstrap estimates (BP and BT).
Now, to compute the Bayesian estimates, the prior distributions of the parameters are needed to specify. Since we have no prior information, we assume that the noninformative gamma priors for , , and , that is, the hyperparameters, are and , . To conduct the MCMC algorithm which is described in Section 4, the initial values for parameters , , and were taken to be their MLEs. In addition, 11,000 MCMC samples were generated. To avoid the effect of the initial values, we expunge the first 1000 samples as “burn-in.” Bayesian estimates as well as CRIs for , , , , and are reported in Tables 12 and 13.
7. Concluding Remarks
The Gompertz distribution is an important lifetime model for representing the unimodal behavior of the failure rate function. In addition, it occupies an important position in mortality modeling, reliability analysis, and life testing. In this article, we have developed different methods to estimate and construct the confidence intervals for the unknown parameters and reliability characteristics of GD under the random censoring model. The MLEs as well as asymptotic confidence intervals for the parameters using expected FIM are obtained. In addition, two bootstrap confidence intervals are constructed. The Bayes estimates as well as the highest posterior density credible intervals using the M-H algorithm within Gibbs sampling are also obtained. The Bayes estimates have been obtained under both SE and LINEX loss functions. A simulation study was implemented to compare the performance of various estimates.
Data Availability
The datasets generated during the current study are available from the corresponding author on reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.