# EM Algorithm for Estimating the Parameters of Weibull Competing Risk Model

Authors: Mohamed Kayid
Journal: Applied Bionics and Biomechanics
 (2021)

---

## Abstract

One of the most commonly used models in survival analysis is the additive Weibull model and its generalizations. They are well suited for modeling bathtub-shaped hazard rates that are a natural form of the hazard rate. Although they have some advantages, the maximum likelihood and the least square estimators are biased and have poor performance when the data set contains a large number of parameters. As an alternative, the expectation-maximization (EM) algorithm was applied to estimate the parameters of the additive Weibull model. The accuracy of the parameter estimates and the simulation study confirmed the advantages of the EM algorithm.


---

## Body

## 1. Introduction

There are many situations in which the hazard rate function shows a bathtub shape (BTS) with three life periods, early decreasing hazard rate period, useful period (where the hazard rate is approximately constant), and eventually increasing hazard rate period. Lai et al. [1] and Nadarajah [2] presented a list of distributions with BTS hazard rate. Many authors assumed a BTS hazard rate model in their research, and among them, we can point to Block et al. [3], Glaser [4], Leemis and Beneke [5], Mi [6, 7], Mitra and Basu [8], and Noughabi et al. [9].Xie and Lai [10] and Jiang and Murthy [11, 12] studied some aspects of the additive Weibull model with the hazard rate function:
(1)hx=α1λ1α1xα1−1+α2λ2α2xα2−1,x≥0,as a good candidate for describing the bathtub-shaped failure rate function. They urged that the statistical inference about this model is complex due to the number of the parameters, and as a remedy, they suggested the reduced model by considering λ1=λ2 and α2=1/α1 which still accommodates the bathtub-shaped hazard rate.Lai et al. [13] added one constant magnitude to the additive Weibull hazard rate (1) to provide more realistic model:
(2)hx=α1λ1α1xα1−1+α2λ2α2xα2−1+λ3,x≥0.In addition, Bebbington et al. [14] considered the additive Weibull models (1) and (2) in their research to express the concept of the useful period of life of a bathtub-shaped hazard rate distribution. They concluded that the additive Weibull model is sufficiently flexible to describe a bathtub-shaped hazard rate. The common estimator of the parameters of these models is the maximum likelihood estimator (MLE).The EM algorithm is an iterative algorithm for estimating parameters of models involving latent variables, e.g., when data is derived from a mixture or competing risk model. It was used by Dempster et al. [15], Balakrishnan et al. [16], Davies et al. [17], Yang et al. [18], and Okamura and Dohi [19] to estimate the parameters in their models. In this paper, we use the EM algorithm to estimate the parameters of the additive model (2) and show by a simulation study that this algorithm gives a better estimate than the MLE and the least square estimator (LSE).The paper has been organized as the following. In Section2, we present a short representation of the MLE and LSE of the parameters. Then, the EM algorithm has been discussed. Section 3 provides a simulation study for comparing the results related to MLE, LSE, and EM estimator. In Section 4, a data set has been analyzed to show the applicability of the proposed estimators.


## 2. The MLE, Least Square, and EM

Letx1,x2,⋯,xn represents a realization of an iid random sample of the additive Weibull hazard rate distribution with the hazard rate (2). The log-likelihood function has five parameters and is of the form
(3)lα1,λ1,α2,λ2,λ3=∑i=1nlnα1λ1α1xiα1−1+α2λ2α2xiα2−1+λ3−∑i=1nλ1xiα1+λ2xiα2+λ3xi.To find the MLE, we should find admissible values ofα1,λ1,α2,λ2,λ3 which maximize the log-likelihood function.To find the LSE of the parameters, we apply the reliability function related to the hazard rate model (2) which is(4)Rx=e−λ1x1α−λ2x2α−λ3x,x≥0.The empirical reliability function is defined to be(5)R~t=1n∑i=1nIt<xi,in which the indicator function It<xi equals 1 when t<xi and otherwise is 0. So, the LSE of the parameters can be computed by minimizing the following sum of squares of errors in terms of α1,λ1,α2,λ2,λ3.
(6)S2=∑i=1nR~xi−e−λ1xiα1−λ2xiα2−λ3xi2.


