TEXT SIZE

CrossRef (0)
A new extended Birnbaum-Saunders model with cure fraction: classical and Bayesian approach

Edwin M. M. Ortegaa, Gauss M. Cordeirob, Adriano K. Suzuki1,c, and Thiago G. Ramiresa

aDepartamento de Ciências Exatas, Universidade de São Paulo, Brazil, bDepartamento de Estatística, Universidade Federal de Pernambuco, Brazil, cDepartamento de Matemática Aplicada e Estatística, Universidade de São Paulo, Brazil
Correspondence to: Departamento de Matemática Aplicada e Estatística, Universidade de São Paulo, Avenida Trabalhador São-carlense, 400 - Centro CEP: 13566-59, São Carlos-SP, Brasil. E-mail: suzuki@icmc.usp.br
Received May 10, 2017; Revised July 6, 2017; Accepted July 7, 2017.
Abstract

A four-parameter extended fatigue lifetime model called the odd Birnbaum-Saunders geometric distribution is proposed. This model extends the odd Birnbaum-Saunders and Birnbaum-Saunders distributions. We derive some properties of the new distribution that include expressions for the ordinary moments and generating and quantile functions. The method of maximum likelihood and a Bayesian approach are adopted to estimate the model parameters; in addition, various simulations are performed for different parameter settings and sample sizes. We propose two new models with a cure rate called the odd Birnbaum-Saunders mixture and odd Birnbaum-Saunders geometric models by assuming that the number of competing causes for the event of interest has a geometric distribution. The applicability of the new models are illustrated by means of ethylene data and melanoma data with cure fraction.

Keywords : Bayesian estimation, Birnbaum-Saunders distribution, lifetime data distribution, maximum likelihood estimation, simulations
1. Introduction

The two-parameter Birnbaum-Saunders (BS) distribution is a popular distribution for modeling lifetime data and phenomenon with monotone failure rates that has been applied in several research areas such as engineering, hydrology and survival analysis. However, it does not provide a reasonable parametric fit for non-monotone failure rates such as bathtub shaped and unimodal failure rates, which are common in reliability and biological studies. Bathtub hazard curves have nearly flat middle portions and the associated densities have a positive anti-mode. Unimodal failure rates can be observed in the course of a disease whose mortality reaches a peak after a finite period and then declines gradually.

Many lifetime distributions have been constructed with a view for applications in several areas such as survival analysis, reliability engineering, demography and actuarial studies, and hydrology. The BS model, also known as the fatigue life distribution, is a very popular model which has been extensively used for modeling the failure times of fatiguing materials and lifetime data in the fields cited above. Many generalizations of the BS model have been recently proposed by introducing shape parameters to better explore the tails and other properties of the generated distributions. For example, Cordeiro and Lemonte (2011) discussed the βBS model, Castillo et al. (2011) introduced a new extension of the BS distribution based on the family of the epsilon-skew-symmetric distributions, Ortega et al. (2012) developed the censored log-BS regression model, Cordeiro et al. (2013) proposed the McDonald-BS class of distributions, Lemonte (2013) defined a new extension of the BS distribution and Cordeiro and Lemonte (2014) studied exponentiated generalized BS distribution.

A random variable W having the BS(λ, β) distribution with shape parameter λ > 0 and scale parameter β > 0 is defined by

$W=β[λ U2+{(λ U2)2+1}12]2,$

where U is a standard normal random variable. The cumulative distribution function (cdf) of W is given by

$Gλβ(w)=Φ(v), w>0,$

where v = λ−1ρ(w/β), ρ(u) = u1/2u−1/2 and Φ(·) is the standard normal cumulative function. The parameter β is the median of the distribution: G(β) = Φ(0) = 1/2. For any k > 0, kW ~ BS(λ, ). Kundu et al. (2008) discussed the shape of its hazard function. The probability density function (pdf) corresponding to (1.1) is

$gλ,β(w)=κ(λ,β)w-32(w+β)exp {-τ(w/β)2λ2}, w>0,$

where $κ(λ,β)=exp(λ-2)/(2λ2πβ)$ and τ(u) = u + u−1.

The cdf of the odd Birnbaum-Saunders (“OBS” for short) distribution with an additional shape parameter α > 0 is defined by

$Πα,λ,β(w)=∫0Gλ,β(w)G¯λ,β(w)α xα-1(1+xα)2dx=Gλ,β(w)αGλ,β(w)α+G¯λ,β(w)α,$

where λ,β(w) = 1 − Gλ,β(w).

The BS cdf Gλ,β(w) is clearly a special case of (1.3) when α = 1. We note that there is no complicated function in equation (1.3) in contrast with the βBS model (Cordeiro and Lemonte, 2011), which involves the beta incomplete function and two extra parameters. The OBS density is given by

$πα,λ,β(w)=α gλ,β(w) {Gλ,β(w) [1-Gλ,β(w)]}α-1{Gλ,β(w)α+ [1-Gλ,β(w)]α}2.$

Therefore, the parameter α represents the quotient of the log odds ratio for the generated and baseline distributions.

We now consider that the random variable W has the OBS cdf (1.3) with parameters (α, λ, β)T, say W ~ OBS(α, λ, β).

The cdf and pdf of W are given by

$Πα,λ,β(w)=Φ(v)αΦ(v)α+[1-Φ(v)]α, w>0,$

and

$πα,λ,β(w)=ακ(λ,β)w-32 (w+β) exp {-τ(w/β)2λ2}{Φ(v)[1-Φ(v)]}α-1{Φ(v)α+[1-Φ(v)]α}2,$

respectively.

The rest of the paper is organized as follows. In Section 2, we define the odd Birnbaum Saunders geometric (OBSG) model by compounding OBS and geometric distributions. Parameter inference is investigated in Section 3, where two methods (maximum likelihood and Bayesian approach) for parameter estimation are discussed. Further, various simulations are performed for different parameter settings and sample sizes. In Section 4, we introduce two models with long-term survivals and estimate their parameters by maximum likelihood and Bayesian approaches. Applications to two real data sets are addressed in Section 5. The paper provides concludes remarks in Section 6.

2. The odd Birnbaum-Saunders geometric model under activation scheme

Generalizing continuous univariate distributions by introducing additional shape parameters in a distribution is a useful practice to obtain multiple density and hazard rate shapes. Suppose that ${Wi}i=1Z$ are independent and identically distributed (iid) random variables having the OBS distribution with cdf and pdf given by (1.5) and (1.6), respectively.

Let Z be a geometric random variable with probability mass function P(z; p) = (1 − p) pz−1 for zN and p ∈ (0, 1). We define the random variable $T=min({Wi}i=1Z)$ having the OBSG distribution. We provide three motivations for the proposed distribution, where Z has a geometric distribution, which can be applied in some interesting situations such as: (i) Reliability: the new distribution can arise in series and parallel systems with identical components, which appear in many industrial applications and biological organisms; (ii) Times to the first and last failure: suppose the failure of a device occurs due to the presence of an unknown number Z of initial defects of the same kind, which can be identifiable only after causing failure and are repaired perfectly. Let Wi be the time to the failure of the device due to the ith defect, for i ≥ 1. Under the assumptions that the Wi’s are iid OBS random variables independent of Z, the OBSG distribution is appropriate for modeling the times to the first and last failures, respectively; (iii) Time to relapse of cancer under the first-activation scheme: suppose that an individual in the population is susceptible to a certain type of cancer. Let Z be the number of carcinogenic cells for an individual left active after the initial treatment and denote by Wi the time spent for the ith carcinogenic cell to produce a detectable cancer mass, for i ≥ 1. Under the assumptions that {Wi}i≥1 is a sequence of iid OBS random variables independent of Z, we conclude that the time to relapse of cancer of a susceptible individual has the OBSG distribution.

The conditional density function of $T=min({Wi}i=1Z)$ given Z = z is

