In this paper a new long-term survival distribution is proposed. The so called long term Fréchet distribution allows us to fit data where a part of the population is not susceptible to the event of interest. This model may be used, for example, in clinical studies where a portion of the population can be cured during a treatment. It is shown an account of mathematical properties of the new distribution such as its moments and survival properties. As well is presented the maximum likelihood estimators (MLEs) for the parameters. A numerical simulation is carried out in order to verify the performance of the MLEs. Finally, an important application related to the leukemia free-survival times for transplant patients are discussed to illustrates our proposed distribution

Keywords: Fréchet distribution; Long-term survival distribution; Survival model

Extreme value models play an important role in statistic. The generalized extreme value (GEV) distribution [1] and its sub-models are widely used in application involving extreme events. These sub-models are the well known Weibull, Fréchet and Gumbel distributions. The Fréchet distribution can be seen as the inverse Weibull distribution which gives a probability density function (PDF) such as

$f\left(t,\lambda ,\alpha \right)=\frac{\alpha }{\lambda }{\left(\frac{t}{\lambda }\right)}^{-\left(\alpha +1\right)}{e}^{-{\left(\frac{t}{\lambda }\right)}^{-\alpha }}$ (1)

The survival function is given by

$S\left(t,\lambda ,\alpha \right)=1-e^{-{\left(\frac{t}{\lambda }\right)}^{-\alpha }}$ (2)

Although the GEV distribution is the most used generalization of the Fréchet model, other distributions has been proposed in the literature. De Gusmão [2] proposed a three parameter generalized inverse Weibull distribution in which includes the Fréchet distribution. Krishna et al. [3] proposed the Marshall-Olkin Fréchet distribution. Barreto-Souza et al. [4] discussed some results for beta Fréchet distribution. However, in survival studies habitually the researches may consider a portion of the population as cured during a given treatment, this type of distribution is called long-term (LT) survival models.

In this study, a long-term survival novel proposing a mixture model introduced by Berkson and Gage [5], hereafter we shall call it the long-term Fréchet distribution or simplistically the LF distribution. Some mathematical properties about the LF distribution were provided such as moments, survival properties and hazard function. The maximum likelihood estimators of the parameters and its asymptotic properties are discussed likewise. Similar studies were presented by Roman et al. [6] for the geometric exponential distribution and by Louzada and Ramos [7] for the weighted Lindley distribution. It was performed a numerical simulation towards to examine the performance of the MLEs. Finally, our proposed methodology is illustrated in a real data set related to the leukemia free-survival times (in years) for the 50 autologous transplant patients.

The paper is organized as follows. Section 2 presents the long term Fréchet distribution and its mathematical properties. Section 3 discusses the parameter estimation under the maximum likelihood approach. Section 4 presents a simulation study under different values of the parameters and different levels of censorship. The proposed methodology is also fully illustrated in a real data set. Lastly, Section 6 summarizes the founds in this study and its potential contribution.

Long-term survivors are an important feature to incorporate in the modeling process, since a portion of the population may no longer be eligible to the event of interest (according to Maller and Zhou, [8]; or Perdona and Louzada, [9]). Hence the population can be segregate as a not eligible to the event of interest with probability $p$ and as eligible (in risk) to the event of interest with probability$\left(1-p\right)$ . The long-term survivor is expressed as

$S\left(t;p,\theta \right)=p+\left(1-p\right)S_0\left(t;\theta \right)$ ,  (3)

where$p\in \left(0,1\right)$ and$S_0\left(t;\theta \right)$ is the survival function related to the eligible group. The obtained survival function (not conditional) is improper and its limit corresponds to the individual proportion cure. From the survival function one can easily derive the PDF (improper) given by

$f\left(t;p,\theta \right)=-\frac{\partial }{\partial t}S\left(t;p,\theta \right)=\left(1-\pi \right)f_0\left(t;\theta \right)$ ,  (4)

where$f_0\left(t;\theta \right)$ is the PDF related to the susceptible group.

Figure 1: It shows some cases about the PDF and the survival function shapes applied to LF distribution.
In Left panel: Probability density function of the LF distribution. Right panel: Survival function of the LF distribution.

Considering that$f_0\left(t;\theta \right)$ follows a Fréchet distribution, then the PDF of the Long Term Fréchet (LF) distribution is given by

