Semiparametric estimation for accelerated failure time mixture cure model allowing non-curable competing risk

Abstract

The mixture cure model is the most popular model used to analyse the major event with a potential cure fraction. But in the real world there may exist a potential risk from other non-curable competing events. In this paper, we study the accelerated failure time model with mixture cure model via kernel-based nonparametric maximum likelihood estimation allowing non-curable competing risk. An EM algorithm is developed to calculate the estimates for both the regression parameters and the unknown error densities, in which a kernel-smoothed conditional profile likelihood is maximised in the M-step, and the resulting estimates are consistent. Its performance is demonstrated through comprehensive simulation studies. Finally, the proposed method is applied to the colorectal clinical trial data.

Keywords:

• AFT mixture cure model
• competing risk
• EM algorithm

1. Introduction

In some medical studies, like prostate cancer or breast cancer, it is often observed that a substantial proportion of study patients never experience the event of interest, which are treated as cured or non-susceptible subjects. A number of mixture cure models have been proposed in the literature for analysing such survival data. The two-component mixture cure model can capture cured and uncured patients simultaneously. The early approach of univariate mixture cure model was used to model the survival data by Boag (1949), which assumed that a fraction of the patients are cured from cancer and will never experience the event. Since then, parametric and semiparametric mixture cure models have been proposed (Li & Taylor, 2002; Peng & Dear, 2000; Peng, Dear, Denham, 1998; Sy & Taylor, 2000). The proportional hazards (PHMC) model was further discussed by Sy Taylor (2000) and Peng (2003). Besides, Li Taylor (2002), Zhang Peng (2007), Xu Zhang (2009) and Lu (2010) developed the accelerated failure time (AFTMC) model.

Although PHMC model specifies that the effects of the covariates act multipicatively on the hazard function, the AFTMC model regresses the logarithm of the failure time over the covariates, assuming a direct relationship between failure time over the covariates. Li and Taylor (2002) and Zhang and Peng (2007) estimated parameters through EM for the AFTMC model with semiparametric methods, and the large sample properties was discussed in Xu and Zhang (2009). Lu (2010) developed the AFTMC model and the error density was estimated by the kernel method.

The competing risk model is the main tool in explaining more than one event in survival analysis (Crowder, 2001; David & Moeschberge, 1978; Kalbfleisch & Prentice, 2002; Kleinbaum & Klein, 2006). There two common approaches for competing risk data, one is the cumulative incidence function (CIF) (Fine & Gray, 1999; Kalbfleisch & Prentice, 2002; Tai, Machin, White, & Gebski, 2001), the other one is cause-specified hazard or subhazard (Fusaro, Bacchetti, & Jewell, 1996; Gaynor et al., 1993; Klein, 2006; Ohneberg & Schumacher, 2017; Pintilie, 2007). Usually, the marginal hazards equals to the cause-specified hazards when the independence among events of interest are assumed (Gamel, Weller, Wesley, & Feuer, 2000; Kuk, 1992; Larson & Dinse, 1985; Ng & McLachlan, 2003). In mixture model, the marginal hazards are the hazards of single risk (Lambert, Thompson, Weston, & Dickman, 2006; Yu & Tiwari, 2007; Yu, Tiwari, Cronin, & Feuer, 2004).

Due to the scarcity of efficient and reliable computational methods, there is little literature dealing with AFTMC model and competing risk data. In this paper, we assume all patients may experience either event 1 or event 2, where event 1 is the primary event with a potential cure rate and event 2 is all other possible events of interest, which are non-curable. Different from the standard competing risk model, there is a cured fraction existed in event 1. For the patients being cured from event 1, they will only experience event 2; while for the patients being uncured from event 1, they will experience either event 1 or event 2. The rest of the paper is organised as follows. In Section 2, the accelerated failure time mixture cure model allowing non-curable competing risk is developed. An EM algorithm is developed to implement the estimation in Section 3, we maximise a kernel-smoothed conditional profile likelihood in the M-step, which is motivated by Zeng and Lin (2007) in efficient estimation for the accelerated failure time model without cure fraction. We show that the estimates are consistent and efficient in the appendix. In addition, the proposed model and method are evaluated by the comprehensive simulation study in Section 4 and illustrated by the real data analysis in Section 5. Conclusions are made in Section 6.

2. Accelerated failure time mixture cure model allowing non-curable competing risk

Consider a study with n subjects. For subject i, let $T_{1,i}$ denote the event time undergoing the risk of primary interest and $T_{2,i}$ the event time under other risks. In addition, let $X_i$ denote the covariates and $C_i$ the censoring time for subject i. ${\epsilon }_i$ is used as event indicator (${\epsilon }_i=1$ for event of interest; ${\epsilon }_i=2$ for other risk), and define $T_i=T_{j,i}$ when ${\epsilon }_i=j,j=1,2$, ${\tilde{T}}_i=min(T_i,C_i)$, ${\delta }_i=\mathrm{I}({\mathrm{T}}_{\mathrm{i}}\le {\mathrm{C}}_{\mathrm{i}})$. The observed data then consist of $\mathbf{O}=\left\{\right({\tilde{t}}_i,{\delta }_i,{\epsilon }_i,X_i),i=1,\dots ,n\}$. $Y_i$ is the indicator of the uncured patient i subject to the event of primary interest. That is, $Y_i=1$ means that subject i is subjected to the event of primary interest with uncured probability $P(Y_i=1)$. Note that $Y_i$ is not directly observable. As usual, we assume that $C_i$ is given or independent of $(T_{1,i},T_{2,i},Y_i)$ given $X_i$.

As commonly done in the mixture cure models, we assume a logistic regression for $Y_i$, i.e. ${\pi }_{\alpha }\left(X_i\right)\equiv P(Y_i=1|X_i)=\frac{\exp \left({\alpha }^{\mathrm{\prime }}{\tilde{X}}_i\right)}{1+\exp \left({\alpha }^{\mathrm{\prime }}{\tilde{X}}_i\right)},$ where ${\tilde{X}}_i=(1,X_i^{\prime }{)}^{\prime }$.

To model the competing risk data, we adopt a similar mixture regression approach as studied in Lu and Peng (2008). Note that when $Y_i=0$, ${\epsilon }_i\equiv 2$. When $Y_i=1$, we assume a logistic regression model ${\pi }_{\theta }\left(X_i\right)\equiv P({\epsilon }_i=1|X_i,Y_i=1)=\frac{\exp \left({\theta }^{\mathrm{\prime }}{\tilde{X}}_i\right)}{1+\exp \left({\theta }^{\mathrm{\prime }}{\tilde{X}}_i\right)}.$ Then, we have $P({\epsilon }_i=1|X_i)={\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)$ and $P({\epsilon }_i=2|X_i)=1-{\pi }_{\alpha }\left(X_i\right)+{\pi }_{\alpha }\left(X_i\right)\{1-{\pi }_{\theta }(X_i\left)\right\}$. Here, we use a logistic link for both ${\pi }_{\alpha }$ and ${\pi }_{\theta }$, but other link functions, such as the log-log and probit links, can also be used to model ${\pi }_{\alpha }$ and ${\pi }_{\theta }$.

In addition, we consider a accelerated failure time regression model for $T_i$ given $X_i$ and ${ϵ}_{j,i}=j$, i.e. (1) $\begin{array}{rl}\log T_{j,i}& ={\beta }_j^T{X}_i+{ϵ}_{j,i},j=1,2,\end{array}$(1) (2) $\begin{array}{rl}{f}_j\left(t|X_i\right)& =\exp (-{\beta }_j^T{X}_i)f_{{ϵ}_j}\left(e^{\log \left(t\right)-{\beta }_j^T{X}_i}\right)\mathrm{or}{S}_j\left(t|X_i\right)\\ & =S_{{ϵ}_j}\left(e^{\log \left(t\right)-{\beta }_j^T{X}_i}\right),j=1,2,\end{array}$(2) where ${\beta }_j$ is a row vector of unknown parameters, ${ϵ}_{j,i}$ denotes the error term for the event j with an unknown survival function, $f_{{ϵ}_j}$ and $S_{{ϵ}_j}$ are, respectively, the density and survival function of $e^{{ϵ}_{j,i}}$.

Define the cumulative incidence functions (CIFs) as (3) $F_j(t,X_i)=P(T_i\le t,{\epsilon }_i=j|X_i),j=1,2.$(3) Then, we have $F_1(t,X_i)={\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)\{1-S_1(t|X_i\left)\right\}$ and $F_2(t,X_i)=[1-{\pi }_{\alpha }(X_i\left){\pi }_{\theta }\right(X_i\left)\right]\{1-S_2(t|X_i\left)\right\}.$ Moreover, (4) $\begin{array}{rl}& S\left(t|X_i\right)\equiv P(T_i>t|X_i)\\ & =\sum _{j=1}^{2}P(T_i>t|{\epsilon }_i=j,X_i)P({\epsilon }_i=j|X_i)\\ & =\{1-{\pi }_{\alpha }(X_i\left)\right\}{S}_2\left(t|X_i\right)+{\pi }_{\alpha }\left(X_i\right)\left[{\pi }_{\theta }\left(X_i\right)S_1\left(t|X_i\right)\right.\\ & \left.+\{1-{\pi }_{\theta }(X_i\left)\right\}{S}_2\left(t|X_i\right)\right]\\ & ={\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)S_1\left(t|X_i\right)\\ & +\left[1-{\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)\right]S_2\left(t|X_i\right).\end{array}$(4) We see that it is just a mixture model with weights ${\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)$ and $1-{\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)$ on the model of event 1 (primary) and event 2 (others). When only one risk was considered, i.e. ${\pi }_{\theta }\left(X_i\right)=1$ and $S_2\left(t|X_i\right)=1$, it becomes the simple mixture cure model $S\left(t|X_i\right)=1-{\pi }_{\alpha }\left(X_i\right)+{\pi }_{\alpha }\left(X_i\right)S_1\left(t|X_i\right)$.