$f(t∣z;α,λ,β)=zα κ(λ,β)t-32 (t+β) exp {-τ(t/β)2λ2}{Φ(v)[1-Φ(v)]}α-1 [1-Φ(v)]α(z-1){Φ(v)α+[1-Φ(v)]α}2 {Φ(v)α+[1-Φ(v)]α}z-1,$

where ν = λ−1ρ(t/β) and ρ(u) = u1/2u−1/2.

Then, the four-parameter OBSG density function (for t > 0) becomes

$f(t;α,λ,β,p)=(1-p) α κ(λ,β)t-32 (t+β) exp {-τ(t/β)2λ2}{Φ(v)[1-Φ(v)]}α-1{Φ(v)α+[1-Φ(v)]α}2×{1-p[1-Φ(v)]αΦ(v)α+[1-Φ(v)]α}-2,$

where β > 0 is a scale parameter and α > 0, λ > 0 and p ∈ (0, 1) are shape parameters. The random variable T having density function (2.1) is denoted by T ~ OBSG(α, λ, β, p). We can omit sometimes the dependence on the vector of parameters and write f (t) = f (t; α, λ, β, p), etc. The cdf of T (for t > 0) is

$F(t)=Φ(v)αΦ(v)α+(1-p)[1-Φ(v)]α.$

Clearly, the OBS distribution is a special case of (2.2) when p → 0+. The BS distribution can also be obtained when p → 0+ and α = 1. The Birnbaum-Saunders geometric (BSG) distribution is obtained when α = 1.

The hazard rate function (hrf) of T is given by

$S(t)=(1-p)[1-Φ(v)]αΦα(v)+(1-p)[1-Φ(v)]α.$

Plots of the OBSG density, survival and hazard functions for selected parameter values are displayed in Figures 13, respectively.

We show some structural properties of the OBSG distribution in the Appendix.

3. Inference

### 3.1. Maximum likelihood estimation

Several approaches for parameter estimation were proposed in the literature but the maximum likelihood method is the most commonly employed. The maximum likelihood estimators (MLEs) enjoy desirable properties and can be used when constructing confidence intervals and regions as well as in test statistics. The normal approximation for these estimators in asymptotic theory is easily handled either analytically or numerically. Therefore, we consider the estimation of the unknown parameters for the OBSG model by maximum likelihood. Let Ti be a random variable following (2.1) with the vector of parameters θ = (α, λ, β, p)T . Data encountered in survival analysis and reliability studies are often censored. A very simple random censoring mechanism that is often realistic is one in which each individual i is assumed to have a lifetime Ti and a censoring time Ci, where Ti and Ci are independent random variables. Suppose the data consist of n independent observations ti = min(Ti,Ci) for i = 1, …, n. The distribution of Ci does not depend on any unknown parameters of Ti. Parametric inference for such data is usually based on likelihood methods and their asymptotic theory. The censored log-likelihood l(θ) for the model parameters is given by

$l(θ)=n log(1-p)+r log[α κ(λ,β)]-32∑i∈Flog(ti)+∑i∈Flog(ti+β)-12λ2∑i∈Fτ (tiβ)+(α-1)∑i∈Flog{Φ(vi)[1-Φ(vi)]}-2∑i∈Flog {Φ(vi)α+(1-p)[1-Φ(vi)]α}+α∑i∈Clog[1-Φ(vi)]-∑i∈Clog {Φ(vi)α+(1-p)[1-Φ(vi)]α},$

where v = λ−1ρ(ti), ρ(z) = z1/2z−1/2, r is the number of failures and F and C denote uncensored and censored sets of observations, respectively.

