Proportional Odds under Conway-Maxwell-Poisson Cure Rate Model and Associated Likelihood Inference

Cure rate models are useful while modelling lifetime data involving long time survivors. In this work, we discuss a flexible cure rate model by assuming the number of competing causes for the event of interest to follow the ConwayMaxwell Poisson distribution and the lifetimes of the non-cured individuals to follow a proportional odds model. The baseline distribution is considered to be either Weibull or log-logistic distribution. Under right censoring, we develop the maximum likelihood estimators using EM algorithm. Model discrimination among some well-known special cases are discussed under both likelihoodand information-based criteria. An extensive simulation study is carried out to examine the performance of the proposed model and the inferential methods. Finally, a cutaneous melanoma dataset is analyzed for illustrative purpose.


Introduction
Statistical models accommodating a surviving fraction are known as cure rate models.The cure rate model was first introduced by [8] and [7], and have been subsequently studied by many authors.Applications of cure rate models are not limited to biomedical studies, and can also be seen in industrial reliability, finance, manufacturing, demography, and criminology.The basic cure rate model can be seen as a two-component mixture model.If S p (t) is the population survival function, it can be expressed as where p 0 is the probability of cure and S s (t) is the survival function of the non-cured or susceptible individuals in the population.Primarily, cure data have been analyzed in the literature by the structure of the underling survival model of the non-cured individuals S s (t) as proportional hazards (PH) mixture cure model [20] [15] [24], accelerated failure time (AFT) mixture cure rate model [25] [17] [14], accelerated hazards (AH) mixture cure rate model [26], and proportional odds (PO) mixture cure rate model [10] [18].More generally, a cure model can be approached through a competing risks set up as follows.Suppose M is an unobservable random variable denoting the number of competing causes related to the occurrence of an event of interest.Let W j , j = 1, . . ., m, be the random variables denoting the time-to-event for the jth competing cause.Given M = m, W 1 , . . ., W m are assumed to be independent and identically distributed (i.i.d.) with a common cumulative distribution function where A M (.) is the probability generating function (p.g.f.) of M .
In the present work, we assume a proportional odds model for the distribution of W j , with a parametric assumption on the baseline odds function.To be more specific, the odds function of W j is taken as where O(w) = S(w)/F (w) is the odds of survival up to w, the term θ is linked to covariates as e x ′ γ γ γ2 with x = (x 1 , . . . ,x p ) ′ being a vector of p covariates, γ γ γ 2 = (γ 21 , . . ., γ 2p ) ′ is the proportional odds regression coefficients, and O 0 (w) is the baseline odds function.We can further obtain the survival function of W j as with the corresponding probability density function (p.d.f.) The rest of this paper proceeds as follows.In Section 2, we present briefly the COM-Poisson cure model and its three special cases.Section 3 describes the data and the likelihood, while the estimation of the cure rate and associated inferential issues are discussed in Section 4. In Section 5, an extensive simulation study is carried out.In Section 6, we discuss model discrimination using information-and likelihood-based methods.A data on cutaneous melanoma is analyzed in Section 7 for illustrative purpose.Some concluding comments are finally made in Section 8.

