Negative Binomial Kumaraswamy-G Cure Rate Regression Model

In survival analysis, the presence of elements not susceptible to the event of interest is very common. These elements lead to what is called a fraction cure, cure rate, or even long-term survivors. In this paper, we propose a unified approach using the negative binomial distribution for modeling cure rates under the Kumaraswamy family of distributions. The estimation is made by maximum likelihood. We checked the maximum likelihood asymptotic properties through some simulation setups. Furthermore, we propose an estimation strategy based on the Negative Binomial Kumaraswamy-G generalized linear model. Finally, we illustrate the distributions proposed using a real data set related to health risk.


Introduction
In survival analysis, the study is based on data relating to the time until the occurrence of a particular event of interest, also known as time to failure.This data can come from the time until there is a failure in an electronic component; time until a particular disease occurs in a patient; time for a particular drug to have the desired effect, among others.The behavior of such data can be verified empirically; this approach is said to not be parametric.If the data follows a probability distribution, then this approach is called parametric; this is the most used form in this work.
The survival and hazard functions, the objects of greatest interest in survival analysis, allow the study of such behavior.The survival function is the probability of an individual or component surviving after a preset time and the hazard function is the instantaneous failure rate, which graphically can take various forms, such as constant, increasing, decreasing, unimodal or bathtub shaped.However, when the behavior of the hazard function is not monotonous, the most commonly known distributions, such as exponential and Weibull, cannot accommodate this kind of behavior.This is because a disadvantage of these models is that they are very limited due to the small number of parameters and therefore the conclusions drawn from the models cannot be sufficiently robust to accommodate deviations from the data.There are some distributions that accommodate the non-monotonic hazard function, but they are usually very complicated and with many parameters.
We can model real survival data using almost any continuous distribution and with positive values; the simplest and most common models, such as exponential and Weibull distribution, may not be appropriate.Therefore, to find a distribution that accommodates non-monotone hazard functions is a known issue in survival analysis.Therefore, it is desirable to consider other approaches to achieve greater flexibility, and this is what has motivated studies to find distributions that accommodate these types of function.Kumaraswamy (1980) proposed the Kumaraswamy distribution, which was widely used in hydrology and, based on this, Cordeiro and de Castro (2011) proposed a new family of generalized distributions, called Kumaraswamy generalized (Kum-G).It is flexible and contains distributions with unimodal and bathtub-shaped hazard functions, as shown by De Pascoa et al. (2011), and has, as special cases, any distribution that is normal, exponential, Weibull, Gamma, Gumbel and inverse Gaussian.The domain of the distribution is the range in which the particular cases are set.Other examples of generalized distributions are the Generalized exponential distribution (Gupta and Kundu 1999) and the Stoppa (or Generalized Pareto) distribution (Stoppa 1990) and (Calderín-Ojeda and Kwok 2016).
In a population, there may be individuals who have not experienced the event of interest until the end of the study; this is called censorship.When there are a large number of censored individuals, we have an indication that in this population there are individuals who are not subjected to the event; they are considered immune, cured or not susceptible to the event of interest.
From the traditional models of survival, it is not possible to estimate the cured fraction of the population, or the percentage of individuals who are considered cured.Thus, statistical models are needed to incorporate such fractions and these are termed long-term or cure rate models.Because of this capability, different fit methods have been proposed in several areas such as biomedical studies, financial, criminology, demography and industrial reliability, among others.For example, in biomedical data, an event of interest may be the death of the patient due to tumor recurrence, but there may be patients who are cured and do not die due to cancer.When the financial data is studied, an event of interest may be the customer's closing of a bank account by default, but there may be customers who will never close their account.In criminology data, the event of interest may be a repeat offence and there may be people who do not repeat an offence.In industrial reliability, long duration models are used to verify the proportion of components that are not tested at zero time and are exposed to various voltage regimes or uses.In market research areas, individuals who will never buy a particular product are considered immune.See, for example, (Anscombe 1961;Farewell 1977;Goldman 1984;Broadhurst and Maller 1991;Meeker and Escobar 1998).
Many authors have contributed to the theory of long-term models.Boag (1949) was the pioneer; the maximum likelihood method was used to estimate the proportion of survivors in a population of 121 women with breast cancer, in an experiment that lasted 14 years.Based on Boag's idea, Berkson and Gage (1952) proposed a mixture model in order to estimate the proportion of cured patients in a population subjected to a treatment of stomach cancer.More complex long-term models, such as Yakovlev and Tsodikov (1996), Chen et al. (1999) among others, have emerged in order to better explain the biological effects involved.More recently, Rodrigues et al. (2009) proposed a unified theory of long duration, considering different competitive causes.In this context, most long-term models make use of this proposal, among which are (Sy and Taylor 2000;Castro et al. 2009;Cancho et al. 2011;Gu et al. 2011), besides (Ibrahim et al. 2005;Cooner et al. 2007;Ortega et al. 2008Ortega et al. , 2009;;Cancho et al. 2009).
A very important point in survival analysis is the study of covariates, because many factors can influence the survival time of an individual.Therefore, incorporate covariates enable us to have a much more complete model, full of valuable information.For example, if we are interested in studying the life time of patients with a particular disease who are receiving a certain treatment, other factors may influence the patient's healing, so we can find new ways to treat the disease from covariates.One real situation is the study of patients that were observed for recurrence after the removal of a malignant melanoma; it is desired to know if the nodule category or the age of the patient may influence the recurrence of melanoma.We will analyze this clinical study latter.
This paper presents the unified long-term model using, as a distribution of the number of competing causes, the negative binomial distribution, as studied in Cordeiro et al. (2015), where the authors use the Birnbaum-Saunders distribution of times.However, our contribution is proposing the use of a different distribution of times, i.e., the family of Kum-G distributions, which were studied only in the usual models of survival analysis, as in De Pascoa et al. (2011), De Santana et al. (2012) and Bourguignon et al. (2013).In this new model, we propose the incorporation of covariates influencing the survival time.In addition, we performed a simulation study to see how this model would behave with different sample sizes, as well as an application to a data set to demonstrate the applicability of this model.
The paper is organized as follows: in Section 2, we have the methodology, in which we present the family of Kumaraswamy generalized distributions, the unified cure rate model and a distribution used in this model, i.e., the negative binomial distribution; then, we propose a unified model Kumaraswamy-G cure rate as well as a regression approach, using the distribution Kumaraswamy exponential and its inferential methods.Section 2.7 presents some simulation studies.Application to a real data set is presented in Section 3. Finally, in Section 4, we conclude the paper with some final remarks.