Define $\mathbf{\Theta }=(\alpha ,\theta ,{\beta }_1,{\beta }_2,f_{{ϵ}_1}(\cdot ),f_{{ϵ}_2}(\cdot \left)\right)$. With known ${\epsilon }_i$ and ${\delta }_i$, $i=1,2,\dots ,n$, the likelihood function can be written as (5) $\begin{array}{rl}{L}_o& =\prod _{i=1}^n\left\{{\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)\right.\\ & {\left.\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)\right\}}^{\mathrm{I}({\epsilon }_i=1){\delta }_i}\\ & \times \left(\phantom{\exp (-{\beta }_2^T{X}_i)f_{{ϵ}_2}\left(e^{R_i\left({\beta }_2\right)}\right)}[1-{\pi }_{\alpha }(X_i)+{\pi }_{\alpha }(X_i\left)\right\{1-{\pi }_{\theta }\left(X_i\right)\left\}\right]\right.\\ & {\left.\exp (-{\beta }_2^T{X}_i)f_{{ϵ}_2}\left(e^{R_i\left({\beta }_2\right)}\right)\right)}^{\mathrm{I}({\epsilon }_i=2){\delta }_i}\left\{S\right({\tilde{t}}_i|X_i){\}}^{1-{\delta }_i},\end{array}$(5) where $R_i\left({\beta }_j\right)=\log \left({\tilde{t}}_i\right)-{\beta }_j^T{X}_i,j=1,2$. Direct maximisation of the above observed likelihood function is very challenging due to its complex structure. In the next Section, we derive an EM algorithm to maximise the complete likelihood function.

3. EM algorithm

Similar to the standard mixture cure mode, if $Y_i$ could be observed, then the complete likelihood function is written as (6) $\begin{array}{rl}{L}_c& =\prod _{i=1}^n\left\{\phantom{\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)}{\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)\right.\\ & {\left.\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)\right\}}^{\mathrm{I}({\epsilon }_i=1){\delta }_i}\\ & \times \left(\phantom{\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)}\left[\left\{1-{\pi }_{\alpha }\left(X_i\right)\right\}+{\pi }_{\alpha }\left(X_i\right)\left\{1-{\pi }_{\theta }\left(X_i\right)\right\}\right]\right.\\ & {\left.\exp (-{\beta }_2^T{X}_i)f_{{ϵ}_2}\left(e^{R_i\left({\beta }_2\right)}\right)\right)}^{\mathrm{I}({\epsilon }_i=2){\delta }_i}\\ & \times {\left\{{\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)S_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)\right\}}^{\mathrm{I}({\epsilon }_i=1)(1-{\delta }_i)}\\ & \times \left(\phantom{\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)}\left[\left\{1-{\pi }_{\alpha }\left(X_i\right)\right\}+{\pi }_{\alpha }\left(X_i\right)\left\{1-{\pi }_{\theta }\left(X_i\right)\right\}\right]\right.\\ & {\left.S_{{ϵ}_2}\left(e^{R_i\left({\beta }_2\right)}\right)\right)}^{\mathrm{I}({\epsilon }_i=2)(1-{\delta }_i)}\\ & =\prod _{i=1}^n\left[\phantom{\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)}{\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)\right.\\ & {\left.\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)\right]}^{\mathrm{I}({\epsilon }_i=1){\delta }_i}\\ & \times \left(\phantom{\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)}{\left\{1-{\pi }_{\alpha }\left(X_i\right)\right\}}^{1-y_i}{\left[{\pi }_{\alpha }\left(X_i\right)\left\{1-{\pi }_{\theta }\left(X_i\right)\right\}\right]}^{y_i}\right.\\ & {\left.\exp (-{\beta }_2^T{X}_i)f_{{ϵ}_2}\left(e^{R_i\left({\beta }_2\right)}\right)\}\right)}^{\mathrm{I}({\epsilon }_i=2){\delta }_i}\\ & \times {\left\{{\pi }_{\alpha }\left(X_i\right){\pi }_{\theta }\left(X_i\right)S_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)\right\}}^{\mathrm{I}({\epsilon }_i=1)(1-{\delta }_i)}\\ & \times \left(\phantom{\exp (-{\beta }_1^T{X}_i)f_{{ϵ}_1}\left(e^{R_i\left({\beta }_1\right)}\right)}{\left\{1-{\pi }_{\alpha }\left(X_i\right)\right\}}^{1-y_i}{\left[{\pi }_{\alpha }\left(X_i\right)\left\{1-{\pi }_{\theta }\left(X_i\right)\right\}\right]}^{y_i}\right.\\ & {\left.S_{{ϵ}_2}\left(e^{R_i\left({\beta }_2\right)}\right)\right)}^{\mathrm{I}({\epsilon }_i=2)(1-{\delta }_i)}.\end{array}$(6) It is easy to write the logarithm of the complete likelihood function into four components regarding to their unknown parameters. $\begin{array}{rl}& l_c(\alpha ,\theta ,{\beta }_1,{\beta }_2|\mathbf{O}\mathbf{,}\mathbf{y})=l_{c_1}\left(\alpha |\mathbf{O}\mathbf{,}\mathbf{y}\right)+l_{c_2}\left(\theta |\mathbf{O}\mathbf{,}\mathbf{y}\right)\\ & +l_{c_3}\left({\beta }_1|\mathbf{O}\mathbf{,}\mathbf{y}\right)+l_{c_4}\left({\beta }_2|\mathbf{O}\mathbf{,}\mathbf{y}\right),\end{array}$ where (7) $\begin{array}{rl}& l_{c1}\left(\alpha |\mathbf{O}\mathbf{,}\mathbf{y}\right)\\ & =\sum _{i=1}^n\left[\mathrm{I}\right({\epsilon }_i=1)+\mathrm{I}({\epsilon }_i=2\left)y_i\right]\log {\pi }_{\alpha }\left(X_i\right)\\ & +\sum _{i=1}^n\mathrm{I}({\epsilon }_i=2)(1-y_i)\\ & \log \left\{1-{\pi }_{\alpha }\left(X_i\right)\right\},\end{array}$(7) (8) $\begin{array}{rl}{l}_{c2}\left(\theta |\mathbf{O}\mathbf{,}\mathbf{y}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\log {\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right)\\ & +\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\mathrm{y}}_{\mathrm{i}}\log \left\{1-{\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right)\right\},\end{array}$(8) (9) $\begin{array}{rl}& l_{c3}\left({\beta }_1|\mathbf{O}\mathbf{,}\mathbf{y}\right)\\ & =\sum _{i=1}^n\left[\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\delta }_{\mathrm{i}}\left[\log {\mathrm{h}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)\right.\right.\\ & \left.-{\beta }_1^T{X}_i\right]\\ & \left.-\mathrm{I}({\epsilon }_{\mathrm{i}}=1)H_{{ϵ}_1}\left(e^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)\right],\end{array}$(9) (10) $\begin{array}{rl}& l_{c4}\left({\beta }_2|\mathbf{O}\mathbf{,}\mathbf{y}\right)\\ & =\sum _{i=1}^n\left[\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\delta }_{\mathrm{i}}\left[\log {\mathrm{h}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right)\right.\right.\\ & \left.\left.-{\beta }_2^T{X}_i\right]-\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\mathrm{H}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right)\right],\end{array}$(10) and $f_{{ϵ}_j}=h_{{ϵ}_j}\cdot S_{{ϵ}_j}$, here $h_{{ϵ}_j}$ and $H_{{ϵ}_j}$ are the hazard function and cumulative functions of $e^{{ϵ}_j}$, j=1,2 respectively.