Data and the likelihood
Suppose the time-to-event is not completely observed and is subject to non-informative right censoring, which means that the data above a certain value are not observed.Therefore, the observation time T i , for the ith subject, would be the minimum of the censoring time C i and actual lifetime Y i , i.e., We define an indicator function δ i = I(Y i ≤ C i ) for the i-th subject such that δ i = 0 if the lifetime is observed while δ i = 1 if the lifetime is right censored, ∆ 0 and ∆ 1 are sets with all the i's equal to 0 and 1, respectively, and set ∆ * contains all the i's.It is to be noted that the cure rate p 0 = Z(η, ϕ) −1 is purely a function of η for a fixed value of ϕ.The range of 1/p 0 is from 1 to infinity and it is monotone in η.Therefore, it is natural to use a logistic regression model H ϕ (η) = 1 + e x x x ′ i β β β to link the covariate x 1 , . . ., x p to the cured proportion p 0i , i.e., where p 0i is the cured proportion for the ith category, x x x i = (1, x x x ′ ic ) ′ = (1, x i1 , . . ., x ip ) ′ is a vector of p + 1 covariates, and β β β is the vector of regression coefficients.Under this link, η would equal H −1 (1 + e x x x ′ i β β β ), i.e., η can be calculated from the inverse function of H ϕ (.) analytically for the Geometric, Poisson and Bernoulli distributions or by using numerical method for the general COM-Poisson distribution.
For n pairs of observations (t t t, δ δ δ)={(t 1 , δ 1 ),. . .,(t n , δ n )}, the observed data likelihood function under the noninformative censoring is then given by L(θ θ θ; t t t, δ δ δ) where θ θ θ is the set of parameters (ϕ, β β β ′ , γ γ γ ′ ), which is equivalent to Here, we consider two baseline distributions for the proportional odds survival model corresponding to the timeto-event random variable, namely, Weibull and log-logistic distributions.It should also be noted that log-logistic distribution in fact processes the proportional odds property, while the Weibull distribution does not.The survival function and p.d.f. of W under a Weibull baseline, for example, are where γ 0 > 0 and γ 1 > 0 are the shape and scale parameters of the baseline Weibull distribution, respectively.If we assume the baseline distribution to be a log-logistic distribution with γ 0 > 0 and γ 1 > 0 as the scale and shape parameters, respectively, then the corresponding odds function of W i is given by We observe that W i still follows a two-parameter log-logistic distribution (γ 0 , γ 1 > 0) with shape parameter γ 1 and scale parameter γ 0 e −x x x ′ c ,γ γ γ2/γ1 , and with corresponding survival function Note that the mean does not exist if γ 1 < 1 and the variance does not exist if γ 1 < 2.

Estimation of parameters
We propose an Expectation-Maximization (EM) algorithm for obtaining the MLE of θ θ θ, and a profile likelihood approach for the estimation of the dispersion patameter ϕ.It is well-known that EM is an effective technique for finding the MLEs of unknown parameters of a model involving unobserved variables (for further discussion, refer to [19]).In our model, the random variable I i 's are observed for i in the set ∆ 1 and unobserved for i in the set ∆ 0 , where I i = 1 if the individual is susceptible and I i = 0 if the individual is cured.Let us denote the set of complete data by (t t t, δ δ δ, x x x, I I I)={(t 1 , δ 1 , x x x 1 , I 1 ), . . ., (t n , δ n , x x x n , I n )}.The complete data likelihood function is then where

E-step
The expectation step is achieved by calculating the expected value of the complete data log-likelihood function with respect to the conditional distribution of the unobserved and the current estimates of the parameters θ θ θ (k) = (β β β ′ , γ γ γ ′ ) ′ for a fixed value of ϕ.Let us denote the function as at the k-th iteration step.In our model, I i 's are Bernoulli random variables and we can easily find the conditional expectation if the ith individual being susceptible is Stat., Optim.Inf.Comput.Vol.
which can be further simplified as where

M-step
The M-step is achieved by maximizing the Q(θ θ θ, π π π (k) ) function in (26) in order to obtain the improved estimate of θ θ θ, i.e., The MLEs of β β β and γ γ γ do not have explicit expressions.In this paper, the numerical maximization is carried out by Newton-Raphson method.
For a fixed value of ϕ, the E-step and M-step are alternated until the parameter estimate converges to a desired level of accuracy.The parameter ϕ is determined by using the profile likelihood technique.We consider a range of ϕ with small increment, and then for each value of ϕ, the MLEs of other parameters are found, and the estimates with the largest likelihood is chosen as the final estimate.The following subsections present explicit forms of the firstand second-order derivatives of the Q function as well as update function for the case of COM-Poisson distribution.

Standard errors and asymptotic confidence intervals
We may approximate the asymptotic variance-covariance matrix of the MLEs ( β β β ′ , γ γ γ ′ ) ′ by inverting the observed Fisher information matrix of β β β and γ γ γ, for a fixed value of ϕ.The components of the observed Fisher information matrix can be calculated from the negative of the second-order derivatives of the complete data likelihood function with respect to β β β and γ γ γ (for detailed information, refer to [16]).Thus, we can obtain the standard errors of the estimates and then construct corresponding asymptotic confidence intervals for the parameters.