$f\left(t;\lambda ,\alpha ,p\right)=\frac{\alpha \left(1-p\right)}{\lambda }{\left(\frac{t}{\lambda }\right)}^{-\left(\alpha +1\right)}{e}^{-{\left(\frac{t}{\lambda }\right)}^{-\alpha }}$ ,  (5)

where $\lambda >0$ , $\alpha >0$ and $p\in \left(0,1\right)$ . The cumulative distribution function is given by

$F\left(t;\lambda ,\alpha ,p\right)=\left(1-p\right)e^{-{\left(\frac{t}{\lambda }\right)}^{-\alpha }}$ .  (6)

In this case, the LF has the quantile function in closed-form and is given by

$t_u=\lambda \log {\left(\frac{\left(1-p\right)}{u}\right)}^{-\frac{1}{\alpha }}$    (7)

where$0\le u<1$ . The $r$  -th moments of T about the origin is

$E\left(T^r;\lambda ,\alpha ,p\right)=\left(1-p\right){\lambda }^{\text{r}}\Gamma \left(1-\frac{r}{\alpha }\right) , \alpha >r$ , (8)

For $r\in \mathbb{N}$ and $\Gamma \left(x\right)=\underset{0}{\overset{\infty }{{\displaystyle \int }}}{e}^{-y}{y}^{x-1}dy$ is called gamma function. Along with some algebraic manipulation the mean and variance of the LF distribution are given, respectively, by

$E\left(T;\lambda ,\alpha ,p\right)=\left(1-p\right)\lambda \Gamma \left(1-\frac{1}{\alpha }\right) , \alpha >1$

and

$V\left(T;\lambda ,\alpha ,p\right)=\left(1-p\right){\lambda }^2\left(\Gamma \left(1-\frac{2}{\alpha }\right)-\left(1-p\right)\Gamma {\left(1-\frac{1}{\alpha }\right)}^2\right), \alpha >2$ .

The survival and hazard functions of $LF\left(\lambda ,\alpha ,p\right)$ distribution is given by

$S\left(t;\lambda ,\alpha ,p\right)=p+\left(1-p\right)\left(1-e^{-{\left(\frac{t}{\lambda }\right)}^{-\alpha }}\right)$    (9)

and

$h\left(t;\lambda ,\alpha ,p\right)=\frac{\frac{\alpha }{\lambda }{\left(\frac{t}{\lambda }\right)}^{-\left(\alpha +1\right)}{e}^{-{\left(\frac{t}{\lambda }\right)}^{-\alpha }}}{p+\left(1-p\right)\left(1-\exp \left(-{\left(\frac{t}{\lambda }\right)}^{-\alpha }\right)\right)}$ . (10)

For each failure time related to the i-th individual, it may not be perceived or subject by the right censoring. Furthermore, the random censoring times $C_i$ s are independent of $T_i$ s (non-censored time) and their distribution does not depend on the parameters. In a scenario of a $n$  sample of size, the data set will be describe by$\mathcal{D}=\left(t_i,{\delta }_i\right)$ , where $t_i=\min \left(T_i,C_i\right)$ and ${\delta }_i=I\left(T_i\le C_i\right)$ . This general random censoring scheme has as special case type I and II censoring mechanism. The likelihood function is given by

$L\left(\theta ;\mathcal{D}\right)={\displaystyle \prod }_{i=1}^{n}f{\left(t_i;\theta \right)}^{{\delta }_i}S{\left(t_i;\theta \right)}^{1-{\delta }_i}$

Let$T_1, \dots ,T_n$ be a random sample of LF distribution, the likelihood function considering data with random censoring is given by

$L\left(\lambda ,\alpha ,p;\mathcal{D}\right)=\frac{{\alpha }^d{\left(1-p\right)}^d}{{\lambda }^{-d\alpha }}{\displaystyle \prod }_{i=1}^n{t}_i^{-{\delta }_i\left(\alpha +1\right)}\exp \left(-{\displaystyle \sum }_{i=1}^n{\delta }_i{\left(\frac{t_i}{\lambda }\right)}^{-\alpha }\right)\times {\left(p+\left(1-p\right)\left(1-e^{-{\left(\frac{t}{\lambda }\right)}^{-\alpha }}\right)\right)}^{1-{\delta }_i}$