Take $Y_i$ as auxiliary variables, and use EM algorithm to find the maximum likelihood estimates of the unknown parameters. The E-step in the EM algorithm computes the conditional expectation of the complete log-likelihood function with respect to three unobserved probabilities $P(Y_i=1,{\epsilon }_i=1|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})$, the probability that uncured patients die from the events of primary interest; $P(Y_i=1,{\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})$, the probability that uncured patients die from other risks; and $P(Y_i=0,{\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})$, the probability that cured patients die from other risks, respectively, where m indicates the m-th step in the EM algorithm. These three probabilities sum to 1 and can, respectively, be given by (11) $\begin{array}{rl}& P(Y_i=1,{\epsilon }_i=1|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=1)P({\mathrm{Y}}_{\mathrm{i}}=1|{\epsilon }_{\mathrm{i}}=1,{\delta }_{\mathrm{i}}=1)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)\mathrm{P}(Y_i=1,{\epsilon }_i=1|{\delta }_i=0)\right|}_{{\mathbf{\Theta }}^{\left(\mathrm{m}\right)}}\\ & =\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=1)+{\left.(1-{\delta }_{\mathrm{i}})\frac{{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right){\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right){\mathrm{S}}_1\left({\tilde{\mathrm{t}}}_{\mathrm{i}}|{\mathrm{X}}_{\mathrm{i}}\right)}{\mathrm{S}\left({\tilde{\mathrm{t}}}_{\mathrm{i}}|{\mathrm{X}}_{\mathrm{i}}\right)}\right|}_{{\mathbf{\Theta }}^{\left(\mathrm{m}\right)}},\end{array}$(11) (12) $\begin{array}{rl}& P(Y_i=1,{\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=2)\mathrm{P}({\mathrm{Y}}_{\mathrm{i}}=1|{\epsilon }_{\mathrm{i}}=2,{\delta }_{\mathrm{i}}=1)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)P(Y_i=1,{\epsilon }_i=2|{\delta }_i=0)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=2)\frac{{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right)[1-{\pi }_{\theta }({\mathrm{X}}_{\mathrm{i}}\left)\right]}{1-{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right){\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right)}\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)\frac{{\pi }_{\alpha }\left(X_i\right)\left\{1-{\pi }_{\theta }\left(X_i\right)\right\}{S}_2\left({\tilde{t}}_i|X_i\right)}{S\left({\tilde{t}}_i|X_i\right)}\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}},\end{array}$(12) (13) $\begin{array}{rl}& P(Y_i=0,{\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=2)\mathrm{P}({\mathrm{Y}}_{\mathrm{i}}=0|{\epsilon }_{\mathrm{i}}=2,{\delta }_{\mathrm{i}}=1)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)P(Y_i=0,{\epsilon }_i=2|{\delta }_i=0)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=2)\frac{1-{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right)}{1-{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right){\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right)}\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)\frac{\left\{1-{\pi }_{\alpha }\left(X_i\right)\right\}{S}_2\left({\tilde{t}}_i|X_i\right)}{S\left({\tilde{t}}_i|X_i\right)}\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}.\end{array}$(13)

Let ${\hat{p}}_{11,i}^{\left(m\right)}=P(Y_i=1,{\epsilon }_i=1|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})$ and ${\hat{p}}_{\epsilon ,1i}^{\left(m\right)}=P({\epsilon }_i=1|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})$, then ${\hat{p}}_{\epsilon ,1i}^{\left(m\right)}={\hat{p}}_{11,i}^{\left(m\right)}$. Let ${\hat{p}}_{12,i}^{\left(m\right)}=P(Y_i=1,{\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})$ and ${\hat{p}}_{02,i}^{\left(m\right)}=P(Y_i=0,{\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})$ and ${\hat{p}}_{\epsilon ,2i}^{\left(m\right)}=P({\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})$, then ${\hat{p}}_{\epsilon ,2i}^{\left(m\right)}={\hat{p}}_{12,i}^{\left(m\right)}+{\hat{p}}_{02,i}^{\left(m\right)}$. The expectations of (7(7) $\begin{array}{rl}& l_{c1}\left(\alpha |\mathbf{O}\mathbf{,}\mathbf{y}\right)\\ & =\sum _{i=1}^n\left[\mathrm{I}\right({\epsilon }_i=1)+\mathrm{I}({\epsilon }_i=2\left)y_i\right]\log {\pi }_{\alpha }\left(X_i\right)\\ & +\sum _{i=1}^n\mathrm{I}({\epsilon }_i=2)(1-y_i)\\ & \log \left\{1-{\pi }_{\alpha }\left(X_i\right)\right\},\end{array}$(7) ), (8(8) $\begin{array}{rl}{l}_{c2}\left(\theta |\mathbf{O}\mathbf{,}\mathbf{y}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\log {\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right)\\ & +\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\mathrm{y}}_{\mathrm{i}}\log \left\{1-{\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right)\right\},\end{array}$(8) ), (9(9) $\begin{array}{rl}& l_{c3}\left({\beta }_1|\mathbf{O}\mathbf{,}\mathbf{y}\right)\\ & =\sum _{i=1}^n\left[\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\delta }_{\mathrm{i}}\left[\log {\mathrm{h}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)\right.\right.\\ & \left.-{\beta }_1^T{X}_i\right]\\ & \left.-\mathrm{I}({\epsilon }_{\mathrm{i}}=1)H_{{ϵ}_1}\left(e^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)\right],\end{array}$(9) ) and (10(10) $\begin{array}{rl}& l_{c4}\left({\beta }_2|\mathbf{O}\mathbf{,}\mathbf{y}\right)\\ & =\sum _{i=1}^n\left[\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\delta }_{\mathrm{i}}\left[\log {\mathrm{h}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right)\right.\right.\\ & \left.\left.-{\beta }_2^T{X}_i\right]-\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\mathrm{H}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right)\right],\end{array}$(10) ) can be written as (14) $\begin{array}{rl}E\left(l_{c_1}\right)& =\sum _{i=1}^n([{\hat{p}}_{\epsilon ,1i}^{\left(m\right)}+{\hat{p}}_{12,i}^{\left(m\right)}]\log \left[{\pi }_{\alpha }\right(X_i\left)\right]\\ & +{\hat{p}}_{02,i}^{\left(m\right)}\log [1-{\pi }_{\alpha }(X_i\left)\right]),\end{array}$(14) (15) $\begin{array}{rl}E\left(l_{c_2}\right)& =\sum _{i=1}^n({\hat{p}}_{\epsilon ,1i}^{\left(m\right)}\log \left[{\pi }_{\theta }\right(X_i\left)\right]\\ & +{\hat{p}}_{12,i}^{\left(m\right)}\log [1-{\pi }_{\theta }(X_i\left)\right]),\end{array}$(15) (16) $\begin{array}{rl}E\left(l_{c_3}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\delta }_{\mathrm{i}}[\log {\mathrm{h}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)-{\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}]\\ & +{\hat{\mathrm{p}}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}{\mathrm{H}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right),\end{array}$(16) (17) $\begin{array}{rl}E\left(l_{c_4}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\delta }_{\mathrm{i}}[\log {\mathrm{h}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right)-{\beta }_2^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}]\\ & +{\hat{\mathrm{p}}}_{\epsilon ,2\mathrm{i}}^{\left(\mathrm{m}\right)}{\mathrm{H}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right).\end{array}$(17) The M-step in the EM algorithm is to maximise (14(14) $\begin{array}{rl}E\left(l_{c_1}\right)& =\sum _{i=1}^n([{\hat{p}}_{\epsilon ,1i}^{\left(m\right)}+{\hat{p}}_{12,i}^{\left(m\right)}]\log \left[{\pi }_{\alpha }\right(X_i\left)\right]\\ & +{\hat{p}}_{02,i}^{\left(m\right)}\log [1-{\pi }_{\alpha }(X_i\left)\right]),\end{array}$(14) ), (15(15) $\begin{array}{rl}E\left(l_{c_2}\right)& =\sum _{i=1}^n({\hat{p}}_{\epsilon ,1i}^{\left(m\right)}\log \left[{\pi }_{\theta }\right(X_i\left)\right]\\ & +{\hat{p}}_{12,i}^{\left(m\right)}\log [1-{\pi }_{\theta }(X_i\left)\right]),\end{array}$(15) ), (16(16) $\begin{array}{rl}E\left(l_{c_3}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\delta }_{\mathrm{i}}[\log {\mathrm{h}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)-{\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}]\\ & +{\hat{\mathrm{p}}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}{\mathrm{H}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right),\end{array}$(16) ) and (17(17) $\begin{array}{rl}E\left(l_{c_4}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\delta }_{\mathrm{i}}[\log {\mathrm{h}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right)-{\beta }_2^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}]\\ & +{\hat{\mathrm{p}}}_{\epsilon ,2\mathrm{i}}^{\left(\mathrm{m}\right)}{\mathrm{H}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right).\end{array}$(17) ) with respect to the unknown parameters α, θ, ${\beta }_1$, ${\beta }_2$, $h_{{ϵ}_1}$ and $h_{{ϵ}_2}$. Maximising (14(14) $\begin{array}{rl}E\left(l_{c_1}\right)& =\sum _{i=1}^n([{\hat{p}}_{\epsilon ,1i}^{\left(m\right)}+{\hat{p}}_{12,i}^{\left(m\right)}]\log \left[{\pi }_{\alpha }\right(X_i\left)\right]\\ & +{\hat{p}}_{02,i}^{\left(m\right)}\log [1-{\pi }_{\alpha }(X_i\left)\right]),\end{array}$(14) ) and (15(15) $\begin{array}{rl}E\left(l_{c_2}\right)& =\sum _{i=1}^n({\hat{p}}_{\epsilon ,1i}^{\left(m\right)}\log \left[{\pi }_{\theta }\right(X_i\left)\right]\\ & +{\hat{p}}_{12,i}^{\left(m\right)}\log [1-{\pi }_{\theta }(X_i\left)\right]),\end{array}$(15) ) with respect to α and θ can be easily carried out using the Newton–Raphson algorithm, and maximising $l_{c_3}$ and $l_{c_4}$ in (16(16) $\begin{array}{rl}E\left(l_{c_3}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\delta }_{\mathrm{i}}[\log {\mathrm{h}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)-{\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}]\\ & +{\hat{\mathrm{p}}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}{\mathrm{H}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right),\end{array}$(16) ) and (17(17) $\begin{array}{rl}E\left(l_{c_4}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\delta }_{\mathrm{i}}[\log {\mathrm{h}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right)-{\beta }_2^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}]\\ & +{\hat{\mathrm{p}}}_{\epsilon ,2\mathrm{i}}^{\left(\mathrm{m}\right)}{\mathrm{H}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right).\end{array}$(17) )with respect to ${\beta }_1$, ${\beta }_2$, $h_{{ϵ}_1}$ and $h_{{ϵ}_2}$ can be carried out utilising the approach proposed by Zeng and Lin (2007).

