On a Class of Transmuted Distributions to Model Survival Data with a Cure Fraction

Abstract

Modelling time to event data, when there is always a proportion of the individuals, commonly referred to as immunes who do not experience the event of interest, is of importance in many biomedical studies. Improper distributions are used to model these situations and they are generally referred to as cure rate models. In the literature, two main families of cure rate models have been proposed, namely the mixture cure models and promotion time cure models. Here we propose a new model by extending the mixture model via a generating function by considering a shifted Bernoulli distribution. This gives rise to a new class of popular distributions called the transmuted class of distributions to model survival data with a cure fraction. The properties of the proposed model are investigated and parameters estimated. The Bayesian approach to the estimation of parameters is also adopted. The complexity of the likelihood function is handled through the Metropolis-Hasting algorithm. The proposed method is illustrated with few examples using different baseline distributions. A real life data set is also analysed.

AMS subject classifications: 62N02, 62F15

Keywords

Improper distributions Cure rate Metropolis–Hasting algorithm posterior distributions bootstrap

1 Introduction

The mixture model^{[1, 2]} and promotion time cure rate model^[3] are the two general cure rate models that have received a lot of attention and have been worked on by many authors under different perspectives. A unified approach to both these models is provided in Yin and Ibrahim,^[4] where the model has a link parameter. Particular values of this parameter yield the mixture cure rate model and the promotion time cure rate model. In both cases, the cure rate model has an improper distribution function. For the mixture cure rate model, the population consist of two groups. The immunes and the susceptible. If $T$ be the failure time, then the probability of survival of the total population, $\overset{̅}{G} (t) = P [T > t]$ for $0 < p < 1$ is given as

\overset{̅}{G} (t) = p + (1 - p) \overset{̅}{F_{0}} (t)

(1)

where

\overset{̅}{F_{0}} (t)

is any baseline survival function. Then

\overset{̅}{G} (t)

is an improper survival function and

lim_{t \to \infty} \overset{̅}{G} (t) \to p_{0}

, is referred to as the cure rate. Improper survival functions are popular as cure rate models. (See Rocha et al.^[5] and references therein). Given a binary random variable

M

such that

P [M = m] = p_{m}

m = 0, 1

any improper survival function can be represented as

\overset{̅}{G} (t) = E [\overset{̅}{F_{0}} (t)^{M}] = p_{0} + (1 - p_{0}) \overset{̅}{F_{0}} (t)

In Cox and Oakes,^[6] a competing risk perspective of this model is presented. Recently in Rodrigues et al.^{[7, 8]} this was further studied by considering $M$ to be an unobservable random variable denoting the number of competing causes related to the occurrence of an event of interest, with the probability mass function

p_{m} = P [M = m] m = 0, 1, 2, . . . .

Also, for $m > 0$ , $W_{k}$ , $k = 1, . ., m$ denote the time to event due to the $k th$ competing cause and $W_{k}$ are independent and identical random variables with survival function $\overset{̅}{F_{0}} (t)$ that does not depend on $M$ . The number of competing causes $M$ and the lifetime $W_{k}$ associated with particular cause is not observable but only the minimum lifetime among all causes is observed. It is defined as

T = \min (W_{1}, W_{2}, . . ., W_{m}) with P [W_{0} = \infty] = 1

(2)

which leads to a proportion

p_{0}

(say) of the population not susceptible to the event occurrence, termed the cure fraction. Then for a given

\overset{̅}{F_{0}} (t)

the long-term survival function of the random variable

T

is given as

\overset{̅}{G} (t) = P [T > t] = \sum_{m = 0}^{\infty} P [M = m] {\overset{̅}{F_{0}} (t)^{m}} = E [\overset{̅}{F_{0}} (t)^{M}]

(3)

where (3) is the probability mass function of

M

Rodrigues et al.^[7] studied the Model (3) by assuming a Weibull baseline distribution for $\overset{̅}{G_{0}} (t)$ and $M$ having a Conway-Maxwell (COM)- Poisson distribution. The survival function of COM-Poisson cure rate model is given as

\overset{̅}{G} (t) = \frac{Z (η \overset{̅}{F_{0}} (t), ϕ)}{Z (η, ϕ)}

(4)

where

Z (η, ϕ) = \sum_{j = 0}^{\infty} \frac{η^{j}}{(j!)^{ϕ}}

η = e^{θ}

and

- \infty < ϕ < \infty

In this paper, we propose a new model by extending the Rodrigues et al.^[7] model via a generating function by considering a shifted Bernoulli distribution for $M$ , given by

P (M = 1) = 1 - p_{1} and P (M = 2) = p_{1}, for 0 < p_{1} < 1

(5)

In Section 2 we propose the shifted Bernoulli Cure Rate (SBC) model. Several properties of the model along with identifiability conditions are presented. In Section 3, the proposed model is illustrated with generalized Gompertz distribution and Destructive Negative Binomial cure rate model.^[9] The properties of these examples and their estimation procedures of parameters are also discussed in Section 3. In Section 4, a Bayesian approach to the estimation of parameters of SBC is explored in detail. The complexity of the likelihood function is handled through the Metropolis- Hasting algorithm and simulation studies are carried out to establish the effectiveness of the Bayesian estimation procedure. In Section 5, an application to a real data set is detailed by using the well analysed melanoma data.^[10] We conclude the paper with a discussion and summary on our findings in Section 6.

2 Shifted Bernoulli Cure Rate Model

It is now assumed that $M,$ the unobservable random variable denoting the number of competing causes related to the occurrence of the event, represent one or two causes is a binary random variable with shifted Bernoulli distribution,

P (M = 1) = 1 - p_{1} and P (M = 2) = p_{1} .

(6)

The covariates are introduced through the parameter $p_{1}$ as $p_{1} (Z, β)$ , using the logistic function^[11] as

1 - p_{1} (Z, β) = \frac{1}{1 + e^{β Z}} and p_{1} (Z, β) = \frac{e^{β Z}}{1 + e^{β Z}}

(7)

where