Suppose β
β β is the MLE of the regression coefficient β β β.The estimated cure rate for the corresponding group i is then p0i = (1 + e x x x ′ i β β β ) −1 , for i = 1, ..., τ .The standard error of p0i can be found through delta method as

Equivalent models
Proposition: The population survival function under Geometric and Bernoulli cure rate model are equivalent through re-parametrization if baseline odds follow log-logistic distribution.
Proof: Cure rates p 0 under Geometric and Bernoulli cure rate models are 1 − η 1 , and 1 1+η3 , respectively.If we equate these cure rates, we obtain a relationship between η 1 and η 3 as Suppose S 1 = γ γ 1 10 e x x xγ 12 γ γ 1 10 e x x xγ 12 +t γ 1 and S 3 = γ γ 1 30 e x x xγ 32 γ γ 1 30 e x x xγ 32 +t γ 1 are the survival functions of the susceptible group under Geometric and Bernoulli cure rate models, respectively.Then, the population survival function for the under the Bernoulli case is If we fix the relationship between γ 30 , γ 10 , γ 32 and γ 12 as then we obtain from (32) that which is the population survival function for the Geometric cure rate model under proportional odds assumption.Thus, these two models in this case are re-parametrizations of each other.

Simulation study
An extensive Monte Carlo simulation study is carried out in this paper to evaluate the performance of the proposed methodology by varying the sample size, cure rate, censoring proportion, lifetime parameter, and underling lifetime distribution.We try to mimic the cutaneous melanoma data, and considered 4 possible categories for the individuals, namely, x = 0, 1, 2, 3. Three different sample sizes are considered in the study: n = 200 (50, 42, 53, 55), n = 400 (95, 102, 97, 106), n = 800 (200,168,212,220) to reflect small, medium and large sample sizes.Moreover, if we assume that β β β = (β 0 , β 1 ) has two parameters, fixing the cure rates for the first and fourth categories would be enough to cover all cases as the cure rates for the second and third categories can then be obtained from β β β.Here, we take (p 00 , p 03 ) = (0.4,0.2) and (p 00 , p 03 ) = (0.6, 0.25) with respect to categories one and four for low and high cure rates, respectively.Also, the cure rate would be in a decreasing order in this way.The β β β's are We thus obtain the true value of β β β as (0.405,0.321) and (-0.405,0.501),respectively.In addition, we consider light and heavy censored data in the simulation.The light and heavy censoring rates are (0.52, 0.45, 0.37, 0.3) with (0.65, 0.49, 0.4, 0.35) as the corresponding cure rates, (0.7, 0.57, 0.45, 0.34) and (0.8, 0.64, 0.5, 0.38) for the corresponding cure rates.It is natural to assume that the probability of censored population for the susceptible group equal to the difference between the probability of getting censored and cured; i.e., If we assume the censoring time C x follows an exponential distribution with rate λ x on x = 0, 1, 2, 3, Eq.(36) can be re-written as The choice of (γ 0 , γ 1 ) in the underlying distribution of the proportional odds survival model are (0.571, 0.307) and (1.75, 3.25) for Weibull and log-logistic distributions, respectively.The odds parameter is specified by γ 2 = −0.75 to ensure a decreasing lifetime for the four nodule categories.We consider an inverse transform sampling method to simulate the actual survival lifetime Y i for each individual under different competing risks, i.e., ] γ0 and w i = γ 0 ( u 1−u e xiγ2 ) 1/γ1 , i = 1, . . ., n, under the proportional odds model with Weibull and log-logistic baseline distribution, respectively, where u follows an uniform distribution with range 0 to 1.
Under the above setting, the procedure to generate the data form different cure rate models is as follows.Geometric cure rate model: For each individual, we simulate the number of competing risk M i from Geometric distributions with probability P (M i = 0) = p 0x ; and we simulate the censoring time C x from exponential distribution with rate λ x , the parameter λ x can be found from (37).If M i does not equal zero, we simulate M i number of actual lifetimes {Y i , .., Y Mi } from proportional odds survival model, and the observed lifetime T i is taken as the minimum of all the actual lifetime and the censoring time, i.e., min{Y i , .., Y Mi , C i }.If Y i = C i , we make the censoring indicator δ i = 1, otherwise δ i = 0. On the other hand, M i = 0 means the individual is censored, we assign C i to the actual lifetime T i , and the censoring indicator is taken to be δ i = 0. Poisson cure rate model: In this case, the procedure is the same as the Geometric cure rate model except that M i is simulated from Poisson distribution with parameter −log(p 0x ).Bernoulli cure rate model: There are two ways to do the data generation in this case.One is the same as Geometric cure rate model except that M i is simulated from Bernoulli distribution with probability of success as 1 − p 0x .Another way is a little bit simpler since M i can only be taken as 0 or 1 in this case.For each individual, we simulate the censoring time C x from exponential distribution with rate λ x .Then, we simulate an uniform random variable U i and if U i ≤ p 0x , the observed lifetime T i is set to C x ; otherwise, we generate the observed lifetime T i from the proportional odds survival model.
In our simulation study, 1000 Monte Carlo runs were considered in each scenario.The estimates were calculated through EM method.We stopped our estimation if the absolute difference between two consecutive estimates was less than 10 −5 .We calculated the empirical Bias, standard errors(SE), root Mean Square Error (RMSE), and 95% coverage probabilities (CPs) for the estimates of the parameters.In addition, we computed the cure rate, SE and 95% CPs.Here, the initial values of the parameters (β β β, γ γ γ) were taken from a grid of parameters with a range from 80% to 120% of the true value, and those estimates having the maximum likelihood were chosen as the initial value.
Tables 1-6 present the bias, SE, RMSE, and coverage probabilities for the three special cases.We can see that the estimates are quite accurate under different cure rate models.The Bias, standard error along with RMSE get reduced as the sample size increases.The same follows when the censoring is light or the cure rate is high.The standard errors and RMSE of β 0 are always larger than other parameters.The coverage probabilities of the confidence intervals based on the asymptotic normality of the MLEs are quite close to the nominal level in most of the cases.To summarize, a larger sample size, smaller censoring proportion, and lower cure rate would result in more accurate estimates.