To find a smooth estimator for $H_{{ϵ}_j}$, or $h_{{ϵ}_j}$, for simplify, we take $H_{{ϵ}_1}$, or $h_{{ϵ}_1}$ as example, $H_{{ϵ}_2}$ and $h_{{ϵ}_2}$ are similar. We start with the simplest case of a piecewise constant $h_{{ϵ}_1}$. To be specific, we partition an interval containing all $e^{\log {\tilde{t}}_i-{\beta }_1^T{x}_i}$ into $J_n$ equally spaced intervals, $0\equiv t_0<t_1<\cdots <t_{J_n}\equiv M$, where M denotes an upper bound for the $e^{\log {\tilde{t}}_i-{\beta }_1^T{x}_i}$ over all possible ${\beta }_1$'s in a bounded set. A piecewise constant $h_{{ϵ}_1}$ takes the form (18) $h_{{ϵ}_1}\left(t\right)=\sum _{k=1}^{J_n}{c}_k\mathrm{I}\left({\tilde{\mathrm{t}}}_{\mathrm{i}}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right).$(18) Then, for any t, (19) $\begin{array}{rl}{H}_{{ϵ}_1}\left(t\right)& =\sum _{k=1}^{J_n}{c}_k({\tilde{t}}_i-t_k)\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\tilde{\mathrm{t}}}_{\mathrm{i}}<t_k\right)\\ & +\frac{\mathrm{M}}{{\mathrm{J}}_{\mathrm{n}}}\sum _{\mathrm{k}=1}^{{\mathrm{J}}_{\mathrm{n}}}{\mathrm{c}}_{\mathrm{k}}\mathrm{I}\left({\tilde{\mathrm{t}}}_{\mathrm{i}}\ge {\mathrm{t}}_{\mathrm{k}}\right).\end{array}$(19) As $\log \left\{\sum _{k=1}^{J_n}{c}_k\mathrm{I}\right(e^{\log {\tilde{t}}_i-{\beta }_1^T{x}_i}\in [t_{k-1},t_k)\left)\right\}=\sum _{k=1}^{J_n}\log c_k\mathrm{I}(e^{\log {\tilde{t}}_i-{\beta }_1^T{x}_i}\in [t_{k-1},t_k\left)\right)$, then (9(9) $\begin{array}{rl}& l_{c3}\left({\beta }_1|\mathbf{O}\mathbf{,}\mathbf{y}\right)\\ & =\sum _{i=1}^n\left[\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\delta }_{\mathrm{i}}\left[\log {\mathrm{h}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)\right.\right.\\ & \left.-{\beta }_1^T{X}_i\right]\\ & \left.-\mathrm{I}({\epsilon }_{\mathrm{i}}=1)H_{{ϵ}_1}\left(e^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)\right],\end{array}$(9) ) can be rewritten as (20) $\begin{array}{rl}& \sum _{i=1}^n-{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}\\ & +\sum _{k=1}^{J_n}\log c_k\left\{\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\right.\\ & \left.\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\right\}\\ & -\sum _{k=1}^{J_n}{c}_k\left\{\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}(e^{R_i\left({\beta }_1\right)}-t_k)\right.\\ & \mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\\ & \left.+\frac{\mathrm{M}}{{\mathrm{J}}_{\mathrm{n}}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right\}.\end{array}$(20) By differentiating with respect to $c_k$, we see that the solution to the score equation of $c_k$ is (21) $c_k=\frac{\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)}{\begin{array}{c}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}(e^{R_i\left({\beta }_1\right)}-t_k)\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\\ +\frac{\mathrm{M}}{{\mathrm{J}}_{\mathrm{n}}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\end{array}}.$(21) After plugging the equations for the $c_k$ into (20(20) $\begin{array}{rl}& \sum _{i=1}^n-{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}\\ & +\sum _{k=1}^{J_n}\log c_k\left\{\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\right.\\ & \left.\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\right\}\\ & -\sum _{k=1}^{J_n}{c}_k\left\{\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}(e^{R_i\left({\beta }_1\right)}-t_k)\right.\\ & \mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\\ & \left.+\frac{\mathrm{M}}{{\mathrm{J}}_{\mathrm{n}}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right\}.\end{array}$(20) ) and discarding the irrelevant component, we obtain the following sieve profile function: (22) $\begin{array}{rl}& l_{c_3}^p\left({\beta }_1\right)\\ & =\sum _{i=1}^n-{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}\\ & +\sum _{k=1}^{J_n}\log \left\{\frac{\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_i=1)\mathrm{I}\left(e^{R_i\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)}{\begin{array}{c}\sum _{i=1}^n{p}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}-{\mathrm{t}}_{\mathrm{k}})\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\\ +\frac{\mathrm{M}}{{\mathrm{J}}_{\mathrm{n}}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\end{array}}\right\}\\ & \times \left\{\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\right\}\\ & -\sum _{k=1}^{J_n}\frac{\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)}{\begin{array}{c}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}(e^{R_i\left({\beta }_1\right)}-t_k)\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\\ +\frac{\mathrm{M}}{{\mathrm{J}}_{\mathrm{n}}}\sum _{\mathrm{i}=1}^n{p}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\end{array}}\\ & \times \left\{\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}(e^{R_i\left({\beta }_1\right)}-t_k)\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\right.\\ & \left.+\frac{\mathrm{M}}{{\mathrm{J}}_{\mathrm{n}}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right\}\\ & =\sum _{i=1}^n-{\delta }_{\mathrm{i}}\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}\\ & +\sum _{\mathrm{k}=1}^{{\mathrm{J}}_{\mathrm{n}}}\left\{\sum _{\mathrm{i}=1}^{\mathrm{n}}{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\right\}\\ & \times \log \left\{\frac{J_n}{nM}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)I\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\right\}\\ & -\sum _{k=1}^{J_n}\left\{\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\right\}\\ & \times \log \left\{\frac{J_n}{nM}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}(e^{R_i\left({\beta }_1\right)}-t_k)\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\right.\\ & \left.+\frac{1}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right\}.\end{array}$(22) Note $l_{c_3}^p\left({\beta }_1\right)$ is not smooth and may have multiple local maxima. Following Zeng and Lin (2007), we further seek a smooth approximation of $l_{c_3}^p\left({\beta }_1\right)$ by the empirical measure. The kernel smoothed approximation of $l_{c_3}^p\left({\beta }_1\right)$ is (23) $\begin{array}{rl}& l^s{c}_3\left({\beta }_1\right)=\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\\ & \log \left\{\frac{1}{{\mathrm{na}}_{\mathrm{n}}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{j}}=1)K\left(\frac{{\mathrm{R}}_{\mathrm{j}}\left({\beta }_1\right)-{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}{{\mathrm{a}}_{\mathrm{n}}}\right)\right\}\\ & -\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\\ & \log \left\{\frac{1}{\mathrm{n}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}{\int }_{-\mathrm{\infty }}^{{\mathrm{R}}_{\mathrm{j}}\left({\beta }_1\right)-{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)/{\mathrm{a}}_{\mathrm{n}}}\mathrm{K}\left(\mathrm{s}\right)\mathrm{ds}\right\}.\end{array}$(23) where $K(\cdot )$ is the kernel function and $a_n$ is the bandwidth. The detail of this derivation can be found in the appendix. The selection of the kernel function and bandwidth can be found in Zeng and Lin (2007). We propose to maximise $l^s{c}_3\left({\beta }_1\right)$ over ${\beta }_1$ and denote the resulting estimator as $\widehat{{\beta }_1}$. Because $K(\cdot )$ is a smooth kernel function, we can use the Newton–Raphson algorithm or other gradient-based search algorithms to calculate $\widehat{{\beta }_1}$. Given $\widehat{{\beta }_1}$, we estimate $h_{{ϵ}_1}\left(t\right)$ by the following kernel-smoothed estimator: (24) ${\hat{h}}_{{ϵ}_1}\left(t\right)=\frac{\frac{1}{n{a}_{n}t}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)K\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left(\widehat{{\beta }_1}\right)-\log \mathrm{t}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_i\left(\widehat{{\beta }_1}\right)-\log t/a_n}K\left(\mu \right)\mathrm{d}\mu }.$(24) The corresponding estimator of $H_{{ϵ}_1}\left(t\right)$ is (25) $\begin{array}{rl}& {\hat{H}}_{{ϵ}_1}\left(t\right)\\ & ={\int }_{-\mathrm{\infty }}^{\log t}\frac{\frac{1}{n{a}_n}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left(\widehat{{\beta }_1}\right)-\mathrm{s}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_i\left(\widehat{{\beta }_1}\right)-\log t/a_n}K\left(\mu \right)\mathrm{d}\mu }\mathrm{d}s.\end{array}$(25) Similar to $h_{{ϵ}_1}$ and $H_{{ϵ}_1}$, the estimator of $h_{{ϵ}_2}$ and $H_{{ϵ}_2}$ are (26) ${\hat{h}}_{{ϵ}_2}\left(t\right)=\frac{\frac{1}{n{a}_{n}t}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=2)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left(\widehat{{\beta }_2}\right)-\log \mathrm{t}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,2i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_i\left(\widehat{{\beta }_2}\right)-\log t/a_n}K\left(\mu \right)\mathrm{d}\mu },$(26) and (27) $\begin{array}{rl}& {\hat{H}}_{{ϵ}_2}\left(t\right)\\ & ={\int }_{-\mathrm{\infty }}^{\log t}\frac{\frac{1}{n{a}_n}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=2)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left(\widehat{{\beta }_2}\right)-\mathrm{s}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,2i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_i\left(\widehat{{\beta }_2}\right)-\log t/a_n}K\left(\mu \right)\mathrm{d}\mu }\mathrm{d}s.\end{array}$(27) The M-step in the EM algorithm is to maximise (14(14) $\begin{array}{rl}E\left(l_{c_1}\right)& =\sum _{i=1}^n([{\hat{p}}_{\epsilon ,1i}^{\left(m\right)}+{\hat{p}}_{12,i}^{\left(m\right)}]\log \left[{\pi }_{\alpha }\right(X_i\left)\right]\\ & +{\hat{p}}_{02,i}^{\left(m\right)}\log [1-{\pi }_{\alpha }(X_i\left)\right]),\end{array}$(14) ), (15(15) $\begin{array}{rl}E\left(l_{c_2}\right)& =\sum _{i=1}^n({\hat{p}}_{\epsilon ,1i}^{\left(m\right)}\log \left[{\pi }_{\theta }\right(X_i\left)\right]\\ & +{\hat{p}}_{12,i}^{\left(m\right)}\log [1-{\pi }_{\theta }(X_i\left)\right]),\end{array}$(15) ), (16(16) $\begin{array}{rl}E\left(l_{c_3}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\delta }_{\mathrm{i}}[\log {\mathrm{h}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right)-{\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}]\\ & +{\hat{\mathrm{p}}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}{\mathrm{H}}_{{ϵ}_1}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\right),\end{array}$(16) ) and (17(17) $\begin{array}{rl}E\left(l_{c_4}\right)& =\sum _{i=1}^n\mathrm{I}({\epsilon }_{\mathrm{i}}=2){\delta }_{\mathrm{i}}[\log {\mathrm{h}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right)-{\beta }_2^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}]\\ & +{\hat{\mathrm{p}}}_{\epsilon ,2\mathrm{i}}^{\left(\mathrm{m}\right)}{\mathrm{H}}_{{ϵ}_2}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_2\right)}\right).\end{array}$(17) ) with respect to the unknown parameters $\mathbf{\Theta }=(\alpha ,\theta ,{\beta }_1,{\beta }_2,H_{{ϵ}_1}(\cdot ),H_{{ϵ}_2}(\cdot \left)\right)$, which can be easily estimated by the Newton–Raphson method using ‘optim’ function in R.

