# Asymmetric Randomly Censored Mortality Distribution: Bayesian Framework and Parametric Bootstrap with Application to COVID-19 Data

Authors: Rashad M. EL-Sagheer; Mohamed S. Eliwa; Khaled M. Alqahtani; Mahmoud EL-Morshedy
Journal: Journal of Mathematics
 (2022)

---

## 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.


---

## Body

## 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 rx, and hazard rate function hx of the GD are given, respectively, by(1)fx|λ,β=λβeβxe−λeβx−1,λ>0,β>0;0<x<+∞,(2)Fx|λ,β=1−e−λeβx−1,0<x<+∞,(3)rx=e−λeβx−1,0<x<+∞,(4)hx=λβeβx,0<x<+∞,where β 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 n items are put on experiment with their lifetimes as X1,X2,…,Xn which are independent and identically distributed (IID) random variables with PDF fXx and CDF FXx. Also, assume that T1,T2,…,Tn be their random censoring times with respective PDF and CDF fTt and FTt. Let Xi’s and Ti’s be mutually independent so that one observes IID random pairs Y1,D1,…,Yn,Dn, where Yi=minXi,Ti, i=1,…,n and Di is the indicator variable can be defined as Di=1 if Xi≤Ti and otherwise Di=0. Thus, the joint PDF of Y and D is(5)fY,Dy,d=fXy1−FTydfTy1−FXy1−d,where y>0 and d=0,1. The random variables X and T satisfy the proportional hazards model with proportionality constant ρ>0 if(6)1−FTt=1−FXtρ.Figure 1
The failure rate function of the GD.From equation (5) and (6), the joint PDF of Y and D can be written as(7)fY,Dy,d=fXy1−FXyρρ1−d.Hence, from (1), (2) and (7), we obtain(8)fY,Dy,d;λ,β,ρ=λβeβxe−λeβx−11+ρρ1−d.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: Section2 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

Lety,d=y1,d1,y2,d2,…,yn,dn be a randomly censored sample from GD as discussed in Section 1; from (8), the likelihood function can be written as(9)Lλ,β,ρ|y,d=λnβneβ∑i=1nxi−λ1+ρ∑i=1neβxi−1ρn−∑i=1ndi.Thus, the log-likelihood functionℓλ,β,ρ=logLα,β,δ,λ;x¯ becomes(10)ℓλ,β,ρ=nlogλ+nlogβ+β∑i=1nxi−λ1+ρ∑i=1neβxi−1+n−∑i=1ndilogρ.Taking partial derivatives ofℓλ,β,ρ with respect to λ, β, and ρ, respectively, and equating each result to zero, the following likelihood equations can be obtained:(11)nλ−1+ρ∑i=1neβxi−1=0,(12)nβ+∑i=1nxi−λ1+ρ∑i=1nxieβxi=0,(13)n−∑i=1ndiρ−λ∑i=1neβxi−1=0.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 λ0,β0,ρ0, and set k=0.(2)
Calculate∂ℓ/∂λ,∂ℓ/∂β,∂ℓ/∂ρλk,βk,ρk and the observed Fisher information matrix I−1λ,β,ρ, given in Section 2.1.(3)
Updateλ,β,ρ as(14)λk+1,βk+1,ρk+1=λk,βk,ρk+∂ℓ∂λ,∂ℓ∂β,∂ℓ∂ρλk,βk,ρk×I−1λ,β,ρ.(4)
Setk=k+1 and then go back to Step 1.(5)
Continue the iterative steps untilλk+1,βk+1,ρk+1−λk,βk,ρk is smaller than a threshold value. The final estimates of λ, β, and ρ are the MLE of the parameters, denoted as λ^, β^, and ρ^.Figure 2
Profile log-likelihood function ofλ, β, and ρ.Once we get desired MLEs ofλ, β, and ρ using the invariant property of MLEs, see Zehna, the MLEs of rt and ht for given mission time t can be obtained as(15)r^t=e−λ^eβ^t−1andh^t=λ^β^eβ^t,t>0.