Model discrimination
The COM-Poisson distribution contains many commonly used discrete distributions under different selection of ϕ.It would be of interest to select the suitable ϕ and make full use of the COM-Poisson distribution to get the best fit for the data.So, we focus here on a model discrimination among the three special cases of the COM-Poisson distribution.

Likelihood-based method
We consider a likelihood ratio test for the null hypothesis H 0 that the competing risk follows one of the three special cases of COM-Poisson distribution, namely, Geometric ϕ = 0, Poisson ϕ = 1, and Bernoulli ϕ −→ ∞ versus the alternative hypothesis H a that the competing risk follows the COM-Poisson distribution.The test statistic is taken as Λ = −2( l0 − l), where l0 and l are the values of the maximized log-likelihood function under the null and alternative hypotheses, respectively.The asymptotic distribution of the test statistic Λ, under H 0 : ϕ = 1 follows a χ 2 distribution with one degree of freedom.However, the boundary distribution of the test statistic Λ when ϕ = 0 (Geometric) and ϕ −→ ∞ (Poisson) has a mixture distribution of χ 2 0 and χ 2 1 distributions such that , where χ 2 0 is chi-square distribution with 0 degrees of freedom and χ 2 1 is the chi-square distribution with one degree of freedom.
The values of ϕ used in the profile likelihood approach for the COM-Poisson distribution are {0,0.25,0.5,2/3,1,1.5,2,4,∞}.Figure 1 provides the histograms of the test statistics Λ on the Poisson cure rate model with proportional odds assumption under log-logistic baseline when sample size equals to 400 and 4000 over 1000 generated datasets.These plots also display the probability density of chi-square distribution with one degree of freedom, and 90% quantiles.The histogram of Λ is not close to the asymptotic distribution when sample size is small while they become close as the sample size increases.This suggests that a parametric bootstrap method would be better when the sample size is small, and we will describe this method in detail in the illustrative example section.Incidentally, the same was observed under the proportional odds assumption with Weibull baseline.