The EM algorithm is described as follows:

Step 0: Given initial values of ${\alpha }^{\left(0\right)}$, ${\theta }^{\left(0\right)}$, ${\beta }_1^{\left(0\right)}$, ${\beta }_2^{\left(0\right)}$ and ${\hat{p}}_{\epsilon ,ji}^{\left(0\right)}=1$. $h_{{ϵ}_1}^{\left(0\right)}\left({\tilde{t}}_i|X_i\right)$ and $h_{{ϵ}_2}^{\left(0\right)}\left({\tilde{t}}_i|X_i\right)$ can be estimated by (24(24) ${\hat{h}}_{{ϵ}_1}\left(t\right)=\frac{\frac{1}{n{a}_{n}t}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)K\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left(\widehat{{\beta }_1}\right)-\log \mathrm{t}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_i\left(\widehat{{\beta }_1}\right)-\log t/a_n}K\left(\mu \right)\mathrm{d}\mu }.$(24) ) and (26(26) ${\hat{h}}_{{ϵ}_2}\left(t\right)=\frac{\frac{1}{n{a}_{n}t}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=2)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left(\widehat{{\beta }_2}\right)-\log \mathrm{t}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,2i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_i\left(\widehat{{\beta }_2}\right)-\log t/a_n}K\left(\mu \right)\mathrm{d}\mu },$(26) ) based on the value of ${\beta }_1^{\left(0\right)}$ and ${\beta }_2^{\left(0\right)}$.

Step 1: Update ${\alpha }^{(m+1)}$ and ${\theta }^{(m+1)}$ by maximising (14(14) $\begin{array}{rl}E\left(l_{c_1}\right)& =\sum _{i=1}^n([{\hat{p}}_{\epsilon ,1i}^{\left(m\right)}+{\hat{p}}_{12,i}^{\left(m\right)}]\log \left[{\pi }_{\alpha }\right(X_i\left)\right]\\ & +{\hat{p}}_{02,i}^{\left(m\right)}\log [1-{\pi }_{\alpha }(X_i\left)\right]),\end{array}$(14) ) and (15(15) $\begin{array}{rl}E\left(l_{c_2}\right)& =\sum _{i=1}^n({\hat{p}}_{\epsilon ,1i}^{\left(m\right)}\log \left[{\pi }_{\theta }\right(X_i\left)\right]\\ & +{\hat{p}}_{12,i}^{\left(m\right)}\log [1-{\pi }_{\theta }(X_i\left)\right]),\end{array}$(15) ). Update ${\beta }_1^{(m+1)}$ and ${\beta }_2^{(m+1)}$ by maxi-mising (23(23) $\begin{array}{rl}& l^s{c}_3\left({\beta }_1\right)=\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\\ & \log \left\{\frac{1}{{\mathrm{na}}_{\mathrm{n}}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{j}}=1)K\left(\frac{{\mathrm{R}}_{\mathrm{j}}\left({\beta }_1\right)-{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}{{\mathrm{a}}_{\mathrm{n}}}\right)\right\}\\ & -\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\\ & \log \left\{\frac{1}{\mathrm{n}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}{\int }_{-\mathrm{\infty }}^{{\mathrm{R}}_{\mathrm{j}}\left({\beta }_1\right)-{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)/{\mathrm{a}}_{\mathrm{n}}}\mathrm{K}\left(\mathrm{s}\right)\mathrm{ds}\right\}.\end{array}$(23) ).

Step 2: Update $h_{{ϵ}_1}^{(m+1)}(\cdot )$ and $h_{{ϵ}_2}^{(m+1)}(\cdot )$ via (24(24) ${\hat{h}}_{{ϵ}_1}\left(t\right)=\frac{\frac{1}{n{a}_{n}t}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)K\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left(\widehat{{\beta }_1}\right)-\log \mathrm{t}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_i\left(\widehat{{\beta }_1}\right)-\log t/a_n}K\left(\mu \right)\mathrm{d}\mu }.$(24) ) and (26(26) ${\hat{h}}_{{ϵ}_2}\left(t\right)=\frac{\frac{1}{n{a}_{n}t}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=2)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left(\widehat{{\beta }_2}\right)-\log \mathrm{t}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,2i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_i\left(\widehat{{\beta }_2}\right)-\log t/a_n}K\left(\mu \right)\mathrm{d}\mu },$(26) ).

