The Slash Half-Normal Distribution Applied to a Cure Rate Model with Application to Bone Marrow Transplantation

Abstract

This paper proposes, for the first time, the use of an asymmetric positive and heavy-tailed distribution in a cure rate model context. In particular, it introduces a cure-rate survival model by assuming that the time-to-event of interest follows a slash half-normal distribution and that the number of competing causes of the event of interest follows a power series distribution, which defines six new cure rate models. Several properties of the model are derived and an alternative expression for the cumulative distribution function of the model is presented, which is very useful for the computational implementation of the model. A procedure based on the expectation–maximization algorithm is proposed for the parameter estimation. Two simulation studies are performed to assess some properties of the estimators, showing the good performance of the proposed estimators in finite samples. Finally, an application to a bone marrow transplant data set is presented.

1. Introduction and Motivation

Survival models that incorporate a proportion of cured subjects are becoming increasingly popular in the analysis of data from clinical trials and other studies, mainly due to advances in medicine. In a cure rate model context, specifically using a competitive risk model framework, it is assumed that the event of interest is produced due to $M\in {\mathbb{N}}_0$ carcinogenic cells (a latent variable), each of them with an associated survival function $S(·;\mathbf{\lambda })$, where $\mathbf{\lambda }$ is a vector of unknown parameters from a distribution defined on the positive line, and $M=0$ represents the cure individuals. Different models are assumed for the number of cells (M) and the survival function $S(·;\mathbf{\lambda })$. Specifically for the latter, we can mention Weibull, log-normal, generalized gamma, Birnbaum–Saunders distributions, among others. In particular, we focus on the proposed extension in Olmos et al. [1]. This model (parsimonious, since it only has two parameters) is called a slash half-normal (SHN) distribution. This distribution with a heavy tail has not yet been used in a cure rate model context. The objective of the proposed model is to increase the kurtosis with respect to its half-normal basal distribution, being more useful for modeling data sets that may have heavy tails. A version of the SHN reparameterized in terms of the mean with covariates was presented in Gómez et al. [2].

On the other hand, bone marrow transplantation is an important treatment for more than 50 fatal diseases such as leukemia, lymphoma, solid organ neoplasia, autoimmune diseases, among others. Bone marrow transplantation has unique biological characteristics and associated problems that markedly differ from solid organ transplantation, such as kidney, liver, or heart (Hardy et al. [3]). The problems are associated with the fact that immunocompetent cells can reject each other, which can cause graft rejection and generate graft-versus-host disease (GVHD). Furthermore, successful engraftment requires a strict match with the major histocompatibility complex (MHC) class II antigen of the donor and the recipient. Added to this is that patients must receive prior treatment with cytotoxic and myeloablative agents to avoid graft rejection (Chinen and Buckley [4], Simpson and Dazzi [5]). For these reasons, bone marrow transplantation is associated with considerable morbidity and mortality.

From a biological point of view, it makes sense to assume a heavy-tailed distribution when following patients over a period of time and following patient survival to assess transplant efficacy. These models can describe in detail the progression of the disease and allow us to make predictions for future clinical trials.

The motivation to introduce the SHN distribution in a cure rate model is given by a real data set application from the European Society for Blood and Marrow Transplantation (EBMT). The EBMT is a nonprofit medical, educational and scientific organization founded more than forty years ago, focused in the field of hematopoietic stem cell transplantation. The EBMT has more than 5000 members from 85 countries. This organization is focused on innovation, research and the advancement in the field of hematopoietic stem cell transplantation. The EBMT hosts a unique patient registry providing a pool of data to researchers all over the world to perform studies and assess new improvements in the field (EBTM [6]). Acute myeloid leukemia (AML) in children is rare in comparison with acute lymphoid leukemia and has overall survival of almost 75% at five years. Patients with high-risk characteristics need to achieve remission and undergo bone marrow transplantation with an overall survival of 70% (de Rooij et al. [7]). Post-transplant relapse is common and involves almost 80% of treated AML, of which 75% occur in the first year after transplantation (Sander et al. [8]).

The data were presented in Fiocco et al. [9] and are available in the mstate R package [10], labeled as ebmt4. The data set includes 2279 patients transplanted between 1985 and 1998. We considered the time in years from transplantation to relapse as the response variable, which had 84% of censored times. Figure 1 shows the Kaplan–Meier (KM) estimator for this data set and the histogram for the failure times.

Note that there is a plateau in the KM estimator, suggesting the presence of cure individuals for this problem. In addition, the histogram of failure times suggest a heavy-tailed distribution for the noncured individuals. Some descriptive statistics for such failure times are mean = 1.14, median = 0.52, standard deviation = 1.84 and kurtosis = 18.71. The high kurtosis for these data, again, suggests the use of a heavy-tailed model for this problem.

Allogeneic stem cell transplantation (allo-HSCT) is a therapeutic approach used for a wide array of malignant and nonmalignant hematologic disorders. Its effectiveness is related to mortality and morbidity. Several complications can occur, including hemorrhages, infections, GVHD and endothelial-damage-related toxicities. There are new immunological and genetic therapies to avoid these complications and reduce mortality (See Horio et al. [11], Radujkovic et al. [12] and Tsai et al. [13]). Other immunodeficiency diseases such as ataxia-telangiectasia can be treated by bone marrow transplantation as mentioned in Sabino Pinho de Oliveira et al. [14]. Further studies with an increased cohort of patients and a longer period of follow-ups are required to validate and monitor the associated risks of transplantation. For this reason, it is important to use a heavy-tailed distribution such as the SHN distribution to study the survival in this type of bone marrow transplant and therapies. We emphasize that, up to this moment, this distribution has not been used in a cure rate model context. In this context, we define the power series–SHN (PS–SHN) class of univariate cure rate models obtained by compounding the SHN and power series distributions. The proposed model defines six new cure rate models as special cases.

The contents in this paper are organized in the following manner. In Section 2, we introduce the SHN distribution and derive some of its distributional properties. In Section 3, we introduce a general cure rate model. In Section 4, we describe a parameter estimation by the maximum likelihood method using the expectation–maximization (EM) algorithm. In Section 5, we present two simulation studies, where the first one is related to assessing the performance of the maximum likelihood estimation (MLE) using the EM algorithm and the second one is the misspecification for the time-to-event for the concurrent causes. In Section 6, we present a real data set on bone marrow transplantation. Finally, some conclusions are included in Section 7.

2. The SHN Distribution

The SHN distribution [1] with parameters $\sigma >0$ and $\alpha >0$ (denoted by SHN$(\sigma ,\alpha )$) is a model with positive support. Its probability density function (pdf) is given by