### 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 function(16)Iθ=E−∂2ℓθ/∂θi∂θji,j=1,2,3=Iij3×3,with(17)∂2ℓθ∂λ2=−nλ2,∂2ℓθ∂λ∂β=∂2ℓθ∂β∂λ=−1+ρ∑i=1nxieβxi,∂2ℓθ∂β2=−nβ2−λ1+ρ∑i=1nxi2eβxi,∂2ℓθ∂β∂ρ=∂2ℓθ∂ρ∂β=−λ∑i=1nxieβxi,(18)∂2ℓθ∂ρ2=−n−∑i=1ndiρ2,∂2ℓθ∂λ∂ρ=∂2ℓθ∂ρ∂λ=−∑i=1neβxi−1.Since MLE has asymptotic normality property under certain regularity conditions, the estimatorθ^=λ^,β^,ρ^ has asymptotic distribution θ^−θ⟶N0,I−1θ, see Lawless [22], and the inverse matrix of Iθ is(19)I−1θ=I−1θ^=varλ^covλ^,β^covλ^,ρ^covβ^,λ^varβ^covβ^,ρ^covρ^,λ^covρ^,β^varρ^.Therefore, two-sided equal tail1001−γ% ACIs for the parameters θ=λ,β,ρ are given by(20)θ^−zγ/2varθ^,θ^+zγ/2varθ^,where zγ/2 is the percentile of the standard normal distribution with right-tail probability γ/2.


### 2.2. Delta Method

In order to construct the ACIs ofrt and ht, we need to find the variances of them. So, Greene [23] implemented the delta method for this purpose as follows.Assume thatG=∂ϑ/∂λ,∂ϑ/∂β,∂ϑ/∂ρ be the first partial derivative of ϑ=rt,ht, where(21)∂rt∂λ=−eβt−1e−λeβt−1,∂rt∂β=−λteβte−λeβt−1,∂rt∂ρ=0,(22)∂ht∂λ=βeβt,∂ht∂β=λβ2+1eβt,∂ht∂ρ=0.Then, the asymptotic variances ofϑ^=r^t,h^t are given by(23)varϑ^≃GTI−1θ^G↓θ=θ^,where GT is the transpose matrix of G. Thus, 1−γ100% ACIs for ϑ^=r^t,h^t can be written as(24)ϑ^−zγ/2varϑ^,ϑ^+zγ/2varϑ^.


## 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 function(16)Iθ=E−∂2ℓθ/∂θi∂θji,j=1,2,3=Iij3×3,with(17)∂2ℓθ∂λ2=−nλ2,∂2ℓθ∂λ∂β=∂2ℓθ∂β∂λ=−1+ρ∑i=1nxieβxi,∂2ℓθ∂β2=−nβ2−λ1+ρ∑i=1nxi2eβxi,∂2ℓθ∂β∂ρ=∂2ℓθ∂ρ∂β=−λ∑i=1nxieβxi,(18)∂2ℓθ∂ρ2=−n−∑i=1ndiρ2,∂2ℓθ∂λ∂ρ=∂2ℓθ∂ρ∂λ=−∑i=1neβxi−1.Since MLE has asymptotic normality property under certain regularity conditions, the estimatorθ^=λ^,β^,ρ^ has asymptotic distribution θ^−θ⟶N0,I−1θ, see Lawless [22], and the inverse matrix of Iθ is(19)I−1θ=I−1θ^=varλ^covλ^,β^covλ^,ρ^covβ^,λ^varβ^covβ^,ρ^covρ^,λ^covρ^,β^varρ^.Therefore, two-sided equal tail1001−γ% ACIs for the parameters θ=λ,β,ρ are given by(20)θ^−zγ/2varθ^,θ^+zγ/2varθ^,where zγ/2 is the percentile of the standard normal distribution with right-tail probability γ/2.


## 2.2. Delta Method

In order to construct the ACIs ofrt and ht, we need to find the variances of them. So, Greene [23] implemented the delta method for this purpose as follows.Assume thatG=∂ϑ/∂λ,∂ϑ/∂β,∂ϑ/∂ρ be the first partial derivative of ϑ=rt,ht, where(21)∂rt∂λ=−eβt−1e−λeβt−1,∂rt∂β=−λteβte−λeβt−1,∂rt∂ρ=0,(22)∂ht∂λ=βeβt,∂ht∂β=λβ2+1eβt,∂ht∂ρ=0.Then, the asymptotic variances ofϑ^=r^t,h^t are given by(23)varϑ^≃GTI−1θ^G↓θ=θ^,where GT is the transpose matrix of G. Thus, 1−γ100% ACIs for ϑ^=r^t,h^t can be written as(24)ϑ^−zγ/2varϑ^,ϑ^+zγ/2varϑ^.