Kumaraswamy Family of Distributions
The time until the occurrence of some event of interest can be generally accommodated by a probability distribution.In the literature, various distributions have been used to describe survival times but most commonly used distributions do not have the flexibility to model non-monotone hazard functions, such as unimodal and bathtub-shaped hazard functions, which are very common in biological studies.Thus, in this section, we will study the Kumaraswamy generalized distribution because it is a flexible but simple distribution.
The Kumaraswamy generalized distribution (Kum-G) presented by Cordeiro and de Castro (2011) has the flexibility to accommodate different shapes for the hazard function, which can be used in a variety of problems for modeling survival data.It is a generalization of the Kumaraswamy distribution with the addition of a distribution function G(t) of a family of continuous distributions.
Definition.Let G(t) be a cumulative distribution function (cdf) of any continuous random variable.The cdf of the Kum-G distribution is given by where λ > 0 and ϕ > 0. Let g(t) = dG(t) dt be the probability density function (pdf) of the distribution of G(t), then the pdf of Kum-G is Thus, we obtain the survival and hazard functions, given respectively by In the literature, there are different generalized distributions, one of which is the beta distribution.The pdf of beta generalizations uses the beta function, which is difficult to handle.On the other hand, the Kum-G distribution is a generalization that shows no complicated function in its pdf, and it is more advantageous than many generalizations.
As the Kum-G distribution depends on a G(t) distribution function, for each continuous distribution, we have a case of Kum-G with the number of parameters of G(t) over the two parameters λ and ϕ.For example, if we take the cumulative distribution function of an exponential distribution as G(t), then in this case we have the Kumaraswamy exponential distribution.In the literature, many cases of this distribution were studied, some of which are Kumaraswamy normal (Correa et al. 2012) (Shahbaz et al. 2012) and Kumaraswamy Rayleigh inverse (Hussian and A Amin 2014).