### 2.1. The EM Algorithm

LetX1i, X2i, and X3i, i=1,2,⋯,n follows from the Weibull distributions with parameters α1,λ1 and α2,λ2 and the exponential distribution with mean 1/λ3, respectively. Assume that in a lifetime experiment, the observations are realizations of the competing risk random variable Xi=minX1i,X2i,X3i. This means that the lifetime event may be due to one of three competing causes. Let the latent random variable Zi with the support 1,2,3 such that
(7)Zi=1whenX1i≤minX2i,X3i,2whenX2i≤minX1i,X3i,3whenX3i≤minX1i,X2i.With these notations, the likelihood function is(8)Lθ;x,z=∏i=1n∏j=13fZi=j∣Xi=xi,θIZi=jfXi=xi∣θ=∏i=1n∏j=13λjxi∣θIZi=jλ1xi∣θ+λ2xi∣θ+λ3xi∣θ×f1xi∣θR2xi∣θR3xi∣θ+f2xi∣θR1xi∣θR3xi∣θ+f3xi∣θR1xi∣θR2xi∣θ=∏i=1n∏j=13λjxi∣θIZi=jR1xi∣θR2xi∣θR3xi∣θ,where Rjx∣θ, fjx∣θ, and λjx∣θ show the corresponding reliability function, the density function, and the hazard rate function of Xji, j=1,2,3. Then, the log-likelihood function is
(9)lθ;x,z=lnLθ;x,z=∑i=1n∑j=13IZi=jlnλjxi∣θ+∑i=1n∑j=13lnRjxi∣θ.The EM algorithm is an iterative algorithm, and every iteration of it consists of two consecutive steps, namely, the E step and the M step. In the E step, the expectation of the log-likelihood with respect to the estimate of the conditional probabilities of the latent variables has been constructed. Then, in the M step, the constructed expectation of the E step is maximized to compute the estimate of the parameters in the current iteration.


### 2.2. The E Step

Suppose that the estimate ofθ at iteration t be denoted by θt, then the conditional distribution of Zi is
(10)pi1,t=PZi=1∣Xi=xi,θt=f1xi∣θtR2xi∣θtR3xi∣θtfxi∣θt,(11)pi2,t=PZi=2∣Xi=xi,θt=f2xi∣θtR1xi∣θtR3xi∣θtfxi∣θt,and pi3,t=1−pi1,t−pi2,t where
(12)fxi∣θt=f1xi∣θtR2xi∣θtR3xi∣θt+f2xi∣θtR1xi∣θtR3xi∣θt+f3xi∣θtR1xi∣θtR2xi∣θt.The probabilitiespi1,t, pi2,t, and pi3,t are called the membership probabilities.Now we define the expectation of the log-likelihood function (8) with respect to the conditional distribution ofZi.
(13)Qθ∣θt=EZ∣X,θtlθ;x,Z=∑i=1np1i,tlnα1λ1α1xiα1−1−λ1xiα1−λ2xiα2−λ3xi+∑i=1np2i,tlnα2λ2α2xiα2−1−λ1xiα1−λ2xiα2−λ3xi+∑i=1np3i,tlnλ3−λ1xiα1−λ2xiα2−λ3xi=∑i=1np1i,tlnα1λ1α1xiα1−1−∑i=1nλ1xiα1+∑i=1np2i,tlnα2λ2α2xiα2−1−∑i=1nλ2xiα2+∑i=1np3i,tlnλ3−∑i=1nλ3xi.So theQθ∣θt can be written as sum of three distinct expressions Q1θ∣θt, Q2θ∣θt, and Q3θ∣θt, where
(14)Q1θ∣θt∑i=1np1i,tlnα1λ1α1xiα1−1−∑i=1nλ1xiα1,(15)Q2θ∣θt=∑i=1np2i,tlnα2λ2α2xiα2−1−∑i=1nλ2xiα2,(16)Q3θ∣θt=∑i=1np3i,tlnλ3−∑i=1nλ3xi.These statements will be applied in the M step, to compute the estimates of the parameters.


### 2.3. The M Step

