|Home | About | Journals | Submit | Contact Us | Français|
In this article, we propose a new generalization of the Weibull distribution, which incorporates the exponentiated Weibull distribution introduced by Mudholkar and Srivastava  as a special case. We refer to the new family of distributions as the beta-Weibull distribution. We investigate the potential usefulness of the beta-Weibull distribution for modeling censored survival data from biomedical studies. Several other generalizations of the standard two-parameter Weibull distribution are compared with regards to maximum likelihood inference of the cumulative incidence function, under the setting of competing risks. These Weibull-based parametric models are fit to a breast cancer dataset from the National Surgical Adjuvant Breast and Bowel Project (NSABP). In terms of statistical significance of the treatment effect and model adequacy, all generalized models lead to similar conclusions, suggesting that the beta-Weibull family is a reasonable candidate for modeling survival data.
A suitable parametric model is often of interest in the analysis of survival data, as it provides insight into characteristics of failure times and hazard functions that may not be available with non-parametric methods. The Weibull distribution  is one of the most commonly used families for modeling such data. However, only monotonically increasing and decreasing hazard functions can be generated from the classic two-parameter Weibull distribution; as such this two-parameter model is inadequate when the true hazard shape is of unimodal or bathtub nature.
The widely used Kaplan-Meier product-limit estimator is a flexible method to model survival data but is noted to often be inefficient . Other semi-parametric methods, such as proportional hazards modelling, require assumptions that may not be plausible in many situations . Meanwhile, various parametric techniques have been developed to incorporate the wide variety of patterns in survival data. Some of the proposed parametric models have incorporated a shape parameter into the classic Weibull distribution to account for additional possible hazard shapes. One of the first models of this type was proposed by Kalbfleisch , but this model may be impractical in the presence of censored data as it often requires the evaluation of an incomplete gamma integral or beta ratio.
Mudholkar and Srivastava  introduced an exponentiated version of the Weibull model that included an additional shape parameter. The distribution has a closed form of probability density, survival, and hazard functions, that are flexible and able to generate a wide variety of frequently observed hazard shapes, including unimodal and bathtub. It can be shown that this extension of the Weibull can be achieved through a simple application of the probability integral transform, using the densities of Beta and two-parameter Weibull distributions. Mudholkar, Srivastava, and Kollia  proposed another generalization of the Weibull model, which is able to generate similar types of hazard shapes as the exponentiated model; however irregularities may arise as the support of the distribution becomes dependent on the parameter space. A four-parameter generalization of the Weibull distribution was introduced by Jeong  based on stable distributions proposed by Hougaard , induced under the semiparametric frailty model. This model can also be viewed as a generalization of the model of Mudholkar et al. .
A decade has passed since the exponentiated Weibull was proposed, however a limited amount of literature is available on its application to survival data. This article introduces the family of beta-Weibull distributions encompassing regular two-parameter Weibull and the exponentiated Weibull distributions. It focuses on the potential use of the beta-Weibull distribution as a model for survival data in the presence of competing risks, in comparison with other generalizations of the Weibull model.
In the analysis of breast cancer data, it is often of interest to investigate only a subset of any first events. For example, investigators may be interested in investigating local or regional recurrence of the original breast cancer, considering other events such as distant recurrence, new primary cancers other than breast, or deaths prior to any disease . In this case, a competing risks situation arises, i.e. local or regional events against other events. It is well-known that the cumulative incidence function estimates the proportion of local or regional recurrences correctly as a function of time in the presence of other events [5, 10]. In this paper, the cumulative incidence function will be parameterized and compared by using various existing Weibull distributions along with the proposed beta-Weibull distribution. These parametric approaches will be illustrated in a data set from one of the phase III breast cancer clinical trials performed by National Surgical Adjuvant Breast and Bowel Project (NSABP).
The article is organized as follows. In Section 2, the beta-Weibull family of distributions is introduced along with an overview of existing Weibull-based families. In Section 3, maximum likelihood estimation and corresponding inferential procedures are described for the generalized Weibull models in the setting of competing risks. In Section 4, these models are fitted to a previously analyzed breast cancer data set, and performance of the Weibull-based approaches are evaluated by using conventional model assessment methods. We conclude with a few remarks in Section 5.
The classic two-parameter Weibull distribution has the probability density function (pdf)
with its cumulative distribution function being
While this convenient distribution can be easily incorporated into practical tools for survival analysis such as proportional hazards and accelerated failure time models, only monotonically increasing or decreasing hazard shapes can be properly estimated. To overcome this limitation, several generalized versions of the classic two-parameter Weibull model have been proposed. In this section we describe two such models that require at most two extra parameters to be estimated.
A three-parameter generalization  of the Weibull distribution is defined by the survival function
and the corresponding hazard function
The regular case of this generalized Weibull distribution occurs when λ ≥ 0 with the density having support (0, ∞), generating monotonically decreasing or unimodal hazard functions. In particular, the distribution approaches the two-parameter Weibull as λ → ∞. When λ < 0, monotonically increasing or bathtub hazard shapes are generated, but the support of the density becomes parameter-dependent, within the range (0, β(−λ)α). Similar to other non-regular densities, in these situations the likelihood may become unbounded and thus maximum likelihood estimates may not exist, or alternately, although the maximum likelihood estimates exist, their asymptotic properties may not hold for the classical asymptotic theory .
A four-parameter Weibull model was proposed by Jeong . This generalization of the Weibull distribution incorporates an additional parameter τ to the three-parameter version proposed by Mudholkar  that was previously discussed. This generalization can be characterized by the survival function
with the corresponding hazard function
This distribution reduces to Mudholkar’s generalized Weibull distribution as τ approaches to zero. As one would expect, it covers a wider range of shapes than the three-parameter extension and is regular, regardless of the parameter values.
It has been demonstrated by Wahed  and Ferreira and Steel  that any parametric family of distribution can be incorporated into larger families through an application of the probability integral transform. Specifically, given two valid probability density functions (pdf) g1(·) and g2(·), with the latter having support on the unit interval, a third valid density function f(·) may be obtained by applying the equation
where G1(·) is the cumulative distribution function (cdf) corresponding to g1(·). Using this technique, we obtain yet another generalization of the two-parameter Weibull distribution as described below.
To obtain the appropriate extension of the Weibull distribution, we considered g1(·) as the two-parameter Weibull pdf given by (1) and g2(·) as a two-parameter beta distribution with pdf
where 0 < x < ∞, 0 < α1, α2, γ, β < ∞ and B(α1, α2) = Γ(α1)Γ(α2)/Γ(α1+α2). We will refer to this distribution as the beta-Weibull distribution throughout the article. The density given by (6) is a proper density and when α1 = α2 = 1, it reduces to the two-parameter Weibull density with scale parameter β and the shape parameter γ. This model is generated entirely based on two commonly used probability distributions, beta and Weibull, and therefore can be easily implemented using conventional statistical software packages. This process of using the beta distribution to generate a new class of distribution has been previously considered by Jones . The survival function of the beta-Weibull distribution cannot be expressed in a closed form; however it may be estimated by evaluating the regularized beta function at the cumulative distribution function of the Weibull distribution. Specifically, the survival function of the beta-Weibull distribution (6) is given by
is the incomplete beta function. The corresponding hazard function is obtained as
The survival and hazard functions are plotted in Figures 1 and and2,2, respectively. The hazard plots show that this model can generate unimodal, v-shaped and bathtub shaped hazards. The parameters, α1 and α2, are related to the beta-component of the distribution; as a result, the mean and variance of the distribution increases as α1 increases, as for the beta distribution. Similarly, the survival probability is higher for larger α1. The shape of the hazard function is influenced by the parameter α1 and γ, while β remains as a scale parameter. The second parameter from the beta distribution, α2, also primarily acts as a scale parameter as it is absorbed into the Weibull function through the exponent term e−(x/β)γ.
The skewing mechanism based on the beta distribution has also been considered in other settings. For example, Ferreira and Steel  used a restricted version of the beta distribution to obtain skewed versions of multivariate distributions where they imposed the restriction α2 = 1/α1 on the beta-parameters. We impose the same constraint on (6) to obtain what we will refer to as a restricted beta-Weibull distribution.
Notice that when α2 = 1, the beta-Weibull distribution defined by (6) reduces to
Setting α = α1 for notational convenience, the derivative of the hazard function with respect to x can be expressed as
Notice that m(x; α, β, γ) is non-negative regardless of any values of x or parameters.
As discussed in Mudholkar and Srivastava , for α < 1, the hazard shapes are unimodal if α < 1/γ and monotonically decreasing otherwise; likewise when α > 1 bathtub shapes are generated for α > 1/γ, with monotonically increasing shapes otherwise. The hazard function is monotonically increasing when α ≥ 1, γ > 1, while α ≤ 1, γ < 1 implies a monotonically decreasing function; when both parameters are equal to 1, the distribution reduces to the exponential distribution (constant hazard). However, when one of the parameters is greater than 1 and the other is less than 1, more shapes can be generated. If α > 1 and γ < 1, a unimodal hazard shape is generated if α > 1/γ and is monotonically decreasing otherwise; if α < 1 and γ > 1, bathtub shapes are generated when α < 1/γ < 1 and are monotonically increasing otherwise. Behavior of the hazard function with respect to the parameters is depicted in Figure 3 for values of α less than 1 (left) and greater than 1 (right). Similarly, the hazard shapes for different values of the shape parameter γ for α < 1 (left) and α > 1 (right) are plotted in Figure 4. Different types of hazard shapes may be obtained by altering this parameter while keeping the other two parameters fixed.
While both exponentiated Weibull distribution and three-parameter generalized Weibull (Section 2.1) reduce to the two-parameter Weibull distribution when λ approaches ∞ and λ approaches 1, there is no clear link between the parameters that comprise these distribution functions. This is apparent from the behavior of respective hazard functions; λ < 0 implies either a monotonically increasing or bathtub shape for the generalized Weibull, while α > 1 may result in a monotone hazard function that is increasing or decreasing, or a unimodal function. Similarly, λ > 0 produces a monotonically decreasing or unimodal shape for the generalized Weibull, whereas if 0 < α < 1 for the exponentiated Weibull, a bathtub shape is possible, along with monotonically decreasing or increasing shapes. The exponentiated Weibull is a generalization of the Burr Type X distribution S(x) = [1 − e−(x/β)2]α, while the family of the three-parameter generalized Weibull includes that of the Burr Type XII family of distributions  as a special case.
When multiple cause-specific events, or competing risks, are present, the cumulative incidence function correctly assesses the cumulative probability of a particular cause-specific event while taking into account other types of events . This prevents the cumulative probability of a sub-distribution of cause specific events of interest from being overestimated, which would occur if failures due to other causes are considered as censored . For data consisting of only two types of cause-specific events, the cumulative incidence function for cause-specific events of type 1, under the assumption of cumulative hazards, is defined as
where S2(·, ϕ2) and f1(·, ϕ1) are survival function for events of type 2 and probability density function for events of type 1, respectively. In this paper, we will focus on parameterizing S2(·, ϕ2) and f1(·, ϕ1) via Weibull-based models to estimate the cumulative incidence function.
The set of parameter vector ϕ = (ϕ1, ϕ2) of a given Weibull-based model for a given data set can be estimated through the maximum likelihood method. For the cause-specific event of type k(k = 1, 2), the observed data for the ith individual is denoted by (xi, δki), where xi represents the event time (failure/censoring) and δki is an indicator variable, taking value 1 when the kth cause-specific event occurs. The likelihood function for the kth cause-specific event can be expressed as
For all the Weibull-based models described in the previous section, the likelihood function (16) cannot be maximized analytically. Thus, an iterative procedure such as the Newton-Raphson method needs to be used to obtain maximum likelihood estimates for the parameter vector ϕk for kth cause-specific event. We have used the S-Plus function nlminb to minimize −lk(ϕk) to obtain the maximum likelihood estimates of ϕk. The initial values of the parameters for more complicated models were chosen from the results of fitting a classical two-parameter Weibull model to the data, while the initial values for both skewing beta parameters were set to 1 (equivalent to the classical Weibull). When the Newton-Raphson procedure is used to obtain the MLEs, in order to facilitate convergence it is suggested to set the options of the minimization procedure such that the iterative step length is of similar magnitude between each parameter . Once the MLE of ϕk is obtained, the variance of these estimates are calculated by inverting the observed information matrix.
The corresponding estimators for the cumulative incidence and survival functions are obtained by substituting the parameter estimates into the respective functions. For instance, the cumulative incidence function (15) is estimated by
The standard errors of the estimated cumulative incidence function and the survival probability is then obtained through an application of the multivariate delta method. For instance, an estimate of the approximate variance of the estimated cumulative incidence function F1(·) is given by (see  for derivation)
Equality of survival or cumulative incidence functions between independent treatment groups can then be tested using the Wald test. For testing the equality of two parametric functions θ(1)(t) and θ(2)(t) from two different treatments at a specific time point t, the appropriate test statistic is
where Z ~ N(0, 1). For instance, the appropriate Wald test statistic to compare two cumulative incidence functions at time t will have the form
We fitted the models to a data set from a breast cancer study (Protocol B-20) performed by the National Surgical Adjuvant Breast and Bowel Project Study (NSABP) . The study was designed to test if the addition of chemotherapy to tamoxifen would result in increased disease-free survival compared to tamoxifen alone. Here the disease-free survival events include breast cancer recurrences in local, regional, or distant sites, second primary cancers other than breast, or death prior to previously mentioned events, whichever occurs first. In this section, the cause-specific events of interest will be the local or regional recurrences, and the other competing events consist of distant recurrences, second primaries, and death without evidence of disease (see Figure 5). A total of 2,363 patients were randomized to tamoxifen, CMF (a chemotherapy containing the alkylating agent cyclephosphamide) or CMFT (CMF plus tamoxifen). The sample consists of 770 of the patients in the tamoxifen group and 766 in the CMFT group. We present our data only up to first 10 years as the risk set become sparse afterwards.
Maximum likelihood estimates for the parameters from different model fits are displayed in Table I for the cause-specific models, and Table II for the disease-free model. The Akaike Informational Criterion (AIC) is also provided. For the 2-parameter, 3-parameter, 4-parameter and exponentiated versions of the Weibull model, the standard errors were obtained from inverting the matrix of minus the second-derivatives of the log-likelihoods, evaluated at the MLEs. Bootstrapping, using 1000 replicates, was used to estimate the standard errors for the restricted and full beta-Weibull models, for which the second derivatives of the log-likelihood were not available in a closed form.
Within the tamoxifen group, the 4-parameter Weibull model had the best fit as assessed by the AIC followed by the restricted and unrestricted beta-Weibull model, although there were hardly any difference between the AIC’s for the latter two models. The two-parameter Weibull model had the largest AIC. By contrast, in the CMFT group there were no significant differences among the fits, regardless of the number of parameters in the model.
The estimated overall hazard function from each of the six model fits is displayed in Figure 6 for each type of the competing events and in Figure 7 for all disease-free events. The non-parametric hazard estimates are life table estimates of the annual hazard rate at the midpoint of each year. Model fits seem generally better for the other events, as local or regional recurrences were relatively infrequent and had more variable hazard shapes (see more on model adequacy in Section 4.2).
The generalized Weibull models, including the beta-Weibull versions, were able to capture more complex unimodal hazard shapes (local or regional recurrences) as well as the basic monotone shapes (other events). Within the tamoxifen group, according to the non-parametric fit there was an apparent peak at 2–3 years which leveled off after 10 years. This peak was inadequately captured by the two-parameter model, which only produced monotonic shapes; the Weibull extensions produced unimodal shapes that had captured this peak. The curves from the exponentiated model were closest to those from the two-parameter model, tending to produce a more conservative peak at 2–3 years than those from the other four extensions of the Weibull. The hazards from the restricted version of beta-Weibull tended to follow more closely those from the full beta-Weibull model. It is interesting to observe that, as the number of parameters increases, a model tends to pick up the peak earlier in this dataset.
There did not appear to be a similar single peak for the observed events in the CMFT group. As such, monotonically increasing hazard shapes from the two-parameter Weibull model provided an adequate fit to the observed non-parametric curves. The three and four-parameter models very closely followed the two-parameter fit.
Model adequacy was investigated using signed deviance residuals with censored observations  to assess goodness of fit (Tables III and andIV).IV). These are obtained by calculating the expected number of events ey during each year y as Nyhy, where Ny is the number of patients at risk at the beginning of each year, and hy is the hazard function integrated over the time interval. The goodness-of-fit between the observed number of events Sy and expected events ey for each of the models can be assessed through signed deviance residuals (Ry), the squared sum of which approximately follows a χ2 distribution, where
The models tended to provide a relatively poor fit to the data for the endpoint of local or regional recurrences, in large part due to their relative scarcity and the irregular hazard shapes; the models provided better fits for the other events and when all events were grouped together to estimate the overall survival curve.
The cumulative incidence function was estimated using all six Weibull models for each type of competing events, i.e. local or regional recurrences and other events.
The non-parametric cumulative incidence was estimated through Gray’s method , while the parametric cumulative incidence functions were obtained through the invariance property of the maximum likelihood estimates from all six Weibull models. To demonstrate how closely the parametric cumulative incidence curves followed the non-parametric estimates, we only presented the estimated cumulative incidence curves from the four-parameter Jeong model and the restricted beta-Weibull model in Figure 8 for each event type. We have chosen these two Weibull models as they had lower AIC values compared to other models. Figure 9 shows the corresponding disease-free survival curves. There were hardly any differences in cumulative incidence or survival estimates between the two models. The fits from the Weibull models generally followed the non-parametric cumulative incidence curves closely. There is a significant reduction in local or regional recurrences in the CMFT group, but no apparent difference between the two treatment groups in other types of events (Table V). A formal Wald test also concluded that there was a significant improvement in disease-free survival in the CMFT group, which was mainly due to reduction in local or regional recurrences. Both Jeong model and the restricted beta-Weibull model resulted in the same conclusions.
In this article, we have presented a new generalization of Weibull distribution, the beta-Weibull family, by applying the beta skewing mechanism on the classical two-parameter Weibull distribution. This generalization is easily obtained by transforming the two-parameter Weibull model through a beta distribution and hence can be implemented easily using standard statistical software packages.
The beta-Weibull, two-parameter, exponentiated and other existing generalized models were used to model competing risks and disease-free survival in data from a previously published breast cancer study. The two-parameter Weibull distribution produced simple monotone hazard shapes, as expected, that did not reflect pattern of the hazard shape in the observed data. On the other hand, the extensions of the Weibull model were able to more accurately capture the observed hazard pattern. In terms of statistical significance of the treatment effect and model adequacy, all extended models lead to similar conclusions, suggesting that beta-Weibull family could play a reasonable role as a candidate for modeling survival data, including competing risks. One limitation of the beta-Weibull family is that the survival and the hazard functions cannot be expressed in closed form and hence numerical integration techniques were needed to evaluate them.