Step 3: In the $(m+1)$th iteration, calculate ${\hat{p}}_{\epsilon ,1i}^{\left(m\right)}$, ${\hat{p}}_{12,i}^{\left(m\right)}$, ${\hat{p}}_{02,i}^{\left(m\right)}$ based on ${\alpha }^{\left(m\right)}$, ${\theta }^{\left(m\right)}$, ${\beta }_1^{\left(m\right)}$, ${\beta }_2^{\left(m\right)}$, $h_{{ϵ}_1}^{\left(m\right)}\left({\tilde{t}}_i|X_i\right)$ and $h_{{ϵ}_2}^{\left(m\right)}\left({\tilde{t}}_i|X_i\right)$ from (11(11) $\begin{array}{rl}& P(Y_i=1,{\epsilon }_i=1|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=1)P({\mathrm{Y}}_{\mathrm{i}}=1|{\epsilon }_{\mathrm{i}}=1,{\delta }_{\mathrm{i}}=1)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)\mathrm{P}(Y_i=1,{\epsilon }_i=1|{\delta }_i=0)\right|}_{{\mathbf{\Theta }}^{\left(\mathrm{m}\right)}}\\ & =\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=1)+{\left.(1-{\delta }_{\mathrm{i}})\frac{{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right){\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right){\mathrm{S}}_1\left({\tilde{\mathrm{t}}}_{\mathrm{i}}|{\mathrm{X}}_{\mathrm{i}}\right)}{\mathrm{S}\left({\tilde{\mathrm{t}}}_{\mathrm{i}}|{\mathrm{X}}_{\mathrm{i}}\right)}\right|}_{{\mathbf{\Theta }}^{\left(\mathrm{m}\right)}},\end{array}$(11) ), (12(12) $\begin{array}{rl}& P(Y_i=1,{\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=2)\mathrm{P}({\mathrm{Y}}_{\mathrm{i}}=1|{\epsilon }_{\mathrm{i}}=2,{\delta }_{\mathrm{i}}=1)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)P(Y_i=1,{\epsilon }_i=2|{\delta }_i=0)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=2)\frac{{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right)[1-{\pi }_{\theta }({\mathrm{X}}_{\mathrm{i}}\left)\right]}{1-{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right){\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right)}\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)\frac{{\pi }_{\alpha }\left(X_i\right)\left\{1-{\pi }_{\theta }\left(X_i\right)\right\}{S}_2\left({\tilde{t}}_i|X_i\right)}{S\left({\tilde{t}}_i|X_i\right)}\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}},\end{array}$(12) ), (13(13) $\begin{array}{rl}& P(Y_i=0,{\epsilon }_i=2|\mathbf{O},{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}})\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=2)\mathrm{P}({\mathrm{Y}}_{\mathrm{i}}=0|{\epsilon }_{\mathrm{i}}=2,{\delta }_{\mathrm{i}}=1)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)P(Y_i=0,{\epsilon }_i=2|{\delta }_i=0)\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & ={\left.\delta \mathrm{I}({\epsilon }_{\mathrm{i}}=2)\frac{1-{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right)}{1-{\pi }_{\alpha }\left({\mathrm{X}}_{\mathrm{i}}\right){\pi }_{\theta }\left({\mathrm{X}}_{\mathrm{i}}\right)}\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}\\ & +{\left.(1-{\delta }_i)\frac{\left\{1-{\pi }_{\alpha }\left(X_i\right)\right\}{S}_2\left({\tilde{t}}_i|X_i\right)}{S\left({\tilde{t}}_i|X_i\right)}\right|}_{{\mathbf{\Theta }}^{\mathbf{\left(}\mathbf{m}\mathbf{\right)}}}.\end{array}$(13) ).

Step 4: Repeat Steps 1 and 3 until convergence attained. For convergence, we use the criterion, $max\left\{\right({\alpha }^{\left(m\right)}-{\alpha }^{(m-1)}{)}^2,({\theta }^{\left(m\right)}-{\theta }^{(m-1)}{)}^2,({\beta }_1^{\left(m\right)}-{\beta }_1^{(m-1)}{)}^2,({\beta }_2^m-{\beta }_2^{m-1}{)}^2\}<0.0001$.

The algorithm always converges in our simulation.

4. Simulation

We conducted numerous simulation studies to examine the proposed inference procedures. We generated failure times from the following model: (28) $\log T_j={\beta }_{j,i}{X}_1+{\beta }_{j,i}{X}_2+{ϵ}_{j,i},j=1,2,$(28) where $X_1$ is standard normal distribution, and $X_2$ is Bernoulli with 0.5 success probability. We considered three distribution of $e^{{ϵ}_{j,i}}$: Weibull distribution with different shape parameters, denoted by Weibull$(0.5,1)$ and Weibull$(0.1,1)$; lognormal distribution with different parameters, denoted by lognormal$(0,1)$ and lognormal$(1,1)$; Weibull distribution for the event of interest and lognormal distribution for the other event, which is Weibull$(0.5,1)$ and lognormal $(1,1)$. ${\beta }_1=(1,1)$ for the event of interest for uncured patient and ${\beta }_2=(2,2)$ for other risk event. We assume $\alpha =(2,1,1)$ in the uncure component ${\pi }_{\alpha }$, $\theta =(0.5,-0.5,0.5)$ for the component in the competing risk ${\pi }_{\theta }$, We generated censoring times from the uniform$[0,\tau ]$ distribution, where τ was chosen to produce a $25\mathrm{\%}$ censoring rate. We set n to 100, 500 and 1000.

Following Zeng and Lin (2007), we choose the kernel function $K(\cdot )$ to be the standard normal density for convenience and tractability. We used the optimal bandwidths $4^{1/3}\sigma n^{-1/3}$, where σ is the sample standard deviations of $\log T$.

Table 1 reports the biases, mean squared errors (MSE), empirical standard errors (SE), average estimated standard deviations (ESD) with bootstrap sample size 500 and $95\mathrm{\%}$ coverage probabilities (CP). From Table 1, we can see that all biases are relatively small, the SEs and ESDs are close to each other, and the $95\mathrm{\%}$ CPs are close to their nominal levels. With sample size increases, the biases, MSE and SE become smaller.

Table 1. Estimates of parameters for three baseline survivals.

	n=100					n=500					n=1000
Par.	Bias	MSE	SE	ESD	CP	Bias	MSE	SE	ESD	CP	Bias	MSE	SE	ESD	CP
Baseline-Weibull AFT mixture cure model
${\alpha }_1$	0.0553	0.0488	0.2017	0.2140	0.960	0.0186	0.0091	0.0912	0.0938	0.957	0.0012	0.0047	0.0627	0.0688	0.956
${\alpha }_2$	0.0798	0.0556	0.2194	0.2219	0.916	0.0593	0.0136	0.1000	0.1003	0.976	0.0580	0.0086	0.0708	0.0723	0.962
${\alpha }_3$	0.0319	0.1681	0.3957	0.4088	0.962	−0.0480	0.0332	0.1776	0.1757	0.915	−0.0479	0.0180	0.1117	0.1253	0.956
${\theta }_1$	0.0184	0.0789	0.2694	0.2803	0.968	0.0128	0.0144	0.1201	0.1194	0.971	−0.0039	0.0071	0.0829	0.0843	0.953
${\theta }_2$	0.0626	0.0537	0.2248	0.2231	0.940	0.0657	0.0125	0.0913	0.0906	0.910	0.0688	0.0087	0.0616	0.0633	0.957
${\theta }_3$	0.0154	0.1744	0.4066	0.4174	0.958	−0.0549	0.0343	0.1750	0.1769	0.927	−0.0566	0.0183	0.1216	0.1227	0.933
${\beta }_{11}$	−0.0030	0.0077	0.0786	0.0878	0.968	−0.0061	0.0007	0.0239	0.0258	0.969	−0.0041	0.0003	0.0164	0.0160	0.931
${\beta }_{12}$	−0.0023	0.0270	0.1429	0.1644	0.974	0.0008	0.0026	0.0503	0.0511	0.971	−0.0034	0.0010	0.0326	0.0317	0.938
${\beta }_{21}$	0.0359	0.0214	0.0766	0.1418	0.982	0.0099	0.0004	0.0161	0.0175	0.922	0.0087	0.0002	0.0108	0.0105	0.944
${\beta }_{22}$	0.0348	0.0538	0.1349	0.2293	0.988	0.0186	0.0078	0.0843	0.0863	0.947	0.0190	0.0020	0.0409	0.0411	0.936
Baseline-Lognormal AFT mixture cure model
${\alpha }_1$	0.0173	0.0473	0.2041	0.2167	0.960	-0.0571	0.0123	0.0810	0.0951	0.932	−0.0713	0.0101	0.0580	0.0711	0.952
${\alpha }_2$	0.0411	0.0493	0.2113	0.2182	0.940	0.0451	0.0108	0.0950	0.0934	0.910	0.0463	0.0067	0.0656	0.0677	0.948
${\alpha }_3$	0.0216	0.1643	0.3750	0.4047	0.970	−0.0360	0.0290	0.1511	0.1664	0.950	−0.0506	0.0173	0.1050	0.1213	0.950
${\theta }_1$	−0.0378	0.0905	0.2877	0.2985	0.964	−0.0962	0.0260	0.1168	0.1293	0.914	−0.1170	0.0222	0.0825	0.0924	0.984
${\theta }_2$	0.0041	0.0519	0.2106	0.2278	0.978	0.0382	0.0095	0.0953	0.0899	0.922	0.0424	0.0058	0.0641	0.0630	0.988
${\theta }_3$	0.0061	0.1820	0.4046	0.4266	0.960	−0.0361	0.0315	0.1631	0.1737	0.968	−0.0382	0.0164	0.1172	0.1223	0.950
${\beta }_{11}$	0.0076	0.0268	0.1478	0.1636	0.966	−0.0114	0.0029	0.0506	0.0525	0.954	−0.0160	0.0015	0.0328	0.0346	0.942
${\beta }_{12}$	0.0185	0.0848	0.2864	0.2907	0.938	−0.0128	0.0105	0.0998	0.1016	0.946	−0.0159	0.0047	0.0646	0.0667	0.944
${\beta }_{21}$	0.0270	0.0636	0.1766	0.2507	0.982	0.0109	0.0036	0.0581	0.0589	0.946	0.0104	0.0015	0.0358	0.0377	0.958
${\beta }_{22}$	0.0095	0.1519	0.3073	0.3896	0.966	0.0152	0.0148	0.1097	0.1206	0.954	0.0078	0.0066	0.0704	0.0809	0.958
Baseline-Weilbull-Lognormal AFT mixture cure model
${\alpha }_1$	0.0049	0.0467	0.2039	0.2161	0.922	−0.0403	0.0106	0.0879	0.0946	0.931	−0.0485	0.0073	0.0617	0.0701	0.906
${\alpha }_2$	0.0924	0.0528	0.2097	0.2104	0.922	0.0443	0.0111	0.0949	0.0957	0.916	0.0471	0.0072	0.0677	0.0703	0.914
${\alpha }_3$	0.0536	0.1509	0.4017	0.3847	0.922	−0.0510	0.0311	0.1659	0.1687	0.945	−0.0610	0.0188	0.1149	0.1226	0.936
${\theta }_1$	−0.0471	0.0925	0.2816	0.3004	0.922	−0.0710	0.0211	0.1195	0.1268	0.934	−0.0855	0.0154	0.0868	0.0901	0.954
${\theta }_2$	0.0611	0.0515	0.2230	0.2185	0.910	0.0482	0.0106	0.0885	0.0912	0.916	0.0478	0.0063	0.0603	0.0635	0.978
${\theta }_3$	0.0432	0.1656	0.4533	0.4047	0.912	−0.0535	0.0335	0.1735	0.1751	0.942	−0.0484	0.0174	0.1250	0.1227	0.912
${\beta }_{11}$	0.0058	0.0069	0.0729	0.0830	0.923	−0.0036	0.0007	0.0247	0.0254	0.965	−0.0051	0.0003	0.0156	0.0157	0.942
${\beta }_{12}$	−0.0004	0.0222	0.1365	0.1491	0.923	0.0011	0.0024	0.0489	0.0494	0.958	−0.0036	0.0010	0.0313	0.0310	0.956
${\beta }_{21}$	0.0689	0.1691	0.4272	0.4054	0.923	−0.0002	0.0035	0.0589	0.0594	0.947	−0.0023	0.0015	0.0403	0.0382	0.928
${\beta }_{22}$	0.0400	0.3310	0.3687	0.5740	0.923	0.0020	0.0143	0.1099	0.1196	0.947	−0.0016	0.0049	0.0705	0.0702	0.940