where $d={\displaystyle \sum }_{i=1}^n{\delta }_i$ . The log-likelihood function is given as

$l\left(\lambda ,\alpha ,p;\mathcal{D}\right)=d\log \left(\alpha \right)+d\log \left(1-p\right)+d\alpha \log \lambda -\left(\alpha +1\right){\displaystyle \sum }_{i=1}^n{\delta }_i\log t_i-{\displaystyle \sum }_{i=1}^n{\delta }_i{\left(\frac{\lambda }{t_i}\right)}^{\alpha }+{\displaystyle \sum }_{i=1}^n\left(1-{\delta }_i\right)\log \left(p+\left(1-p\right)\left(1-e^{-{\left(\frac{t}{\lambda }\right)}^{-\alpha }}\right)\right).$ (11)

The maximum likelihood estimators (MLEs) are widely explored as statistical inferential methodology due its many desirable properties, in which includes consistency, asymptotic efficiency and invariance. The MLEs are obtained from the maximization of the log-likelihood function (11). Before we derive the MLEs of the LF, let us define the following function

${\eta }_j\left(\lambda ,\alpha ,p;\mathcal{D}\right)={\displaystyle \sum }_{i=1}^n\left(1-{\delta }_i\right)\frac{\log S\left(t_i;\theta \right)}{\partial {\theta }_j}, j=1,2,3.$

Then, the likelihood equations are given by

$\frac{d\alpha }{\lambda }-\alpha {\lambda }^{\alpha -1}{\displaystyle \sum }_{i=1}^n{\delta }_i{t}_i^{-\alpha }+{\eta }_1\left(\lambda ,\alpha ,p;\mathcal{D}\right)=0$ ,

$\frac{d}{\alpha }+d\log \lambda -{\displaystyle \sum }_{i=1}^n{\delta }_i{\left(\frac{\lambda }{t_i}\right)}^{\alpha }\log \left(\frac{\lambda }{t_i}\right)-{\displaystyle \sum }_{i=1}^n{\delta }_i\log t_i+{\eta }_2\left(\lambda ,\alpha ,p;\mathcal{D}\right)=0$ and

$\frac{d}{p-1}+{\eta }_3\left(\lambda ,\alpha ,p;\mathcal{D}\right)=0$ .

The maximization of the log-likelihood function can be performed directly by using existing statistical packages. Further information about the numerical procedures will be discussed in the next section.

According to Migon et al. [10], under mild conditions the obtained estimators are consistent and efficient with an asymptotically normal joint distribution given by

$\left(\stackrel{}{\widehat{\lambda }},\stackrel{}{\widehat{\alpha }},\stackrel{}{\widehat{p}}\right)~N_3\left(\left(\lambda ,\alpha ,p\right),H^{-1}\left(\lambda ,\alpha ,p\right)\right)$ ,

where${\rm H}\left(\lambda ,\alpha ,P\right)$ , is the 3×3 observed Fisher information matrix and $H_{ij}\left(\lambda ,\alpha ,p\right)$ is the Fisher information given by

$H_{ij}\left(\theta \right)=-\frac{\partial }{\partial {\theta }_i\partial {\theta }_j}l\left(\theta ;\mathcal{D}\right), i,j=1,2,3$ .

Note that, the observed Fisher information matrix was used since it is not possible to compute the expected Fisher information matrix due its lack of closed form expression. For large samples, confidence intervals approximation can be constructed for the individual parameters ${\theta }_i$ i=1,2,3, assuming a confidence coefficient $100\left(1-\gamma \right)\%$ the marginal distributions are given by

${\stackrel{}{\widehat{\theta }}}_i~N\left({\theta }_i,H_{ii}^{-1}\left(\theta \right)\right), i=1,2,3.$

The maximum likelihood method efficiency was analyzed through a simulation study on the LF distribution. This procedure was conducted by computing the mean relative errors (MRE) and the mean square errors (MSE) given by

$MR{E}_i=\frac{1}{N}{\displaystyle \sum }_{j=1}^N\frac{{\widehat{\theta }}_{i,j}}{{\theta }_i}, MS{E}_i=\frac{1}{N}{\displaystyle \sum }_{j=1}^N{\left({\widehat{\theta }}_{i,j}-{\theta }_i\right)}^2, \text{for} i=1,2,3.$

