EM Algorithm for Estimating the Parameters of Weibull Competing Risk Model

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.


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: 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: 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 bathtubshaped 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 Section 2, 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 Let x 1 , x 2 , ⋯, x n 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 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 The empirical reliability function is defined to bẽ in which the indicator function Iðt < x i Þ equals 1 when t < x i 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 Þ.
2.1. The EM Algorithm. Let X 1i , X 2i , and X 3i , 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 X i = min fX 1i , X 2i , X 3i g. This means that the lifetime event may be due to one of three competing causes. Let the latent random variable Z i with the support f1, 2, 3g such that With these notations, the likelihood function is where R j ðx | θÞ, f j ðx | θÞ, and λ j ðx | θÞ show the corresponding reliability function, the density function, and the hazard rate function of X ji , j = 1, 2, 3. Then, the log-likelihood function is 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.

The E
Step. Suppose that the estimate of θ at iteration t be denoted by θ t , then the conditional distribution of Z i is and p i3, Applied Bionics and Biomechanics The probabilities p i1,t , p i2,t , and p i3,t are called the membership probabilities. Now we define the expectation of the log-likelihood function (8) with respect to the conditional distribution of So the Qðθ | θ t Þ can be written as sum of three distinct These statements will be applied in the M step, to compute the estimates of the parameters.

The M
Step. The estimate of the parameters at iteration t + 1 can be obtained by b λ 3,t+1 = arg max We should optimize Q 1 ðθ | θ t Þ and Q 2 ðθ | θ t Þ numerically since they have not closed form for their critical points. But, the point which maximizes Q 3 ðθ | θ t Þ has a closed form, and by solving the equation ð∂/∂λ 3 The iterative process can be concluded if Qðθ t+1 | θ t+1 Þ < Qðθ t | θ t Þ + ε for some small predefined ε.

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, X 11 , one random instance of Weibull with parameters ðα 2 , λ 2 Þ, namely, X 12 , and one random instance of the exponential distribution with parameter λ 3 , namely, X 13 . Then, the random variable X 1 = min fX 11 , X 12 , X 13 g follows from the desired competing risk model.
In every run of the simulation study, we drive r = 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.
where b α 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.

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 3 Applied Bionics and Biomechanics competing risk model (2). Also, the reduced model with the reliability function 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. Moreover, some results of fit have been abstracted in Table 3. Based on the Kolmogorov-Smirnov (K-S) statistics Table 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.    Table 2.  Applied Bionics and Biomechanics 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.

Conclusion
The competing risk model of the baseline distribution Weibull plays a vital role in describing nonmonotone hazard rate      34  59  61  69  80  123  142  165  210  381  479  556  574  839  917  969  991  1064  1088  1091  1270  1275  1355  1397  1477  1578  1649  1702  1893  1932  2161  2292  2326  2337  2628  2785  2811  2886  2993  3122  3715  3790  3857  3912  4100  4106  4116  4315  4510  4584  5299  5583  6065  9701 5 Applied Bionics and Biomechanics 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: in which Γðα, tÞ = Ð ∞ t y α−1 e −y dy is the upper incomplete gamma function and α 1 > 0, λ 1 > 0, α 2 > 0, and λ 2 > 0.

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

Conflicts of Interest
There is no any conflict of interest.