## 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 datay1,d1,y2,d2,…,yn,dn, 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 r^∗t=e−λ^∗eβ^∗t−1 and h^∗t=λ^∗β^∗eβ^∗t, t>0.(5)
Repeat Steps (2)–(4)N times to obtain η^1∗,η^2∗,…,η^N∗ based on N bootstrap samples, where η^i∗=λ^∗,β^∗,ρ^∗,r^∗t,h^∗t.(6)
Arrange all components in the ascending order, the bootstrap estimates areη^1∗,η^2∗,…,η^N∗.(7)
Letφx=Pη^∗≤x be the CDF of η^∗. Define η^BPx=φ−1x for given x. Then, a two-side 1001−γ% confidence intervals of η are given by(25)η^BP∗Nγ2,η^BP∗N1−γ2.


### 3.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 statisticsTη=η^∗−η^/varη^∗, where η^∗=λ^∗,β^∗,ρ^∗,r^∗t,h^∗t and varη^∗ obtained using the FIM for λ, β, ρ, and delta method for rt and ht.(6)
Repeat Steps (2)–(5)N times; then, we have T1η,T2η,…,TNη.(7)
ArrangeT1η,T2η,…,TNη in the ascending order to obtain T1η,T2η,…,TNη.(8)
Letϕx=PTη≤x be the CDF of Tη. Define η^BTx=η^+ϕ−1xvarη^∗ for given x. Thus, two-side 1001−γ% bootstrap-t confidence intervals of η=λ,β,ρ,rt or ht are given by(26)η^BTNγ2,η^BTN1−γ2


## 3.1. Percentile Bootstrap Algorithm

(1)
Based on the original datay1,d1,y2,d2,…,yn,dn, 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 r^∗t=e−λ^∗eβ^∗t−1 and h^∗t=λ^∗β^∗eβ^∗t, t>0.(5)
Repeat Steps (2)–(4)N times to obtain η^1∗,η^2∗,…,η^N∗ based on N bootstrap samples, where η^i∗=λ^∗,β^∗,ρ^∗,r^∗t,h^∗t.(6)
Arrange all components in the ascending order, the bootstrap estimates areη^1∗,η^2∗,…,η^N∗.(7)
Letφx=Pη^∗≤x be the CDF of η^∗. Define η^BPx=φ−1x for given x. Then, a two-side 1001−γ% confidence intervals of η are given by(25)η^BP∗Nγ2,η^BP∗N1−γ2.


## 3.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 statisticsTη=η^∗−η^/varη^∗, where η^∗=λ^∗,β^∗,ρ^∗,r^∗t,h^∗t and varη^∗ obtained using the FIM for λ, β, ρ, and delta method for rt and ht.(6)
Repeat Steps (2)–(5)N times; then, we have T1η,T2η,…,TNη.(7)
ArrangeT1η,T2η,…,TNη in the ascending order to obtain T1η,T2η,…,TNη.(8)
Letϕx=PTη≤x be the CDF of Tη. Define η^BTx=η^+ϕ−1xvarη^∗ for given x. Thus, two-side 1001−γ% bootstrap-t confidence intervals of η=λ,β,ρ,rt or ht are given by(26)η^BTNγ2,η^BTN1−γ2