as N is the number of estimates obtained through the MLE approach. The 95% coverage probabilities of the asymptotic confidence intervals were also evaluated. The adopted approach prioritize that the expected MLEs returns the MREs closer to one with smaller MSEs. Additionally, by considering a 95% confidence level, the interval covers the true values of θ closer to 95%.Considering scenarios with sample sizes n=(10,25,50,100, 200) and N=100,000 for the simulation study, two situations are presented by considering the proportion of cure in the population of $0.3$ and$0.5$ . In these cases, the censored proportions are observed in different levels.

In pursuance to find the maximization of the log-likelihood function, described in the equation (11), the package called maxLik available in R developed by Henningsen and Toomet [11] was used. The numerical results are well-behaved since was not found numerical problems using the SANN method (Simulated-annealing), such as failure evidence of convergence or end on multiple maxima. The programs can be obtained, upon request.

The estimates obtained from Tables 1-4 for $\alpha , \lambda$  and $p$ are asymptotically unbiased, implying that MREs tend to one when n increases and the MSEs decrease to zero for n large. Analyzing the MLEs performance, with a coverage probabilities tending to 0.95, good coverage properties may be deliberated for the parameter estimators. In practical applications, those estimation procedures will be relevant as shown in the next section.

In this section, we considered the data set presented by Kersey et al. [12]. The results were collected in a group of 46 patients, per years, upon the recurrence of leukemia whom received autologous marrow. Table 5 shows the full data set (+ indicates censored observations).

The proposed model is compared with some usual long-term survival models, such as the LT Weibull and LT weighted Lindley (Louzada and Ramos, [13]). Different discrimination criterion methods are considered:  the negative of the maximum value of the likelihood function$l\left(\stackrel{}{\widehat{\theta }};t\right)$ , the Akaike information criterion $\left(AIC=-2l\left(\stackrel{}{\widehat{\theta }};t\right)+2k\right)$ and the corrected AIC$\left(AIC+2k\left(k+1\right)/ \left(n-k-1\right)\right)$ , where k is the number of parameters to be fitted. The best model is the one which provides the minimum criterion method values.

Figure 2 presents the empirical survival function adjusted by the Kaplan-Meier estimator and different LT survival distributions.

Figure 2: Survival function adjusted by the empirical survival function (Kaplan-Meier estimator), LT Fréchet, LT Weibull and LT WL distribution.

Table 6 presents the results of the different discrimination criteria for different probability distributions. Comparing the results of the different discrimination methods, we observed that the LT Fréchet distribution has better fit then the LT models under the Weibull and weighted Lindley baseline distribution.

The MLEs were obtained through the same procedure as described in Section 3. The standard error (SE) and the confidence intervals, considering a 95% confidence level for $\alpha , \lambda$ and $p$  are displays in Table 7.

Note that, in Kersey et al. [12] they use the non-parametric KM estimate of the cure fraction in which was $0.20$ where $\left(0.08;0.32\right)$ is the 95% confidence interval. Therefore, results showed to be consistence with Kersey et al. [12] results while our estimate was contained in the non-parametric interval. By using our parametric model the estimate obtained for $p$ was $0.125$ showing an overestimation of the long term survival patients. As it can be seen, through our proposed methodology the data related to the leukemia free-survival times (in months) for the 50 autologous transplant patients can be described by the LF distribution.

In this paper, we have proposed a new long-term survival distribution called long term Fréchet distribution and its mathematical properties were studied. It was presented results towards the maximum likelihood parameters’ estimators and their asymptotic properties. The estimators’ efficient were present in the simulation study as the MLEs for the three unknown parameters obtained acceptable results even for small sample sizes. As such of the real dataset problem, related to the leukemia, free-survival times (in months) for the 50 autologous transplant patients. Many extensions from this present work can be considered, for instance, the parameters estimation may also be studied under an objective Bayesian analysis (Ramos et al., [14,15]) or using different classical methods (Louzada et al., Bakouch et al. [16]). Other approach should be to include covariates under the assumption of Cox model, i.e., proportional hazards. In conclusion, this regression model can be extended for the Bayesian approach as well.

The authors are thankful to the Editorial Board and to the reviewers for their valuable comments and suggestions which led to this improved version.