This log-likelihood can be maximized either directly by using the R (optim function; https://www.rproject.org/), SAS (NLMIXED procedure; SAS Institute Inc., Cary, NC, USA), or Ox (MaxBFGS function; http://www.doornik.com/ox/), or by solving the nonlinear likelihood equations by its differentiation. The score functions for the parameters α, λ, β, and p are given by

$Uα(θ)=rα+∑i∈Flog {Φ(vi) [1-Φ(vi)]}+∑i∈Clog [1-Φ(vi)]-2∑i∈FΦ(vi)α log [Φ(vi)]+(1-p) [1-Φ(vi)]α log [1-Φ(vi)]Φ(vi)α+(1-p) [1-Φ(vi)]α-∑i∈CΦ(vi)α log [Φ(vi)]+(1-p) [1-Φ(vi)]α log [1-Φ(vi)]Φ(vi)α+(1-p) [1-Φ(vi)]α,Uλ(θ)=rκλ(λ,β)κ(λ,β)+1λ3∑i∈Fτ (tiβ)+(α-1)∑i∈Fϕ(vi)vi(λ) {1-2Φ(vi)Φ(vi) [1-Φ(vi)]}-α∑i∈Cϕ(vi)vi(λ)1-Φ(vi)-2α∑i∈Fϕ(vi)vi(λ) {Φ(vi)α-1-(1-p) [1-Φ(vi)]α-1Φ(vi)α+(1-p) [1-Φ(vi)]α}-α∑i∈Cϕ(vi)vi(λ) {Φ(vi)α-1-(1-p) [1-Φ(vi)]α-1Φ(vi)α+(1-p) [1-Φ(vi)]α},Uβ(θ)=rκβ(λ,β)κ(λ,β)+∑i∈F1ti+β-12λ2∑i∈F{1ti-tiβ2}+(α-1)∑i∈Fϕ(vi)vi(β) {1-2Φ(vi)Φ(vi) [1-Φ(vi)]}-2α∑i∈Fvi(β)ϕ(vi){Φα-1(vi)-(1-p) [1-Φ(vi)]α-1Φα(vi)+(1-p) [1-Φ(vi)]α}-α∑i∈Cvi(β)ϕ(vi)1-Φ(vi)-α∑i∈Cvi(β)ϕ(vi) {Φα-1(vi)-(1-p) [1-Φ(vi)]α-1Φα(vi)+(1-p) [1-Φ(vi)]α},Up(θ)=-n1-p+2∑i∈F[1-Φ(vi)]αΦ(vi)α+(1-p)[1-Φ(vi)]α+∑i∈C[1-Φ(vi)]αΦα(vi)+(1-p)[1-Φ(vi)]α,$

where

$κλ(λ,β)=exp (λ-2) (λ-2-2)2λ42πβ,κβ(λ,β)=-exp (λ-2)2λ2π(β-322),vi(λ)=1λ2 [(tiβ)-12-(tiβ)12],vi(β)=-λ-12 [(tiβ3)12+(1βti)12].$

The MLE θ̂ of θ can be obtained by solving the nonlinear likelihood equations Uα(θ) = 0, Uλ(θ) = 0, Uβ(θ) = 0, and Up(θ) = 0. We employ the NLMIXED procedure in SAS because these equations cannot be solved analytically.

For interval estimation of (α, λ, β, p) and tests of hypotheses on these parameters we obtain the observed information matrix since the expected information matrix requires numerical integration. The elements of the 4 × 4 observed information matrix J(θ), namely Jαα, Jαλ, Jαβ, Jαp, Jλλ, Jλβ, Jλp, Jββ, Jβp, and Jpp, can be calculated numerically.

Under standard regularity conditions, the asymptotic distribution of (θ̂θ) is N4(0, I(θ)−1), where I(θ) is the expected information matrix. This asymptotic behavior is valid if I(θ) is replaced by J(θ̂), i.e., the observed information matrix evaluated atθ̂. The multivariate normal N4(0, J(θ̂)−1) distribution can then be used to construct approximate confidence intervals for the individual parameters.

Further, we can compute the maximum values of the unrestricted and restricted log-likelihoods to obtain the likelihood ratio (LR) statistics to test some sub-models of the OBSG distribution. For example, the test of H0 : α = 1 versus H1 : α ≠ 0 is equivalent to compare the BSG distribution with the OBSG distribution. In this case, the LR statistic is

$w=2 {l (α^,λ^,β^,p^)-l (1,λ˜,β˜,p˜)},$

where α̂, λ̂, β̂, and are the MLEs under H1 and λ̃, β̃, and are the estimates under H0.

### 3.2. Bayesian estimation

An alternative statistical approach is to use the Bayesian method that allows the incorporation of previous knowledge of the parameters through an informative prior distribution. A non-informative prior structure may be considered when previous knowledge is unavailable. In this work we adopt proper prior distributions according to variations of the parametric space, but ensuring non-informativeness according to the fixed hyper-parameters that lead to such a situation. We assume α, λ, β, and p are mutually independent, that is, we consider a prior distribution of θ = (α, λ, β, p) given by π(θ) = π(α)π(λ)π(β)π(p), where α ~ Gamma(0.1, 0.01), λ ~ Gamma(0.1, 0.01), β ~ Gamma(0.1, 0.01) and p ~ Beta(1, 1), i.e., all the hyper-parameters have been specified to express non-informative priors. In order to generate our samples we then implement an Metropolis-Hasting algorithm within Gibbs iterations (Chib and Greenberg, 1995). Simulations were performed using the OpenBUGS software (Spiegelhalter et al., 2007). In all computations we simulate two chains of size 60,000 for each parameter, disregarding the first 20,000 iterations to eliminate the effect of the initial values and to avoid autocorrelation problems. We then consider a spacing of size 10 and obtain an effective sample of size 8,000, upon which the posterior inference is based on. The MCMC convergence was monitored according to the methods recommended by Cowles and Carlin (1996), as well as trace plots.

3.2.1. Model comparison criteria

In the literature, there are various methodologies which intend to analyze the suitability of distribution, as well as select the best fit from among a collection of distributions. One of the most used approaches is derived from the conditional predictive ordinate (CPO) statistic. For a detailed discussion on the CPO statistic and its applications to model selection, see Gelfand et al. (1992).

Let denote the full data and denote the data with the deleted ith observation. We denote the posterior density of θ given by π(θ| ), i = 1, …, n. For the ith observation, CPOi can be written as

$CPOi=∫θ∈Θf(ti∣θ)π (θ∣D(-i)) dθ={∫θπ(θ∣D)f(ti∣θ)dθ}-1.$

The CPOi can be interpreted as the height of the marginal density of the time for an event at ti.

Therefore, high CPOi implies a better fit of the model. For model comparisons we use the log pseudo marginal likelihood (LPML) defined by $LPML=Σi=1nlog(CPO^i)$. The higher the LPML value, the better the fit of the model.

Other criteria, such as the deviance information criterion (DIC) proposed by Spiegelhalter et al. (2002), the expected Akaike information criterion (EAIC)-Brooks (2002), and the expected Bayesian (or Schwarz) information criterion (EBIC)-Carlin and Louis (2001) can also be used. These other criteria are based on the posterior mean of deviance, which can be approximated by $d¯=Σq=1Qd(θq)/Q$, where $d(θ)=-2Σi=1n log [f(ti∣θ)]$. The DIC criterion can be estimated using the MCMC output by $DIC^=d¯+ρd^=2d¯-d^$, where ρD is the effective number of parameters defined as E{d(θ)} − d{E(θ)}, where d{E(θ)} is the deviance evaluated at the posterior mean. Similarly, EAIC and EBIC criteria can be estimated by means of $EAIC^=d¯+2#(θ)$ and $EAIC^=d¯+#(θ) log(n)$, where #(θ) is the number of model parameters. Comparing alternative models, the preferred model is the one with the smallest criteria values.

### 3.3. Simulation

We simulate the OBSG distribution (for α = 0.7, λ = 0.8, β = 2, p = 0.1) from equation (A.3) by using a random variable U with a uniform distribution in (0, 1). We simulate n = 100, 200, 300, and 500 variates and, for each replication, we compute the MLEs and the posterior mean of the parameter α, λ, β, and p. We repeat this process 1,000 times for classical inference and 200 times for the Bayesian procedure. We determine the average estimates (mean), biases, and means squared errors (MSEs). The results are given in Table 1 provides the results.

The results of the Monte Carlo study in Table 1 indicate that the MSEs of the mean of α, λ, β and p decay toward zero as the sample size increases, as expected under standard asymptotic theory. As n increases, the mean of the parameters tend to be closer to the true parameter values. These results support a asymptotic normal distribution that provides an adequate approximation to the finite sample distribution of the MLEs. The normal approximation can often be improved using bias adjustments for these estimators. Approximations to the their biases in simple models may be obtained analytically. Bias correction typically does a very good job correct the MLEs; however, it may either increase the MSEs. Whether bias correction is useful in practice depends on the shape of the bias function and the variance of the MLE. Several cumulants of log likelihood derivatives need to be obtained in order to improve the accuracy of the estimators using analytical bias reduction; however, these derivatives are notoriously cumbersome for the proposed model. Also, the results under a Bayesian approach (Table 1) are similar to those obtained via MLEs.

4. The odd Birnbaum-Saunders model for survival data with long-term survivors

For censored survival times, the presence of an immune proportion of individuals who are not subject to death, failure or relapse may be indicated by a relatively high number of individuals with large censored survival times. Here, the OBS model is modified in two different ways for the possible presence of long-term survivors in the data.

### 4.1. The odd Birnbaum-Saunders mixture model

In population-based cancer studies, cure is said to occur when mortality in the group of cancer patients returns to the same level as that expected in the general population. The cure fraction is a useful measure that is of interest to patients analyzing trends in cancer patient survival. Models for survival analysis typically assume that every subject in the study population is susceptible to the event under study and will eventually experience such event if follow-up is sufficiently long. However, there are situations when a fraction of individuals are not expected to experience the event of interest such as when those individuals are cured or not susceptible. Cure rate models have been used to model time-to-event data for various types of cancers, including breast cancer, non-Hodgkin’s lymphoma, leukemia, prostate cancer, and melanoma. The most popular type of cure rate models are often the mixture models (MMs) introduced by Boag (1949), Berkson and Gage (1952), and Farewell (1982). Additionally, MMs allow both the cure fraction and the survival function of uncured patients (latency distribution) to depend on covariates. Longini and Halloran (1996) and Price and Manatunga (2001) have further introduced frailty to MMs for individual survival data. Further, Peng and Dear (2000) investigated a nonparametricMMfor cure estimation, Sy and Taylor (2000) considered estimation in a proportional hazards cure model and Yu and Peng (2008) extended MMs to bivariate survival data by modeling marginal distributions. Fachini et al. (2014) recently proposed a scale model for bivariate survival times based on the copula that model the dependence of bivariate survival data with cure fraction, Hashimoto et al. (2015) introduced the new long-term survival model with interval-censored data, Ortega et al. (2015) proposed the power series beta Weibull regression model to predict breast carcinoma, Lanjoni et al. (2016) conducted extended Burr XII regression models and Ortega et al. (2017) proposed regression models generated by gamma random variables with long-term survivors.

To formulate the odd Birnbaum-Saunders mixture (OBSM) model, we consider that the population under study is a mixture of susceptible (uncured) individuals, who may experience the event of interest, and non-susceptible (cured) individuals, who never will experience the event of interest (Maller and Zhou, 1996). This approach allows to estimate simultaneously whether the event of interest will occur, which is called incidence, and if it does occur then when it will occur, which is called latency. Let Ni (for i = 1, …, n) be the indicator denoting that the ith individual is susceptible (Ni = 1) or non-susceptible (Ni = 0). The MM is given by

$Spop(ti)=π+(1-π)S(ti∣Ni=1),$

where Spop(ti) is the unconditional survival function of Ti for the entire population, S (ti|Ni = 1) is the survival function for susceptible individuals and π = P(Ni = 0) is the probability of cure for an individual. The OBSM model is defined by assuming that the survival function for susceptible individuals in (4.1) is given by

$S(ti∣Ni=1)=[1-Φ(vi)]αΦ(vi)α+[1-Φ(vi)]α.$

Let θ = (π, α, λ, β) be the vector of parameters of interest for the OBSM model. For π = 0, the OBSM model reduces to the OBS model.

Suppose we have data in the form ti, i = 1, …, n, where ti denotes the observed logarithm of time for the ith subject, i.e., ti = min {Ti,Ci}, Ti is the lifetime for the ith individual and Ci is the censoring time for the ith individual. Under this assumption, the contribution of an individual that failed at ti to the likelihood function is given by

$(1-π) α κ(λ,β) ti-32 (ti+β) exp {-τ(ti/β)2λ2}{Φ(vi)[1-Φ(vi)]}α-1{Φ(vi)α+[1-Φ(vi)]α}2$

and the contribution of an individual that is at risk at ti is

$π+(1-π)[1-Φ(vi)]αΦ(vi)α+[1-Φ(vi)]α.$

The full log-likelihood function for the parameter vector θ = (π, α, λ, β) is given by

$l(θ)=r log[(1-π)α κ(λ,β)]-32∑i∈Flog(ti)+∑i∈Flog(ti+β)-12λ2∑i∈Fτ (tiβ)+(α-1)∑i∈Flog{Φ(vi)[1-Φ(vi)]}-2∑i∈Flog {Φα(vi)+[1-Φ(vi)]α}+∑i∈Clog {π+(1-π)[1-Φ(vi)]αΦ(vi)α+[1-Φ(vi)]α},$

where F and C denote the sets of individuals who failed and are censoring, respectively. The MLEθ̂ of θ can be obtained by maximization of (4.3). The asymptotic covariance matrix of θ̂ can be estimated by the Hessian matrix evaluated at θ̂.

### 4.2. The odd Birnbaum-Saunders geometric model with cure fraction

Models with a cure fraction play an important role in survival analysis. Recently, some of these models have been published in the literature. For example, Cooner et al. (2007) proposed a flexible cure rate under latent activation schemes, Ortega et al. (2009) presented the generalized log-gamma regression model with cure fraction and Cancho et al. (2011) defined a flexible cure rate survival model by considering that the number of competing causes of the event of interest follows the Conway-Maxwell distribution and the time for the event has the generalized gamma distribution. Recently, Hashimoto et al. (2014) proposed the Poisson BS model with long-term survivors and Suzuki et al. (2016) studied the Poisson inverse-Gaussian regression cure rate model. In the same context, we present the OBSG model with a cure fraction.

For an individual in the population, let N denote the unobservable number of causes of the event of interest for this individual. We assume that N has a geometric distribution. The time for the jth cause to produce the event of interest is denoted by Wj, j = 1, …, N. We consider that, conditional on N, the Wj’s are iid random variables with the OBS cdf. Further, we assume that N is independent of W1,W2, …. The observable time to the event is defined by T = min{W1, …,WN}, if N ≥ 1 and T = ∞ if N = 0, with P (T = ∞|N = 0) = 1. Under this setup, the survival function for the population is given by

$Spop(t)=P(N=0)+P(W1>t,…,WN>t∣N≥1)P(N≥1).$

Tsodikov et al. (2003), among others, demonstrated that Spop(t) = Ap[S (t)], where Ap(·) is the probability generating function (pgf) of the number of competing causes (N). We consider that N has the geometric distribution with parameter 0 ≤ p < 1. The pgf of N is given by $Ap(w)=∑m=0∞pm wm=(1-p)/(1-pw)$, 0 ≤ w ≤ 1, and then the survival function for the population becomes

$Spop(t)=Ap [S(t)]=(1-p) {1-p[1-Φ(v)]αΦ(v)α+[1-Φ(v)]α}-1.$

Equation for Spop(t) is called improper survival function because it has two main properties:

$limt→0 Spop(t)=1 and limt→0 Spop(t)=p0,$

where p0 denoted the cured fraction, in this model the cured fraction is given by p0 = 1 − p. The pdf corresponding to (4.5) reduces to

$fpop(t)=p(1-p)α κ(λ,β)t-32 (t+β)exp {-τ(t/β)2λ2} {Φ(v)[1-Φ(v)]}α-1×{Φ(v)α+[1-Φ(v)]α}-2 {1-p [1-Φ(v)]αΦ(v)α+ [1-Φ(v)]α}-2.$

Equation (4.6) refers to the OBSG model with long-term survivors in a competitive risk structure. The associated population hrf is

$hpop(t)=pα κ(λ,β)t-32(t+β)exp {-τ(t/β)2λ2} {Φ(v)[1-Φ(v)]}α-1×{Φ(v)α+[1-Φ(v)]α}-2 {1-p [1-Φ(v)]αΦ(v)α+ [1-Φ(v)]α}-1.$

We note that fpop(t) and hpop(t) are improper functions, since Spop(t) is not a proper survival function.

We consider a situation where the time to the event of interest is not completely observed but subjected to right censoring. Let Ci denote the censoring time. We observe ti = min{Ti,Ci}, for i = 1, …, n. Let θ = (p, α, λ, β)T be the parameter vector, from n pairs of times and censoring indicators t1, …, tn. The full log-likelihood function under non-informative censoring can be expressed as

$l(θ)=r log[p(1-p)α κ(λ,β)]+(n-r) log(1-p)-32∑i∈Flog(ti)+∑i∈Flog(ti+β)-12λ2∑i∈Fτ (tiβ)+(α-1)∑i∈Flog{Φ(vi)[1-Φ(vi)]}-2∑i∈Flog{Φα(vi)+[1-Φ(vi)]α}-2∑i∈Flog {1-p [1-Φ(vi)]αΦ(vi)α+[1-Φ(vi)]α}-∑i∈Clog {1-p [1-Φ(vi)]αΦ(vi)α+[1-Φ(vi)]α},$

where r is the number of failures and F and C denote the uncensored and censored sets of observations, respectively.

The MLE θ̂ of the parameter vector θ can be determined by maximizing (4.7). We use the NLMIXED procedure in SAS to compute θ̂. The inference method for θ = (p, α, λ, β)T can be based on the normal approximation

$(p^,α^,λ^,β^)T~N4 {(p,α,λ,β)T,J-1 (θ^)},$

where J−1(θ) = {−2l(θ)/∂θθT } is the 4 × 4 observed information matrix.

The Bayesian procedure was used analogously to Section 3.2. We simulate the OBSG cure fraction distribution (for α = 1.5, λ = 1.5, β = 0.5, p = 0.25) and the OBSM distribution (for α = 2.0, λ = 0.5, β = 1.5, p = 0.2). Censoring times were then sampled from the uniform distribution in the (0, τ) interval, where τ controlled the censoring proportion of the uncured population. Here, the proportions of censored observations were approximately 42% and 38%, respectively.

Analogously to the previous simulation study, 1,000 data sets for classical inference and 200 data sets for Bayesian procedure were generated each one with n = 100, 200, 300, and 500 variates. For each replication, we compute the MLEs and the posterior mean of the parameter α, λ, β, and p. We determine the average estimates (mean), biases, and MSEs.

Tables 2 and 3 show the results for the OBSGcr model and OBSM model, respectively. Similar results to the previously simulation study can also be observed.

5. Applications

In this section, we provide two applications to real data to demonstrate the flexibility of the OBS, OBSM, and OBSG models. The computations are performed using the NLMixed subroutine in the SAS package.

### 5.1. Application 1: Ethylene data

These data are taken from a study by the University of São Paulo, ESALQ (Laboratory of Physiology and Post-colheita Biochemistry), which evaluated the effects of mechanical damage on banana fruits (genus Musa spp.); see Saavedra del Aguila et al. (2010) for more details. The major problem affecting bananas during and after harvest is the susceptibility of the mature fruit to physical damage caused by transport and marketing. Ethylene is a plant hormone important in post-harvest fruit. A high ethylene production can generate a fast senescence of fruit. We use n = 630 points on ethylene data and assume a parent normal distribution.

Recently, Alizadeh et al. (2015) proved that the Kumaraswamy odd log-logistic-normal (KwOLLN) distribution yields a better fit to these data than the McDonald normal (McN) (Cordeiro et al., 2012), beta normal (BN) (Alexander et al., 2012), and Kumaraswamy normal (KwN) (Cordeiro et al., 2010). We now compare the results of the fits of the KwOLLN, McN, BN, and KwN models with the OBSG model as well as compare two sub-models of the OBSG model, namely the OBS and BS models.

Table 4 lists the MLEs (and the corresponding standard errors (SEs) in parentheses) of the model parameters and the values for the fitted models of the statistics: Akaike information criterion (AIC), Bayesian information criterion (BIC), and consistent Akaike information criterion (CAIC). The figures in this table indicate that the OBSG model could be selected as the best model because it has the lowest AIC, BIC, and CAIC values among the values for the fitted models. In addition, SEs of the estimates for all fitted models are also quite small.

Formal tests for the extra skewness parameters in the OBSG model can be based on LR statistics described in Section 3. Table 5 list the values of applying the LR statistics to the ethylene data. We then reject the null hypotheses of the two LR tests in favor of the OBSG distribution (Table 5). More information is provided by a visual comparison of the histogram of the data with the fitted density functions. Figure 4 displays the plots of the fitted OBSG, OBS, KwOLLN, and McN densities.

The empirical scaled total-time-on-test (TTT) transform can be used to identify the shape of the hazard function. The TTT plot for the ethylene data given in Figure 5(a) shows a unimodal and bathtub-shaped hrf and therefore indicates the appropriateness of the OBSG distribution as indicated by Figure 3. Further, Figure 5(b) reveals that the estimated hrf of the OBSG distribution can capture both behaviors (unimodal and bathtub-shaped) are presented in the ethylene data. The OBSG distribution provides a closer fit to the histogram of the current data than the other models.

Under a Bayesian approach we also fitted the OBSG, OBS, and BS models for the ethylene data as well as obtained the values DIC, EAIC, EBIC, and LPML criteria to compare these models. Table 6 shows the results that indicate estimates close to those obtained via MLEs; therefore, the OBSG model presented the best fit when using criteria found in

### 5.2. Application 2: Melanoma data with cure fraction

Here, we discuss thee application of the OBSM and OBSG models with cure fraction to cancer recurrence with long-term survivors. These data are part of a study on cutaneous melanoma (a type of malignant cancer) for the evaluation of postoperative treatment performance with a high dose of a certain drug (interferon alfa-2b) in order to prevent recurrence. Patients from the years 1991 to 1995 were included in the study with follow-up conducted until 1998. The data were collected by Ibrahim et al. (2001). The survival time T is defined as the time until the patient’s death. The original sample size was n = 427 patients, 10 of whom did not present a value for the explanatory variable tumor thickness. When such cases were removed, a sample of size n = 417 patients was retained. A sample of size n = 417 patients was retained after the removal of these cases. The percentage of censored observations was 56%. Tables 7 and 8 lists the MLEs and Bayesian estimates of the model parameters.

The values of AIC, CAIC, and BIC statistics (Table 7) and of the DIC, EAIC, EBIC, and LPML statistics (Table 8) are smaller for the OBSM and OBSG models when compared to the values of the Birnbaum-Saunders mixture (BSM) and BSG models thus indicating that the first two models could be chosen as the best models to fit these data.

The LR statistics given in Table 9 indicate that the OBSM and OBSG models are superior to BSM and BSG models in terms of model fitting.

In order to assess if the model is appropriate, Figures 6 and 7 display the empirical survival functions and the estimated OBS, OBSM, BSG, and OBSG survival functions, respectively, from which a significant fraction of survivors can be observed. It can be verified that the OBSM and OBSG models with cure fraction provide good fits to these data.

6. Conclusions

Many lifetime distributions have been constructed with a view for applications in several areas, in particular, survival analysis, reliability engineering, demography, and actuarial studies, and hydrology. Generalizing continuous univariate distributions by introducing additional shape parameters is an important way to provide very different shapes for the density and hrf of the generated model. We introduce and study a new lifetime model called the OBSG distribution. We derive some general mathematical properties based on a simple mixture representations for its density function. Cure rate models have been used to model time-to-event data for various types of cancers. Based on OBSG distribution, we define two new models with cure rate called the OBSM and OBSG models, where the unobservable number N of causes of the event of interest has a geometric distribution and the time for the jth cause to produce this event is modeled by a random variable Wj such that, conditional on N, the Wj’s are considered iid random variables with the OBS distribution independent of N. Several simulation studies are performed for different parameter settings, sample sizes, and censoring percentages. Maximum likelihood and Bayesian approaches are described to estimate the model parameters. Two empirical applications to real data are presented in order to illustrate the flexibility of the proposed models.

Acknowledgments

We are very grateful to a referee and an associate editor for helpful comments that improved the paper. We gratefully acknowledge financial support from CAPES and CNPq.

Appendix
Structural properties of the OBSG distribution

The necessity and importance of the quantile function, moments and generating function is always important in any statistical analysis such as applied work.

### Quantile function

Let QOBS(u) = Πα,λ,β(t)−1(u) be the qf of the OBS model (for 0 < u < 1). Inverting F(t) = u in (2.2), we obtain the qf of T as

$t=F-1(u)=QOBS [u(1-p)1-p u].$

Hence, equation (A.1) reveals that the OBSG qf can be expressed in terms of the OBS qf. Using (2.2), we obtain

$Φ(v)=au,α,p=[u (1-p)]1α(1-u)1α+[u (1-p)]1α,$

where v = λ−1ρ(t/β) and ρ(u) = u1/2u−1/2. After some algebra, we can write from (A.1)

$t-λ(β t)12 QN(au,α,p)-β=0,$

where QN(·) denotes the standard normal qf. Solving equation (A.2), the OBSG qf becomes

$t=0.5 β λ QN(au,α,p){λ QN(au,α,p)+[λ QN(au,α,p)]2+4}+β.$

We can write the standard normal qf in terms of the inverse error function, namely

$QN(au,α,p)=2 erf-1(2 au,α,p-1).$

An expansion of the inverse error function erf−1(x) can be obtained from Wolfram website ( http://mathworld.wolfram.com/InverseErf.html) as

$erf-1(x)=π(12x+124π x3+7960π2 x5+12780640π3 x7+⋯).$

Quantiles of interest can be determined from (A.1) or (A.3) by substituting appropriate values for u. In particular, the median of T follows when u = 1/2. By using (A.3), we have

$t=β λ2QN(a0.5,α,p){λ QN(a0.5,α,p)+[λ QN(a0.5,α,p)]2+4}+β,$

where

$QN(a0.5,α,p)=(1-p)1α1+(1-p)1α.$

We can use (A.3) for simulating OBSG random variables by setting u as a uniform random variable in the unit interval (0, 1).

### Mixture representation

For any real α > 0, we can expand Φ(v)α as

$Φ(v)α=∑k=0∞vk Φ(v)k,$

where

$vk=vk(α)=∑j=k∞(-1)k+j(αj)(jk).$

Further, the generalized binomial expansion holds

$[1-Φ(v)]α=∑k=0∞(-1)k(αk)Φ(v)k.$

Inserting (A.4) and (A.5) in equation (2.2), we obtain

$F(t)=∑k=0∞vk Φ(v)k∑k=0∞wk Φ(v)k,$

where $wk=wk(α)=vk+(-1)k(1-p) (αk)$.

The ratio of the two power series can be expressed as

$F(t)=∑k=0∞dk Φ(v)k,$

where d0 = v0/w0 and the coefficients dk’s (for k ≥ 1) are determined from the recurrence equation

$dk=w0-1 (vk-w0-1∑r=1kwr vk-r).$

By differentiating (A.6), we obtain

$f(t)=∑k=0∞dk+1 hk+1(t),$

where hk+1(t) = (k + 1)Φ(v)kgλ,β(t) is the exponentiated Birnbaum-Saunders (exp-BS) model with power parameter k + 1.

Cordeiro and Lemonte (2011) recently proposed the βBS distribution which has the exp-BS distribution as a special case when b = 1. Therefore, some properties of the OBSG distribution can be obtained from the exp-BS distribution.

The formulae derived from equation (A.7) can be easily handled in most symbolic computation software platforms such as Maple, Mathematica, and Matlab, which have the ability to deal with analytic expressions of formidable size and complexity. In addition, we can change infinity by twenty or thirty for the majority of these formulae

### Moments

Some of the most important features and characteristics of a distribution can be studied through moments. Cordeiro and Lemonte (2011) determined the probability weighted moments (PWMs) of the BS distribution since they are required for the ordinary moments of the βBS distribution. Here, we also use the PWMs of the BS distribution to demonstrate the moments of the OBSG distribution. The PWMs (for p ≥ 1, r ≥ 0) of the BS distribution are defined by τp,r = E[TpGλ,β(T)r]. We obtain

$τp,r=βp2r∑j=0r(rj)∑k1,…,kj=0∞A(k1,…,kj)∑m=02sj+j(-1)m(2sj+jm)β-2sj+j-2m2I(p+2sj+j-2m2,λ),$

where

$I(p,λ)=Kp+12 (λ-2)+Kp-12 (λ-2)2K12 (λ-2)$

and Kν(z) denotes the modified Bessel function of the third kind with ν representing its order and z the argument. Its integral representation is $Kν(z)=0.5∫-∞∞exp{-z cosh(t)-νt} dt$ (Watson, 1995).

From now on, let Ym+1 have the exp-BS distribution with power parameter m + 1. A first formula for the sth moment of T follows from (A.7) as

$E(Ts)=∑m=0∞dm+1 E (Ym+1s).$

Based on the results by Cordeiro and Lemonte (2011), the moments of T are given by

$E(Ts)=∑m=0∞(m+1) dm+1τs,m,$

Where τs,m comes from (A.8).

A second formula for E(Ts) follows from (A.7) in terms of the OBS qf as

$E(Ts)=∑k=0∞(k+1) dk+1 κ(s,k),$

where

$κ(s,k)=∫0∞ts Gλ,βk(t) gλ,β(t)dx=∫01QOBS(u)s ukdu.$

### Generating function

The moment generating function of T can be determined from (A.7) using a result by Cordeiro and Lemonte (2011) that can be expressed as

$M(s)=∑m,k=0∞(m+1) dm+1 τk,mk!sk,$

where τk,m is given by (A.8).

Figures
Fig. 1. The odd Birnbaum-Saunders geometric density function for some parameter values: (a) For values λ = 0.1 and β = 0.5, (b) For values λ = 0.1, β = 5.0, and α = 0.1.
Fig. 2. The odd Birnbaum-Saunders geometric survival function for some parameter values: (a) For values β = 3.0 and p = 0.8, (b) For values λ = 0.15, β = 3.0, and α = 0.3.
Fig. 3. The odd Birnbaum-Saunders geometric hazard rate function for some parameter values: (a) For values β = 3.0 and p = 0.8, (b) For values λ = 0.15, β = 3.0, and α = 0.3.
Fig. 4. Estimated densities for the: (a) OBSG, OBS, and BS models; (b) OBSG, KwOLLN, and McN models fitted to the ethylene data. OBSG = odd Birnbaum-Saunders geometric; OBS = odd Birnbaum-Saunders; BS = Birnbaum-Saunders; KwOLLN = Kumaraswamy odd log-logistic-normal; McN = McDonald normal.
Fig. 5. TTT plot (a) and estimated hazard function (b) for the OBSG distribution for the ethylene data. TTT = total-time-on-test; OBSG = odd Birnbaum-Saunders geometric.
Fig. 6. (a) Kaplan-Meier curves and BSM model and the estimated cure fraction for the cutaneous melanoma data and (b) Kaplan-Meier curves and estimated survival function for the OBSM model. BSM = Birnbaum-Saunders mixture; OBSM = odd Birnbaum-Saunders mixture; CL = confidence level.
Fig. 7. (a) Kaplan-Meier curves and BSG model and the estimated cure fraction for the cutaneous melanoma data and (b) Kaplan-Meier curves and estimated survival function for the OBSG model. BSG = Birnbaum Saunders geometric; OBSG = odd Birnbaum-Saunders geometric.
TABLES

### Table 1

The mean, biases, and MSEs of the OBSG distribution for α = 0.7, λ = 0.8, β = 2, and p = 0.1, and n = 100, 200, 300, and 500

Sample sizeParameter (true value)Classical estimationBayesian estimation

MeanBiasMSEMeanBiasMSE
N = 100p (0.1)0.12380.02380.02200.12270.02270.0270
α (0.7)0.6562−0.04380.07130.6949−0.00510.0511
λ (0.8)0.7553−0.04470.05880.83230.03230.0453
β (2.0)2.04100.04100.05432.03770.03770.0324

N = 200p (0.1)0.11690.01690.01630.11940.01940.0261
α (0.7)0.6920−0.00800.05130.6970−0.00300.0399
λ (0.8)0.7953−0.00470.04190.82340.02340.0340
β (2.0)2.05620.05620.04122.03390.03390.0198

N = 300p (0.1)0.10970.00970.01290.12070.02070.0256
α (0.7)0.6952−0.00480.03950.6874−0.01260.0353
λ (0.8)0.7988−0.00120.03190.81050.01050.0295
β (2.0)2.03760.03760.02932.02850.02850.0140

N = 500p (0.1)0.10800.00800.00970.12180.02180.0250
α (0.7)0.70170.00170.02380.70270.00270.0205
λ (0.8)0.80760.00760.01910.82110.02110.0179
β (2.0)2.03280.03280.02052.02280.02280.0098

MSEs = means squared errors; OBSG = odd Birnbaum-Saunders geometric.

### Table 2

Mean, biases, and MSEs of the OBSGcr distribution for α = 1.5, λ = 0.5, β = 2.0, and p = 0.25 and n = 100, 200, 300, and 500

Sample sizeParameter (true value)Classical estimationBayesian estimation

MeanBiasMSEMeanBiasMSE
N = 100p (0.25)0.25950.00950.00260.2197−0.03030.0043
α (1.5)1.3903−0.10970.08451.50380.00380.0816
λ (0.5)0.4737−0.02630.01560.63790.13790.0382
β (2.0)2.01580.01580.04542.39420.39420.3605

N = 200p (0.25)0.25440.00440.00150.2379−0.01210.0016
α (1.5)1.3790−0.12100.09481.51960.01960.0763
λ (0.5)0.4737−0.02630.01280.55730.05730.0160
β (2.0)1.9872−0.01280.01942.11350.11350.0688

N = 300p (0.25)0.25330.00330.00090.2443−0.00570.0012
α (1.5)1.4113−0.08870.08721.58600.08600.0878
λ (0.5)0.4854−0.01460.01040.55470.05470.0119
β (2.0)1.9966−0.00340.01322.04830.04830.0191

N = 500p (0.25)0.2476−0.00240.00050.2463−0.00370.0005
α (1.5)1.4343−0.06570.07671.52090.02090.0824
λ (0.5)0.4844−0.01560.00670.52520.02520.0087
β (2.0)2.01760.01760.00582.03160.03160.0058

MSEs = means squared errors; OBSGcr = odd Birnbaum-Saunders geometric cure.

### Table 3

The mean, biases, and MSE of the OBSG mixture distribution for α = 2.0, λ = 0.5, β = 1.5, and p = 0.2 and n = 100, 200, 300, and 500

Sample sizeParameter (true value)Classical estimationBayesian estimation

MeanBiasMSEMeanBiasMSE
N = 100p (0.2)0.1960−0.00400.00190.1838−0.01620.0027
α (2.0)1.8903−0.10970.08372.27310.27310.1187
λ (0.5)0.4820−0.018 00.00740.65730.15730.0291
β (1.5)1.4378−0.06220.00631.4437−0.05630.0061

N = 200p (0.2)0.20340.00340.00120.1943−0.00570.0011
α (2.0)1.9346−0.06540.09462.19710.19710.0844
λ (0.5)0.4891−0.01090.00660.59260.09260.0122
β (1.5)1.4331−0.06690.00551.4352−0.06480.0053

N = 300p (0.2)0.20470.00470.00080.1987−0.00130.0008
α (2.0)1.9164−0.08360.09212.11510.11510.0820
λ (0.5)0.4914−0.00860.00620.55880.05880.0082
β (1.5)1.4342−0.06580.00511.4357−0.06430.0048

N = 500p (0.2)0.20240.00240.00040.1992−0.00080.0004
α (2.0)1.8891−0.11090.09122.0740.07400.0838
λ (0.5)0.4815−0.01850.00500.54250.04250.0067
β (1.5)1.4371−0.06290.00451.4370−0.06300.0044

MSE = means squared error; OBSG = odd Birnbaum-Saunders geometric.

### Table 4

MLEs of the model parameters for the ethylene data, the corresponding standard errors, and the AIC, CAIC, and BIC statistics

ModelEstimatesAICCAICBIC
KwN(μ,σ,a,b)−32.77 (2.55)29.40 (0.81)13.47 (1.42)0.45 (0.03)5775.15775.25792.9
BN(μ,σ,a,b)−56.17 (2.16)32.24 (0.96)50.93 (2.57)0.41 (0.02)5709.95710.05727.7
McN(μ,σ, a, b, c)−186.04 (7.92)47.99 (1.77)10021.00 (8.85)0.46 (0.03)4.63 (0.63)5638.35638.45660.5
KwOLLN(μ,σ, a, b, α)6.53 (2.77)113.18 (14.60)2.26 (0.33)0.27 (0.01)11.29 (2.31)5547.05547.15569.3
BS(λ, β)0.71 (0.02)27.47 (0.73)5424.05424.05432.9
OBS(α, λ, β)0.21 (0.01)0.28 (0.01)27.81 (0.87)5425.85425.95439.2
OBSG(α, λ, β, p)0.249 (0.013)0.269 (0.008)30.641 (0.920)0.309 (0.072)5254.85254.85272.5

MLEs = maximum likelihood estimators; AIC = Akaike information criterion; CAIC = consistent Akaike information criterion; BIC = Bayesian information criterion; KwN = Kumaraswamy normal; BN = beta normal; McN = McDonald normal; KwOLLN = Kumaraswamy odd log-logistic-normal; BS = Birnbaum-Saunders; OBS = odd Birnbaum-Saunders; OBSG = odd Birnbaum-Saunders geometric.

### Table 5

Likelihood ratio tests

EthyleneHypothesesStatisticwp-value
OBSG vs. OBSH0 : p = 0 vs. H1 : H0 is false173.0<0.0001
OBSG vs. BSH0 : α = 1 and p = 0 vs. H1 : H0 is false173.2<0.0001

OBSG = odd Birnbaum-Saunders geometric; OBS = odd Birnbaum-Saunders; BS = Birnbaum-Saunders.

### Table 6

Posterior mean, SD, and 95% CI of the parameters for the ethylene data and the DIC, EAIC, EBIC, and LPML statistics

ParameterMeanSD95% CIDICEAICEBICLPML

2.50%97.50%
OBSGα0.2670.0100.2490.2875220.9315227.2695245.052−2610.026
β31.8540.14131.48131.996
λ0.2660.0030.2620.274
p0.3050.0120.2760.320

OBSα0.2480.0020.2440.2505347.1095353.035366.368−2673.559
β27.9430.05727.78627.998
λ0.2870.0010.2860.290

BSβ27.5210.72426.11028.9625423.9645425.9925434.883−2711.748
λ0.7160.0200.6770.757

SD = standard deviation; CI = credible interval; DIC = deviance information criterion; EAIC = expected Akaike information Scriterion; EBIC = expected Bayesian information criterion; LPML = log pseudo marginal likelihood;OBSG = odd Birnbaum- Saunders geometric; OBS = odd Birnbaum-Saunders; BS = Birnbaum-Saunders.

### Table 7

MLEs of the model parameters for the melanoma data, the corresponding standard errors (given in parentheses), and the AIC, CAIC, and BIC statistics

ModelπαλβAICCAICBIC
OBSM0.4672 (0.0384)7.3671 (3.8938)5.9259 (1.1123)1.8335 (0.1705)1057.21057.31073.3
BSM0.4475 (0.0512)10.9378 (0.0998)1.8843 (0.2582)1062.51062.61074.6
pαλβ
OBSG0.5511 (0.0495)504.560 (29.2503)441.570 (33.4233)2.9234 (0.6284)1058.11058.21074.2
BSG0.6201 (0.0877)11.3107 (0.3799)4.6027 (2.6620)1069.51069.61081.6

MLEs = maximum likelihood estimators; AIC = Akaike information criterion; CAIC = consistent Akaike information criterion; BIC = Bayesian information criterion; OBSM = odd Birnbaum-Saunders mixture; BSM = Birnbaum-Saunders mixture; OBSG = odd Birnbaum-Saunders geometric; BSG = Birnbaum Saunders geometric.

### Table 8

Posterior mean, SD, and 95% CI of the parameters for the melanoma data and the DIC, EAIC, EBIC, and LPML statistics

ParameterMeanSD95% CIDICEAICEBICLPML

2.50%97.50%
OBSMπ0.4420.0570.3100.5261055.3651060.5971076.73−527.929
α7.0632.0762.5979.886
λ6.0571.8562.2759.253
β1.9890.3241.6052.758

BSMπ0.4500.0320.4020.5201061.4771065.0421077.141−531.093
λ0.9540.0810.8081.123
β1.9200.2061.5672.369

OBSGp0.5690.0320.5100.6331053.9031060.0631076.195−527.008
α491.03512.154470.983510.809
λ450.02726.685403.082497.458
β3.1870.3772.5303.911

BSGπ0.6150.0230.5690.662
λ1.3050.0731.1681.4511067.4511071.4991083.598−534.168
β4.5060.2884.0204.978

SD = standard deviation; CI = credible interval; DIC = deviance information criterion; EAIC = expected Akaike information criterion; EBIC = expected Bayesian information criterion; LPML = log pseudo marginal likelihood;OBSM= odd Birnbaum- Saunders mixture; BSM = Birnbaum-Saunders mixture; OBSG = odd Birnbaum-Saunders geometric; BSG = Birnbaum Saunders geometric.

### Table 9

Likelihood ratio tests

MelanomaHypothesesStatistic wp-value
OBSM vs. BSMH0 : α = 1 vs. H1 : H0 is false7.30.0069
OBSG vs. BSGH0 : α = 1 vs. H1 : H0 is false13.00.0003

OBSM = odd Birnbaum-Saunders mixture; BSM = Birnbaum-Saunders mixture; OBSG = odd Birnbaum-Saunders geometric; BSG = Birnbaum Saunders geometric.

References
1. Alexander, C, Cordeiro, GM, Ortega, EMM, and Sarabia, JM (2012). Generalized beta-generated distributions. Computational Statistics and Data Analysis. 56, 1880-1897.
2. Alizadeh, M, Emadi, M, Doostparast, M, Cordeiro, GM, Ortega, EMM, and Pescim, RR (2015). A new family of distributions: the Kumaraswamy odd log-logistic, properties and applications. Hacettepe Journal of Mathematics and Statistics. 44, 1491-1512.
3. Berkson, J, and Gage, RP (1952). Survival curve for cancer patients following treatment. Journal of the American Statistical Association. 47, 505-515.
4. Boag, JW (1949). Maximum likelihood estimates of the proportion of patients cured by cancer therapy. Journal of the Royal Statistical Statistical Society Series B (Methodological). 11, 15-53.
5. Brooks, SP (2002). Discussion on the paper by Spiegelhalter, Best, Carlin, and van der Linde. Journal of the Royal Statistical Society Series B (Methodological). 64, 616-618.
6. Cancho, VG, Ortega, EMM, Barriga, GDC, and Hashimoto, EM (2011). The Conway-Maxwell-Poisson-generalized gamma regression model with long-term survivors. Journal of Statistical Computation and Simulation. 81, 1461-1481.
7. Carlin, BP, and Louis, TA (2001). Bayes and Empirical Bayes Methods for Data Analysis. Boca Raton: Chapman and Hall
8. Castillo, NO, Gómez, HW, and Bolfarine, H (2011). Epsilon Birnbaum-Saunders distribution family: properties and inference. Statistical Papers. 52, 871-883.
9. Chib, S, and Greenberg, E (1995). Understanding the metropolis-Hastings algorithm. The American Statistician. 49, 327-335.
10. Cooner, F, Banerjee, S, Carlin, BP, and Sinha, D (2007). Flexible cure rate modeling under latent activation schemes. Journal of the American Statistical Association. 102, 560-572.
11. Cordeiro, GM, Cintra, RJ, Rego, LC, and Ortega, EMM (2012). The McDonald normal distribution. Pakistan Journal of Statistics and Operation Research. 8, 301-329.
12. Cordeiro, GM, and Lemonte, AJ (2011). The β-Birnbaum–Saunders distribution: an improved distribution for fatigue life modeling. Computational Statistics and Data Analysis. 55, 1445-1461.
13. Cordeiro, GM, and Lemonte, AJ (2014). The exponentiated generalized Birnbaum-Saunders distribution. Applied Mathematics and Computation. 247, 762-779.
14. Cordeiro, GM, Lemonte, AJ, and Ortega, EMM (2013). An extended fatigue life distribution. Statistics. 47, 626-653.
15. Cordeiro, GM, Ortega, EMM, and Nadarajah, S (2010). The Kumaraswamy Weibull distribution with application to failure data. Journal of the Franklin Institute. 347, 1399-1429.
16. Cowles, MK, and Carlin, BP (1996). Markov chain Monte Carlo convergence diagnostics: a comparative review. Journal of the American Statistical Association. 91, 883-904.
17. Fachini, JB, Ortega, EMM, and Cordeiro, GM (2014). A bivariate regression model with cure fraction. Journal of Statistical Computation and Simulation. 84, 1580-1595.
18. Farewell, VT (1982). The use of mixture models for the analysis of survival data with long-term survivors. Biometrics, 1041-1046.
19. Gelfand, AE, Dey, DK, and Chang, H (1992). Model determination using predictive distributions with implementation via sampling based methods (with discussion). Bayesian Statistics, Bernardo, JM, Berger, JO, Dawid, AP, and Smith, AFM, ed. New York: Oxford University Press, pp. 7-167
20. Hashimoto, EM, Ortega, EMM, Cordeiro, GM, and Cancho, VG (2014). The Poisson Birnbaum-Saunders model with long-term survivors. Statistics. 48, 1394-1413.
21. Hashimoto, EM, Ortega, EMM, Cordeiro, GM, and Cancho, VG (2015). A new long-term survival model with interval-censored data. Sankhya B. 77, 207-239.
22. Ibrahim, JG, Chen, MH, and Sinha, D (2001). Bayesian Survival Analysis. New York: Springer-Verlag
23. Kundu, D, Kannan, N, and Balakrishnan, N (2008). On the hazard function of Birnbaum–Saunders distribution and associated inference. Computational Statistics and Data Analysis. 52, 2692-2702.
24. Lanjoni, BR, Ortega, EMM, and Cordeiro, GM (2016). Extended Burr XII regression models: theory and applications. Journal of Agricultural Biological and Environmental Statistics. 21, 203-224.
25. Lemonte, AJ (2013). A new extension of the Birnbaum-Saunders distribution. Brazilian Journal of Probability and Statistics. 27, 133-149.
26. Longini, IM, and Halloran, ME (1996). A frailty mixture model for estimating vaccine efficacy. Journal of the Royal Statistical Society Series C Applied Statistics. 45, 165-173.
27. Maller, RA, and Zhou, X (1996). Survival Analysis with Long-Term Survivors. New York: Wiley
28. Ortega, EMM, Cancho, VG, and Paula, GA (2009). Generalized log-gamma regression models with cure fraction. Lifetime Data Analysis. 15, 79-106.
29. Ortega, EMM, Cordeiro, GM, Campelo, AK, Kattan, MW, and Cancho, VG (2015). A power series beta Weibull regression model for predicting breast carcinoma. Statistics in Medicine. 34, 1366-1388.
30. Ortega, EMM, Cordeiro, GM, Hashimoto, EM, and Suzuki, AK (2017). Regression models generated by gamma random variables with long-term survivors. Communications for Statistical Applications and Methods. 24, 43-65.
31. Ortega, EMM, Cordeiro, GM, and Lemonte, AJ (2012). A log-linear regression model for the β-Birnbaum-Saunders distribution with censored data. Computational Statistics Data Analysis. 56, 698-718.
32. Peng, Y, and Dear, KB (2000). A nonparametric mixture model for cure rate estimation. Biometrics. 56, 237-243.
33. Price, DL, and Manatunga, AK (2001). Modelling survival data with a cured fraction using frailty models. Statistics in Medicine. 20, 1515-1527.
34. Saavedra del Aguila, J, Heiffig-del Aguila, LS, Sasaki, FF, Tsumanuma, GM, das Graças Ongarelli, M, Spoto, MHF, Jacomino, AP, Ortega, EMM, and Kluge, RA (2010). Postharvest modifications of mechanically injured bananas. Revista Iberoamericana de Tecnologia Postcosecha. 10, 73-85.
35. Spiegelhalter, DJ, Best, NG, Carlin, BP, and van der Linde, A (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society Series B (Statistical Methodology). 64, 583-639.
36. Spiegelhalter, DJ, Thomas, A, Best, N, and Lunn, D (2007). OpenBUGS user manual version 3.0.2.Retrieved July 2, 2017, from http://mathstat.helsinki.fi/openbugs
37. Suzuki, AK, Cancho, VG, and Louzada, F (2016). The Poisson-Inverse-Gaussian regression model with cure rate: a Bayesian approach and its case influence diagnostics. Statistical Papers. 57, 133-159.
38. Sy, JP, and Taylor, JM (2000). Estimation in a Cox proportional hazards cure model. Biometrics. 56, 227-236.
39. Tsodikov, AD, Ibrahim, JG, and Yakovlev, AY (2003). Estimating cure rates from survival data: an alternative to two-component mixture models. Journal of the American Statistical Association. 98, 1063-1078.
40. Watson, GN (1995). A Treatise on the Theory of Bessel Functions. Cambridge: Cambridge University Press
41. Yu, B, and Peng, Y (2008). Mixture cure models for multivariate survival data. Computational Statistics and Data Analysis. 52, 1524-1532.