## 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λ, β, ρ, rt, ht, 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:(27)ϕθ^,θ=eqθ^−θ−qθ^−θ−1,where q≠0 is constant. Then, we have(28)Λ^k=−1qlnEe−qΛk,q≠0,k=1,…,5,where Λ1=λ, Λ2=β, Λ3=ρ, Λ4=rt, and Λ5=ht. The sign and magnitude of q represent the direction and degree of symmetry, respectively. When q>0, 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 q, the LINEX loss function is approximated as a squared loss function. Thus, selecting appropriate values of q, one can assign desired unequal weights to the positive and negative estimation errors.Regarding variation in values ofq, 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 q that satisfies(29)ϕθ+δ,θϕθ−δ,θ=r,δ>0,where r 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(30)π1λ=b1a1Γa1λa1−1e−b1λ,λ>0,a1>0,b1>0,π2β=b2a2Γa1βa2−1e−b2β,β>0,a2>0,b2>0,π3ρ=b3a3Γa3ρa3−1e−b3ρ,ρ>0,a3>0,b3>0.Here, the hyperparametersai and bi, i=1,2,3, 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 S number of samples available from Gλ,β distribution; furthermore, let the associated MLEs of λ,β,ρ are λj,βj,ρj, j=1,2,…,S. Now, on equating the mean and variance of λj,βj,ρj, j=1,2,…,S with the mean and variance of considered priors, the hyperparameter values can be obtained. In this paper, we have considered the gamma prior πζ∝ζh1−1e−h1ζ, for which mean ζ=h1/h2 and variance ζ=h1/h22, ζ=h1/h22; therefore, for ζ=λ, we have h1=a1andh2=b1, h1=a1andh2=b1; for ζ=β, we have h1=a2andh2=b2; and for ζ=ρ, we have h1=a3andh2=b3. Therefore, on equating the mean and variance of ζ^j, j=1,2,…,S, with the mean and variance of gamma priors, we get(31)1S∑j=1Sζ^j=h1h2,1S−1ζ^j−1S∑j=1Sζ^j2=h1h22.Now, on solving the above equations, the estimated hyperparameters turn out to have in the following form:(32)h1=1/S∑j=1Sζ^j21/S−1ζ^j−1/S∑j=1Sζ^j2,h2=1/S∑j=1Sζ^j1/S−1ζ^j−1/S∑j=1Sζ^j2.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(33)πλ,β,ρ∝λa1−1βa2−1ρa3−1e−b1λ+b2β+b3ρ.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:(34)π∗λ,β,ρ|y,d=Lλ,β,ρ|y,dπλ,β,ρ∫0∞∫0∞∫0∞Lλ,β,ρ|y,dπλ,β,ρdλdβdρ∝λn+a1−1βn+a2−1ρn−∑i=1ndi+a3−1e−b1λ+b2β+b3ρ×eβ∑i=1nxi−λ1+ρ∑i=1neβxi−1.Therefore, the Bayes estimator for any function ofλ, β, and ρ, say gλ,β,ρ under SE loss function is the posterior expectation of gλ,β,ρ and is given by(35)Egλ,β,ρ|y,d=∫0∞∫0∞∫0∞gλ,β,ρLλ,β,ρ|y,dπλ,β,ρdλdβdρ∫0∞∫0∞∫0∞Lλ,β,ρ|y,dπλ,β,ρdλdβdρ.The full conditional posterior densities ofλ, β, and ρ can be obtained from equation (34) as follows:(36)π1∗λ|β,ρ,y,d∝λn+a1−1e−λ1+ρ∑i=1neβxi−1+b1,(37)π2∗β|λ,ρ,y,d∝βn+a2−1e−βb2−∑i=1nxi−λ1+ρ∑i=1neβxi,(38)π3∗ρ|λ,β,y,d∝ρn−∑i=1ndi+a3−1e−ρλ∑i=1neβxi−1+b3.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λ0,β0,ρ0=λ^,β^,ρ^ and set k=1(2)
Generateλk from gamma distribution π1∗λ|β,ρ,y,d(3)
Generateρk from gamma distribution π3∗ρ|λ,β,y,d(4)
Using M-H algorithm, generateβk from π2∗β|λ,ρ,y,d with the Nβk−1,varβ^ proposal distribution, where varβ^ represents the variance of β^ in a variance-covariance matrix(i)
Generateβ∗ from Nβk−1,varβ^.(ii)
Evaluate the acceptance probability(39)Ωβ=min1,π2∗β∗|λk,ρk,y,dπ2∗βk−1|λk,ρk,y,d.(iii)
Generate aυ from a Uniform 0,1 distribution.(iv)
Ifυ<Ωβ, accept the proposal and set βk=β∗, else set βk=βk−1.(5)
Calculatert and ht as(40)rkt=e−λkeβkt−1,hkt=λkβkeβkt,t>0.(6)
Setk=k+1.(7)
Repeat Steps2−5N times.(8)
Under SE and LINEX Loss functions the Bayes estimate ofζ, where ζ=λ, β, ρ, rt and ht can be obtained by(41)ζ^BS=1N−M∑k=M+1Nζk,ζ^BL=−1qlog1N−M∑k=M+1Ne−qζk,q≠0,whereM is the burn-in period, i.e., a number of iterations in MCMC before the stationary distribution is achieved and ζk=λk, βk,ρk, rkt, and hkt.To compute the highest posterior density CRIs ofζ, order ζM+1,ζM+2,…,ζN as ζ1,ζ2,…,ζN−M. Then, 1−γ100% symmetric CRIs of ζ are(42)ζN−Mγ/2,ζN−M1−γ/2.