The Unified Cure Rate Model
The unified model of the cured fraction of Rodrigues et al. ( 2009) is a statistical model capable of estimating the proportion of a cured population, that is, in data sets in which many individuals do not experience the event of interest, even if observed over a long period of time, part of the population is cured or immune to the event of interest; we can estimate the cured fraction.Several authors have worked with this modeling, for example, (Rodrigues et al. 2009, Peng andXu 2012;Balakrishnan and Pal 2012, 2013a, 2013b, 2015), and others.
In general, the basic idea of the unified model of the cured fraction is based on the notion of occurrence of the event of interest in a process in two stages: Initiation stage.Let N be a random variable representing the number of causes or competitive risks of occurrence of an event of interest.The cause of the occurrence of the event is unknown, and the variable N is not observed, with probability distribution p n and its tail given respectively by .., n, continuous random variables (non-negative), independent of a cumulative distribution function F(z) = 1 − S(z) and independent of N, represent the time of occurrence of an event of interest because of the k-th cause.
In order to include individuals who are not susceptible to the event of interest, its time of occurrence is defined as where P(Z 0 = ∞) = 1, admitting the possibility that a proportion p 0 of the population lacks the occurrence of an event of interest, T is an observable or censored random variable, and Z j and N are latent variables.
Let {a n } be a sequence of real numbers.A(s) is defined as a function of the sequence {a n } as follows The survival function of the random variable T (population survival function) will be indicated by where A(•) corresponds to a genuine generating function of the sequence p n .That is, in the survival function of the random variable T, corresponding to a long-term model in two stages, the composition is a probability-generating function and survival function.The long-term survival function, in two stages S pop (t), is not a survival function.
Note that for the survival function, lim t→0 S(t) = 1 and lim t→∞ S(t) = 0.As for the improper survival function, lim t→0 S(t) = 1 and lim t→∞ S pop (t) = P(N = 0) = p 0 .Thus, p 0 is the proportion of non-event occurrences in a population of interest, that is, the cured fraction.
The population survival function has the following properties: The density and hazard functions associated with long-term survival function are given respectively by and Any discrete distribution can be used to model N, such as Bernoulli, binomial, Poisson, negative binomial and Geometric.What follows is the negative binomial distribution, which will be used because it is a very flexible distribution with various special cases, including those resulting in the standard model mix.
The probability generating function is given by Thus, the long-term survival function for the negative binomial model is given by where F(t) is the cumulative distribution function of the random variable T and the cured fraction of the population is The density function of the model (1) is where f (t) = −S (t).Furthermore, the corresponding hazard function is given by We observed some particular cases in this model: from the Equation (1), when η → 0, we obtain the density function of the Poisson distribution, resulting in the promotion time model; if η = −1, we fall into the Bernoulli distribution, where we have the model of the standard mixture; if η = 1, we have the geometric distribution; when η = 1/m (m integer), we have a binomial distribution (m, θ/m), where 0 ≤ θ/m ≤ 1.We also observed, from expressions of expectation and the variance of the model, that the variance of the number of competing causes is very flexible.If −1/θ < η < 0, there is an underdispersion relative to the Poisson distribution; if η > 0, there is an overdispersion.
Table 1 presents the long-term survival function, improper density and cure rate corresponding to negative binomials and their particular cases.

Distribution
S pop (t)

Negative Binomial Kumaraswamy-G Cure Rate Model
Considering the negative binomial distribution for the number of competing causes and the time following the Kumaraswamy-G distribution, we obtain a family of long-term distributions, wherein the population survival function of the model is given by with the cured fraction of the population given by So, by replacing the function G(t) by the cumulative distribution function of some distribution, one obtains a negative binomial Kumaraswamy-G model of long-term survival.
The population density function is and the population hazard function is given by Table 2 shows the particular cases of this model.It is noteworthy that for every G(t), we will have different distributions.
Negative Binomial Kumaraswamy Exponential Cure Rate Model Considering G(t), following an Exponential(α) distribution and substituting in (2), we have the NegBinKumExp(α, λ, ϕ, η, θ), i.e., a family of cure rate models where their population survival function is given by The population density and hazard function of this model are, respectively, and

Negative Binomial Kumaraswamy-G Regression Cure Rate Model
The use of covariate information is essential when analyzing survival data.Here, we discuss an approach to including covariate information for the proposed models.
Suppose that x = (1, x 1 , . . . ,x k ) is a vector of covariates from a data set and β = (β 0 , β 1 , . . . ,β k ) is a vector of regression coefficients.We are going to set θ(x) = exp β x to link the cure rate to the covariates.The only two parameters that link the cure rate to the covariates are θ and η.Since θ has a positive domain, we can use it to simply model the covariates through the exponential function.
This way, the Negative Binomial Kumaraswamy-G generalized linear model is given by for t > 0. The cure rate p is given by This way, the cured fraction does not depend on the parameters of the Kumaraswamy family or the baseline distribution, but on the parameters η and θ.They are estimated differently for each baseline distribution and then they are incorporated into the cure rate.
Other particular cure rate models can be obtained.The Bernoulli Kumaraswamy-G generalized linear model and its respective cure rate is given by The Poisson Kumaraswamy-G generalized linear model is given by The Geometric Kumaraswamy-G generalized linear model is given by In Section 2.6, we discuss the estimation procedures of the NegBinKum-G cure rate generalized linear model.An application of these models is presented in Section 3.