$$\begin{array}{c}\hfill f_{SHN}(t;\mathbf{\lambda })=\alpha \sqrt{\frac{2^{\alpha }}{\pi }}{\sigma }^{\alpha }{t}^{-(\alpha +1)}\gamma \left(\frac{t^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right),t>0,\end{array}$$

(1)

where $\lambda =(\alpha ,\sigma )$ and $\gamma (x,a)={\int }_0^x{u}^{a-1}{e}^{-u}du$ is the lower incomplete gamma function. The cumulative distribution function (cdf) of the SHN model is given by

$$F_{SHN}(t;\mathbf{\lambda })={\int }_0^t{f}_{SHN}(w;\mathbf{\lambda })dw=\alpha \sqrt{\frac{2^{\alpha }}{\pi }}{\sigma }^{\alpha }{\int }_0^t{w}^{-(\alpha +1)}\gamma \left(\frac{w^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)dw.$$

In the next Proposition, we present an alternative expression for the cdf.

Proposition 1.

The cdf of the SHN distribution is given by

$$F_{SHN}(t;\mathbf{\lambda })=2\mathsf{\Phi }\left(\frac{t}{\sigma }\right)-1-\sqrt{\frac{2^{\alpha }}{\pi }}{\left(\frac{\sigma }{t}\right)}^{\alpha }\Gamma \left(\frac{\alpha +1}{2}\right)G\left(\frac{w^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right),t>0,$$

where $\mathsf{\Phi }(·)$ denotes the cdf of the standard normal model and $G(·,a)$ is the cdf of the gamma model with shape and rate parameters equal to a and 1, respectively.

Proof. 

By definition

$$\begin{array}{cc}\hfill F_{SHN}(t;\mathbf{\lambda })& =\alpha \sqrt{\frac{2^{\alpha }}{\pi }}{\sigma }^{\alpha }{\int }_0^t{w}^{-(\alpha +1)}\Gamma \left(\frac{\alpha +1}{2}\right)G\left(\frac{w^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)dw.\hfill \end{array}$$

Using techniques of integration by parts, we have that

$$\begin{array}{cc}\hfill u=G\left(\frac{w^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)\Rightarrow du& =\frac{w}{{\sigma }^2}g\left(\frac{w^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)dw=\frac{w}{{\sigma }^2\Gamma \left(\frac{\alpha +1}{2}\right)}{\left(\frac{w^2}{2{\sigma }^2}\right)}^{\frac{\alpha +1}{2}-1}exp\left\{-\frac{w^2}{2{\sigma }^2}\right\}dw\hfill \\ \hfill dv=w^{-(\alpha +1)}\Rightarrow v& =-\frac{w^{-\alpha }}{\alpha },\hfill \end{array}$$

where $g(·,·)$ denotes the density function of the gamma distribution. Therefore,

$$\begin{array}{cc}\hfill F_{SHN}(t;\mathbf{\lambda })& =\alpha \sqrt{\frac{2^{\alpha }}{\pi }}{\sigma }^{\alpha }\Gamma \left(\frac{\alpha +1}{2}\right)\hfill \\ \hfill & \times \left[-G\left(\frac{t^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)\frac{t^{-\alpha }}{\alpha }+\frac{1}{\Gamma \left(\frac{\alpha +1}{2}\right)\alpha {\sigma }^{\alpha }}\sqrt{\frac{\pi }{2^{\alpha }}}{\int }_0^t\frac{1}{\sqrt{2\pi }\sigma }exp\left\{-\frac{w^2}{2{\sigma }^2}\right\}dw\right].\hfill \end{array}$$

and the result follows after identifying the pdf for the normal distribution within the integral. □

Corollary 1.

The survival function for the SHN model is given by

$$S_{SHN}(t;\mathbf{\lambda })=2\mathsf{\Phi }\left(-\frac{t}{\sigma }\right)+\sqrt{\frac{2^{\alpha }}{\pi }}{\left(\frac{\sigma }{t}\right)}^{\alpha }\Gamma \left(\frac{\alpha +1}{2}\right)G\left(\frac{t^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right),t>0.$$

(2)

Note that Corollary 1 provides an alternative way to program the survival function, especially from a computational point of view, because many software tools include the cdf of the standard normal and gamma models and the lower incomplete gamma function.

The SHN distribution has a convenient stochastic representation. If $X\sim HN(\sigma )$ and $Y\sim U(0,1)$ are independent (i.e., the half-normal and uniform distributions, respectively), then

$$T=\frac{X}{Y^{1/\alpha }}\sim SHN(\sigma ,\alpha ).$$

(3)

The mean, variance, asymmetry ($\sqrt{{\beta }_1}$) and kurtosis (${\beta }_2$) coefficients from this distribution are

• $\mathbb{E}(T)=\sqrt{\frac{2}{\pi }}\frac{\alpha \sigma }{\alpha -1},\alpha >1$.
• Var$(T)=\frac{\alpha {\sigma }^2\left[(\pi -2)\alpha (\alpha -2)+\pi \right]}{\pi {(\alpha -1)}^2(\alpha -2)},\alpha >2$.
• $\sqrt{{\beta }_1}=\frac{\pi \sqrt{2(\alpha -2)}\left[\frac{4}{\pi }{\alpha }^2(\alpha -2)(\alpha -3)-{(\alpha -1)}^2(\alpha -4)(\alpha +1)\right]}{\sqrt{\alpha }(\alpha -3){\left[\alpha (\pi -2)(\alpha -2)+\pi \right]}^{3/2}},\alpha >3$.
• ${\beta }_2=\frac{3\alpha {(\alpha -2)}^2(\alpha -3)\left[{\pi }^2{(\alpha -1)}^4-4{\alpha }^3(\alpha -4)\right]-4\pi {\alpha }^2{(\alpha -1)}^2(\alpha -2)(\alpha -4)({\alpha }^2-3\alpha +8)}{{\alpha }^2(\alpha -3)(\alpha -4){\left[\alpha (\pi -2)(\alpha -2)+\pi \right]}^2},\alpha >4$.

Note that $\sigma$ is a scale parameter and $\alpha$ is related mainly with the kurtosis of the model. In fact, as discussed in Olmos et al. [1], low values for $\alpha$ allow us to accommodate distributions with greater kurtosis coefficients. On the other hand, when $\alpha \to \infty$, the HN distribution is recovered. Besides being a parsimonious model, an interesting property from the SHN distribution is that its tails are heavier than common distributions with positive support, as we enunciate in the following proposition.

Proposition 2.

Suppose a distribution with positive support satisfying the following conditions:

(A1) 

Its density function can be written in the form $f(t;\mathbf{\lambda })=C{t}^{a}exp\{-g(t;\mathbf{\lambda })\},t>0,$ where C does not depend on t, $a\in \mathbb{R}$ and $g(t;\mathbf{\lambda })\ge 0$, $\forall t>0$.

(A2) 

${\displaystyle \underset{t\to \infty }{lim}{t}^{a^{\prime }}exp\{-g(t;\mathbf{\lambda })\}=0,\forall a^{\prime }>a+1}$.

Then, the SHN model has heavier tails than such a distribution.

Proof. 

To show that the SHN model has heavier tails than a model with a density function given as in assumption A1, it is enough to prove that

$$\begin{array}{c}\hfill \underset{t\to +\infty }{lim}\frac{f(t;\mathbf{\lambda })}{f_{SHN}(t;\sigma ,\alpha )}=0.\end{array}$$

Therefore,

$$\begin{array}{cc}\hfill \underset{t\to +\infty }{lim}\frac{f(t;\mathbf{\lambda })}{f_{SHN}(t;\sigma ,\alpha )}& =\underset{t\to +\infty }{lim}\frac{C{t}^{a}exp\{-g(t;\mathbf{\lambda })\}}{\alpha {\sigma }^{\alpha }\sqrt{\frac{2^{\alpha }}{\pi }}{t}^{-(\alpha +1)}\gamma \left(\frac{t^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)}\hfill \\ \hfill & =\frac{C}{\alpha {\sigma }^{\alpha }\sqrt{\frac{2^{\alpha }}{\pi }}}\underset{t\to +\infty }{lim}\frac{t^{a+\alpha +1}exp\{-g(t;\mathbf{\lambda })\}}{\gamma \left(\frac{t^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)}.\hfill \end{array}$$

As ${lim}_{t\to +\infty }\gamma \left(t,c^{\prime }\right)=\Gamma (c^{\prime })$, $\forall c^{\prime }>0$, $\alpha >0$, by assumption A2 and by the continuity of the limit, we have the result. □

Corollary 2.

If $X_1$ and $X_2$ are two random variables with density functions satisfying conditions A1 and A2 in Proposition 2 and $Y=p{X}_1+(1-p)X_2$, $p\in (0,1)$, then the SHN model also has heavier tails than the distribution related to Y.

Remark 1.

Common distributions used in a cure rate model context satisfy the conditions A1 and A2 in Proposition 2. For instance:

(a) 

If $C=\alpha /{\sigma }^{\alpha }$, $a=\alpha -1$ and $g(t;\mathbf{\lambda })={(t/\sigma )}^{\alpha }$, the Weibull model is recovered (Chen et al. [15], Rodrigues et al. [16], Pal and Balakrishnan [17]).

(b) 

If $C={[\Gamma (\alpha ){\sigma }^{\alpha }]}^{-1}$, $a=\alpha -1$ and $g(t;\mathbf{\lambda })=t/\sigma$, the gamma distribution is recovered (Balakrishnan and Pal [18], Wiangnak and Pal. [19], Ortega et al. [20]).

(c) 

If $C={\left[\sigma \sqrt{2\pi }\right]}^{-1},a=-1$ and $g(t;\mathbf{\lambda })={(logt-\alpha )}^2/(2{\sigma }^2)$, we obtain the log-normal model (Balakrishnan and Pal [21], Gallardo et al. [22]).

(d) 

If $C=\sqrt{\sigma /(2\pi )}exp(\sigma /\alpha )$, $a=-3/2$ and $g(t;\mathbf{\lambda })=(\sigma /2)\left[t/{\alpha }^2+1/t\right]$, the inverse Gaussian model is recovered. We denote it as $IG(\alpha ,\sigma )$.

(e) 

If $C=\sqrt{\sigma /(2\pi )}exp(\sigma \alpha )$, $a=-1/2$ and $g(t;\mathbf{\lambda })=(\sigma /2)\left[t/{\alpha }^2+1/t\right]$, the density function of the reciprocal of an inverse Gaussian random variable is recovered. Up to this moment, this model has not been used in a cure rate model context. We denote it as $RIG(\alpha ,\sigma )$.

(f) 

If $C={\left[B(\alpha ,\sigma )\right]}^{-1}$, $a=\alpha -1$ and $g(t;\mathbf{\lambda })=(\alpha +\sigma )log(1+t)$, the density of the beta prime distribution is recovered (Leao et al. [23]).

Remark 2.

The Birnbaum–Saunders (BS, Birnbaum and Saunders [24,25]) model with density function given by

$$f(t;\alpha ,\sigma )=\frac{1}{\sqrt{2\pi }}exp\left(-\frac{1}{2{\alpha }^2}\left[\frac{t}{\sigma }+\frac{\sigma }{t}-2\right]\right)\frac{(t+\sigma )}{2\alpha \sqrt{\sigma t^3}},t>0,$$

can be written as a mixture between $IG(\beta ,\beta /{\alpha }^2)$ and $RIG(\beta ,1/(\beta {\alpha }^2))$ with equal probabilities (see Desmond [26] for details). Therefore, by Corollary 2, the SHN distribution has heavier tails than the BS distribution.

In the next Section, we discuss the use of the SHN model in the presence of long-term survivors.

3. A General Cure Rate Model and the SHN Distribution

To exemplify the use of the SHN distribution in a cure rate model context, we consider a general class of models based on a competitive risk scheme. The reader is referred to Chen et al. [15] for more details. Let M be a random variable denoting the initial number of carcinogenic cells of an individual. Several models have been considered for this scheme, among others, Bernoulli [27], Poisson [15] and negative binomial [16].

Let $W_a$, $a=1,\dots ,M$ be the random variable representing the time at which the ath concurrent cause produces the event of interest. For instance, in a cancer context, $W_a$ represents the progressive time for the ath tumor cell. As shown in the Appendix A, a general class of model is produced when it is assumed that M follows a power series distribution (PS, Cancho et al. [28]). As is common in this context, it is also assumed that, conditional on M, the $W_a$’s are independent identically distributed with a common cumulative distribution function $F(·;\mathbf{\lambda })$ and survival function $S(·;\mathbf{\lambda })=1-F(·;\mathbf{\lambda })$, where $\mathbf{\lambda }$ is a vector of unknown parameters. Moreover, it is assumed that the sequence $W_1,W_2,\dots ,$ is independent of M. Assuming that only one concurrent cause can produce the of interest, the observable time of the occurrence of the event of interest is given by

$$T=\left\{\begin{array}{cc}+\infty \hfill & ,\mathrm{if}M=0,\hfill \\ min(W_1,\dots ,W_M)\hfill & ,\mathrm{if}M>0.\hfill \end{array}\right.$$

Under this setting, the population survival function is given by

$$\begin{array}{c}\hfill S_{pop}(t;\theta ,\mathbf{\lambda })=\frac{A(\theta S(t;\mathbf{\lambda }))}{A(\theta )},\end{array}$$

(4)

where $A(·)$ is related to the PS distribution. Details about the function $A(·)$ for different particular models are presented in Appendix A. From Equation (4), it is immediate that the cure rate of the model is given by $q_0=a_0/A(\theta )$ and the survival function for the susceptible individuals is

$$\begin{array}{c}\hfill S_{sus}(t;\theta ,\mathbf{\lambda })=\frac{a_0-A(\theta S(t;\mathbf{\lambda }))}{a_0-A(\theta )}.\end{array}$$

(5)

For heterogeneous populations, we consider a set of covariates of dimension $r+1$ for the ith individual, say ${\mathbf{x}}_i={(1,x_{i1},x_{i2},\dots ,x_{ir})}^{\top }$ (i.e., considering the intercept term), related to the cure rate term of the model. For the Poisson, binomial and COM-Poisson, such covariates can be introduced as

$$\begin{array}{c}\hfill {\theta }_i=exp\left({\mathbf{x}}_i^{\top }\mathbf{\beta }\right),i=1,\dots ,n,\end{array}$$

(6)

where n is the sample size and $\mathbf{\beta }=({\beta }_0,{\beta }_1,\dots ,{\beta }_r)$ is the vector of unknown regression coefficients. For the rest of the discrete models, the covariates can be included as

$$\begin{array}{c}\hfill {\theta }_i=\frac{exp\left({\mathbf{x}}_i^{\top }\mathbf{\beta }\right)}{1+exp\left({\mathbf{x}}_i^{\top }\mathbf{\beta }\right)},i=1,\dots ,n.\end{array}$$

(7)

In the literature, many models have been considered for $S(·;\mathbf{\lambda })$. As discussed in Remark 1, Weibull, gamma and log-normal models are common alternatives, to name a few. However, we remark that, up to this moment, no distribution with heavy tails has been used in this context and in particular, the SHN distribution has not been considered to model the survival function for the concurrent causes. Henceforth, we refer this model as PS–SHN. It is interesting to explore the use of the SHN distribution in a cure rate model context because it is usual that most failures occur until a certain time and from that moment mainly censoring times are observed. However, in diseases, from a certain time, it is considered that the patient is “cured”. In practice, and especially in clinical trials with large sample sizes, it may occur that failure times are observed in unusual times, higher than expected.

4. Estimation

In order to estimate the model parameters, we utilize the maximum likelihood (ML) method. Let us consider the situation when the time-to-event is not completely observed and is subject to right-censoring. Let $c_i$ denote the censoring time for the ith individual. We observe $y_i=min(t_i,c_i)$ and ${\delta }_i=\mathbb{I}(t_i\le c_i)$, where ${\delta }_i=1$ if $t_i$ is a time-to-event and ${\delta }_i=0$ if right-censored, $i=1,\dots ,n$. From n observed vectors $(y_1,{\delta }_1,{\mathit{x}}_1^{\top }),\dots ,(y_n,{\delta }_n,{\mathit{x}}_n^{\top })$, the corresponding likelihood function, under noninformative censoring, can be expressed as

$$\begin{array}{cc}\hfill \ell (\mathbf{\psi })& =\sum _{i=1}^n\left\{{\delta }_{i}log\left[f_{pop}(y_i;{\theta }_i,\mathbf{\lambda })\right]+(1-{\delta }_i)log\left[S_{pop}(y_i;{\theta }_i,\mathbf{\lambda })\right]\right\}\hfill \\ \hfill & =\sum _{i=1}^n\{(1-{\delta }_i)log\left[A({\theta }_i{S}_{SHN}(y_i;\mathbf{\lambda }))\right]-log\left[A({\theta }_i)\right]\hfill \\ \hfill & +{\delta }_i\left(log({\theta }_i)+log\left[f_{SHN}(y_i;\mathbf{\lambda })\right]+log\left[A^{\prime }({\theta }_i{S}_{SHN}(y_i;\mathbf{\lambda }))\right]\right)\},\hfill \end{array}$$

(8)

where $\mathbf{\psi }={({\mathbf{\beta }}^{\top },{\mathbf{\lambda }}^{\top })}^{\top }$, and $f_{SHN}(·;·)$ and $S_{SHN}(·;·)$ are defined in (1) and (2). Details for the construction of Equation (8) can be seen in Cancho et al. [28] and Gallardo et al. [29].

The MLE of $\mathbf{\psi }$, denoted by $\widehat{\mathbf{\psi }}$, can be obtained by maximizing Equation (8) in relation to $\mathbf{\psi }$. However, this procedure can be computationally expensive and it would not necessarily provide a global maximum, especially when the dimension of $\mathbf{\beta }$ is high. To facilitate the obtainment of the MLE for $\mathbf{\psi }$, we can use the EM algorithm developed in Gallardo et al. [29] for this class of models. The kth iteration of such procedure consists of the following steps:

• E-step: For $i=1,\dots ,n$, compute

  $$\begin{array}{cc}\hfill {\tilde{M}}_i^{(k)}& =(1-{\delta }_i)\mathbb{E}\left(M_i;{\theta }_i^{(k-1)}{S}_i^{(k-1)}\right)+{\delta }_i\frac{\mathbb{E}\left(M_i^2;{\theta }_i^{(k-1)}{S}_i^{(k-1)}\right)}{\mathbb{E}\left(M_i;{\theta }_i^{(k-1)}{S}_i^{(k-1)}\right)},\hfill \end{array}$$

  where $\mathbb{E}\left(M_i;{\mu }_i\right)$ and $\mathbb{E}\left(M_i^2;{\mu }_i\right)$ denote the two first moments of the PS distribution with parameter ${\mu }_i$ given by

  $$\mathbb{E}\left(M_i;{\mu }_i\right)={\mu }_i\frac{\partial logA({\mu }_i)}{\partial {\mu }_i}\mathrm{and}\mathbb{E}\left(M_i^2;{\mu }_i\right)={\mu }_i\frac{\partial \mathbb{E}\left(M_i;{\mu }_i\right)}{\partial {\mu }_i}+{\mathbb{E}}^2\left(M_i;{\mu }_i\right).$$
• M1-step: Given ${\mathit{M}}^{(k)}=({\tilde{M}}_1^{(k)},\dots ,{\tilde{M}}_n^{(k)})$, find $\mathbf{\beta }$ that maximizes $Q_1(\mathbf{\beta }\mid {\psi }^{(k)})$ in relation to ${\mathbf{\beta }}^{(k)}$, where

  $$Q_1(\mathbf{\beta }\mid {\psi }^{(k)})=\sum _{i=1}^n\left\{{\tilde{M}}_i^{(k)}log({\theta }_i)-log\left[A({\theta }_i)\right]\right\}.$$
• M2-step: Given ${\mathit{M}}^{(k)}$, find ${\mathbf{\lambda }}^{(k)}$ that maximizes $Q_2(\mathbf{\lambda }\mid {\psi }^{(k)})$ in relation to $\mathbf{\lambda }$, where

  $$\begin{array}{cc}\hfill Q_2(\mathbf{\lambda }\mid {\psi }^{(k)})& =\sum _{i=1}^n\{({\tilde{M}}_i^{(k)}-{\delta }_i)log[2\mathsf{\Phi }\left(-\frac{y_i}{\sigma }\right)\hfill \\ \hfill & +\sqrt{\frac{2^{\alpha }}{\pi }}{\left(\frac{\sigma }{y_i}\right)}^{\alpha }\Gamma \left(\frac{\alpha +1}{2}\right)G\left(\frac{y_i^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)]\hfill \\ \hfill & +{\delta }_i[log(\alpha )+\frac{\alpha }{2}log(2)+\alpha log(\sigma )-(\alpha +1)log(y_i)\hfill \\ \hfill & +log\left[\gamma \left(\frac{y_i^2}{2{\sigma }^2},\frac{\alpha +1}{2}\right)\right]\left]\right\}.\hfill \end{array}$$

The E-, M1- and M2-steps are repeatedly alternated until a suitable convergence rule is satisfied, e.g., the difference in successive values of the estimates is less than a prespecified tolerance. Note that this procedure requires two independent maximization procedures in relation to $\mathbf{\beta }$ and $\mathbf{\lambda }$ of dimensions $r+1$ and 2, respectively, instead of a unique maximization procedure in relation to $\mathbf{\psi }$, which has a dimension $(r+3)$. Standard errors for $\widehat{\mathbf{\psi }}$ can be computed based on the Hessian matrix, i.e., $\mathbf{H}(\widehat{\mathbf{\psi }})={\partial }^2\ell (\mathbf{\psi })/\partial {\mathbf{\psi }}^{\top }\partial \mathbf{\psi }{|}_{\mathbf{\psi }={\mathbf{\psi }}^{\top }}$. This matrix can be estimated numerically using, for instance, the package pracma [30] of R [31]. As the estimators obtained via the EM algorithm are MLE, they satisfy the asymptotic distribution

$$\sqrt{n}\left(\widehat{\mathbf{\psi }}-\mathbf{\psi }\right)\stackrel{\mathcal{D}}{\to }N\left(\mathbf{0},{\left[\mathbf{H}(\widehat{\mathbf{\psi }})\right]}^{-1}\right),\mathrm{as}n\to +\infty ,$$

(9)

where $\stackrel{\mathcal{D}}{\to }$ denotes convergence in distribution. This result allows us to compute an approximated confidence interval for the parameters. Moreover, based on the delta method (see Sen et al. [32] for details), we can also build a confidence interval for any scalar function of $\mathbf{\psi }$, say $g(\mathbf{\psi })$, using

$$\sqrt{n}\left[g\left(\widehat{\mathbf{\psi }}\right)-g\left(\mathbf{\psi }\right)\right]\stackrel{\mathcal{D}}{\to }N\left(\mathbf{0},{\left[\nabla g\left(\mathbf{\psi }\right)\right]}^{\top }{\left[\mathbf{H}(\widehat{\mathbf{\psi }})\right]}^{-1}\nabla g\left(\mathbf{\psi }\right)\right),\mathrm{as}n\to +\infty ,$$

(10)

where $\nabla g\left(\mathbf{\psi }\right)=\partial g\left(\mathbf{\psi }\right)/\partial \mathbf{\psi }$ denotes the gradient of $g\left(\mathbf{\psi }\right)$. In particular, a very useful case is $g\left(\mathbf{\psi }\right)=a_0/A({\theta }_i)$, i.e., the cure rate of the corresponding model. The described estimation procedure was incorporated in the package PScr [33] of R [31].

5. Simulation Studies

In this Section, we present two simulation studies. The first one is related to assessing the performance of the MLE using the EM algorithm discussed in Section 3. The second study is devoted to assessing the performance of the PS–SHN compared with the PS–Weibull model when data are drawn from the PS–Weibull model with some modified scheme. The software tool used for this study was R [31].

5.1. Assessing the Performance of MLE in a Finite Sample Size

In this subsection, we present a simulation study to evaluate the performance of the MLE estimators of the PS–SHN in a finite sample size. We considered three particular distributions belonging to the PS model for the number of concurrent causes: Poisson, logarithmic and geometric. We also considered three sample sizes: 100, 200 and 400. For each individual, we drew two covariates, say $x_1$ and $x_2$, where $x_1$ was drawn from the Bernoulli distribution with a success probability of 0.5 and $x_2$ from the standard normal model. For this reason, $\mathbf{\beta }=({\beta }_0,{\beta }_1,{\beta }_2)$ had dimension three. For each particular distribution for the concurrent causes, we considered two scenarios, named high and low cure rates. For the high cure rate, we considered that for ${\beta }_2$ fixed at zero, the values of the two levels for $x_1$ for the cure rate were 0.7 and 0.45. With such information, the values of ${\beta }_0$ and ${\beta }_1$ could be determined from the equations

$$exp\left(-exp({\beta }_0)\right)=0.7\mathrm{and}exp\left(-exp({\beta }_0+{\beta }_1)\right)=0.45$$

for the Poisson case and

$$\frac{1}{A\left(\frac{exp\left({\beta }_0\right)}{1+exp\left({\beta }_0\right)}\right)}=0.7\mathrm{and}\frac{1}{A\left(\frac{exp\left({\beta }_0+{\beta }_1\right)}{1+exp\left({\beta }_0+{\beta }_1\right)}\right)}=0.45,$$

for the logarithmic and geometric cases. A similar procedure could be applied for the low cure rate case, where we considered that, for ${\beta }_2$ fixed at zero, the values of the two levels for $x_1$ for the cure rate were 0.5 and 0.25. Summarizing, for the Poisson model, we had ${\beta }_0=-1.031$, ${\beta }_1=0.806$ and ${\beta }_2=0$ for the high cure rate case and ${\beta }_0=-0.367$, ${\beta }_1=0.693$ and ${\beta }_2=0$ for the lower cure rate case; for the logarithmic distribution we had ${\beta }_0=0.132$, ${\beta }_1=1.588$ and ${\beta }_2=0$ for the high cure rate case and ${\beta }_0=1.366$, ${\beta }_1=2.534$ and ${\beta }_2=0$ for the lower cure rate case; and for the geometric model, we had ${\beta }_0=-0.847$, ${\beta }_1=1.048$ and ${\beta }_2=0$ for the high cure rate case and ${\beta }_0=0$, ${\beta }_1=1.099$ and ${\beta }_2=0$ for the lower cure rate case. Finally, we also considered two values for $\alpha$, 0.75 and 2.5 and two values for $\sigma$, 1 and 3. The combinations of distributions for the concurrent causes, sample size, type of cure rate (high or low), $\alpha$ and $\sigma$ produced 72 different scenarios. However, we only present the cases for the Poisson and geometric models, although conclusions are similar for the logarithmic distribution. For each, we considered 1000 replicates and computed the MLE based on the EM algorithm described in Section 3. Such results are summarized in

Table 1. Simulation study to assess the performance of MLE in the Poisson–SHN model.

Cure				$\mathit{n}=100$				$\mathit{n}=200$				$\mathit{n}=400$
Rate	$\mathit{\sigma }$	$\mathit{\alpha }$	Parameter	bias	s.e.	RMSE	CP	bias	s.e.	RMSE	CP	bias	s.e.	RMSE	CP
High	1	0.75	${\widehat{\beta }}_0$	0.023	0.370	0.255	0.978	0.012	0.248	0.200	0.974	0.010	0.177	0.147	0.965
			${\widehat{\beta }}_1$	−0.019	0.335	0.302	0.961	−0.006	0.234	0.210	0.953	−0.001	0.164	0.145	0.952
			${\widehat{\beta }}_2$	0.005	0.165	0.144	0.963	0.003	0.116	0.104	0.958	−0.002	0.081	0.074	0.949
			$\widehat{\alpha }$	0.294	0.979	0.751	0.975	0.172	0.600	0.505	0.960	0.055	0.376	0.335	0.951
			$\widehat{\sigma }$	0.082	0.452	0.372	0.965	0.052	0.321	0.276	0.955	0.015	0.227	0.206	0.952
		2.5	${\widehat{\beta }}_0$	0.002	0.208	0.180	0.962	0.002	0.145	0.124	0.959	0.001	0.102	0.089	0.956
			${\widehat{\beta }}_1$	−0.014	0.316	0.281	0.959	−0.013	0.221	0.192	0.958	−0.008	0.155	0.140	0.953
			${\widehat{\beta }}_2$	−0.011	0.155	0.140	0.942	−0.005	0.109	0.103	0.943	−0.003	0.077	0.071	0.951
			$\widehat{\alpha }$	0.287	1.637	1.416	0.929	0.140	1.342	1.091	0.937	0.127	0.995	0.854	0.946
			$\widehat{\sigma }$	−0.039	0.302	0.229	0.933	−0.011	0.211	0.171	0.941	−0.010	0.146	0.131	0.943
	3	0.75	${\widehat{\beta }}_0$	−0.140	0.922	0.307	0.962	−0.125	0.738	0.250	0.952	−0.101	0.614	0.212	0.950
			${\widehat{\beta }}_1$	0.013	0.371	0.327	0.960	0.009	0.258	0.226	0.954	0.003	0.181	0.170	0.953
			${\widehat{\beta }}_2$	0.006	0.183	0.164	0.948	0.002	0.127	0.111	0.948	0.000	0.089	0.080	0.950
			$\widehat{\alpha }$	0.687	1.658	1.304	0.972	0.657	1.389	1.269	0.961	0.584	1.193	1.169	0.955
			$\widehat{\sigma }$	−0.269	1.073	0.846	0.957	−0.088	0.892	0.719	0.956	−0.040	0.652	0.602	0.953
		2.5	${\widehat{\beta }}_0$	0.100	0.505	0.263	0.969	0.085	0.339	0.201	0.963	0.078	0.223	0.160	0.958
			${\widehat{\beta }}_1$	−0.011	0.325	0.295	0.958	−0.009	0.226	0.201	0.952	−0.008	0.159	0.144	0.951
			${\widehat{\beta }}_2$	0.005	0.161	0.151	0.939	0.005	0.112	0.097	0.942	0.003	0.078	0.068	0.949
			$\widehat{\alpha }$	−0.791	1.976	1.583	0.938	−0.549	1.797	1.496	0.941	−0.409	1.490	1.442	0.949
			$\widehat{\sigma }$	−0.641	1.428	0.893	0.919	−0.409	1.237	0.691	0.928	−0.319	0.939	0.592	0.930
Low	1	0.75	${\widehat{\beta }}_0$	0.013	0.441	0.296	0.977	0.009	0.285	0.226	0.960	0.005	0.206	0.172	0.954
			${\widehat{\beta }}_1$	−0.015	0.413	0.369	0.960	−0.013	0.289	0.270	0.952	−0.001	0.202	0.181	0.951
			${\widehat{\beta }}_2$	0.013	0.198	0.180	0.951	0.008	0.138	0.125	0.950	0.003	0.096	0.088	0.950
			$\widehat{\alpha }$	0.293	1.127	0.747	0.972	0.262	0.778	0.632	0.965	0.094	0.459	0.391	0.959
			$\widehat{\sigma }$	0.084	0.531	0.385	0.944	0.063	0.392	0.324	0.947	0.037	0.272	0.243	0.949
		2.5	${\widehat{\beta }}_0$	0.004	0.238	0.210	0.961	0.002	0.163	0.144	0.953	0.001	0.114	0.104	0.951
			${\widehat{\beta }}_1$	−0.010	0.384	0.329	0.965	−0.009	0.270	0.239	0.957	−0.004	0.189	0.178	0.953
			${\widehat{\beta }}_2$	0.012	0.185	0.166	0.960	0.008	0.129	0.118	0.959	0.004	0.091	0.081	0.951
			$\widehat{\alpha }$	0.175	1.975	1.271	0.924	0.134	1.722	1.144	0.938	0.087	1.198	0.968	0.944
			$\widehat{\sigma }$	−0.057	0.353	0.258	0.912	−0.016	0.244	0.190	0.940	−0.002	0.170	0.146	0.943
	3	0.75	${\widehat{\beta }}_0$	−0.151	0.999	0.346	0.941	−0.139	0.805	0.282	0.948	−0.105	0.679	0.227	0.950
			${\widehat{\beta }}_1$	−0.012	0.464	0.402	0.969	−0.010	0.323	0.279	0.962	−0.007	0.225	0.202	0.958
			${\widehat{\beta }}_2$	0.012	0.221	0.206	0.943	0.010	0.153	0.140	0.946	0.004	0.107	0.098	0.950
			$\widehat{\alpha }$	0.526	1.660	1.336	0.971	0.371	1.378	1.251	0.965	0.294	1.276	1.154	0.959
			$\widehat{\sigma }$	−0.467	1.285	0.938	0.929	−0.176	2.077	0.719	0.936	−0.025	1.785	0.619	0.940
		2.5	${\widehat{\beta }}_0$	0.107	0.586	0.295	0.973	0.093	0.405	0.225	0.962	0.086	0.261	0.177	0.959
			${\widehat{\beta }}_1$	0.013	0.398	0.365	0.955	0.008	0.276	0.239	0.954	0.000	0.194	0.179	0.951
			${\widehat{\beta }}_2$	−0.017	0.192	0.173	0.965	−0.002	0.133	0.119	0.961	0.000	0.093	0.081	0.959
			$\widehat{\alpha }$	−0.912	4.181	1.545	0.934	−0.682	3.910	1.489	0.936	−0.493	3.601	1.462	0.942
			$\widehat{\sigma }$	−0.743	1.490	0.969	0.961	−0.526	1.287	0.764	0.960	−0.379	1.038	0.647	0.958

and

Table 2. Simulation study to assess the performance of MLE in the geometric–SHN model.

Cure				$\mathit{n}=100$				$\mathit{n}=200$				$\mathit{n}=400$
Rate	$\mathit{\sigma }$	$\mathit{\alpha }$	Estimator	bias	s.e.	MSE	CP	bias	s.e.	MSE	CP	bias	s.e.	MSE	CP
High	1	0.75	${\widehat{\beta }}_0$	0.025	0.556	0.391	0.970	0.008	0.374	0.278	0.983	0.006	0.259	0.219	0.966
			${\widehat{\beta }}_1$	0.023	0.469	0.424	0.955	0.021	0.328	0.289	0.956	0.001	0.230	0.204	0.955
			${\widehat{\beta }}_2$	0.004	0.234	0.213	0.954	0.003	0.162	0.150	0.942	0.001	0.114	0.106	0.936
			$\widehat{\alpha }$	0.284	1.203	0.741	0.996	0.218	0.796	0.627	0.982	0.109	0.495	0.431	0.961
			$\widehat{\sigma }$	0.078	0.546	0.379	0.949	0.070	0.394	0.342	0.950	0.033	0.280	0.251	0.948
		2.5	${\widehat{\beta }}_0$	−0.014	0.353	0.302	0.959	−0.003	0.241	0.218	0.947	0.003	0.168	0.152	0.949
			${\widehat{\beta }}_1$	0.016	0.453	0.403	0.959	0.012	0.315	0.294	0.956	0.008	0.221	0.199	0.941
			${\widehat{\beta }}_2$	0.006	0.225	0.200	0.959	0.004	0.156	0.143	0.944	0.002	0.109	0.100	0.943
			$\widehat{\alpha }$	0.242	1.871	1.434	0.854	0.163	1.752	1.270	0.910	0.081	1.332	1.018	0.946
			$\widehat{\sigma }$	−0.078	0.375	0.272	0.911	−0.017	0.266	0.206	0.944	−0.014	0.184	0.153	0.964
	3	0.75	${\widehat{\beta }}_0$	−0.153	0.609	0.417	0.989	−0.135	0.488	0.329	0.960	−0.111	0.323	0.262	0.964
			${\widehat{\beta }}_1$	−0.012	0.510	0.448	0.959	−0.003	0.353	0.310	0.961	−0.001	0.247	0.225	0.946
			${\widehat{\beta }}_2$	−0.006	0.253	0.226	0.948	−0.004	0.174	0.155	0.959	−0.003	0.121	0.105	0.961
			$\widehat{\alpha }$	0.534	1.577	1.389	1.000	0.381	1.402	1.313	1.000	0.194	1.356	1.309	1.000
			$\widehat{\sigma }$	−0.490	1.302	0.974	0.910	−0.197	1.153	0.747	0.972	−0.029	0.825	0.638	0.987
		2.5	${\widehat{\beta }}_0$	0.121	0.505	0.385	0.989	0.095	0.380	0.289	0.990	0.083	0.232	0.215	0.985
			${\widehat{\beta }}_1$	−0.005	0.460	0.427	0.961	−0.003	0.321	0.279	0.961	−0.002	0.225	0.199	0.957
			${\widehat{\beta }}_2$	0.012	0.228	0.196	0.960	0.002	0.159	0.137	0.955	−0.001	0.111	0.099	0.949
			$\widehat{\alpha }$	−0.921	2.029	1.755	0.841	−0.709	1.903	1.570	0.835	−0.569	1.661	1.535	0.819
			$\widehat{\sigma }$	−0.788	1.333	1.004	0.851	−0.555	1.052	0.825	0.886	−0.407	0.899	0.672	0.914
Low	1	0.75	${\widehat{\beta }}_0$	0.008	0.477	0.344	0.981	0.003	0.338	0.254	0.969	0.001	0.226	0.194	0.971
			${\widehat{\beta }}_1$	0.004	0.422	0.384	0.954	0.003	0.294	0.256	0.957	0.003	0.207	0.189	0.951
			${\widehat{\beta }}_2$	−0.011	0.212	0.197	0.940	−0.006	0.147	0.133	0.953	−0.003	0.103	0.095	0.948
			$\widehat{\alpha }$	0.280	1.043	0.702	0.998	0.190	0.693	0.583	0.981	0.097	0.441	0.387	0.957
			$\widehat{\sigma }$	0.071	0.480	0.349	0.960	0.059	0.352	0.300	0.960	0.037	0.251	0.217	0.952
		2.5	${\widehat{\beta }}_0$	0.012	0.309	0.274	0.954	0.010	0.213	0.186	0.956	0.004	0.149	0.136	0.953
			${\widehat{\beta }}_1$	0.010	0.409	0.360	0.957	0.007	0.286	0.247	0.964	0.003	0.201	0.180	0.942
			${\widehat{\beta }}_2$	−0.007	0.206	0.190	0.961	−0.005	0.143	0.125	0.957	−0.002	0.100	0.084	0.962
			$\widehat{\alpha }$	0.353	1.622	1.301	0.868	0.274	1.438	1.258	0.932	0.269	1.025	0.968	0.939
			$\widehat{\sigma }$	−0.078	0.340	0.260	0.900	−0.040	0.244	0.184	0.956	−0.015	0.167	0.144	0.952
	3	0.75	${\widehat{\beta }}_0$	−0.155	0.524	0.379	0.970	−0.128	0.415	0.284	0.973	−0.109	0.266	0.240	0.961
			${\widehat{\beta }}_1$	0.009	0.448	0.399	0.963	0.005	0.311	0.271	0.962	0.003	0.219	0.187	0.956
			${\widehat{\beta }}_2$	−0.016	0.223	0.196	0.965	−0.006	0.155	0.136	0.959	−0.003	0.109	0.099	0.950
			$\widehat{\alpha }$	0.605	1.893	1.506	1.000	0.481	1.501	1.388	1.000	0.276	1.372	1.292	1.000
			$\widehat{\sigma }$	−0.340	1.372	0.894	0.929	−0.089	1.116	0.703	0.966	−0.053	0.773	0.612	0.993
		2.5	${\widehat{\beta }}_0$	0.090	0.607	0.347	0.987	0.074	0.439	0.268	0.990	0.051	0.282	0.201	0.982
			${\widehat{\beta }}_1$	0.011	0.415	0.378	0.954	0.005	0.289	0.270	0.946	0.001	0.203	0.180	0.949
			${\widehat{\beta }}_2$	−0.005	0.208	0.187	0.952	−0.001	0.145	0.135	0.937	−0.001	0.101	0.094	0.944
			$\widehat{\alpha }$	−0.865	1.780	1.612	0.852	−0.630	1.708	1.517	0.837	−0.449	1.629	1.462	0.838
			$\widehat{\sigma }$	−0.695	1.283	0.927	0.887	−0.492	1.097	0.762	0.905	−0.341	0.732	0.605	0.924

through the average of the estimated bias (bias), the average of the estimated standard error (s.e.) based on the Hessian matrix numerically evaluated at $\mathbf{\psi }=\widehat{\mathbf{\psi }}$, the root of the estimated mean squared error (RMSE) and the 95% coverage probability (CP). As expected, the bias of the estimators was reduced and the terms s.e. and RMSE were closer when the sample size was increased, suggesting the asymptotic consistency of the estimators. Additionally, the terms CP were closer to the nominal values used for their construction when n increased.

5.2. Misspecification for the Time-to-Event for the Concurrent Causes

In this subsection, we assessed the performance of the MLE for the Poisson–SHN and geometric–SHN models when data were drawn from an artificial model with heavy tails based on the Weibull distribution.

The generation of the data was very similar to that in the last subsection, considering the following changes:

• The times related to the concurrent causes were drawn from the Weibull model with mean $\mathbb{E}(W)=\mu$ and variance 1. We considered three values for $\mu$: 2.5, 5.0 and 7.5. For the parameterization considered in part a) from Remark 1, this corresponded to $\alpha =2.70$ and $\sigma =7.90$, $\alpha =5.80$ and $\sigma =29.16$ and $\alpha =8.97$ and $\sigma =62.75$, respectively.
• A determined number of failure times (2% of the sample size, i.e., two observations for $n=100$, four observations for $n=200$ and eight observations for $n=400$), were imputed as $max\left(t^{fail}\right)+U{\sigma }_t^{fail}$, where $max\left(t^{fail}\right)$ and ${\sigma }_t^{fail}$ denote the maximum and the standard error of the failure times and $U\sim U(1,2)$. This scheme artificially simulated a distribution for the time-to-event of the concurrent causes with a heavier tail than the Weibull distribution, but not corresponding to the SHN.

The main goal of this study was to compare the coefficients related to the cure (i.e., the components of the vector $\mathbf{\beta }$) obtained with the Poisson–SHN and Poisson–Weibull models (when the true model for M is the Poisson distribution) and geometric–SHN and geometric–Weibull models (when the true model is the geometric distribution). All the considered combinations of $\mathbb{E}(W)$, sample size, type of cure (high/low) and model for M (Poisson/geometric), totalizing 36 cases. Each was replicated 1000 times.

Table 3. Misspecification simulation study for Poisson–SHN and geometric–SHN models.

			$\mathit{n}=100$				$\mathit{n}=200$				$\mathit{n}=400$
Poisson Case
Cure			bias		RMSE		bias		RMSE		bias		RMSE
Rate	$\mathbb{E}(\mathit{W})$	Estimator	SHN	Weibull	SHN	Weibull	SHN	Weibull	SHN	Weibull	SHN	Weibull	SHN	Weibull
High	2.5	${\widehat{\beta }}_0$	0.059	0.089	0.394	0.431	0.053	0.084	0.248	0.309	0.051	0.073	0.168	0.221
		${\widehat{\beta }}_1$	−0.018	−0.066	0.356	0.430	−0.006	−0.049	0.248	0.316	−0.003	−0.048	0.175	0.227
		${\widehat{\beta }}_2$	−0.002	−0.003	0.176	0.255	−0.002	−0.003	0.122	0.187	−0.001	−0.002	0.086	0.131
	5.0	${\widehat{\beta }}_0$	−0.109	−0.159	0.344	0.347	−0.077	−0.126	0.225	0.239	−0.058	−0.123	0.134	0.169
		${\widehat{\beta }}_1$	−0.014	−0.020	0.402	0.773	−0.007	−0.018	0.280	0.601	−0.005	−0.017	0.197	0.451
		${\widehat{\beta }}_2$	−0.005	−0.006	0.197	0.288	−0.003	−0.004	0.138	0.209	0.000	0.000	0.097	0.146
	7.5	${\widehat{\beta }}_0$	−0.307	−0.359	0.246	0.297	−0.276	−0.355	0.201	0.210	−0.251	−0.322	0.107	0.159
		${\widehat{\beta }}_1$	−0.017	−0.025	0.443	0.933	−0.012	−0.015	0.311	0.773	−0.011	−0.013	0.219	0.650
		${\widehat{\beta }}_2$	0.011	0.012	0.218	0.317	−0.003	−0.003	0.152	0.230	0.001	0.001	0.107	0.160
Low	2.5	${\widehat{\beta }}_0$	0.087	0.097	0.287	0.314	0.070	0.079	0.181	0.211	0.058	0.063	0.139	0.144
		${\widehat{\beta }}_1$	−0.014	−0.069	0.437	0.448	−0.012	−0.061	0.306	0.334	−0.007	−0.052	0.215	0.228
		${\widehat{\beta }}_2$	−0.001	0.000	0.207	0.297	0.000	0.000	0.146	0.215	−0.001	−0.001	0.102	0.155
	5.0	${\widehat{\beta }}_0$	−0.093	−0.100	0.234	0.255	−0.071	−0.092	0.165	0.176	−0.052	−0.078	0.110	0.121
		${\widehat{\beta }}_1$	−0.029	−0.048	0.489	0.809	−0.022	−0.025	0.344	0.656	−0.008	−0.012	0.244	0.487
		${\widehat{\beta }}_2$	−0.003	−0.004	0.231	0.336	−0.001	−0.001	0.164	0.241	−0.001	−0.001	0.116	0.172
	7.5	${\widehat{\beta }}_0$	−0.204	−0.285	0.197	0.217	−0.194	−0.282	0.133	0.155	−0.192	−0.279	0.102	0.111
		${\widehat{\beta }}_1$	−0.012	−0.018	0.533	0.893	−0.003	−0.014	0.379	0.742	−0.002	−0.004	0.267	0.583
		${\widehat{\beta }}_2$	−0.008	−0.008	0.252	0.372	−0.005	−0.006	0.179	0.265	−0.003	−0.003	0.127	0.187
Geometric Case
High	2.5	${\widehat{\beta }}_0$	0.071	0.092	0.196	0.244	0.053	0.082	0.163	0.190	0.031	0.048	0.124	0.139
		${\widehat{\beta }}_1$	−0.105	−0.127	0.209	0.281	−0.083	−0.094	0.153	0.218	−0.054	−0.056	0.148	0.164
		${\widehat{\beta }}_2$	−0.003	−0.005	0.254	0.435	0.001	0.002	0.176	0.299	0.001	0.001	0.123	0.207
	5.0	${\widehat{\beta }}_0$	−0.091	−0.129	0.135	0.154	−0.088	−0.119	0.098	0.103	−0.071	−0.106	0.042	0.071
		${\widehat{\beta }}_1$	−0.089	−0.130	0.352	0.427	−0.080	−0.120	0.287	0.329	−0.057	−0.092	0.172	0.243
		${\widehat{\beta }}_2$	0.001	0.000	0.274	0.478	−0.001	−0.001	0.192	0.330	0.001	0.001	0.135	0.229
	7.5	${\widehat{\beta }}_0$	−0.288	−0.333	0.094	0.125	−0.251	−0.311	0.067	0.082	−0.242	−0.271	0.049	0.059
		${\widehat{\beta }}_1$	−0.055	−0.099	0.493	0.512	−0.054	−0.096	0.314	0.389	−0.042	−0.053	0.290	0.301
		${\widehat{\beta }}_2$	0.009	−0.010	0.293	0.524	0.005	−0.008	0.204	0.362	0.001	0.001	0.143	0.248
Low	2.5	${\widehat{\beta }}_0$	−0.075	−0.109	0.245	0.298	−0.062	−0.081	0.219	0.234	−0.046	−0.069	0.147	0.165
		${\widehat{\beta }}_1$	−0.094	−0.130	0.217	0.245	−0.074	−0.083	0.133	0.192	−0.050	−0.061	0.128	0.138
		${\widehat{\beta }}_2$	0.012	0.016	0.235	0.397	0.011	0.012	0.163	0.270	0.002	0.003	0.114	0.187
	5.0	${\widehat{\beta }}_0$	−0.224	−0.259	0.197	0.209	−0.216	−0.237	0.101	0.144	−0.153	−0.203	0.061	0.089
		${\widehat{\beta }}_1$	−0.108	−0.126	0.303	0.385	−0.097	−0.106	0.251	0.283	−0.084	−0.099	0.146	0.207
		${\widehat{\beta }}_2$	0.002	0.002	0.253	0.426	0.002	0.002	0.176	0.292	0.001	0.001	0.123	0.201
	7.5	${\widehat{\beta }}_0$	−0.260	−0.459	0.116	0.150	−0.256	−0.424	0.096	0.106	−0.202	−0.392	0.042	0.073
		${\widehat{\beta }}_1$	−0.050	−0.112	0.337	0.471	−0.049	−0.104	0.274	0.359	−0.035	−0.091	0.202	0.268
		${\widehat{\beta }}_2$	−0.002	−0.004	0.269	0.456	−0.001	−0.002	0.187	0.315	−0.001	−0.001	0.131	0.216

summarizes the bias and RMSE for each combination. Note that, in all the considered cases, the bias and RMSE were lower (or at least the same) for the models considering the SHN rather than the Weibull model. Such differences in favor of the SHN were minimized when the sample size was increased.

6. Application to Bone Marrow Transplantation

In this section, we present a real data set application from the EBMT presented in the introduction section. Additionally, the prophylaxis covariate (dichotomous, 1730 patients with “no” and 549 patients with “yes”) was available.

For the analysis, we considered some existent models in the literature, all of them belonging to the power-series cure rate model, such as the Poisson, logarithmic, negative binomial, Bernoulli, COM-Poisson and polylogarithm models. Additionally, we considered three commons models used for the concurrent causes such as the Weibull, gamma and BS models.

Table 4. AIC and BIC criteria for some cure rate models in the literature based on a competing risk scheme for the ebmt4 data set.

Model for M	Weibull	Gamma	BS	SHN
Po	4650.0/4672.9	4683.0/4706.0	4580.8/4603.8	4537.9/4560.8
Lo	4628.7/4651.6	4639.6/4662.6	4628.1/4651.0	4547.3/4570.3
NB	4561.4/4590.1	4547.4/4576.1	4583.0/4611.7	4540.2/4568.8
Bern	4665.1/4688.0	4712.8/4735.8	4577.0/4599.9	4535.8/4558.8
COM-Po	4639.3/4667.9	4659.5/4688.1	4578.9/4607.6	4537.8/4566.5
PL	4612.1/4658.0	4615.9/4661.8	4585.1/4630.9	4545.3/4591.2

shows the Akaike information criterion (AIC) [34] and Bayesian information criterion (BIC) criteria [35] for those models in the ebmt4 data. Note that both criteria favored models where the SHN for concurrent causes was used.

Table 5. Estimates and standard errors (s.e.) for the Bern/SHN cure rate model in the ebmt4 data set.

Parameter	${\mathit{\beta }}_0$	${\mathit{\beta }}_1$	$\mathit{\alpha }$	$\mathit{\sigma }$
Estimate	$-0.4912$	$-0.2224$	$0.3844$	$0.3120$
s.e.	$0.0907$	$0.1030$	$0.0753$	$0.0310$

shows the estimated parameters for the Bern/SHN cure rate model and its standard errors. Note that both parameters related to the cure were significant with 5% of significance. We also considered a residuals analysis in the selected model based on the Cox–Snell residuals [36,37] and the randomized quantile residuals [38]. The first residual was compared with its Nelson–Aalen estimator and the second with the quantile of the standard normal distribution. In both cases, if the model was correct, we expected to observe the identity function for the respective graphs. According to Figure 2, such residuals suggested that the choice model was correct. Finally, we also computed the Kaplan–Meier (K–M) estimator for patients with and without prophylaxis and the respective point estimation and 95% confidence interval for the survival function for each time. Such confidence interval was based on the delta method discussed in Equation (10). Figure 3 shows those results. Note the closeness of the estimation for the survival function and that provided by the K–M estimator. As expected, patients without prophylaxis had a lower survival function than patients with prophylaxis.

7. Conclusions

In this work, we introduced the PS–SHN cure rate model. We provided a mathematical treatment of the new distribution including an alternative expression for the cumulative distribution function. An EM algorithm to compute maximum likelihood estimation procedures in finite samples was discussed. Two simulation studies were performed to assess some properties of the estimators, showing the good performance of the proposed estimation procedure in finite samples. In the analysis of the bone marrow transplant data, the SHN cure rate model yielded a better fit when compared to other well-known models in the literature, according to the AIC and BIC criteria. We developed an R package to perform inference in the power-series cure rate models. As part of future research, we plan to explore other estimation methods for the proposed cure rate model, for instance, a Bayesian approach. A second line of a future research will be a bias-reduction methodology to estimate the parameters of the proposed regression model.