The estimate of the parameters at iterationt+1 can be obtained by
(17)α^1,t+1,λ^1,t+1=argmaxα,λQ1θ∣θt,(18)α^2,t+1,λ^2,t+1=argmaxα,λQ2θ∣θt,(19)λ^3,t+1=argmaxλQ3θ∣θt.We should optimizeQ1θ∣θt and Q2θ∣θt numerically since they have not closed form for their critical points. But, the point which maximizes Q3θ∣θt has a closed form, and by solving the equation ∂/∂λ3Q3=0, we have
(20)λ^3,t+1=∑i=1np3i,t∑i=1nxi.The iterative process can be concluded ifQθt+1∣θt+1<Qθt∣θt+ε for some small predefined ε.


## 2.1. The EM Algorithm

LetX1i, X2i, and X3i, i=1,2,⋯,n follows from the Weibull distributions with parameters α1,λ1 and α2,λ2 and the exponential distribution with mean 1/λ3, respectively. Assume that in a lifetime experiment, the observations are realizations of the competing risk random variable Xi=minX1i,X2i,X3i. This means that the lifetime event may be due to one of three competing causes. Let the latent random variable Zi with the support 1,2,3 such that
(7)Zi=1whenX1i≤minX2i,X3i,2whenX2i≤minX1i,X3i,3whenX3i≤minX1i,X2i.With these notations, the likelihood function is(8)Lθ;x,z=∏i=1n∏j=13fZi=j∣Xi=xi,θIZi=jfXi=xi∣θ=∏i=1n∏j=13λjxi∣θIZi=jλ1xi∣θ+λ2xi∣θ+λ3xi∣θ×f1xi∣θR2xi∣θR3xi∣θ+f2xi∣θR1xi∣θR3xi∣θ+f3xi∣θR1xi∣θR2xi∣θ=∏i=1n∏j=13λjxi∣θIZi=jR1xi∣θR2xi∣θR3xi∣θ,where Rjx∣θ, fjx∣θ, and λjx∣θ show the corresponding reliability function, the density function, and the hazard rate function of Xji, j=1,2,3. Then, the log-likelihood function is
(9)lθ;x,z=lnLθ;x,z=∑i=1n∑j=13IZi=jlnλjxi∣θ+∑i=1n∑j=13lnRjxi∣θ.The EM algorithm is an iterative algorithm, and every iteration of it consists of two consecutive steps, namely, the E step and the M step. In the E step, the expectation of the log-likelihood with respect to the estimate of the conditional probabilities of the latent variables has been constructed. Then, in the M step, the constructed expectation of the E step is maximized to compute the estimate of the parameters in the current iteration.


## 2.2. The E Step

Suppose that the estimate ofθ at iteration t be denoted by θt, then the conditional distribution of Zi is
(10)pi1,t=PZi=1∣Xi=xi,θt=f1xi∣θtR2xi∣θtR3xi∣θtfxi∣θt,(11)pi2,t=PZi=2∣Xi=xi,θt=f2xi∣θtR1xi∣θtR3xi∣θtfxi∣θt,and pi3,t=1−pi1,t−pi2,t where
(12)fxi∣θt=f1xi∣θtR2xi∣θtR3xi∣θt+f2xi∣θtR1xi∣θtR3xi∣θt+f3xi∣θtR1xi∣θtR2xi∣θt.The probabilitiespi1,t, pi2,t, and pi3,t are called the membership probabilities.Now we define the expectation of the log-likelihood function (8) with respect to the conditional distribution ofZi.
(13)Qθ∣θt=EZ∣X,θtlθ;x,Z=∑i=1np1i,tlnα1λ1α1xiα1−1−λ1xiα1−λ2xiα2−λ3xi+∑i=1np2i,tlnα2λ2α2xiα2−1−λ1xiα1−λ2xiα2−λ3xi+∑i=1np3i,tlnλ3−λ1xiα1−λ2xiα2−λ3xi=∑i=1np1i,tlnα1λ1α1xiα1−1−∑i=1nλ1xiα1+∑i=1np2i,tlnα2λ2α2xiα2−1−∑i=1nλ2xiα2+∑i=1np3i,tlnλ3−∑i=1nλ3xi.So theQθ∣θt can be written as sum of three distinct expressions Q1θ∣θt, Q2θ∣θt, and Q3θ∣θt, where
(14)Q1θ∣θt∑i=1np1i,tlnα1λ1α1xiα1−1−∑i=1nλ1xiα1,(15)Q2θ∣θt=∑i=1np2i,tlnα2λ2α2xiα2−1−∑i=1nλ2xiα2,(16)Q3θ∣θt=∑i=1np3i,tlnλ3−∑i=1nλ3xi.These statements will be applied in the M step, to compute the estimates of the parameters.