Inference
Here, we present a procedure to obtain maximum likelihood estimates for the Negative Binomial Kumaraswamy Exponential generalized linear model.We consider data with right-censored information.Let D = (t, δ, x), where t = (t 1 , . . . ,t n ) are the observed failure times, δ = (δ 1 , . . . ,δ n ) are the right-censored times and x is the covariates information.The δ i is equal to 1 if a failure is observed and 0 otherwise.Suppose that the sample data is independently and identically distributed and comes from a distribution with density and survival functions specified by f (•, ν) and S (•, ν), respectively, where ν = (α, λ, ϕ, η, β) denotes a vector of 4 + (k + 1) parameters, with θ = exp(β x), as described in Section 2.5.By combining (4) and the expression (3), the log-likelihood function of ν for the NegBinKumExp distribution is The maximum likelihood estimates are the simultaneous solutions of ∂l (ν, D) ∂ν i = 0.
The estimates are obtained using the BFGS algorithm of maximization, which is an option for the optim function in R (R Core Team 2013).
If ν denotes the maximum likelihood estimator of ν, then it is well known that the distribution of ν − ν can be approximated by a (k + 5)-variate normal distribution with zero means and a covariance matrix I −1 ( ν), where I (ν) denotes the observed information matrix defined by for i and j in 1, . . ., k + 5.This approximation can be used to deduce confidence intervals and tests of hypotheses.For example, an approximate 100(1 − γ) percent confidence interval for , where I ii denotes the ith diagonal element of the inverse of I and z γ denotes the 100(1 − γ) percentile of a standard normal random variable.
Evidence of the existence of cured individuals is given in cases where the Kaplan-Meier curve reaches a plateau between zero and one.In some cases, this is clearer than others, as one can see in our examples.We can assume that some of the censored observations at the end of the study belong to the cured group.If everyone censored at the end were indeed cured, then the plateau reached by the Kaplan-Meier curve is a good estimate of the cured fraction.In general, a lower value of this plateau or a value close to it is an acceptable estimate.This data set collected in the period 1991-1998 is related to a clinical study in which patients were observed for recurrence after the removal of a malignant melanoma.Melanoma is a type of cancer that develops in melanocytes, responsible for skin pigmentation.It is a potentially serious malignant tumor that may arise in the skin, mucous membranes, eyes and the central nervous system, with a great risk of producing metastases and high mortality rates in the latter stages.In total, 417 cases were observed, of which 232 were censored (55.63 percent).The overall survival is 3.18 years.This data set has covariate information, which is used to illustrate the generalized linear model proposed in Section 2.5.The covariates taken represent the nodule category (n 1 = 82, n 2 = 87, n 3 = 137, n 4 = 111) and age (continuous covariate).The overall survival times for the categories are 3.60, 3.27, 3.07, 2.55 years.For more details on this data, see Ibrahim et al. (2001).
Tables 3-6 show the results for the Bernoulli, Poisson, Geometric and Negative Binomial Kumaraswamy Exponential models.The estimated cure rates p 1 , p 2 , p 3 and p 4 for groups 1, 2, 3 and 4, respectively, are calculated by ( 5).The age covariate is taken as their average, 48, for the necessary computations.The estimates of β 0 , β 1 and β 2 are in agreement in all models.For β 0 , the value is around −1.50, for β 1 , the value is around 0.50 and for β 2 , the value is around 0.01.
In Figure 1, the fitted survival curves for each nodule category and each proposed model are given.We can see that the one that best captures the Kaplan-Meier curve is the Negative Binomial Kumaraswamy Exponential distribution.This result is also sustained by the AIC values.The values obtained for the Bernoulli, Poisson, Geometric and Negative Binomial Kumaraswamy Exponential models are 1029.53, 1022.77, 1017.64 and 1016.27.The Negative Binomial Kumaraswamy Exponential achieves a better AIC value even with one extra parameter than the others.
Considering the Negative Binomial Kumaraswamy Exponential model and given the average age of 48 in this study, the estimated cure rate for nodule category 1 is around 0.65.For nodule category 2, it is around 0.54.For nodule category 3, it is around 0.41.For nodule category 4, it is around 0.28.

Table 1 .
Survival function S pop (t), density function f pop (t), and cured fraction of different distributions of latent causes.

Table 2 .
S pop (t), f pop (t) and the cured fraction for different distributions of N. η

Table 3 .
MLEs of the Bernoulli Kumaraswamy Exponential model for the melanoma data set.

Table 4 .
MLEs of the Poisson Kumaraswamy Exponential model for the melanoma data set.

Table 5 .
MLEs of the Geometric Kumaraswamy Exponential model for the melanoma data set.

Table 6 .
MLEs of the Negative Binomial Kumaraswamy Exponential model for the melanoma data set.