Figures 1– 3 show the estimates of CIFs of the proposed model along with the true CIFs. We compare the CIFs of the event of primary interest ($F_1$) and other risks ($F_2$) with their $95\mathrm{\%}$ pointwise confidence intervals (CI). The estimated CIFs are close to the true values, and it becomes closer as sample sizes increase. All these evidences show that the proposed method performs well, even for small and medium size samples.

Figure 1. Estimated CIFs with Weibull distributions (solid lines), their 95% pointwise confidence intervals (dashed and dotted lines), and the true CIFs (dotdash lines). (a) $n=100,F_1$. (b) $n=100,F_2$. (c) $n=500,F_1$. (d) $n=500,F_2$. (e) $n=1000,F_1$. (f) $n=1000,F_2$.

Figure 2. Estimated CIFs with lognormal distributions (solid lines), their 95% pointwise confidence intervals (dashed and dotted lines), and the true CIFs dotdash lines). (a) $n=100,F_1$. (b) $n=100,F_2$. (c) $n=500,F_1$. (d) $n=500,F_2$. (e) $n=1000,F_1$. (f) $n=1000,F_2$.

Figure 3. Estimated CIFs with Weibull and lognormal distributions (solid lines), their 95% pointwise confidence intervals (dashed and dotted lines), and the true CIFs (dotdash lines). (a) $n=100,F_1$. (b) $n=100,F_2$. (c) $n=500,F_1$. (d) $n=500,F_2$. (e) $n=1000,F_1$. (f) $n=1000,F_2$.

5. Colorectal cancer clinical trial data

To illustrate the proposed estimation method for the AFT mixture cure model with competing risks data, we consider the colorectal cancer clinical trial data from González et al. (2005), which contains rehospitalisation and death data of patients diagnosed with colorectal cancer between January 1996 and December 1998. The data includes calendar time (in days) of the successive hospitalisations after surgical procedure. Gray test with P=0.0006 provides the evidence for considering the competing risk model in this data set, and González et al. (2005) considered death to be a competing risk to rehospitalisation. However, some patients may never have rehospitalisation after discharge. Thus, we consider time to rehospitalisation as the primary interest and death as other possible event which is non-curable, and analyse the data via the proposed model.