Z = (1, z_{1}, . . ., z_{p})^{'}

be the vector of covariates and

β = (β_{0}, β_{1}, . . ., β_{p})

be the

(p + 1) \times 1

vector of regression coefficients. Define the class

P

such that

P = {p_{1} (Z, β) : p_{1} (Z, β) = \frac{e^{β Z}}{1 + e^{β Z}}} .

(8)

Let $\overset{̅}{G}$ denote the class of improper distribution $\overset{̅}{G_{0}} (t)$ , and be defined as

\overset{̅}{G_{0}} = {\overset{̅}{G_{0}} (t; θ_{0}) : lim_{t \to \infty} \overset{̅}{G_{0}} (t; θ_{0}) = c_{0} (θ_{0}), 0 < c_{0} (θ_{0}) < 1, θ_{0} \in Θ}

(9)

where the

k

dimensional vector,

θ_{0} = (θ_{01}, . . ., θ_{0 k})^{'}

, is the set of parameters in the parameter space

Θ

\overset{̅}{G_{0}} (\cdot)

. Then for a given baseline

\overset{̅}{G_{0}} (t; θ_{0})

\in

\overset{̅}{G}

, the long-time survival function of

T

with

M

having distribution as in (1) is denoted by

\overset{̅}{G} (t, Z; θ)

and is defined as

\overset{̅}{G} (t, Z; θ) = P [T > t] = E [\overset{̅}{G_{0}} (t)^{M}] = 1 - p_{1} (Z, β) \overset{̅}{G_{0}} (t; θ_{0}) + p_{1} (Z, β) {\overset{̅}{G_{0}}}^{2} (t; θ_{0}) .

(10)

where

θ = (θ_{0}, β)

. We refer to the distribution in (10) as the Shifted Bernoulli Cure Rate (SBC) model. For some baseline,

\overset{̅}{G_{0}} (t; θ_{0}) \in \overset{̅}{G_{0}}

, define the class

H

of survival functions of Shifted Bernoulli Cure Rate (SBC) model as

\begin{matrix} H = {\overset{̅}{G} (t, Z; θ) : \overset{̅}{G} (t, Z; θ) = (1 - p_{1} (Z, β)) \overset{̅}{G_{0}} (t; θ_{0}) + p_{1} (Z, β) ({\overset{̅}{G_{0}}}^{2} (t; θ_{0})), \\ p (\cdot) \in P, \overset{̅}{G_{0}} (\cdot) \in \overset{̅}{G_{0}}} \end{matrix}

(11)

This class of distributions fall into what is now popularly referred to as quadratic rank transmuted distributions which resurfaced with the work of Shaw and Buckley.^[12] These distributions are well researched upon. See Balakrishnan^[13] and Balakrishnan and He^[14] on the genesis of these distributions. The proposed class (11) provides a larger flexible class of cure rate models. The skewness of this distribution is managed through $p_{1}$ , which motivates us to include the effect of covariates through $p_{1}$ only. Here we do not assume any covariate effect on the baseline distribution.

2.1 Properties of SBC Model

We proceed this section by deriving the properties of the SBC Model (11). The density of the SBC model is obtained as

g (t, Z; θ) = g_{0} (t; θ_{0}) (1 + p_{1} (Z, β) - 2 p_{1} (Z, β) G_{0} (t; θ_{0}))

(12)

where

g (t, Z; θ)

and

g_{0} (t; θ_{0})

are the corresponding probability density functions of the survival function

\overset{̅}{G} (t, Z; θ) = 1 - G (t, Z; θ)

and baseline distribution

\overset{̅}{G_{0}} (t; θ_{0}) = 1 - G_{0} (t; θ)

respectively. The hazard rate of SBC model is given by

h (t, Z; θ) = h_{0} (t; θ_{0}) [\frac{1 + p_{1} (Z, β) - 2 p_{1} (Z, β) G_{0} (t; θ_{0})}{1 - p_{1} (Z, β) + p_{1} (Z, β) \overset{̅}{G_{0}} (t; θ_{0})}]

(13)

where

h (t, Z; θ)

and

h_{0} (t; θ_{0})

are hazard rates corresponding to the distribution in the class (11) and (9) respectively.

Theorem 2.1. The SBC model $\overset{̅}{G} (t, Z; θ)$ , has a decreasing failure rate if the baseline distribution $\overset{̅}{G_{0}} (t, θ_{0})$ has a decreasing failure rate.

Proof. Observe,

\begin{array}{l} h^{'} (t, Z; θ) = h_{0}^{'} (t; θ_{0}) (1 + \frac{p_{1} (Z, β) {\bar{G}}_{0} (t; θ)}{1 - p_{1} (Z, β) G_{0} (t; θ_{0})}) \\ + h_{0} (t; θ_{0}) (\frac{p_{1} (Z, β) g_{0} (t; θ_{0}) (p_{1} (Z, β) - 1)}{{(1 - p_{1} (Z, β) G_{0} (t; θ_{0}))}^{2}}) \end{array}

(14)

where

h^{'} (t, Z; θ)

and

h^{'}_{0} (t; θ_{0}

are the first derivatives of

h (t; Z, θ)

and

h_{0} (t; θ_{0})

respectively with respect to

t

. Let assume that

\overset{̅}{G_{0}} (t, θ_{0})

has a decreasing failure rate, then from (14), it is now straight forward that

\overset{̅}{G} (t, Z; θ)

has a decreasing failure rate.

Note that the converse of the Theorem 1 need not be true. For more on behaviour of hazard function in relation to the transmuted class see Bourguignon et al.^[15]

Theorem 2.2 The cure rate of the SBC model is a decreasing function in the mixing parameter $p_{1} (Z, β)$ .

Proof. It follows directly by observing that the cure rate of the SBC model is

c (θ | Z) = c_{0} (θ_{0}) (1 - p_{1} (Z, β) + p_{1} (Z, β) c_{0} (θ_{0}))

(15)

where

c_{0} (θ_{0})

as in (9). □

Remark 2.1. The cure rate of the SBC model is a function of covariates through the mixing parameter $p_{1}$ .

Remark 2.2. The cure rate need not uniquely determine the underlying distribution.

This prompts us to investigate the identifiability of the SBC model. We can notice that the SBC Model (10) can be represented as

U (t, Z; θ) = (1 - p_{1} (Z, β)) + p_{1} (Z, β) \overset{̅}{G_{0}} (t; θ_{0})

(16)

where

U (t, Z; θ) = \frac{\overset{̅}{G} (t, Z; θ)}{\overset{̅}{G_{0}} (t; θ_{0})}

. This is of the form of the standard mixture cure rate model. The identifiability of the standard mixture cure rate model has been extensively discussed.^{[16, 17]} Using their definition,

Definition 2.1. $U (t, Z, θ)$ is identifiable with in the family $P$ and $\overset{̅}{G}$ if the equality

U (t, Z, θ) = U^{*} (t, Z, θ^{*})

(17)

where

U^{*} (t, Z, θ^{*}) = \frac{\overset{̅}{G^{*}} (t, Z; θ^{*})}{\overset{̅}{G_{0}^{*}} (t; θ_{0}^{*})} = (1 - p_{1}^{*} (Z, β^{*})) + p_{1}^{*} (Z, β^{*}) \overset{̅}{G_{0}} (t; θ_{0}^{*})

(18)

for some $p_{1} (Z, β), p_{1}^{*} (Z, β^{*}) \in P$ , $\overset{̅}{G_{0}} (t; θ_{0}), \overset{̅}{G_{0}^{*}} (t; θ_{0}^{*}) \in \overset{̅}{G}$ implies that $p_{1} (Z, β) = p_{1}^{*} (Z, β^{*})$ and $\overset{̅}{G_{0}} (t; θ_{0}) = \overset{̅}{G_{0}^{*}} (t; θ_{0}^{*})$ .

From Theorem 2.3 of Hanin and Huang^[17], we note that $U (t, Z; θ_{0})$ is identifiable for any $p_{1} (Z, θ_{0}) \in P$ and for covariate vector with dimension greater than 2. Appealing to the result we have,

Theorem 2.3. The SBC model $\overset{̅}{G} (t, Z; θ)$ is identifiable if and only if the function $U (t, Z; θ)$ is identifiable.

Proof. If $U (t, Z, θ_{0})$ is identifiable, then

\begin{array}{l} U (t, Z; θ) = U^{*} (t, Z; θ^{*}) \\ \Rightarrow {\bar{G}}_{0} (t, θ_{0}) = G_{0}^{*} (t, θ_{0}^{*}) \end{array}

which in turn implies

\overset{̅}{G} (t, Z, θ) \overset{̅}{G^{*}} (t, Z, θ^{*})

. Hence the proof. □

2.2 The Likelihood Contributions

If $T_{i}$ and $C_{i}$ denote the failure time and censoring time of the $i th$ subject, we observe $Y_{i} = \min (T_{i}, C_{i})$ for $i = 1, 2, . ., n$ . Correspondingly, the failure-time is indicated by $δ = (δ_{1}, . . ., δ_{n})$ where $δ_{i} = I (T_{i} \geq C_{i})$ where $I (.)$ is the indicator function and $δ_{i}$ is equal to $1$ if $T_{i}$ is a time-to-event and is $0$ if it is right-censored. Let $Z = (1, z_{1}, . . ., z_{p})^{'}$ be the vector of covariates and $β = (β_{0}, β_{1}, . . ., β_{p})$ be the $(p + 1) \times 1$ vector of regression coefficient and the covariates are introduced through the parameter $p_{1}$ . If $θ_{0} = (θ_{01}, . . ., θ_{0 k})^{'}$ are the parameters of the baseline distribution, then $θ = (θ_{0}, β)$ forms the vector of parameters for the SBC model. The observed data are of the form $(t_{i}, δ_{i}, Z)$ . Suppose that the data are independently and identically distributed and come from SBC Model with density function $g (t_{i}, Z; θ)$ and survival function $\overset{̅}{G} (t_{i},, Z; θ) \in H$ , the likelihood function is given by

L (θ; t, Z, δ) \prod_{i = 1}^{n} g (t_{i}, Z; θ)^{δ_{i}} \overset{̅}{G} (t_{i}, Z; θ)^{(1 - δ_{i})} .

(19)

Note that this expression is valid for type I censoring, type II censoring, random censoring and when the censoring mechanism is not informative. For the proposed model in (11) the log-likelihood function is

\begin{matrix} l (θ; t, Z, δ) = \sum_{i = 1}^{n} & δ_{i} log (g_{0} (t; θ_{0}) (1 + p_{1} (Z, β) - 2 p_{1} (Z, β) G_{0} (t; θ_{0}))) \\ + (1 - δ_{i}) log ((1 - p_{1} (Z, β)) \overset{̅}{G_{0}} (t; θ_{0}) + p_{1} (Z, β) (\overset{̅}{G_{0}} (t; θ_{0}))^{2}) \end{matrix}

(20)

3 Examples

3.1 Shifted Bernoulli Gompertz Cure Rate Model

It is often observed that in some clinical data on cancer survival, the hazard rate decreases exponentially with time. This motivates us to consider the generalized Gompertz distribution as the baseline distribution in SBC model for our first example. It is specified as,

\overset{̅}{G_{0}} (t; θ_{0}) = e^{\frac{b}{a} (e^{- a t} - 1)}, a \neq 0, b > 0

(21)

where

θ_{0} = (a, b)

. If

a < 0

the survival function (3.1) corresponds to the proper Gompertz distribution and if

a > 0

the distribution is an improper distribution with surviving fraction to

e^{- \frac{b}{a}}

.^[18]

The survival function of Shifted Bernoulli Gompertz Cure Rate (SBGC) Model is

\overset{̅}{G} (t, Z; θ) = e^{\frac{b}{a} (e^{- at} - 1)} [1 - p_{1} (Z, β) + p_{1} (Z, β) e^{\frac{b}{a} (e^{- at} - 1)}]

(22)

where

θ = (a, b, β)

a > 0

b > 0

p_{1} (Z, β) = \frac{e^{β Z}}{1 + e^{β Z}}

Z = (1, z_{1}, . . ., z_{p})^{'}

and

β = (β_{0}, β_{1}, . . ., β_{p})

are the vector of covariates and the vector of regression coefficients respectively.

The cure fraction of SBGC Model is

c (θ | Z) = e^{- \frac{b}{a}} (1 - p_{1} (Z, β) + p_{1} (Z, β) e^{- \frac{b}{a}})

(23)

The probability density function and hazard functions are given by

g (t, Z; θ) = b e^{- at} e^{\frac{b}{a} (e^{- at} - 1)} (1 - p_{1} (Z, β) + 2 p_{1} (Z, β) e^{\frac{b}{a} (e^{- at} - 1)})

(24)

h (t, Z; θ) = {be}^{- at} [\frac{1 - p_{1} (Z, β) + 2 p_{1} (Z, β) e^{\frac{b}{a} (e^{- at} - 1)}}{1 - p_{1} (Z, β) + p_{1} (Z, β) e^{\frac{b}{a} (e^{- at} - 1)}}]

(25)

Figures 1(a) and 1(b) give the hazard rate of the SBGC model for different parametric values. From Theorem (1) and Figure 1, we can observe that $h (t)$ is monotonically decreasing with $t$ . Thus the SBC model is more suitable for the population, where more the time of survival, lesser the probability to experience the event. The distribution function of the SBGC model in (3.2) is given as

G (t, Z; θ) = [1 - e^{\frac{b}{a} (e^{- at} - 1)}] [1 + p_{1} (Z, β) e^{\frac{b}{a} (e^{- at} - 1)}]

(26)

Figure 1:

Source: The authors.

The quantile of the distribution of the SBGC Model (3.2) is derived as

G^{- 1} (u) = - \frac{1}{a} log [1 + \frac{a}{b} log (1 - \frac{1 + p_{1} (Z, β) - \sqrt{(1 + p_{1} (Z, β))^{2} - 4 p_{1} (Z, β) u}}{2 p_{1} (Z, β)})]

(27)

where

u \sim

Uniform

(0, 1)

3.1.1 Likelihood Contributions of SBGC Model

The log likelihood function corresponding to (2.15) for the SBGC is

\begin{matrix} l (θ; t, Z, δ) = & constant + log b \sum_{i = 1}^{n} δ_{i} - a \sum_{i = 1}^{n} δ_{i} t_{i} + \sum_{i = 1}^{n} \frac{b}{a} (e^{- {at}_{i}} - 1) \\ + \sum_{i = 1}^{n} δ_{i} log (1 - p_{1} (Z, β) + 2 p_{1} (Z, β) e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) \\ + \sum_{i = 1}^{n} (1 - δ_{i}) log (1 - p_{1} (Z, β) + p_{1} (Z, β) e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) \end{matrix}

(28)

where

θ = (a, b, β)

and

n

is the sample size.

3.1.2 Estimation

Given the observed data ( $t, δ$ , $Z$ ) we wish to find the value of $θ$ that maximizes (3.8). Here $θ = (a, b, β)$ . The likelihood equations can be obtained by solving first order partial derivatives of the likelihood and equating to zero. These equations are non-linear in nature and difficult to solve. We consider the profile likelihood method. In the first stage, we estimate $a$ by fixing $τ_{1} = (b, β)$ . Then re-substitute the estimate of $a$ from the first stage and estimate $τ_{1}$ in the second stage. The uni-modality of the likelihood function is observed in Figure 2. This process is continued until the convergence is met in both stages. This estimation method is an alternative one when closed form of MLE's is not possible and iterative procedures fail to converge. The asymptotic properties of the estimator $\hat{θ}$ are discussed extensively.^{[19, 20]} The asymptotic distribution of $\hat{θ}$ is multivariate normal with mean vector $θ$ and covariance matrix $\sum (\hat{θ})$ , which can be estimated by

\sum^{̂} (\hat{θ}) = {- \frac{\partial^{2} l (θ; t, Z, δ)}{\partial θ \partial θ^{T}}}^{- 1}

(29)

evaluated at

θ = \hat{θ}

from the observed Fisher information matrix.

Figure 2:

Source: The authors.

3.1.3 Simulation Studies

The algorithm to simulate a sample from SBGC model in (3.2) is detailed below.

Here we consider only one covariate. Generate $n$ samples of the covariate from $Z \sim N (0, 1)$ .

For fixed value of $θ = (β$ , $a, b)$ compute the cure fraction $c (θ | Z)$ as in (3.2).

Generate $M_{i}$ , where $M_{i} \sim$ Bernoulli $(1 - c (θ | Z))$ for each $i = 1, 2, . ., n$

Generate $U$ for Uniform $(0, 1 - c (θ | Z))$

For $M_{i} = 1$ , $t_{i} =^{- 1} G (u) =$ as in (3.7) and $M_{i} = 0$ , $t_{i}^{'} = \infty$ .

For each $i = 1, 2, . ., n$ generate $u_{i}^{'} \sim$ Uniform $(0, \max (t_{i}^{'}))$ , considering only the finite $t_{i}^{'}$ .

Calculate $t_{i} =$ min $(t_{i}^{'}, u_{i}^{'})$ . If $t_{i} < u_{i}^{'}$ set $δ_{i} = 1$ , otherwise set $δ_{i} = 0$ .

This gives us a right censored sample of size $n$ , from the SBGC model with covariate $Z$ .

A simulation study of the performance of maximum likelihood method of estimation was carried out based on simulated samples of

n = 25, 100

. The parameter values was fixed at

a = 3, b = 2 and p_{1} = 0.2

with censoring proportion of the population was equal to 0.46. The results for 1,000 replications in the absence and presence of the covariates are presented in Table 1 and 2 respectively. The results indicate that the maximum likelihood estimation works well. The biases and mean square error(MSE) decrease as the sample sizes increases.

Table 1:

Bias and MSE of Estimates in the Presence of Covariates.

		$n = 25$	$n = 100$
	$\hat{a}$	2.83	3.148
$a = 3$	Bias	$-$ 0.172	0.142
	MSE	1.699	0.680
	Coverage probability	0.756	0.893
	$\hat{b}$	1.882	2.025
$b = 2$	Bias	$-$ 0.118	0.025
	MSE	0.360	0.206
	Coverage probability	0.795	0.915
	$\hat{p}$	0.310	0.248
$p_{1} = 0.2$	Bias	0.110	0.048
	MSE	0.015	0.013
	Coverage probability	0.779	0.928
	$\hat{c}$	0.423	0.457
$c = 0.463$	Bias	$-$ 0.0401	$-$ 0.006
	MSE	0.024	0.006

Source: The authors.

Table 2:

Bias and MSE of Estimates in the Presence of Covariates.

		25	100
	$\hat{a}$	3.241	3.030
$a = 3$	Bias	0.241	0.030
	MSE	2.264	0.710
	Coverage probability	0.889	0.939
	$\hat{b}$	2.202	2.010
$b = 2$	Bias	0.202	0.010
	MSE	0.735	0.146
	Coverage probability	0.820	0.928
	$\hat{β_{0}}$	0.245	0.225
$β_{0} = 0.2$	Bias	0.026	0.045
	MSE	0.022	0.023
	Coverage probability	0.829	0.945
	$\hat{β_{1}}$	0.630	0.608
$β_{1} = 0.6$	Bias	0.029	0.008
	MSE	0.085	0.080
	Coverage probability	0.859	0.924

Source: The authors.

From Table 2 the cure rates of the SBGC model with sample size 25 and 100 are

c (θ | Z) = 0.507 (\frac{1}{1 + e^{0.245 + 0.630 z_{1}}} + 0.507 \frac{e^{0.245 + 0.630 z_{1}}}{1 + e^{0.245 + 0.630 z_{1}}})

(30)

and

c (θ | Z) = 0.515 (\frac{1}{1 + e^{0.225 + 0.608 z_{1}}} + 0.515 \frac{e^{0.225 + 0.608 z_{1}}}{1 + e^{0.225 + 0.608 z_{1}}})

(31)

respectively. Figure 3 gives the plots of the

c (θ | Z)

for

n = 25

and

n = 100

respectively.

Figure 3:

Source: The authors.

3.2 Shifted Bernoulli-Negative Binomial Cure Rate Model

The second example considers the Destructive Negative Binomial cure rate model^[9] as the baseline distribution. The model is transmuted using the mixing parameter $p_{1}$ . The resulting transmuted model is called Shifted Destructive Negative Binomial cure rate model(SBDNBC). The Destructive Negative Binomial Cure rate (DNBC) model given as

\overset{̅}{G_{0}} (t, θ_{0}) = {1 + ϕ η p F (t; θ_{01})}^{- \frac{1}{ϕ}}

(32)

where

F (t; θ_{01})

be the proper baseline distribution of the DNBC model with parameter

θ_{01}

, and

θ_{0} = (ϕ, η, p, θ_{01})

ϕ > 0, η > 0, 0 < p < 1

. Using (3.12) as the baseline in SBC we get, the survival function of Shifted Bernoulli- Destructive Negative Binomial Cure rate (SBDNBC) model as

\overset{̅}{G} (t, Z; θ) = {1 + ϕ η p F (t; θ_{01})}^{- \frac{1}{ϕ}} [1 - p_{1} (Z, β) + p_{1} (Z, β) {1 + ϕ η p F (t; θ_{01})}^{- \frac{1}{ϕ}}]

(33)

where

θ = (θ_{0}, β)

and

p_{1} (Z, β) = \frac{e^{β Z}}{1 + e^{β Z}}

The cure rate of SBDNBC model is

C (θ | Z) = (1 + ϕ η p)^{- 1 / ϕ} (1 - p_{1} (Z, β) + p_{1} (Z, β) (1 + ϕ η p)^{- 1 / ϕ})

(34)

where

0 < C (θ | Z) < 1

. The density of SBDNBC model is

g (t, Z; θ) = \frac{η p}{1 + ϕ η p F (t; θ_{01})} \overset{̅}{G} (t, Z; θ) f (t; θ_{01}) (1 + p_{1} (Z, β) - 2 p_{1} (Z, β) (1 - \overset{̅}{G} (t, Z; θ)))

(35)

where

\overset{̅}{G} (t, Z; θ)

is the survival function of SBDNBC model given in (3.13) and

f (t; θ_{01})

be the pdf of

F (t; θ_{01})

respectively. Note that in DNBC model we assume a Weibull distribution for the lifetime in (3.13) with

F (t; θ_{01}) = 1 - exp (- t^{γ_{1}} e^{γ_{2}})

and

f (t; θ_{01}) = γ_{1} t^{γ_{1} - 1} exp (γ_{2} - t^{γ_{1}} e^{γ_{2})}

for

γ_{1} > 0

γ_{2} \in R

and

θ_{01} = (γ_{1}, γ_{2})

. For more on this distribution we refer to Rodrigues et al.^[9]

Note Any cure rate model can be a baseline for SBC model. In particular any member of broad class of Destructive Weighted Poisson cure rate models^[9] can be used as the baseline in SBC.

3.2.1 Likelihood Contributions of SBDNBC Model

The log likelihood function corresponding to (2.15) for the SBDNBC model is

\begin{matrix} l (θ; t, Z, δ) = & \sum_{i = 1}^{n} {δ_{i} log (ϕ η) - δ_{i} log (1 + ϕ η p F (t_{i}; θ_{01})) + δ_{i} log f (t_{i}; θ_{01}) \\ - \frac{1}{ϕ} log (1 + ϕ η p F (t_{i}; θ_{01})) + δ_{i} log (1 - p_{1} (Z, β) + 2 p_{1} (Z, β) \\ {1 + ϕ η p F (t_{i}; θ_{01})}^{- \frac{1}{ϕ}}) + (1 - δ_{i}) log (1 - p_{1} (Z, β) + p_{1} (Z, β) \\ {1 + ϕ η p F (t_{i}; θ_{01})}^{- \frac{1}{ϕ}})} \end{matrix}

4 The Bayesian Estimators

The Bayesian approach to cure rate modelling allows one to combine prior information about the parameter with observed data and has been comprehensively presented by many authors. See Rodrigues et al.^[21] and Ibrahim et al.^[8] to mention a few.

The posterior distribution of $θ = (a, b, β)$ based on the observed data is given by

π (\hat{θ} | t) L (t | θ) π (a) π (b) π (β) .

(36)

where,

π (a), π (b)

and

π (β)

are the prior distribution of the parameters

a, b

and

β

respectively. With priors of

a, b, β

specified as

Gamma (m_{1}, n_{1})

Gamma (m_{2}, n_{2})

and

N (μ, σ)

respectively, the joint posterior distribution is obtained as

\begin{matrix} π (a, b, β | t, δ, Z) = & log b \sum_{i = 1}^{n} δ_{i} + \frac{b}{a} \sum_{i = 1}^{n} e^{- {at}_{i}} - \frac{n b}{a} + \sum_{i = 1}^{n} δ_{i} (1 + p_{1} - 2 p_{1} e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) \\ - a \sum_{i = 1}^{n} t_{i} δ_{i} + \sum_{i = 1}^{n} (1 - δ_{i}) (1 - p_{1} + p_{1} e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) - log Γ m_{1} \\ - m_{1} log n_{1} + (m_{1} - 1) log a - \frac{a}{n_{1}} - log Γ m_{2} - m_{2} log n_{2} \\ + (m_{2} - 1) log b - \frac{b}{n_{2}} + log (σ \sqrt{2 π}) - \frac{1}{2} {(\frac{β - μ}{σ})}^{2} \end{matrix}

(37)

where

p_{1} = p_{1} (Z, β) = \frac{e^{β Z}}{1 + e^{β Z}}

Posterior summaries of interest will be obtained from simulated samples for the joint posterior distribution using standard MCMC procedure as the posterior is analytically intractable. Gibbs sampling algorithm is used when the complete conditional distributions have known forms. The posterior conditionals of (4.2) for each variable is given below.

\begin{matrix} π (a | b, p_{1}; t) & \sum_{i = 1}^{n} {\frac{b}{a} (e^{- {at}_{i}} - 1) + δ_{i} (1 + p_{1} - 2 p_{1} e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) - a t_{i} δ_{i} \\ + (1 - δ_{i}) (1 - p_{1} + p_{1} e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) + (m_{1} - 1) log a - \frac{a}{n_{1}}} \\ π (b | a, p_{1}; t) & \sum_{i = 1}^{n} {log b δ_{i} + \frac{b}{a} (e^{- {at}_{i}} - 1) + δ_{i} (1 + p_{1} - 2 p_{1} e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) \\ + (1 - δ_{i}) (1 - p_{1} + p_{1} e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) + (m_{2} - 1) log b - \frac{b}{n_{2}}} \\ π (β | a, b; t) & \sum_{i = 1}^{n} {δ_{i} (1 + p_{1} - 2 p_{1} e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) + (1 - δ_{i}) (1 - p_{1} + p_{1} e^{\frac{b}{a} (e^{- {at}_{i}} - 1)}) \\ - \frac{1}{2} {(\frac{β - μ}{σ})}^{2}} \end{matrix}

(38)

The Gibbs sampling algorithm is ruled out as the conditionals are in form (4.3). We resort to the Metropolis-Hasting (M-H) algorithm which enables us to sample from the posterior distribution $π (θ | t)$ . The algorithm is based on a candidate-generating or proposal density $k (θ, ν)$ such that $\int k (θ, ν) d ν = 1$ . The M-H algorithm is

Set an initial value for $θ$ say $θ_{0}$ . Set iteration number $m = 0$

Generate $θ^{*}$ from the proposal density $k (θ_{m}, \cdot)$ . Draw a value $u$ from $U (0, 1)$ .

Calculate a probability of move $α (θ_{m}, θ^{*})$ . The probability of move from $θ$ to $θ^{*}$ is the ratio

α (θ, θ^{*}) = \min (\frac{π (θ^{*} | t) k (θ^{*}, θ)}{π (θ | t) k (θ, θ^{*})}, 1)

If $u \leq α (θ_{m}, θ^{*})$ , accept the new candidate, and set $θ_{m + 1} = θ^{*}$ . If $u > α (θ_{m}, θ^{*})$ , then do not accept the new candidate, and set $θ_{m + 1} = θ_{m}$ . Thus, at iteration $(m + 1)$ , $θ_{m + 1}$ either takes on the value $θ^{*}$ or remains at the previous iterate value $θ_{m}$ .

Set $m = m + 1$ , and repeat.

The iteration is continued till convergence is attained.

Convergence of the samples $θ_{m}$ for $m = 1, 2, . ., M$ to the invariant target posterior distribution occurs only after the Markov chain has passed the transient stage and the effect of the fixed starting value $θ_{0}$ has become negligible and can be ignored and occurs under mild regularity conditions. These regularity conditions are (a) irreducibility and (b) aperiodicity of the Markov chain. We use the package mhadaptive in R software

For a sample size $n$ , the mhadaptive take sample of size $m$ from posterior distribution $π (θ | t)$ . This process is replicated $R$ times. The estimate of $θ$ is $\hat{θ} = \frac{1}{R m} \sum_{i = 1}^{R} \sum_{j = 1}^{m} \hat{θ_{ij}}$ , where $\hat{θ_{ij}}$ be the $j th$ sample from the $i th$ replication and $R$ , $m$ are the total number of replications and the iterations respectively. The bias and MSE of the parameter $θ$ is $Bias = \hat{θ} - θ$ and $MSE = \frac{1}{R} \sum_{i = 1}^{R} (\frac{1}{m} \sum_{j = 1}^{m} \hat{θ_{} ij} - θ)^{2}$ respectively. Posterior summaries of interest will be obtained by using R software.^[22] To compare the fit of different model formulations, Bayesian information criterion defined by $BIC = 2 L (t, Z | θ) + q log (n)$ is used, where $q$ and $n$ are the total number of model parameters and the sample size respectively.

Simulation studies were carried out for various sample sizes including and excluding the covariates. The estimates were obtained for both the informative and flat reference prior. The results from the simulation studies are presented in Tables (3)–(5). The study shows that the biases and MSE's decreases as the sample sizes increases. Also the posterior is not prior sensitive. The prior robustness is evident from Table 3 and 4.

Table 3:

Bayes’ Estimates of $?$ with Reference Prior.

		25	50	100	200	500
True parameter value of	$\hat{a}$	2.567	2.634	2.786	2.899	2.978
$a = 3$	Bias	$-$ 0.432	$-$ 0.364	$-$ 0.215	$-$ 0.101	$-$ 0.0213
	MSE	0.679	0.719	0.515	0.263	0.087
True parameter value of	$\hat{b}$	1.750	1.849	1.858	1.884	1.935
$b = 2$	Bias	$-$ 0.250	$-$ 0.150	$-$ 0.142	$-$ 0.116	$-$ 0.065
	MSE	0.239	0.226	0.142	0.080	0.039
True parameter value of	$\hat{p_{1}}$	0.253	0.251	0.250	0.251	0.248
$p_{1} = 0.2$	Bias	0.531 $\times 10^{- 1}$	0.507 $\times 10^{- 1}$	0.509 $\times 10^{- 1}$	0.505 $\times 10^{- 1}$	0.486 $\times 10^{- 1}$
	MSE	0.507 $\times 10^{- 2}$	0.484 $\times 10^{- 2}$	0.477 $\times 10^{- 2}$	0.494 $\times 10^{- 2}$	0.504 $\times 10^{- 2}$
True parameter value of	$\hat{c}$	0.421	0.414	0.438	0.456	0.461
$c = 0.46$	Bias	$-$ 0.429 $\times 10^{- 1}$	$-$ 0.498 $\times 10^{- 1}$	$-$ 0.251 $\times 10^{- 1}$	$-$ 0.079 $\times 10^{- 1}$	$-$ 0.259 $\times 10^{- 2}$
	MSE	0.134 $\times 10^{- 1}$	0.121 $\times 10^{- 1}$	0.586 $\times 10^{- 2}$	0.240 $\times 10^{- 2}$	0.725 $\times 10^{- 3}$

Source: The authors.

Table 4:

Baye's Estimates of $θ$ with Reference Prior.

		25	50	100	200	500
True parameter value of	$\hat{a}$	2.370	2.532	2.752	2.970	3.057
$a = 3$	Bias	$-$ 0.630	$-$ 0.468	$-$ 0.247	$-$ 0.029	0.057
	MSE	0.597	0.561	0.294	0.190	0.086
True parameter value of	$\hat{b}$	1.806	1.905	1.906	1.968	2.006
$b = 2$	Bias	$-$ 0.194	$-$ 0.094	$-$ 0.094	$-$ 0.031	0.006
	MSE	0.375	0.141	0.082	0.049	0.035
True parameter value of	$\hat{p_{1}}$	0.324	0.243	0.236	0.237	0.225
$p_{1} = 0.2$	Bias	0.124	0.043	0.036	0.037	0.024
	MSE	0.155 $\times 10^{- 1}$	0.429 $\times 10^{- 2}$	0.358 $\times 10^{- 2}$	0.344 $\times 10^{- 2}$	0.365 $\times 10^{- 2}$
True parameter value of	$\hat{c}$	0.372	0.394	0.432	0.453	0.463
$c = 0.46$	Bias	$-$ 0.910 $\times 10^{- 1}$	$-$ 0.719 $\times 10^{- 1}$	$-$ 0.318 $\times 10^{- 1}$	$-$ 0.998 $\times 10^{- 2}$	0.011 $\times 10^{- 2}$
	MSE	0.829 $\times 10^{- 2}$	0.791 $\times 10^{- 2}$	0.633 $\times 10^{- 2}$	0.241 $\times 10^{- 2}$	0.068 $\times 10^{- 2}$

Source: The authors.

Table 5:

Baye's Estimates of $θ$ with Informative Prior.

		25	50	100	200	500
True parameter value of	$\hat{a}$	2.704	2.625	2.768	2.887	2.923
$a = 3$	Bias	$-$ 0.296	$-$ 0.375	$-$ 0.231	$-$ 0.113	$-$ 0.075
	MSE	0.932	0.761	0.518	0.254	0.086
True parameter value of	$\hat{b}$	1.839	1.866	1.928	1.992	1.992
$b = 2$	Bias	$-$ 0.161	$-$ 0.134	$-$ 0.072	$-$ 0.008	$-$ 0.007
	MSE	0.217	0.203	0.130	0.068	0.029
True parameter value of	$\hat{β_{0}}$	0.150	0.150	0.172	0.192	0.197
$β_{0} = 0.2$	Bias	$-$ 0.499 $\times 10^{- 1}$	$-$ 0.479 $\times 10^{- 1}$	$-$ 0.200 $\times 10^{- 1}$	$-$ 0.0800 $\times 10^{- 1}$	$-$ 0.0300 $\times 10^{- 1}$
	MSE	0.250 $\times 10^{- 2}$	0.150 $\times 10^{- 2}$	0.145 $\times 10^{- 2}$	0.050 $\times 10^{- 2}$	0.0025 $\times 10^{- 2}$
True parameter value of	$\hat{β_{1}}$	0.539	0.548	0.567	0.578	0.585
$β_{1} = 0.6$	Bias	$-$ 0.605 $\times 10^{- 1}$	$-$ 0.511 $\times 10^{- 1}$	$-$ 0.328 $\times 10^{- 1}$	$-$ 0.211 $\times 10^{- 1}$	$-$ 0.137 $\times 10^{- 1}$
	MSE	0.402 $\times 10^{- 2}$	0.322 $\times 10^{- 2}$	0.141 $\times 10^{- 2}$	0.094 $\times 10^{- 2}$	0.037 $\times 10^{- 2}$

Source: The authors.

5 An Application

The data from a study that included 205 patients diagnosed and treated for Malignant Melanoma between 1962 and 1997 is studied in Andersen et al.^[10] Each patient had their tumour removed by surgery at the Department of Plastic Surgery, University Hospital of Odense, Denmark. The surgery consisted of complete removal of the tumour together with about 2.5cm of the surrounding skin. Data is available in the timereg package in R.^[23] The observed time (T) ranges from 10 to 5,565 days (0.024 to 15.25 years) and refers to the time until the patient’s death or censoring time. The patients still alive at the end of the study are considered as censored observations ( $65.36$ per cent). This is a well analysed data discussed in detail in Rodrigues et al.^[9]

We fit the data for some existing cure rate models along with the SBC and SBDNBC model. The maximum likelihood estimates and the corresponding standard errors (SE), lower and upper confidence intervals (LCL, UCL) are presented in Table 6. The LCL and UCL of the maximum likelihood estimates of the parameter

θ

\hat{θ} - {SE}_{\hat{θ}} z_{α / 2}

and

\hat{θ} + {SE}_{\hat{θ}} z_{α / 2}

respectively, where

\hat{θ}

is the maximum likelihood estimate and

{SE}_{\hat{θ}}

is the standard error of the parameter

θ

. The Standard error of the parameter

θ

is the square root of the variance, which is obtained from the Fisher information matrix (3.9) by plugging the estimated parameter values of

θ

. The model comparison by the AIC criteria evidences the SBDNBC model to be a better fit than few popular cure rate models. (refer Table 6). Also Figure 4 suggest that the SBDNBC model is able to capture the behaviour of Kaplan–Meier curve relatively better compared to the other models. These findings are consistent with the AIC value presented in Table 6.

Table 6:

Maximum Likelihood Estimation for Different Cure Rate Model

Model	Parameter	Estimates	S.E	LCL	UCL	AIC
	$λ$	2.806	1.819	$-$ 0.758	6.374	474.1
Generalized Exponential cure rate model^[26]	$α$	0.316	0.154	$-$ 0.002	0.602
	$p$	0.693	0.049	0.597	0.792
	$a$	$-$ 0.247	0.064	$-$ 0.373	$-$ 0.120	573.6
Standard Mixture Gompertz cure rate model	$b$	0.104	0.023	0.057	0.150
	$p$	0.751	0.031	0.687	0.812
Defective Gompertz cure rate model	a	0.073	0.051	$-$ 0.031	0.1716	495.5
	b	0.045	0.006	0.026	0.053
	a	0.089	0.064	$-$ 0.038	0.216	631.0
Shifted Bernoulli Gompertz cure rate model	b	0.022	0.003	0.015	0.289
	$p_{1}$	0.279	0.022	$-$ 0.152	0.710
	$γ_{1}$	3.012	0.500	2.028	3.991	464.2
Destructive Negative Binomial cure rate model^[9]	$γ_{1}$	$-$ 7.845	0.815	$-$ 9.447	$-$ 6.252
	$η$	3.961	0.729	2.530	5.389
	$p$	0.551	0.101	0.351	0.748
	$ϕ$	6.425	2.931	0.688	12.167
	$γ_{1}$	3.668	0.439	2.807	4.528	420.9
Shifted Bernouli Destructive	$γ_{2}$	$-$ 6.871	0.870	$-$ 7.73144	$-$ 5.165
Negative Binomial Cure Rate model	$η$	4.839	1.402	2.091	7.586
	$p$	0.251	0.072	0.109	0.392
	$ϕ$	8.121	3.138	1.970	14.271
	$p_{1}$	0.369	0.085	0.202	0.535

Source: The authors.

Figure 4:

Source: The authors.

To illustrate the effect of covariates on the cure rate, we consider tumour thickness (in mm) (

z_{1}

), Gender (Male and Female)

(z_{2})

and ulceration status (present and absent)(

z_{3}

). The covariate Gender (

z_{2}

) did not show a significant effect on the long-term survival; hence we eliminate the covariate

z_{2}

from further analysis. We fit the data in DNBC model,^[8] and SBDNBC model with covariate

z_{1}

and

z_{3}

since these two models were stand out the best ones in analysing data with covariates excluded (see Table 6). Result from the analysis as evidenced in Table 6 indicate that the SBDNBC model fits the data best. This algorithm is cross checked using bootstrap technique. Here we draw

B

bootstrap samples with replacement from the melanoma data and each bootstrap sample consist of 205 observations. We obtain the the bootstrap estimate for each bootstrap sample, say

\hat{θ_{b}}

, where

\hat{θ_{b}}

be the

b th

bootstrap sample estimate of the parameter

θ

b = 1, 2, \dots B

.^[24] The estimated values are presented in Table 7 thereby confirming our claim. The bootstrap technique is used to check the accuracy of algorithm.

Table 7:

Bootstrap Estimates to Check Accuracy of Algorithm

Model	Parameter	Estimates	S.E	LCL	UCL	AIC
	$λ$	2.149	0.557	1.839	3.061	482.4
Generalized Exponential cure rate model^[26]	$α$	0.316	0.110	0.201	0.497
	$p$	0.596	0.216	0.206	0.779
	$a$	$-$ 0.125	0.095	$-$ 0.281	0.106	626.5
Standard Mixture Gompertz cure rate model	$b$	0.067	0.037	0.0182	0.115
	$p$	0.548	0.150	0.437	0.677
Defective Gompertz cure rate model	a	0.056	0.010	0.041	0.073	514.2
	b	0.043	0.031	0.014	0.101
	a	0.054	0.033	0.018	0.098	671.7
Shifted Bernoulli Gompertz cure rate model	b	0.036	0.009	0.023	0.046
	$p_{1}$	0.321	0.073	0.241	0.398
	$γ_{1}$	3.012	0.317	3.341	3.987	477.3
Destructive Negative Binomial cure rate model^[9]	$γ_{1}$	$-$ 7.876	0.317	$-$ 8.069	$-$ 7.093
	$η$	3.252	0.531	2.935	4.047
	$p$	0.718	0.318	0.497	0.938
	$ϕ$	7.350	2.014	6.925	11.973
	$γ_{1}$	3.595	0.200	3.203	3.987	457.0
Shifted Bernoulli Destructive	$γ_{2}$	$-$ 6.389	0.205	$-$ 6.7908	$-$ 5.987
Negative Binomial Cure Rate model	$η$	3.242	1.038	1.207	5.276
	$p$	0.213	0.034	0.146	0.279
	$ϕ$	8.416	2.200	4.104	12.728
	$p_{1}$	0.311	0.054	0.205	0.4168

Source: The authors.

We used the Bayesian inference procedure, and considered Gamma priors for the parameters ( $γ_{1}, ϕ, η$ ), Normal priors for ( $γ_{2}, β_{0}, β_{1}, β_{2 pre}, β_{2 pre}$ ) and uniform prior for the parameter $p$ respectively. Using M.H algorithm, we sample from the posterior distribution and the posterior summaries of interest will be obtained by using R software. The plots of marginal posterior densities are presented in Figure 5 and 6 respectively. The Bayesian estimates for the parameters are given in Table 8, where $β_{2_{abs}}$ and $β_{2_{pre}}$ are the regression coefficients corresponds to ulceration status ‘absent’ and ‘present’ respectively. The standard error of the Bayes estimator is calculated using the Bootstrap algorithm.^[25] The corresponding SE, Higher posterior density (HPD) intervals and BIC value for DNBC and SBDNBC models are presented in Table 8. We observe that the SBDNBC model shows better fit to the data. The BIC values corresponds to the SBDNBC model shows a marked improvement compared to DNBC.^[9] This analysis further reaffirms our claim. Figure 7 displays the combined effect of ulceration status and tumour thickness on the cure rate. We observe that the cure rate of patients who have ulceration status is ‘absent’ is higher than compared with the patients who have ulceration status is present. The cure rate of patients who have tumour thickness greater than 5mm, decreases and reach a plateau.

Table 8:

Bayesian Estimates for Different Cure Rate Models

Model	Parameter	Estimates	S.E	$95 %$ HPD Intervals	K.S Test Statistic	P-Value	BIC
Destructive	$γ_{1}$	3.963	0.030	(3.012, 4.045)
Negative Binomial	$γ_{2}$	$-$ 8.129	0.034	( $-$ 9.207, $-$ 7.970)
cure rate model	$β_{0}$	$-$ 7.012	0.210	( $-$ 7.303, $-$ 6.976)
	$β_{1}$	1.937	0.219	(0.217, 2.051)	0.012	0.781	446.1
	$β_{2_{pre}}$	5.021	0.331	(4.797, 5.303)
	$β_{2_{abs}}$	6.995	0.122	(6.681, 7.012)
	$ϕ$	8.877	0.451	(7.842, 9.951)
Shifted Bernoulli Destructive	$γ_{1}$	1.809	0.443	(0.134, 2.025)
Negative Binomial	$γ_{2}$	$-$ 7.950	0.0394	( $-$ 8.085, $-$ 7.390)
cure rate model	$β_{0}$	$-$ 7.237	0.557	( $-$ 9.067, $-$ 6.981)
	$β_{1}$	2.076	0.116	(1.814, 3.029)
	$β_{2_{pre}}$	4.945	0.060	(4.326, 5.185)	0.321	0.801	412.0
	$β_{2_{abs}}$	6.92	0.135	(5.991, 7.136)
	$ϕ$	8.889	0.189	(7.795, 9.020)
	$η$	3.847	0.290	(2.422, 4.048)
	$p$	0.594	0.060	(0.257, 0.997)

Source: The authors.

Figure 5:

Source: The authors.

Figure 6:

Source: The authors.

Figure 7:

Source: The authors.

6 Discussion and Summary

We propose a general class of cure rate models motivated by rank transmutation of improper distribution functions. The proposed model called Shifted Bernoulli Cure Rate (SBC) model and the identifiability of the model established by appealing to the result in Li et al.^[16] We have assumed that the baseline distribution is independent of covariate effect and the effect of covariates are modelled through the transmuted parameter $p_{1}$ , which manages the skewness of the distribution. The examples of the proposed model is illustrated through the Gompertz and the DNBC^[9] as baseline distributions to give the SBGC and SBDNBC models respectively. The unimodality of the likelihood enables us to use the method of maximum likelihood estimation. The simulation studies show that the MLE is an effective method. Under additional prior information, the Bayesian approach is shown to give more accurate estimates. The complexity of the posterior is handled through the M H algorithm. The Melanoma data is used to illustrate our results. The data includes 205 patients observed after operation for removal of malignant melanoma in the period 1962–1977. The surviving probability of melanoma patient is modelled through existing cure rate models and proposed models. To study the effect of covariates on the long term survival probability of melanoma patients, the parameter $p_{1} (Z, β)$ is linked with the covariates gender, ulceration status and tumour thickness. The covariate gender showed no significant effect on the surviving probability while ulceration status and tumour thickness are substantial. The t-test statistic and the p-value for the regression coefficient corresponding to gender of melanoma patients is obtained as $-$ 0.38775 and 0.4986. This confirms that the long-term surviving probability of melanoma patients is not affected by patient's gender. The results in Table 6 suggests that the DNBC shows a better fit than SBGC model consistent with the finding in Rodrigues et al.^[9] However the SBDNBC gives a marked improvement over DNBC. Our algorithm is verified by the bootstrap technique. From Table 8 the cure rate of the melanoma population is observed as

c (θ) = 0.708 [1 - p_{1} (β, Z) + 0.708 p_{1} (β, Z)]

(39)

where

p_{1} (β, Z) = \frac{e^{- 7.237 + 1.771 z_{1} + 4.945 z_{pre} + 6.992 z_{abs}}}{1 + e^{- 7.237 + 1.771 z_{1} + 4.945 z_{pre} + 6.992 z_{abs}}}

, and

z_{pre}

and

z_{abs}

are ulceration status of melanoma patients defined as

z_{pre} = 1

if ulceration is present and

0

otherwise, similarly

z_{abs} = 1

if ulceration is absent and

0

otherwise. The cure rate (6.1) implies that, higher values of tumour thickness imply smaller estimates for cure rate. If the tumour thickness of a melanoma patient increases 1.771mm then the probability of cure decreases at the rate of 0.01. To conclude, the proposed SBDNBC model shows a better fit compared to the other existing cure rate models for this data.

Footnotes

Acknowledgements

The authors thank Professor N Balakrishnan, McMaster University, for the discussions on many aspects of this article which led to an improvised version. We thank the reviewers for their constructive comments for the earlier version of the article.

Declaration of Conflicting Interests

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors acknowledges the financial assistance from Cochin University of Science & Technology through its UJRF scheme.

References

Berkson

and Gage

Survival curve for cancer patients following treatment. J Am Stat Assoc . 1952; 47(259): 501–515.

Boag

JW.

Maximum likelihood estimates of the proportion of patients cured by cancer therapy. J R Stat Soc . 1949; 11(1): 15–53.

Yakovlev

Asselain

Bardou

Fourquet

Hoang

, Rochefodiere

and Tsodikov

. A simple stochastic model of tumor recurrence and its application to data on premenopausal breast cancer. In: Asselain B, Boniface M, Duby C, Lopez C, Masson JP and Tranchefort J (eds) Biometrie et Analyse de Dormees Spatio-Temporelles . Rennes: ENSA; 1993, pp. 66–82.

Yin

Ibrahim

JG.

Cure rate models: A unified approach. Can J Stat . 2005; 33(4): 559–570.

Rocha

Nadarajah

Tomazella

Louzada

A new class of defective models based on the Marshall–Olkin family of distributions for cure rate modeling. Comput Stat Data Anal . 2017; 107: 48–63.

Cox

, Oakes

Analysis of survival data . London: Chapman and Hall; 1984.

Rodrigues

, de Castro

, Cancho

, Balakrishnan

COM–Poisson cure rate survival models and an application to a cutaneous melanoma data. J Stat Plan Inference . 2009; 139(10): 3605–3611.

Rodrigues

, Cancho

, de Castro

, Balakrishnan

. A Bayesian destructive weighted Poisson cure rate model and an application to a cutaneous melanoma data. Stat Methods Med Res . 2012; 21(6): 585–597.

Rodrigues

, de Castro

, Balakrishnan

, Cancho

. Destructive weighted Poisson cure rate models. Lifetime Data Anal . 2011; 17(3): 333–346.

10.

Andersen

, Borgan

, Gill

, Keiding

Statistical models based on counting processes . New York, NY: Springer Science & Business Media; 2012.

11.

Farewell

VT.

The use of mixture models for the analysis of survival data with longterm survivors. Biometrics . 1982; 38(4): 1041–1046.

12.

Shaw

Buckley

IR.

The alchemy of probability distributions: Beyond Gram-Charlier expansions, and a Skew-Kurtotic-normal distribution from a rank transmutation map. arXiv https://arxiv.org/abs/0901.0434; 2009.

13.

Balakrishnan

My musings on a pioneering work of Erich Lehmann and its rediscoveries. Commun Stat . 2020. https://doi.org/10.1080/03610926.2021.1896003

14.

Balakrishnan

A record-based transmuted family of distributions. In: Indranil Ghosh I, Balakrishnan N and Ng HKT (eds) Advances in statistics–Theory and applications . New York, NY: Springer; 2021.

15.

Bourguignon

Ghosh

, Cordeiro

. general results for the transmuted family of distributions and new models. J Probab Stat . 2016. https://doi.org/10.1155/2016/7208425

16.

Taylor

JP.

Identifiability of cure models. Stat Probab Lett . 2001; 54(4): 389–395.

17.

Hanin

Huang

LS.

Identifiability of cure models revisited. J Multivar Anal . 2014; 130: 261–274.

18.

Cantor

AB.

Projecting the standard error of the Kaplan–Meier estimator. Stat Med . 2001; 20(14): 2091–2097.

19.

Patefield

On the maximized likelihood function. Sankhya: Indian J Stat , Series B. 1977; 39(1): 92–96.

20.

Cox

, Barndorff-Nielsen

. Inference and asymptotics (vol. 52). Boca Raton, FL: CRC Press; 1994.

21.

Ibrahim

Chen

Sinha

Bayesian semiparametric models for survival data with a cure fraction. Biometrics . 2001; 57(2): 383–388.

22.

Chivers

, Chivers

. Package ‘MHadaptive’ . 2011. https://cran.r-project.org/web/packages/MHadaptive/MHadaptive.pdf

23.

Scheike

Package ‘timereg’ . 2009. http://www2.uaem.mx/r-mirror/web/packages/timereg/timereg.pdf

24.

Efron

, Hastie

Computer age statistical inference (vol. 5). Cambridge: Cambridge University Press; 2016.

25.

Lee

Carroll

Sinha

Frequentist standard errors of Bayes estimators. Comput Stat 2017; 32(3): 867–888.

26.

Kannan

Kundu

Nair

Tripathi

The generalized exponential cure rate model with covariates. J Appl Stat . 2010; 37(10): 1625–1636.

On a Class of Transmuted Distributions to Model Survival Data with a Cure Fraction

Abstract

Keywords

1 Introduction

3.1 Shifted Bernoulli Gompertz Cure Rate Model

Table 1:

Bias and MSE of Estimates in the Presence of Covariates.

Bias and MSE of Estimates in the Presence of Covariates.

4 The Bayesian Estimators

Bayes’ Estimates of ? with Reference Prior.

Baye's Estimates of θ with Reference Prior.

Baye's Estimates of θ with Informative Prior.

Table 6:

Maximum Likelihood Estimation for Different Cure Rate Model

Bootstrap Estimates to Check Accuracy of Algorithm

Bayesian Estimates for Different Cure Rate Models

Footnotes

Acknowledgements

Declaration of Conflicting Interests

Funding

References

Bayes’ Estimates of $?$ with Reference Prior.

Baye's Estimates of $θ$ with Reference Prior.

Baye's Estimates of $θ$ with Informative Prior.