## 2.3. The M Step

The estimate of the parameters at iterationt+1 can be obtained by
(17)α^1,t+1,λ^1,t+1=argmaxα,λQ1θ∣θt,(18)α^2,t+1,λ^2,t+1=argmaxα,λQ2θ∣θt,(19)λ^3,t+1=argmaxλQ3θ∣θt.We should optimizeQ1θ∣θt and Q2θ∣θt numerically since they have not closed form for their critical points. But, the point which maximizes Q3θ∣θt has a closed form, and by solving the equation ∂/∂λ3Q3=0, we have
(20)λ^3,t+1=∑i=1np3i,t∑i=1nxi.The iterative process can be concluded ifQθt+1∣θt+1<Qθt∣θt+ε for some small predefined ε.


## 3. Simulation Study

To provide a random instance of the competing risk model with hazard rate (2), we simulate one random instance of Weibull with parametersα1,λ1, namely, X11, one random instance of Weibull with parameters α2,λ2, namely, X12, and one random instance of the exponential distribution with parameter λ3, namely, X13. Then, the random variable X1=minX11,X12,X13 follows from the desired competing risk model.In every run of the simulation study, we driver=500 replicates of samples of sizes n=50 and 100. Then, for each sample, the parameters have been estimated applying the MLE, LSE, or EM algorithm (see Supplementary Materials for all R codes used for simulation study). Every cell of Table 1 shows the results of one run. The results contain the bias (B), the absolute bias (AB), and the mean squared error (MSE) which, for example for α1, have been computed by the following relations.
(21)Bα1=1r∑i=1rα^1i−α1,(22)ABα1=1r∑i=1rα^1i−α1,(23)MSEα1=1r∑i=1rα∧1i−α12,where α^1i is the estimate of α1 based on the ith replication. Some important observations of the simulation results have been pointed out in the following.
(i)
As sample size increases, the AB and MSE decrease in all approaches(ii)
The EM estimator outperforms the LSE and MLE in terms of AB and MSE in all cases(iii)
The LSE outperforms the MLE in terms of AB and MSE in all casesTable 1
Every cell consists of the bias, the absolute bias, and the mean squared error for five parametersα1, λ1, α2, λ2, and λ3 from top to bottom, respectively.
nThe EM algorithmThe LSEThe MLEBABMSEBABMSEBABMSEα1=0.8λ1=0.01α2=1.3λ2=0.02λ3=0.04500.1147880.2094610.0733270.0649840.2664470.2954910.1507920.2782280.475591-0.0014100.0033860.0000170.0159990.0206590.0007070.0269290.0305710.0013990.3548830.4143570.3256450.4504760.5960591.2160751.7102621.82958730.454580.0016440.0042110.0000300.00083240.0152190.0003260.0037050.0139340.0003150.0062230.0076800.000101-0.0185700.0253290.000851-0.0255770.0318220.0011721000.0904200.1610600.0417280.0756340.2086680.0904260.0743980.2217430.092066-0.0018570.0028830.0000110.0141070.0185730.0005580.0226710.0270870.0011320.2698570.3000740.1511750.3513610.5519100.7877561.2668891.40722138.772690.0007110.0028020.0000120.0078000.0130410.0002640.0001710.0124930.0002430.0061170.0068010.000072-0.0173430.0240080.000788-0.0183730.0271220.000936α1=1.1λ1=0.1α2=0.9λ2=0.2λ3=0.3500.0872470.1893550.0671840.3880060.4767040.6566961.6860551.736100193.38740.0109440.0261640.0012470.1311340.1848520.0537990.0577700.1247140.0253030.0491170.1658280.047906-0.0117560.3061560.2738520.0192760.2395810.2410920.0118330.0334450.0019130.0825360.2155830.0702230.1614460.2508800.0846620.0109690.0333250.001867-0.1498150.2122980.056467-0.1650210.2450530.0682121000.0573420.1280940.0273250.3323000.4242070.5837800.9198100.9777018.0220460.0050990.0180840.0005250.1044500.1601000.0402680.0302510.1066440.0175810.0354470.1073490.018491-0.0111640.2277230.1324040.0084670.1799670.3118760.0042820.0230960.0008830.0878480.1958580.0578350.1702130.2487560.0824940.0035640.0255120.001046-0.1518750.2067080.053532-0.1588070.2362920.064980Table 2
Failure time of electrical appliances in terms of 1000 s cycles.
34596169801231421652103814795565748399179699911064108810911270127513551397147715781649170218931932216122922326233726282785281128862993312237153790385739124100410641164315451045845299558360659701