Information-based method
We use Akaike information criterion (AIC) and Bayesian information criterion (BIC) for model selection among Geometric, Poisson, and Bernoulli distribution cure rate models.The AIC and BIC are given by where k is the number of model parameters to be estimated, l is the maximized likelihood value, and n is the sample size.We select the model with the smallest value of AIC or BIC.AIC and BIC would give us the same selection rate since k = 5, n = 400 or 800 are the same in each of the scenarios among Geometric, Poisson, and Bernoulli distributions, i.e., the model with the largest l is the model that fits the data best.We examine the total relative bias (TRB) and total root mean square error (TRMSE) due to misspecification of the cure rate model.TRB is the sum of the absolute bias of the estimated cure rates to that of the true cure rates for each of the four groups.Similarly, TRMSE is the sum of the absolute MSE of the estimated cure rates.TRB and TRMSE due to misspecification is    10 and 11).
Table 12 shows that the selection rates under AIC or BIC for the correct models increase as the sample size increases, while it decrease as censoring rate increases, and the selection rates for the correct models are always the highest among all the cases.If the true model is under Weibull baseline, the rate to select log-logistic as the true baseline is low, and the rate becomes even lower for large sample size, or light censoring.

Illustration with melanoma data
In this section, we consider a cutaneous melanoma (a type of manlignant cancer) data to illustrate the performance of the proposed methodology.The data was first introduced by [13], and subsequently studied by many authors including [2], [4], [5], [6], [1], [21].These data were taken from [11], and were originally used to detect the prospective treatment performance on the high-dose interferon alfa-2b therapy in order to prevent the recurrence of the disease.The study included 427 patients in total from years 1991 to 1995 and follow up until year 1998.Among them, 10 patients were removed in our analysis due to the missingness of the tumor thickness data.The overall percentage of censored observations is 55.6%.The mean and standard deviation of the observed lifetimes are 3.18 and 1.69 in years, respectively.We choose the nodule categories based on the tumor thickness as the only covariate.The subjects were therefore divided into four different categories (x = 0, 1, 2, 3), with corresponding sample sizes n 1 = 111, n 2 = 137, n 3 = 87, n 4 = 82.The percentage of censored observations for each group are 67.57%,61.31%, 52.87%, 32.93%.See Figure 2 for the lifetimes of susceptibles.The initial values for the EM algorithm are chosen in the following way.We consider the censored rate as the over-estimated cured rate of groups one and four, and then calculated β 0 and β 1 .The initial guess for γ γ γ are estimated from the linear relationship between Nelson-Aalen estimates of log odds and log t, i.e., for log-logistic and Weibull baseline distributions, respectively.We then fitted these data by Geometric (ϕ = 0), COM-Poisson (ϕ = 0.5), Poisson (ϕ = 1), COM-Poisson (ϕ = 2) and Bernoulli (ϕ ≈ ∞) cure rate model.These models, along with cure rate models proposed by [6] and [1], are compared on the basis of AIC and BIC.From Table 13, we observe that φ ≈ ∞ and φ = 0 provide the maximized log-likelihood values, which means that the Bernoulli and Geometric cure rate models provide a good fit for the data under log-logistic odds and Weibull odds models, respectively.Moreover, the maximized log-likelihood value increases and decreases as ϕ increases From Table 13, we can see that the maximum difference between l is only 0.09, which means that the difference between the maximized likelihood value is e 0.09 = 1.09 among different ϕ's.Table 14 presents the estimates, standard errors and 95% CI of the cure rates stratified by nodule category.The confidence interval for all the models do not overlap for the first and forth nodule categories.Table 15 presents the MLES and their standard errors for different cure rate models under proportional odds assumption.It is to be noted that l is very close to each other under proportional odds model with log-logistic baseline even-though the estimates are quite different.The same behaviour was also seen in the model discrimination section.
In order to further investigate the effect of ϕ under the COM-Poisson distribution, we fix ϕ from 0 to 5 with an increment of 0.1 and evaluate the maximum log-likelihood value for each ϕ through likelihood approach.And we test the null hypothesis H 0 : ϕ = ∞ vs. H 1 : 0 ≤ ϕ < ∞ using the likelihood ratio test for the log-logistic odds baseline.And H 0 : ϕ = 0 vs. H 1 : ϕ > 0 under the Weibull odds baseline.The test statistic is given by Λ = −2( l0 − l). Figure 3 shows that the likelihood ratio test statistic decreases and increases as ϕ increases for log-logistic and Weibull odds, respectively, which suggest that the maximized likelihoods increase and decrease as ϕ increases.As we mentioned during the model discrimination, the asymptotic distribution is not suitable when the sample size is small.So, we use a bootstrap method to obtain the distribution of the likelihood ratio test statistic Λ.We generated 1000 samples from Geometric, Poisson, Bernoulli cure rate models under proportional odds model with log-logistic and Weibull distributions, respectively.For each of the dataset, we fit the true cure rate model as well as the COM-Poisson cure rate model, then we calculate the values of Λ.The histograms of Λ are given in Figure 4.The p-value is the proportion of times Λ greater than the corresponding value determined from the data.We obtained p-values of 0.142, 0.132 and 0.681 if we test for Geometric, Poisson and Bernoulli cure rate models with log-logistic odds.Also, we obtained p-values of 0.599, 0.001 and 0.000 if we test for Geometric, Poisson and Bernoulli cure rate models under Weibull odds.Moreover, it would be of interest to get an acceptable range of ϕ if we are using the Weibull baseline for the proportional odds model.Figure 3 present the values of Λ against ϕ with  [6].Section * * is taken from [1].We may reject the null hypothesis H 0 : ϕ = 0 with 10% level of significance if Λ is greater than 0.16.This implies that ϕ ∈ [0, 0.2), and the Geometric model under Weibull odds adequately fits the data.We also set up a test on the effect of the proportional odds parameter: H 0 : γ 2 = 0 as null hypothesis vs. H 1 : γ 2 ̸ = 0 as alternative hypothesis for Geometric, Poisson, Bernoulli cure rate models with ϕ = 0, 1, ∞ under Weibull (log-logistic) baseline.Note that the covariate or the nodule categories would not affect the analysis if γ 2 = 0, and the lifetime would just follow a Weibull (log-logistic) distribution.The test statistic turned out to be 3.194 (14.95), 10.414 (19.564), 17.767 (25.99) with corresponding p-values 0.0739(0.00011),0.00125 (9.73 × 10 −6 ), 0.000025 (3.4 × 10 −7 ).Most of the p-values were less than 0.05, which shows that the proportional odds model provides a better fit than a constant lifetime model over the four nodule categories.
Deviance residual is examined to check the error, which is defined as