This study included 523 patients. The first hospital readmission time was considered as the day between date of surgery and the first rehospitalisation after discharge related to colorectal cancer. There are $23\mathrm{\%}$ patients do not experience both events till the end of study and being treated as the right censoring. We considered drug treatment group (chem=1 for thiotepa drug and 0 for control group) and readmissions for comorbidity (Charlson's index:0, 1–2, 3) for possible risk factors for two events of interest. Thiotepa is a member of the class of alkylating agents, which were among the first anticancer drugs used. Alkylating agents are highly reactive and bind to certain chemical groups found in nucleic acids. These compounds inhibit proper synthesis of DNA and RNA, which leads to apoptosis or cell death. However, since alkylating agents cannot discriminate between cancerous and normal cells, both types of cells will be affected by this therapy. For example, normal cells can become cancerous due to alklyating agents. Thus, thiotepa is a highly cytotoxic compound and can potentially have adverse effects. Consequently, the effects of thiotepa on cancer recurrence and death are not obvious.

For the proposed model, all variables were included in the competing risk event with cured and without cured component. The estimated coefficients and their standard deviations for proposed model are listed in Table 2. Based on the uncure part and competing risk part, both Charlson's index (0) and Charlson's index(1−2) have significant impact. Similarly, both treatment and Charlson's index(1−2) show significant impact on whether patients will experience the rehospitalisation or death.

Table 2. Results for colorectal cancer clinical trial data.

	Proposed model
Parameters	Coef	Sd	P
Uncure part
Intercept	0.6968	0.1460	<0.001
Chem	0.1329	0.3021	0.660
Charlson(0)	1.7539	0.3044	<0.001
Charlson(1−2)	1.5675	0.6664	0.019
Competing part
Intercept	0.5339	0.2392	0.026
Chem	−0.2515	0.4433	0.570
Charlson(0)	2.2674	0.4442	<0.001
Charlson(1−2)	1.7019	0.9262	0.066
PH regression for rehospitalisation
Chem	−0.0991	0.0465	0.033
Charlson(0)	−0.0810	0.0479	0.091
Charlson(1−2)	−1.6663	0.3322	<0.001
PH regression for death
Chem	−0.0991	0.0482	0.040
Charlson(0)	0.0548	0.0581	0.346
Charlson(1−2)	2.2774	0.2298	<0.001

After using the (3(3) $F_j(t,X_i)=P(T_i\le t,{\epsilon }_i=j|X_i),j=1,2.$(3) ), we calculated and plotted the cumulative incidence curves for rehospitalisation and death of colorectal cancer in the Figure 4. Figure 4 shows the competing relationship between the rehospitalisaion and death clearly, and the rehospita-lisaion has a higher rate than death.

Figure 4. Proposed model.

6. Discussions and conclusions

In this paper, we developed a new accelerated failure time mixture cure model allowing non-curable competing risk. Comparing with the existing models, this model can better capture the cure rate of disease than the traditional mixture cure model. The semiparametric estimation based on the EM algorithm can be easily got with the help of existing popular R packages. Although the variance estimation is based on the bootstrap due to the complex structure of model, the comprehensive simulation shows that the performance is reasonable even when the resampling size is small.

Acknowledgements

The research is supported by the Natural Science Foundation of China (Nos. 11271136, 81530086) and the 111 Project of China (No. B14019).

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Notes on contributors

Yijun Wang

Yijun Wang is a PhD candidate in the College of Statistics, East China Normal University, Shanghai, China. Her research interests include Bayesian statistics, reliability statistics and biostatistics.

Jiajia Zhang

Jiajia Zhang is a professor of biostatistics in the Department of Epidemiology and Biostatistics, University of South Carolina, USC. She received her PhD degree from Memorial University. Her professional publications and research interests have focused on survival analysis,semiparametric estimation methods and spatial survival analysis.

Yincai Tang

Yincai Tang is professor of statistics in the College of Statistics, East China Normal University, Shanghai, China. He received his PhD degree from East China Normal University. His professional publications and research interests have focused on lifetime data analysis, degradation data analysis, big data analysis and Bayesian inference.

Appendix

Brief description of approximation for the $l^{\ast }{c}_3\left({\beta }_1\right)$. When $n\to \mathrm{\infty }$, $J_n\to \mathrm{\infty }$, and $J_n/n\to 0$, according to the Donsker theorem, we obtain $\begin{array}{rl}& \frac{1}{n}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}\to \mathrm{E}\left\{\delta \mathrm{I}\right(\epsilon =1\left){\beta }_1^{\mathrm{T}}\mathrm{X}\right\},\\ & \frac{1}{n}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\to \\ & \mathrm{P}\left(\delta =1,\epsilon =1,{\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right),\\ & \frac{J_n}{nM}\sum _{i=1}^n{\delta }_{\mathrm{i}}\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\to \\ & {\left.\frac{\mathrm{d}\mathrm{P}(\delta =1,\epsilon =1,{\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\le \mathrm{s})}{\mathrm{d}s}\right|}_{s=t_{k-1}}.\end{array}$ Using the multiplier central limit theorem of Van Der VaartJon and Wellner (1996), we have $\begin{array}{rl}& max\left|\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}(e^{R_i\left({\beta }_1\right)}-t_k)\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\right.\\ & -\left.\mathrm{E}\left[{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}-{\mathrm{t}}_{\mathrm{k}})\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\right]\right|\\ & =\mathrm{Op}\left(\frac{1}{\sqrt{\mathrm{n}}}\right).\\ & max\left|\frac{1}{n}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)-\mathrm{E}\left[{\mathrm{p}}_{\epsilon ,1i}^{\left(m\right)}\mathrm{I}\left({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right]\right|\\ & =\mathrm{Op}\left(\frac{1}{\sqrt{\mathrm{n}}}\right).\end{array}$ Since $J_n/n\to 0$ as $n\to \mathrm{\infty }$, we can obtain $\begin{array}{rl}max& \left|\frac{J_n}{nM}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}(e^{R_i\left({\beta }_1\right)}-t_k)\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\right.\\ & +\frac{1}{\mathrm{n}}\sum _{i=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\\ & -\mathrm{E}\left[\frac{{\mathrm{J}}_{\mathrm{n}}}{\mathrm{M}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}-{\mathrm{t}}_{\mathrm{k}})\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\right.\\ & \left.\left.+{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right]\right|\to 0.\end{array}$ Uniformly in β and $t_k$. Note, for the last term, we have $\begin{array}{rl}& \mathrm{E}\left[\frac{{\mathrm{J}}_{\mathrm{n}}}{\mathrm{M}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}-{\mathrm{t}}_{\mathrm{k}})\mathrm{I}\left({\mathrm{t}}_{\mathrm{k}-1}\le {\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}<{\mathrm{t}}_{\mathrm{k}}\right)\right.\\ & \left.+{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right]\\ & =O\left(J_n^{-1}\right)+\mathrm{E}\left[{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right]=\mathrm{op}\left(1\right)\\ & +\mathrm{E}\left[{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right].\end{array}$ Then, $\begin{array}{rl}& sup\left|-l_{c_3}^p\left({\beta }_1\right)/n-E\left(\delta \mathrm{I}\right(\epsilon =1\left){\beta }_1^{\mathrm{T}}\mathrm{X}\right)\right.\\ & +\sum _{k=1}^{J_n}\mathrm{P}\left(\delta =1,\epsilon =1,{\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\\ & \times \log {\left.\frac{\mathrm{d}\mathrm{P}(\delta =1,\epsilon =1,{\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\le \mathrm{s})}{\mathrm{d}s}\right|}_{\mathrm{s}={\mathrm{t}}_{\mathrm{k}-1}}\\ & -\sum _{k=1}^{J_n}\mathrm{P}\left(\delta =1,\epsilon =1,{\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\in [{\mathrm{t}}_{\mathrm{k}-1},{\mathrm{t}}_{\mathrm{k}})\right)\\ & \left.\times \log \mathrm{E}\left[{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}\left({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\ge {\mathrm{t}}_{\mathrm{k}}\right)\right]\right|\to 0.\end{array}$ Similar to Zeng and Lin (2007), we choose a kernel function $K(\cdot )$ with bandwidth $a_n$. The theory of kernel estimation indicates that under suitable regularity conditions, $\begin{array}{rl}& \frac{1}{n{a}_n}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)-\log \mathrm{t}}{{\mathrm{a}}_{\mathrm{n}}}\right)\to \\ & {\left.\frac{\mathrm{d}\mathrm{P}(\delta =1,\epsilon =1,\mathrm{R}({\beta }_1)\le \mathrm{s})}{\mathrm{ds}}\right|}_{\mathrm{s}=\log \mathrm{t}}\\ & =\frac{\mathrm{d}P(\delta =1,\epsilon =1,{\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\le \mathrm{t})}{\mathrm{d}t}t\end{array}$ and $\begin{array}{rl}& \frac{1}{n{a}_n}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{\log t}K\left(\frac{R_i\left({\beta }_1\right)-s}{a_n}\right)\mathrm{d}s\to \\ & \mathrm{E}({\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)},\mathrm{I}({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\le \mathrm{t}\left)\right).\end{array}$ Thus we approximate $\frac{\mathrm{d}P(\delta =1,\epsilon =1,e^{R\left({\beta }_1\right)}\le t)}{\mathrm{d}t}/\mathrm{E}({\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)},\mathrm{I}({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\ge \mathrm{t}\left)\right)$ by $\frac{1}{t}\frac{\frac{1}{n{a}_n}\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)-\log \mathrm{t}}{{\mathrm{a}}_{\mathrm{n}}}\right)}{\frac{1}{n{a}_n}\sum _{i=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}{\int }_{\log t}^{\mathrm{\infty }}K\left(\frac{R_i\left({\beta }_1\right)-s}{a_n}\right)\mathrm{d}s}.$ The Kernel-smoothed approximation of the likelihood function is $\begin{array}{rl}& l^s{c}_3\left({\beta }_1\right)=n\left\{\phantom{\mathrm{d}P(\delta =1,\epsilon =1,e^{R\left({\beta }_1\right)}\le t)}-E\left[\delta I\right(\epsilon =1\left){\beta }_1^{T}X\right]\right.\\ & +{\int }_0^{\mathrm{\infty }}\log \frac{\mathrm{d}P(\delta =1,\epsilon =1,e^{R\left({\beta }_1\right)}\le t)/\mathrm{d}t}{\mathrm{E}\left({\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}\mathrm{I}({\mathrm{e}}^{\mathrm{R}\left({\beta }_1\right)}\ge \mathrm{t})\right)}\\ & \left.\mathrm{d}P(\delta =1,\epsilon =1,e^{R\left({\beta }_1\right)}\le t)\right\}\\ & =\sum _{i=1}^n-{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\beta }_1^{\mathrm{T}}{\mathrm{X}}_{\mathrm{i}}-\sum _{i=1}^n{\delta }_{\mathrm{i}}\mathrm{I}({\epsilon }_{\mathrm{i}}=1){\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)\end{array}$

$\begin{array}{rl}& +\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\log \left\{\frac{1}{{\mathrm{na}}_{\mathrm{n}}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\delta }_{\mathrm{i}}\mathrm{I}({\epsilon }_{\mathrm{j}}=1)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{j}}\left({\beta }_1\right)-{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}{{\mathrm{a}}_{\mathrm{n}}}\right)\right\}\\ & -\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\log \left\{\frac{1}{\mathrm{n}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{p}}_{\epsilon ,1\mathrm{i}}^{\left(\mathrm{m}\right)}{\int }_{-\mathrm{\infty }}^{{\mathrm{R}}_{\mathrm{j}}\left({\beta }_1\right)-{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)/{\mathrm{a}}_{\mathrm{n}}}\mathrm{K}\left(\mathrm{s}\right)\mathrm{ds}\right\}.\end{array}$ Discarding the constant term, we obtain the following log-likelihood function, $\begin{array}{rl}& l^s{c}_3\left({\beta }_1\right)=\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\\ & \log \left\{\frac{1}{{\mathrm{na}}_{\mathrm{n}}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\delta }_{\mathrm{iI}}({\epsilon }_{\mathrm{j}}=1)\mathrm{K}\left(\frac{{\mathrm{R}}_{\mathrm{j}}\left({\beta }_1\right)-{\mathrm{R}}_{\mathrm{i}}\left({\beta }_1\right)}{{\mathrm{a}}_{\mathrm{n}}}\right)\right\}\end{array}$ $\begin{array}{rl}& -\sum _{i=1}^n{\delta }_i\mathrm{I}({\epsilon }_{\mathrm{i}}=1)\\ & \log \left\{\frac{1}{n}\sum _{j=1}^n{p}_{\epsilon ,1i}^{\left(m\right)}{\int }_{-\mathrm{\infty }}^{R_j\left({\beta }_1\right)-R_i\left({\beta }_1\right)/a_n}K\left(s\right)\mathrm{d}s\right\}.\end{array}$