## 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λ=0.1, β=1, and ρ=0.774. Then, the true values of rt and ht at time t=0.4 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 MSE=1/N∑k=1Nφ^ik−φi2, where N=1000, i=1,2,…,5, φ1=λ, φ2=β, φ3=ρ, φ4=rt, and φ5=ht 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 onN=11000 MCMC samples and discard the first values M=1000 as burn-in. In addition, we assume the informative gamma priors for λ, β, and ρ that is, when the hyperparameters are ai=0.3 and bi=3.7, i=1,2,3. Moreover, 95% CRIs were computed for each simulated sample. In our study, we consider different sample sizes n=15,20,30,40,50,60,70,80,90, and 100. The results of means, MSEs, ALs, and CPs of estimates are displayed in Tables 1–10.Table 1
Means estimates ofλ (first row) with their MSEs (second row).
nMLEBootstrapBayesianBPBTSELINEXq = −1q = 1150.124530.112420.110120.112340.111250.109450.007580.007160.006990.006820.006850.00578200.102140.102450.112410.132650.113540.101240.007120.007090.006520.006430.006470.00553300.123650.113540.128740.108970.109620.099980.007010.006990.006340.006210.006250.00504400.114700.112450.111190.109660.100970.100070.006850.006370.006020.005870.005900.00488500.109990.112410.108940.100080.099970.098890.005570.005240.005010.004880.004890.00442600.110020.109810.099780.102340.110210.110010.004960.004810.004650.004470.004510.00410700.112140.102580.103690.102540.112360.111900.003880.003750.003660.003570.003590.00311800.123500.114520.117820.104890.108990.102250.003330.003270.003150.002990.003020.00279900.099890.102350.099870.100230.100790.099840.002950.002880.002750.002660.002690.002181000.100890.111000.112400.103690.102500.100560.001980.001750.001530.001360.001380.00099Table 2
ALs and CPs of estimates for parameterλ.
nMLEBPBTBayesianACIsACIsACIsCRIsALsCPsALsCPsALsCPsALsCPs150.552120.9280.516320.9290.487510.9380.399840.957200.513620.9390.495610.9560.471230.9440.369910.949300.488720.9410.467810.9380.449630.9510.330120.948400.435680.9520.412350.9540.399870.9480.310010.964500.382650.9450.376540.9510.356140.9510.289990.965600.357910.9370.341120.9520.336510.9290.259440.948700.310070.9510.299460.9470.283690.9580.213480.958800.284560.9440.263640.9390.243680.9450.195470.961900.246810.9530.225890.9410.214780.9390.169990.9521000.199870.9540.184590.9530.179840.9590.120970.963Table 3
Means estimates ofβ (first row) with their MSEs (second row).
nMLEBootstrapBayesianBPBTSELINEXq=−1q=1151.124531.24571.124731.089791.124781.096540.080220.079980.077650.075470.076000.07123201.096541.112381.092451.084511.045120.999880.078630.075360.072390.069740.070240.06784300.102561.036910.999790.103641.100681.00790.074350.071470.068470.065380.065990.06215401.036901.099971.047960.999661.005791.008990.071230.068450.065470.062550.062610.05874501.023681.024571.124701.032891.025471.024580.066320.063240.059450.056380.056470.05124601.100911.023980.999121.100321.006870.997450.061580.057330.053650.051470.051550.04762701.023581.065871.036980.989991.100011.004780.055640.053210.049740.046580.047010.04139801.036251.100781.124701.056891.024780.999840.051000.047690.045680.044120.044230.03761901.035461.024781.089741.104571.089741.096540.043660.041240.038840.037120.037650.029841000.999741.024870.995871.002350.999471.00030.031810.030010.029950.025640.025710.02093Table 4
ALs and CPs of estimates for parameterβ.
nMLEBPBTBayesianACIsACIsACIsCRIsALsCPsALsCPsALsCPsALsCPs152.023550.9391.824510.9281.636540.9291.336690.947201.845620.9371.547660.9291.365470.9381.100240.943301.547210.9291.236540.9391.102350.9290.969890.951401.365810.9471.036990.9421.000950.9480.914780.948501.136920.9410.999880.9510.992540.9510.899780.951601.100020.9450.984570.9540.979840.9470.856340.961700.987410.9510.936840.9490.913640.9480.802380.956800.925870.9480.872350.9520.836470.9470.779940.959900.896720.9520.834720.9510.793690.9510.728910.9491000.823600.9550.783640.9540.736720.9530.696650.962Table 5
Means estimates ofρ (first row) with their MSEs (second row).
nMLEBootstrapBayesianBPBTSELINEXq=−1q=1150.847920.865720.864750.798650.814570.786910.022750.021650.019970.018360.018530.01589200.812450.823410.802230.798450.814730.779990.021360.019940.018360.017540.017850.01423300.798990.812470.805690.817840.785490.795680.019650.017870.016390.015880.015870.01286400.802140.789890.798640.800810.778990.769540.016990.015350.014680.013570.013660.01102500.799840.789640.812000.802450.798420.773690.013650.012580.011890.011010.011030.00998600.802360.801120.798540.776390.802140.789650.011890.010870.010230.009870.009760.00921700.809640.815740.802470.802350.796960.778930.009990.009530.009260.009110.009180.00872800.798990.810470.759990.774580.762980.789210.009470.009160.008940.008850.008860.00823900.802690.779890.765830.802140.801010.774120.008210.007690.007540.007370.007400.006841000.800930.798450.814750.798650.800240.769890.006950.006130.005860.005570.005590.00499Table 6
ALs and CPs of estimates for parameterρ.
nMLEBPBTBayesianACIsACIsACIsCRIsALsCPsALsCPsALsCPsALsCPs151.456720.9281.365440.9311.245730.9290.987450.945201.354780.9361.245780.9411.136540.9510.887460.949301.136920.9391.036470.9401.005870.9450.813650.952401.098470.9451.003640.9450.997840.9510.799950.948500.999870.9370.988470.9510.976870.9380.776580.951600.956540.9450.945820.9490.936520.9520.754410.961700.912350.9510.899450.9520.887450.9530.729360.956800.889360.9480.864510.9390.852470.9450.700110.963900.846920.9540.816820.9510.810010.9570.687450.9591000.813650.9550.799870.9540.785430.9610.661770.964Table 7
Means estimates ofrt (first row) with their MSEs (second row).
nMLEBootstrapBayesianBPBTSELINEXq=−1q=1150.965470.959840.948750.958880.947790.954670.054570.051340.049870.048890.048990.04564200.972230.964520.975800.958790.948910.949960.052230.049650.047650.046870.046910.04265300.964580.958870.953210.957830.947820.942310.048230.046840.044650.043550.043670.04087400.939980.957810.947830.939450.947810.939950.045690.044020.041890.040990.041010.03879500.945510.958210.963210.958740.964570.948510.042670.041880.040570.039250.039270.03768600.973650.968750.954780.961120.945710.951240.039550.038740.037980.036590.036770.03499700.963210.947890.978410.963250.947820.939980.037680.035660.034580.033560.033980.03197800.948720.958740.964750.978450.973330.961120.034620.033670.032580.031560.031740.02982900.957800.965470.973620.945210.939990.946590.031770.029990.028740.027980.027990.026951000.954790.947820.958870.949970.949990.951470.027640.025670.024320.022780.022940.02081Table 8
ALs and CPs of estimates for parameterrt.
nMLEBPBTBayesianACIsACIsACIsCRIsALsCPsALsCPsALsCPsALsCPs150.244560.9250.226580.9250.214580.9310.198410.949200.232450.9350.215430.9290.199450.9390.179680.951300.213560.9390.196580.9550.178920.9450.159470.948400.197140.9390.173620.9410.156780.9380.140210.961500.174520.9510.163890.9390.145670.9520.128990.949600.161150.9470.156970.9470.132280.9470.110900.957700.146580.9520.136210.9370.118940.9450.109950.974800.123690.9480.120110.9470.108580.9520.099970.965900.113590.9490.110010.9510.100020.9610.984550.9671000.105640.9510.101120.9490.099980.9580.975640.971Table 9
Means estimates ofht (first row) with their MSEs (second row).
nMLEBootstrapBayesianBPBTSELINEXq=−1q=1150.164120.157840.149780.158990.149840.146990.021450.020140.019650.018870.018950.01663200.154750.158820.152340.163540.147920.139450.019450.018650.017630.016580.016720.01495300.139940.152340.142310.150020.139450.142210.018240.017360.016870.015780.015990.01302400.152360.142380.138740.143560.138870.149250.016940.015950.014960.013560.013770.01208500.151230.154780.139650.142870.139950.152000.015610.014970.013550.012780.012950.01099600.160110.154260.147300.152410.138740.136580.013870.012990.012010.011860.011930.01001700.148710.152230.160020.148720.152470.138970.011860.011240.010870.010210.010340.00991800.139990.142780.135620.147640.152200.139090.010870.010550.010210.009940.009980.00982900.145580.152230.144490.145660.139980.145610.010030.009970.009910.009830.009850.009731000.148760.152000.145780.138790.142210.137890.009920.009840.009780.009720.009740.00965Table 10
ALs and CPs of estimates for parameterht.
nMLEBPBTBayesianACIsACIsACIsCRIsALsCPsALsCPsALsCPsALsCPs150.324510.9270.299450.9280.284570.9320.253610.949200.293620.9350.274510.9350.269140.9480.236630.951300.263690.9370.258940.9480.246850.9510.216540.962400.236140.9410.225430.9450.214580.9390.198730.948500.215460.9390.205640.9480.199840.9440.169450.952600.186930.9570.176350.9510.165480.9540.136570.959700.165640.9380.158970.9460.148720.9480.125400.961800.153620.9410.145620.9370.135470.9410.115570.963900.136580.9510.124780.9450.120010.9390.100040.9591000.117490.9590.110220.9520.108990.9530.099930.961From the results, we observe the following:(i)
It is clear that, from Tables1–10, as sample size n 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 withq=1 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 choiceq=1 are performed better than their estimates for the choice q=−1 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 choiceq=−1 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 choiceq=1 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 Section5. 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 27 April 2020. The observations with + sign indicate censored times. Data are shown in Table 11.Table 11
The COVID-19 mortality rate data.
4.571+7.2013.6068.479+11.4108.961+10.91910.9086.503+18.47411.01017.33716.561+13.22615.137+8.69715.78713.33311.82214.242+11.27314.330+16.04611.95010.282+11.775+10.1389.037+12.39610.6448.646+8.9058.906+7.407+7.4457.2146.194+4.6405.452+5.073+4.4164.859+4.4084.6393.148+4.040+4.2534.011+3.5643.8273.134+2.780+2.881+3.3412.686+2.814+2.508+2.450+1.518+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 associatedp value is 0.5683. Therefore, according to the p value, we can say that the GD fits quite well to the above data. Empirical, Q−Q and P−P 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 λ, β, ρ, r0.4, and h0.4 are determined to be as in Table 12. Moreover, Table 13 shows the 95% approximate confidence intervals for λ, β, ρ, r0.4, and h0.4 under MLEs and bootstrap estimates (BP and BT).Figure 3
Graphical fitting of the GD.Table 12
ML, bootstrap, and Bayes estimates.
ParameterMLEBootstrapBayesianBPBTSELINEXq=−1q=1λ0.982720.975320.970270.969910.975470.94243β0.051740.521410.049820.047630.49960.04258ρ0.925930.886870.911050.876950.896470.78699rt0.979710.965570.958410.954460.964750.95004ht0.051830.049880.051120.048630.049660.04432Table 13
95% confidence and credible intervals of λ, β, ρ, rt, and ht.
ParameterMLEBPBTBayesianACIsACIsACIsCRIsLowerUpperLowerUpperLowerUpperLowerUpperλ0.42141.54400.51471.76320.34711.35470.52391.2775β0.00611.00030.00720.99570.00480.98650.00320.8763ρ0.31961.27540.45321.42370.41291.23440.28741.0069rt0.95621.00310.94791.01210.93320.99970.93990.9898ht0.00780.47310.00570.39980.00420.35290.00410.3411Now, 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 ai=0.0001 and bi=0.0001, i=1,2,3. 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 95% CRIs for λ, β, ρ, r0.4, and h0.4 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.



---

*Source: 8300753-2022-03-22.xml*