Model diagnosis
In this section, we check if the proportional odds model with log-logistic baseline assumption on the lifetime is met for the cutaneous melanoma data that we have discussed in the previous section.From (40), we can see that log Ô(t (i) ) and logt (i) should have a linear relationship between them.For the cutaneous melanoma data, we calculated the cumulative hazard function Ĥ(t) as the observed hazard through the non-parametric Nelson-Aalen estimator, and then get the log-odds function accordingly as where d i and n i are the number of events and total individuals at risk at time t i , respectively.Figure 9 presents the scatter plot of log Ô(t (i) ) vs. logt (i) for each category.The plot shows almost a linear relationship among the four groups.It is to be noticed that there is an intersection among x = 2 and x = 3 for short survival times which violates our proportional odds assumption.Figure 10 shows the difference of log-odds between nodule categories 1, 2, 3 and baseline 0. In this figure, the log odds for each of the nodule categories are calculated by using the linear interpolation within the range of discrete points from 0.4 to 2. From our proportional odds assumption, logO − logO 0 = x x xγ γ γ 2 , we know that the difference should be a linear horizontal line and does not depend on the time.However, the lines in Figure 10 look parallel but do show give a little curvature.Since our data include the cured individuals and are not independent, linear regression test may not be good for model diagnosis.
We use the parametric bootstrap and Monte Carlo methods to develop a goodness of fit test to check whether the Bernoulli cure rate model with proportional odds assumption under log-logistic baseline is sufficient.The critical region for this test will be to the left.We simulated 1000 data based on Bernoulli cure rate model with proportional odds survival assumption under log-logistic baseline.The parameters are β 0 = −0.413,β 1 = 0.379, γ 0 = 2.461, γ 1 = 2.266, γ 2 = −0.473.The censored proportion for the nodule categories are (0.676, 0.613, 0.529, 0.329), respectively.We calculated the values of maximum likelihood l1 , ..., l1000 for each of these generated datasets, and order them ( l) 1 , ..., l(1000) .Then, we determine the proportion of times the l is smaller than the maximum likelihood we obtained from the data as -506.342.Figure 11 presents the histogram of the log-likelihood over 1000   simulated datasets.We obtained the p-value for this test as 0.895, which indicates that our model is quite good as suitable for these data.In this paper, we develop a flexible COM-Poisson cure rate model under a proportional odds assumption for the lifetime distribution of susceptibles with the baseline function being that of a Weibull distribution or loglogistic distribution.An EM algorithm is developed for the maximum likelihood estimation of the parameters from the proposed cure rate model.We perform an extensive Monte Carlo simulation study by varying sample sizes, censoring proportion, cure rates, parameters in different distributions to evaluate the performance of our proposed methodology.Overall, our methodology provides accurate estimates of the model parameters as well as of the cure rates.Moreover, a real data on cutaneous melanoma is analyzed and model diagnosis is performed for illustrative purpose.There are many potential future works in this direction.One may consider the use of a non-parametric specification of the baseline distribution in the proportional odds model of the lifetimes of susceptible group.In addition, a destructive cure survival rate model, by including a damage or destruction term for the initial risk factors, can also be considered in the context of proportional odds model.