## 4. Applications

Lawless [20] analyzed failure times of some electrical appliances. The scaled TTT transform plot drawn in Figure 1 shows a bathtub shape for the hazard rate function. This gives us some nonparametric information indicating that the data come from a BT hazard rate model. So we tried to fit some distributions accommodating bathtub-shaped hazard rate to this data set. We use the MLE and the EM algorithm to fit the five parameters competing risk model (2). Also, the reduced model with the reliability function
(24)Rx=exp−λxα−λx1/α,has been fitted by computing the MLE of the parameters. In Figure 2, the cumulative distribution function (CDF) for fitted models has been drawn. There is a significant distance between the fitted model (14) and the empirical CDF which shows that the reduced model may be improper in some examples.Figure 1
The scaled TTT transform plot for data sets of Table2.Figure 2
The empirical CDF and some fitted models to data set of Table2.Moreover, some results of fit have been abstracted in Table3. Based on the Kolmogorov-Smirnov (K-S) statistics and p value, the competing risk model (2) which has been fitted by the EM algorithm gives the best description of the data. The MLE has also provided good results, but it is worthy to denote that we applied the estimates of the EM algorithm as the initial values in the likelihood maximization process. All of the fitted models confirm a BT hazard rate model which was firstly recognized by the TTT transform plot. So it may be interesting to investigate the point which maximizes the mean residual life and/or the median residual life functions. These points are referred to burn-in points and show the time at which the component is in its most reliable condition. The left side of Figure 3 draws the mean residual life function along with the median residual life function related to the best fitted model. Also, the burn-in points related to both functions have been determined in the figure. The right side of Figure 2 draws the hazard rate function of this model and shows a BT hazard rate model.Table 3
The results of fitting some suitable models to the data set.
The modelThe methodEstimationsK-S statisticsK-Sp valueAICModel (2)EMα^1=0.62125,λ^1=0.000110.052460.9936-α^2=1.98826,λ^2=0.00021λ^3=0.00012Model (2)MLEα^1=0.62445,λ^1=0.000110.053590.99171049.014α^2=1.99617,λ^2=0.00021λ^3=0.00012Model (20)MLEα^=0.00023,λ^=0.628910.099870.55391044.176Figure 3
(a) The mean residual life and the median residual life of the fitted model. (b) The hazard rate function of this fitted model.


## 5. Conclusion

The competing risk model of the baseline distribution Weibull plays a vital role in describing nonmonotone hazard rate models. One drawback of this model is that it has large number of the parameters which causes the estimation problem harder. Some authors suggested reduced versions to overcome this problem. But, there are many examples showing that the reduced model may not be proper. So, we implemented the EM algorithm for estimating the parameters. The simulation results confirm that this algorithm is better than MLE and LSE. As future works, such EM algorithm may be constructed for similar competing risk models or mixture models, for example, the gamma competing risk model with the following reliability function may be a good candidate:(25)Rx=1Γα1Γα2Γα1,λ1xΓα2,λ2x,x≥0,in which Γα,t=∫t∞yα−1e−ydy is the upper incomplete gamma function and α1>0, λ1>0, α2>0, and λ2>0.



---

*Source: 1179856-2021-10-21.xml*