Acknowledgements
Our sincere thanks go to the reviewers and the guest editors for their useful comments and suggestions on an earlier version of this manuscript which led to this improved version.

Results for the Special cases of COM-Poisson cure rate model
As mentioned earlier, the COM-Poisson distribution includes the Bernoulli, Poisson and Geometric distributions as special cases.Here, we detail the steps of the EM algorithm for these three special cure models.
Bernoulli cure rate model Let the competing cause random variable M follow a Bernoulli distribution with probability of success η/(1 + η).The probability density function for the whole population can then be expressed as The survival function for the susceptible group is just the survival function for the time to event W , i.e., S s (t i ; θ θ θ) = S(t i ; γ γ γ).The inverse of the cure rate under this setting is 1/p 0 = 1 + η.We, therefore, have H ϕ (η) = 1 + η under the logistic link with a fixed value of ϕ, which implies η = e x x x ′ i β β β .The Q(θ θ θ * , π π π (k) ) function is then given by ∑ It is readily seen that some of the terms in the Q function are only corresponding to β β β while the others are only corresponding to γ γ γ.So, it can be split into two parts as follows: with the update step for the ith censored observation.The required first-and second-order derivatives of Q(θ θ θ * , π π π (k) ) with respect to β β β and γ γ γ are as follows: for l, l ′ = 0, . . ., p, j, j ′ = 1, 2, h = 21, . . ., 2p, i = 1, . . ., n.
Poisson cure rate model Let the competing cause random variable M follow a Poisson distribution.The probability density function for the whole population in this case can be expressed as and the survival function for the susceptible group as The cure rate is p 0 = e −η .We would then have H ϕ (η) = e η under the logistic link with a fixed value of ϕ, which implies that η = ln(1 + e with the update step for the ith censored observation.The required first-and second-order derivatives of Q(θ θ θ * , π π π (k) ) with respect to β β β and γ γ γ are as follows: Geometric cure rate model Let the competing cause random variable M follow a Geometric distribution.The probability density function for the whole population in this case can be expressed as and the survival function for the susceptible group as The cure rate under this setting is p 0 = 1 − η, and under the logistic link with a fixed value of ϕ, we would have with the update step for the ith censored observation.The required first-and second-order derivatives of Q(θ θ θ * , π π π (k) ) with respect to β β β and γ γ γ are as follows: x il e x x x ′ β β β R G (t i , θ θ θ) 2 ∂S(t i , θ θ θ) ∂γ j , for l, l ′ = 0, . . ., p, j, j ′ = 1, 2, h = 21, . . ., 2p, i = 1, . . ., n.Hence, the components of the observed information matrix are for l, l ′ = 0, . . ., p, x i0 ≡ 1, h, h ′ = 0, 1, j * , j * = 21, 22, . . ., 2p, i = 1, . . ., n.

Figure 7 ,
Figure 7, 8 and 5, 6 present the deviance residuals as well as qq plot for various fitted cure rate models.It can be seen that the deviance residuals are distributed around 0, and satisfied the normality assumption.

Figure 4 .
Figure 4.The histogram of Λ = −2( l − l0 ) from 1000 generated datasets with respect to MLEs on Geometric (top), Poisson (middle), Bernoulli (bottom) cure rate model with the lifetime distribution as a proportional odds model with log-logistic(left) and Weibull (right) baseline .

Figure 5 .
Figure 5. QQ plot for deviance residual on proportional odds model with log-logistic baseline for cutaneous melanoma data.

Figure 6 .
Figure 6.QQ plot for deviance residual on proportional odds model with Weibull baseline for cutaneous melanoma data.

Figure 7 .
Figure 7. Deviance residual on proportional odds model with log-logistic baseline for cutaneous melanoma data.

Figure 8 .
Figure 8. Deviance residual on proportional odds model with log-logistic baseline for cutaneous melanoma data.

Figure 9 .
Figure 9. Log odds against log t for cutaneous melanoma data based on Nelson-Aalen estimator on patients who died.

Figure 10 .
Figure 10.Difference of log odds between nodule categories against log t for cutaneous melanoma data based on Nelson-Aalen estimator.

Figure 11
Figure 11.Histogram of l

Table 1 .
Bias, SE, RMSE, and CP for the estimates of the parameters of the Geometric cure rate model under proportional odds with log-logistic baseline.

Table 2
. Bias, SE, RMSE, and CP for the estimates of the parameters of the Poisson cure rate model under proportional odds model with log-logistic baseline.

Table 3 .
Bias, SE, RMSE, and CP for the estimates of the parameters of the Bernoulli cure rate model under proportional odds model with log-logistic baseline.

Table 4 .
Bias, SE, RMSE, and CP for the estimates of the parameters of the Geometric cure rate model under proportional odds model with Weibull baseline.

Table 5 .
Bias, SE, RMSE, and CP for the estimates of the parameters of the Poisson cure rate model under proportional odds model with Weibull baseline.

Table 6 .
Bias, SE, RMSE, and CP for the estimates of the parameters of the Bernoulli cure rate model under proportional odds model with Weibull baseline.

Table 7 .
Bias, SE, RMSE, and CP for the cure rates of the Geometric cure rate model under proportional odds with loglogistic (Weibull) baseline.

Table 8 .
Bias, SE, RMSE, and CP for the cure rates of the Poisson cure rate model under proportional odds with log-logistic (Weibull) baseline.

Table 9 .
Bias, SE, RMSE, and CP for the cure rates of the Bernoulli cure rate model under proportional odds with log-logistic (Weibull) baseline.

Table 10 .
TRB (in %) in the estimation of cure proportions when fitting different models for a given true model under proportional odds due to misspecification.

Table 11 .
TRMSE in the estimation of cure proportions when fitting different models for a given true model under proportional odds due to misspecification.

Table 13 .
AIC, BIC and l under different cure rate models * is taken from

Table 14 .
Estimates, standard errors and 95% CI for the cure rates