Likelihood-Based Inference for Semi-Parametric Transformation Cure Models with Interval Censored Data

1. Introduction

Over recent years, cure rate survival models or cure models have gained much popularity in both medical and statistical fields (Othus et al., 2012; Peng & Taylor, 2014; Pal et al., 2018; Pal & Balakrishnan, 2018b; Majakwara & Pal, 2019; Felizzi et al., 2021; Davies et al., 2021; Pal, 2024). With advent of revolutionary treatment techniques in the drug development industry, chances of patients getting cured from various critical and life-threatening conditions have gone higher. In contrast to regular survival models which rely on the assumption that all subjects under study will inevitably face the event of interest (death, relapse etc.) when observed long enough (Overman & Pal, 2022; Hoang et al., 2023), cure models incorporate a cure probability (or, cure rate) to account for the proportion of subjects who may not face the event of interest at all. The notion was introduced in the works of Boag (1949), and later theorized and validated under a formal model fsramework by Berkson & Gage (1952). In general, cure models assume the population to be a mixture of susceptible subjects (or, non-cured) who will eventually face the event and immune subjects (or, cured) who will not face the event even when they are followed-up for a sufficiently long duration of time. For well-designed prospective studies with higher efficacy of novel therapeutic agents or treatments for critical conditions, the end-of-study survival plot (e.g., the Kaplan Meier plot) is likely to level-off to a value significantly higher than zero after showing a plateau, indicating the existence of cure rate (Sy & Taylor, 2000).

A popular way of statistically modeling survival data that indicate the presence of a cure component is by considering time-to-event random variable $Y=(1-\eta )Y_c+\eta Y_s$, where $\eta =1$ if the subject is susceptible and $\eta =0$ if the subject is cured. Therefore, if $\eta =1,Y=Y_s$ is the susceptible time-to-event, and if $\eta =0,Y=Y_c=\infty$ is the cured time-to-event. Note that, $P(\eta \notin \{\mathrm{0,1}\left\}\right)=0$ and $\pi =P(\eta =0)$ is the cure probability whose estimate is of primary interest. If $S(\cdot )$ denotes the survival function corresponding to susceptible subject, the resultant unconditional survival function (population survival function) of any subject in the population can be expressed as

Here, $\pi$ and $S(\cdot )$ are called the incidence and latency parts of a cure model, respectively. Also, note that $\pi =P(\eta =0)={\text{lim}}_{y\to \infty }{S}_p\left(y\right)$. The cure model in (1), known as the mixture cure model (MCM), has been widely studied under various framework and considerations in survival analysis literature. For example, Farewell studied the MCM by assuming a Weibull distribution for time-to-events corresponding to the susceptible individuals, and linked covariates to the cured probability through logistic-link function (Farewell, 1982). Kuk and Chen linked covariates to the susceptible time-to-events by introducing proportional hazards (PH) structure where a marginal-likelihood approach with expectation-maximization (EM) algorithm was applied to estimate baseline hazard function and regression parameters (Kuk & Chen, 1992). By considering a model structure similar to Kuk and Chen, Sy and Taylor studied a Breslow-type estimator while applying EM algorithm to jointly estimate the incidence and latency parameters (Sy & Taylor, 2000; Sy & Yaylor, 2001). The properties of the MCM as a special case of some flexible cure models have been investigated in details by Balakrishnan & Pal (2012, 2013, 2014, 2015a, b); Balakrishnan et al. (2016). Traditionally, a logistic-link function has been applied to link covariates to the incidence part (Farewell, 1982; Balakrishnan & Pal, 2012, 2015a; Balakrishnan et al., 2017, 2022; Balakrishnan & Barui, 2023). More recently, it has been shown by Pal, Peng, & Aselisewine (2023) and Aselisewine & Pal (2023) that the predictive accuracy of the cure probability can be improved by applying suitable machine learning algorithm. The MCM has also been examined under a semi-parametric set-up by considering additive hazards model with study covariates prone to measurement error (Barui & Yi, 2020).

A more realistic approach to modeling time-to-event data in survival studies marked by cured proportion is to consider a competing causes scenario where it is assumed that an unknown number of competing risk factors affect time-to-events corresponding to susceptible subjects Pal & Balakrishnan (2016, 2017c, a); Pal (2021). An estimate of the probability of a subject whose time-to-event is not influenced by any competing causes is considered as the cure rate. This kind of modeling technique has been studied by Yakovlev & Tsodikov (1996); Ibrahim et al. (2001); Tsodikov et al. (2003); Rodrigues et al. (2009); Balakrishnan & Pal (2016); Balakrishnan et al. (2017); Xie & Yu (2021) and Pal & Aselisewine (2023), to name a few. Considering $M$ as the unobserved number of risk factors competing to produce an event, where $M$ is assumed to follow a Poisson distribution with mean $\varphi >0$, gives rise to the promotion-time cure model (PCM) with population survival function given by

In this case, the cured probability is expressed as $\pi ={\text{lim}}_{y\to \infty }{S}_p\left(y\right)=\text{exp}\{-\varphi \}$, where $S(\cdot )$ denotes the common survival function of the time-to-events corresponding to $M$ competing causes.

A convenient and flexible way to unify the MCM and PCM under a single class of models is by considering a Box-Cox transformation cure model (BCTM) where the Box-Cox transformation is applied on $S_p\left(y\right)$ with transformation parameter $\alpha$. The BCTM was introduced by Yin & Ibrahim (2005a) and later studied by Zeng et al. (2006); Koutras & Milienos (2017); Pal & Balakrishnan (2018a); Milienos (2022); Wang & Pal (2022) and Pal & Roy (2023) under various assumptions and model settings, especially for right censored data. However, based on a detailed and systematic literature search, we did not find any published study that deals specifically with the application of BCTM when the observed data is non-informative interval censored in nature. Analysis of interval censored data under cure model set-up, in general, has been carried out by few researchers. For example, Kim & Jhun (2008); Ma (2010); Xiang et al. (2011); Pal & Balakrishnan (2017b); Wiangnak & Pal (2018); Treszoks & Pal (2022, 2023, 2024) considered fully parametric, semi-parametric and distribution-free frameworks. For survival studies where subjects are monitored or followed-up at scheduled discrete time points during the study period, interval censored data are obtained. For subjects who experience the event of interest, the time interval (with finite limits) in which the event occurs is observed and recorded. For subjects who do not experience the event within the study period, the last visit time is considered as the left limit of the corresponding time interval while the right limit is taken as infinite. Therefore, in this article, we carry out an in-depth analysis of the BCTM in the presence of non-informative interval censored data.

In our proposed model, the survival function $S(\cdot )$ of the lifetime corresponding to the latency part is modeled by assuming the corresponding hazard function $\lambda (\cdot )$ to have a PH structure, where the baseline hazard function ${\lambda }_0(\cdot )$ is estimated by piecewise linear approximation (PLA) in a distribution-free or non-parametric way (Balakrishnan et al., 2022). This is a convenient and robust way of modeling time-to-events associated with the susceptible subjects without making any strong assumptions about their lifetime distribution. The PH structure allows us to account for heterogeneity that is present across subjects with respect to covariates. In the context of PH structure, the baseline hazard function has been estimated by piecewise constant function (Larson & Dinse, 1985) and PLA (Lo et al., 1993). For the PCM, $S(\cdot )$ has been approximated by piecewise constant function (Chen & Ibrahim, 2001). Some semi-parametric cure models, e.g., mixture, promotion-time, geometric and Conway-Maxwell (COM) Poisson cure models have been studied where the baseline hazard function is approximated by PLA (Balakrishnan et al., 2016, 2022). A comparison of the model fit and estimation saccuracy between fully parametric MCM with $S(\cdot )$ corresponding to an exponentiated Weibull distribution and semi-parametric MCM with ${\lambda }_0(\cdot )$ estimated by PLA has also been studied recently (Pal et al., 2022). Typically, piecewise linear function in PLA is characterized by the number of cut-points or knots $(B\in {\mathbb{Z}}^{+})$ and choices of the cut-points denoted by ${\tau }_0<,\dots ,<{\tau }_B$, both of which are user-defined. Similar to smoothing splines, $B$ represents the degree of non-parametricity of the model, meaning models with higher $B$ exhibit greater extent of distribution-free features than models with lower $B$ (Balakrishnan et al., 2016). However, models with higher $B$ suffer from over-fitting and estimation complications associated with models having many parameters. The choices of ${\tau }_b,b=0,\dots ,B$, should be made logically to estimate ${\lambda }_0(\cdot )$ with reasonable accuracy. Therefore, our main contributions to the area of cure models, through this article, can be described by the following: For interval censored data, through a detailed simulation study and a real-life example, the behavior of the BCTM is examined numerically by assuming that the hazard function associated with the latency part follows a proportional hazards structure, and the baseline hazard function is estimated by piecewise linear function in a distribution-free manner. The model parameters of the resulting semi-parametric BCTM are estimated by the usual maximum likelihood (ML) estimation by incorporating the EM algorithm. Most literature implements a profile-likelihood technique on the transformation parameter $\alpha$ of the BCTM (Yin & Ibrahim, 2005b; Yin, 2005; Pal & Balakrishnan, 2018a). However, we carry out the estimation by maximizing the conditional expectation of the log-likelihood function with respect to the all parameters simultaneously (also known as simultaneous-maximization technique).

The rest of the article is arranged in the following manner. In Section 2, the theoretical framework of the proposed model is presented in detail. For interval censored data, the mathematical expressions required to implement the EM algorithm are derived in Section 3. In Section 4, the performance (e.g., robustness, accuracy, etc.) of the proposed model and estimation method on interval censored data is demonstrated by carrying out a simulation study under various parameter setting. Section 5 re-establishes the model performance and usefulness via an application on the interval censored data based on a real-life study on smoking cessation. Finally, we provide some concluding remarks and suggestions for future research in Section 6.

2. Transformation cure model framework

A Box-Cox transformation with transformation parameter $\alpha$ on a variable $u$ is given by

Theoretically, $\alpha$ can take any real value. However, the proposed cure model (BCTM) that unifies MCM and PCM under a single framework can be established with $\alpha \in [0,1]$. Considering $G\left(S_p\left(y\right),\alpha \right)=-\pi \left(\alpha \right)(1-S(y\left)\right)$, from eqn.(3), we obtain the following:

where $\pi (\cdot )$ is linked to the covariate vector $\mathit{z}={\left(1,{\mathit{z}}^{*\top }\right)}^{\top }$ of order $q_1+1$ by

Here, $\pi (\cdot )$ is the part of the model that accounts for the incidence while $S(\cdot )$ accounts for the latency. Note that $\alpha$ is the parameter that regulates the structure of the cure model. The MCM and PTCM are, therefore, BCTM with $\alpha =1$ and $\alpha =0$, respectively. Therefore, combining eqns.(4) and (5), the population survival function for the BCTM can be expressed as

The discussion on modeling the latency part is as follows. Being inspired by the works of Koutras & Milienos (2017) and Balakrishnan et al. (2022), $S\left(y\right)$ is modeled by a PH structure with the baseline hazard function ${\lambda }_0\left(y\right)={\lambda }_0(y;\mathit{\psi })$ as

where

for $b=1,\dots ,B$ and a set of user-defined cut-points ${\tau }_0<{\tau }_1<\cdots <{\tau }_B$. Eqn. (7) is considered such that ${\text{lim}}_{y\to {\tau }_b}{\lambda }_0(y;\mathit{\psi })={\psi }_b$. In other words, we approximate the baseline hazard function of any susceptible subject in the population with a piecewise linear function defined in eqn. (7) with continuity at the cut-points in the set $\left\{{\tau }_0,{\tau }_1,\dots ,{\tau }_B\right\}$. Further, to ensure that time-to-event for a susceptible individual varies with respect to covariate vector $\mathit{x}$ of order $q_2$, a PH structure is considered for modeling the hazard function $\lambda (y;\mathit{x},\mathit{\psi },\mathit{\gamma })$ in the usual way, i.e., $\lambda (y;\mathit{x},\mathit{\psi },\mathit{\gamma })={\lambda }_0(y;\mathit{\psi })\text{exp}\left\{{\mathit{x}}^{\mathbf{\top }}\mathit{\gamma }\right\}$. Note that, $\mathit{z}$ and $\mathit{x}$ may share same covariate elements or may be the same. This implies that the survival function $S\left(y\right)$ can be expressed as

where

is the cumulative hazard function for susceptible individual, and

Finally, eqn. (6) takes the following form:

Some approaches to select the cut points in $\left\{{\tau }_0,{\tau }_1,\dots ,{\tau }_B\right\}$ are discussed in Sections 3 and 4. The linear function obtained for the interval $\left[{\tau }_{B-1},{\tau }_B\right]$, i.e., $\left\{{\psi }_B-{\tau }_B\left(\frac{{\psi }_B-{\psi }_{B-1}}{{\tau }_B-{\tau }_{B-1}}\right)\right\}+\left\{\left(\frac{{\psi }_B-{\psi }_{B-1}}{{\tau }_B-{\tau }_{B-1}}\right)t\right\}$ is extended to the interval $\left[{\tau }_B,\infty \right)$ for any $t>{\tau }_B$. The parameter vector of the model is $\mathit{\theta }=(\alpha ,\mathit{\psi },\mathit{\beta },\mathit{\gamma }{)}^{\top }$ with corresponding parameter space as $\Theta \in \mathbb{D}\subset {\mathbb{R}}^{\left(B+q_1+q_2+3\right)}$. The cure rate of the BCTM model is simply given by

3. Simultaneous-maximization under EM framework for interval censored data

In this Section, we discuss the method of parameter estimation by maximum likelihood approach under the EM framework for our proposed model in the presence of interval censored data. The expected value of the complete data log-likelihood function, given the observed data and current value of the parameter estimates, is maximized with respect to the parameter vector $\mathit{\theta }$, that is, considering all parameters simultaneously, also known as the simultaneous-maximization technique.

In the context of right censored data, earlier studies such as Yin & Ibrahim (2005b), Yin (2005) and Pal & Balakrishnan (2018a) have used profile-likelihood technique on $\alpha$ for the BCTM citing model identifiability as a concern which makes parameter estimates unstable with high standard errors. However, few recent studies such as Wang & Pal (2022) and Pal & Roy (2023) have established accuracy and consistency for the simultaneous-maximization technique for right censored data when the search of $\alpha$ is restricted to [0, 1]. To the best of our knowledge, this article is the first one to propose and carry out simultaneous-maximization under the EM-based BCTM framework with PLA of the baseline hazard function for interval censored data. Since one of our objectives is to establish a flexible family of cure model that unifies both the MCM and PCM under the same umbrella, we focus our search for $\alpha$ only in the interval [0, 1].

For a sample of size $n$, the observed interval censored data can be represented as

For subject $i,l_i$ and $r_i$ are the lower and upper limits (with $l_i<r_i$) of the interval in which the actual time-to-event (unobserved) $y_i$ belongs, and ${\delta }_i$ is the censoring indicator such that ${\delta }_i=1$ if both $l_i$ and $r_i$ take finite values (interval censored), else ${\delta }_i=0$ (right censored). Similarly, ${\mathit{x}}_i$ and ${\mathit{z}}_i$ are the observed covariate vectors corresponding to the latency and incidence parts, respectively, for subject $i$. For interval-censored data as defined above, the observed data log-likelihood function is given by

where ${\Delta }_1=\left\{i:{\delta }_i=1\right\}$ and ${\Delta }_0=\left\{i:{\delta }_i=0\right\}$. However, we have partial information on the cure status ${\eta }_i$ for subject $i$, i.e., ${\eta }_i=1$ if $i\in {\Delta }_1$ and ${\eta }_i$ is missing for $i\in {\Delta }_0$. Therefore, the complete interval-censored data can be represented by ${𝒟}_C=\left\{\left(l_i,r_i,{\delta }_i,{\mathit{x}}_i,{\mathit{z}}_i,{\eta }_i\right):i=1,\dots ,n\right\}$, which includes both observed and unobserved ${\eta }_i\mathrm{s}$. An effective way to handle missing information related to ${\eta }_i$ for $i\in {\Delta }_0$ under maximum likelihood estimation set-up is through the application of the EM algorithm. For the BCTM, the complete data log-likelihood function is given by

where $S_p\left(\cdot ;{\mathit{x}}_i,{\mathit{z}}_i,\mathit{\theta }\right)$ and $\pi \left(\mathit{\theta };{\mathit{z}}_i\right)$ are given in eqns. (11) and (12), respectively, and $S_u\left(y_i;{\mathit{x}}_i,{\mathit{z}}_i,\mathit{\theta }\right)$ is the survival function of susceptible subject, which can be expressed as

Following the steps in the works of Pal & Balakrishnan (2017b) and Wiangnak & Pal (2018), we obtain the ML estimate of $\mathit{\theta }$ by implementing EM algorithm to interval censored data.

4. Simulation study

To illustrate the performance (robustness, accuracy, consistency and asymptotic normality) of the proposed model and estimation method under interval censored data, a detailed study is carried out under different parameter settings based on simulated data. The effects on the performance of our method due to variation in the parameter settings are captured by simulating data with transformation parameter $\alpha =0$ (PCM), 0.5 (BCTM) and 1 (MCM), and $n=200$ (smaller sample size) and 400 (moderate sample size). For $i=1,\dots ,n,{\mathit{x}}_i={\left(x_{1i},x_{2i}\right)}^{\top }$ is the observed covariate vector for the $i$th subject where $x_{1i}$ is generated from a Bernoulli distribution with success probability 0.5, and $x_{2i}$ is generated from a uniform distribution in the interval [0.1, 20.0]. Further, ${\mathit{z}}_i={\mathit{x}}_i$ for $i=1,\dots ,n$, and the true values of ${\beta }_0=0.6,{\beta }_1=-1.5,{\beta }_2=0.1,{\gamma }_1=-1.2$ and ${\gamma }_2=0.1$ are chosen conveniently. Assuming the true baseline hazard function to be a real constant $(\zeta >0)$, corresponding to which the simulated data is generated, we obtain the common cumulative hazard function associated with the latency as $\Lambda \left(y_i;{\mathit{x}}_i,\mathit{\gamma }\right)=\zeta y_i{e}^{{\gamma }_1{x}_{1i}+{\gamma }_2{x}_{2i}}$ for some $y_i>0$. Since, $S_u\left(.;{\mathit{x}}_i,{\mathit{z}}_i,{\mathit{\theta }}^{*}\right)$ in eqn. (16) is a proper survival function, it takes values in [0, 1] and is uniformly distributed. Therefore, generating $u_i$ from a uniform distribution [0, 1], we may obtain $y_i$ by solving eqn. (16) for true value of ${\mathit{\theta }}^{*}={\left(\alpha ,\zeta ,{\beta }_0,{\beta }_1,{\beta }_2,{\gamma }_1,{\gamma }_2\right)}^{\top }$. Therefore, the time-to-events $y_i\mathrm{s}$ associated with latency are generated from an exponential distribution with rate $\zeta e^{{\gamma }_1{x}_{1i}+{\gamma }_2{x}_{2i}}$. Further, note that $\mathit{\theta }={\left(\alpha ,{\psi }_0,\dots ,{\psi }_B,{\beta }_0,{\beta }_1,{\beta }_2,{\gamma }_1,{\gamma }_2\right)}^{\top }$ is of order $(B+7)$ corresponding to the fitted model. Censoring times $c_i,i=1,\dots ,n$, are generated from an exponential distribution with rate parameter ${\zeta }^{*}>0$ and $t_i=\text{min}\left\{y_i,c_i\right\}$ is set. If $u_i\le \pi \left(\mathit{\theta };{\mathit{z}}_i\right)$, then $l_i=c_i,r_i=\infty$ and ${\delta }_i=0$. If $u_i>\pi \left(\mathit{\theta };{\mathit{z}}_i\right)$ and $t_i=c_i$, then $l_i=c_i,r_i=\infty$ and ${\delta }_i=0$. Finally, if $u_i>\pi \left(\mathit{\theta };{\mathit{z}}_i\right)$ and $t_i=y_i$, then ${\delta }_i=1$. In this case, $d_1$ and $d_2$ are generated from uniform distribution in [0.2, 0.7] and [0, 1], respectively. If $t_i\in \left(0,d_2\right)$, we set $l_i=0$ and $r_i=d_2$. If $t_i\ge d_2$, we set $l_i=d_2+(\rho -1)d_1$ and $r_i=d_2+\rho d_1$, where $\rho =⌊\frac{t_i-d_2}{d_1}⌋+1$ (Pal & Balakrishnan, 2017b). $\zeta$ is chosen as 0.1 or 0.2 with ${\zeta }^{*}=0.1$ to maintain censoring rate of 20% – 30% for every simulated scenario.

The performance of our method is assessed in terms of measures related to estimation accuracy and consistency. These include parameter estimate (EST), standard error (SE), bias in estimation (BIAS), root mean squared error (RMSE), coverage probability (CP) at a conventional nominal level of 95%. Further, model fit, and hence, performance of the BCTM model by incorporating PLA is assessed in terms of maximized log-likelihood value $(\widehat{l)}$ and Akaike information criterion (AIC). The above-mentioned measures are calculated based on 400 Monte-Carlo simulation runs. An initial guess for $\alpha$ is taken as 0.5 since the search region for $\alpha$ is [0, 1]. The initial guesses for ${\beta }_0,{\beta }_1,{\beta }_2,{\gamma }_1$ and ${\gamma }_2$ are generated randomly from an interval with 10% margin of error from their true values. For the simulation study, the cut-points or knots, ${\tau }_0,\dots ,{\tau }_B$, are considered as the equi-proportion quantiles of the finite lower and upper limits of the intervals. That is, ${\tau }_0=0,{\tau }_B=\text{max}\left\{S_B\right\}$ where $S_B=\left\{l_i,r_i:l_i<\infty ,r_i<\infty ,i=1,\dots ,n\right\}$. For any $b\in \{1,\dots ,B-1\}$, ${\tau }_b$ is taken as the $⌊\frac{b}{B}⌋\text{th}$ quantile of $S_B$. First, all ${\psi }_b,b=0,\dots ,B$, are taken as 0.001. Then, for $b=1,\dots ,B,{\psi }_b$ is set to ${\psi }_b\leftarrow :{\psi }_{b-1}+{\psi }_b+{\psi }_b^2$. These values provide initial estimates of ${\psi }_0,\dots ,{\psi }_B$. The choices are considered keeping in mind that the hazard values may gradually and non-linearly increase with time at every cut-point.

For all six scenarios described above, i.e, $\left(\alpha ,n\right)=\left(0.0,200\right),\left(0.0,400\right),\left(0.5,200\right),\left(0.5,400\right),(1.0,200),(1.0,400)$, for which data are simulated, the BCTM models with PLA-based baseline hazard function with $B=1,2,3$ and 4 are fitted. The simulation study results are presented in Tables 1–6. The Monte-Carlo simulation run based estimates are quite accurate with relatively similar standard error values (relative with respect to similar studies on cure models) for most parameters. The standard errors associated with ${\beta }_0$ are high compared to other parameters, most likely because ${\beta }_0$ is an intercept parameter which is susceptible to small fluctuations in data. The baseline hazard function is estimated by the PLA function in the time interval $\left[{\tau }_0,{\tau }_B\right]$ where often the Euclidean distance between ${\tau }_{B-1}$ and ${\tau }_B$ is large. Therefore, we may notice that standard error associated with ${\psi }_B$ is almost 3 – 6 folds high when compared with respect to other ${\psi }_b,b=0,\dots ,B-1$, owing to lesser accuracy of the PLA estimation at the last cut-point. Further, the simulated data is generated from an exponential baseline distribution corresponding to the latency part of the model. Thus, the simulated data is generated with true time-to-event values in the interval $[0,\infty )$ already influenced by our interval censoring schemes define earlier, and the baseline hazard function is estimated by the PLA function in the partial time interval $\left[{\tau }_0,{\tau }_B\right]$ where ${‖{\tau }_B-{\tau }_{B-1}‖}_2$ is large. These are some of the factors likely to be resulting relatively large value of SE corresponding to ${\psi }_B$. This phenomenon is along the lines of previous findings by Balakrishnan et al. (2022). The SE corresponding to $\alpha$ is high, noting that it can only take values in [0, 1]. Moreover, in some settings, especially for true $\alpha =0.5$, the CPs show under-coverage (70.1863%) even when average of the MLEs of $\alpha$ show closeness to 0. 5 and SE is as high as 0.7806. These can be attributed to the flatness of the likelihood function surface with respect to $\alpha$ resulting in fluctuating and inconsistent ML estimate for $\alpha$ with high variability. However, CPs corresponding to $\alpha$ are closer to the nominal level when data are generated from the PCTM $(\alpha =0.0)$ or the MCM $(\alpha =1.0)$. The SE for $\alpha$ is comparatively smaller when data are generated from $\alpha =0.0$ and higher when data is generated from $\alpha =1.0$. Even though for cases under interval-censored data with true $\alpha =0.5$ the CPs show under-coverage, we would like to emphasize that the BIAS in the estimation of $\alpha$ are not high, the ML estimates related to other parameters show satisfactory results (with respect to all measures considered) and the running time for the entire Monte-Carlo simulation runs (with 400 replications) is much less in comparison to the profile-likelihood approach to estimate $\alpha$. For profile-likelihood based technique, even for a basic grid search with $\alpha =0.0,0.1,\dots ,1.0$, it can take 5 – 10 times greater computation time (running time) than the simultaneous-maximization based approach depending on the values of $\alpha$ and $n$. For better perspective, for interval-censored data generated from true $\alpha =0.5$ and $n=200$, the computation time related to estimation and model fitting is approximately 11.57 minutes by applying the profile-likelihood method with $\alpha \in \{\mathrm{0.0,0.1},\dots ,1.0\}$ while that of the simultaneous-maximization method is 1.49 minutes for a single Monte-Carlo run.

Some obvious trends are observed, e.g., increasing $n$ from 200 to 400 results in the overall decrease in SE, BIAS and RMSE, while the CPs get closer to the nominal level for almost all settings. Finally, since our proposed BCTM with $B=b$ is nested in the BCTM with $B=b^{\prime }$ for $b<b^{\prime }$, the maximized log-likelihood value increases as $b$ increases. As expected, this trend has been captured in all simulation settings. Due to computational challenges, we did not continue with PLAs for $b\ge 5$.

5. Analysis of real data from a smoking cessation study

The performance of our proposed model and estimation method is further demonstrated through an application on interval censored data obtained from the well-known smoking cessation study (Murray et al., 1998). The main objective of the study is to illustrate potency of an aggressive intervention therapy on patients with nicotine addiction. The analysis data set contains intervals (‘Timept1’ and ‘Timept2’ denoting the left and right limits) in which actual times to relapse occur for 223 subjects with history of nicotine addiction. These subjects were randomly assigned to either the smoking intervention (SI) therapy group (denoted by ‘SI.UC’=1) or the usual care (UC) group (‘SI.UC’=0). Available covariate information include duration of smoking in years prior to the study enrollment (‘Duration’) and average numbers of cigarettes smoked per day prior to the study enrollment (‘F10Cigs’). If a subject experienced relapse then the variable ‘Relapse’, present in the data set, takes the value 1, else ‘Relapse’ takes the value 0. ‘Timept2’ is ‘NA’ (which means infinite time to relapse) for ‘Relapse’=0. A descriptive summary of the data with some relevant measures is provided in Table 7.

In our analysis, ‘Duration’ ($x_1$; continuous), ‘F10Cigs’ ($x_2$; continuous) and ‘SI.UC’ ($x_3$; categorical) are considered as covariates. The proposed BCTM model with $B=1,\dots ,6$ are fitted to the smoking cessation data and their model adequacies are compared with respect to the maximized observed log-likelihood value $(\widehat{l)}$ and Akaike information criterion (AIC) (see Table 8). $\text{AIC}=-2\hat{l}+2p_n$ where $p_n$ indicates number of estimated parameters. For user-defined $B$, the initial parameter values, all denoted by notations used earlier with subscript ‘0’, are chosen in the following sequential manner:

1. ${\alpha }_0=0.5$ chosen as the mid-point of the search interval [0, 1].

2. Some logical assumptions are made to obtain ${\mathit{\beta }}_0$. Subjects with higher values of $x_1$ and $x_2$ are less likely to be cured. Similarly, subjects with smaller values of $x_1$ and $x_2$ are more likely to be cured. The observed censoring proportion for the smoking intervention and usual care groups are 0.7337 and 0.6296, respectively. Therefore, a subject with minimum values of $x_1$ and $x_2$ is more likely to have cure probability estimate as 0.7337 (if belongs to ‘SI.UC=1’) or 0.6296 (if belongs to ‘ SI.UC=0’). Note that, censoring proportion is the maximum attainable value for cure probability estimate for a given data. Further, halved cure probability estimates, i.e., 0.3669 (‘SI.UC=1’) or 0.3148 (‘SI.UC=0’) can be assigned to a subject with average values of $x_1$ and $x_2$, respectively and heuristically. Equating (5), with respective covariate input, to the four values of cure probabilities mentioned above and solving lead to ${\mathit{\beta }}_0=(0.5819,0.3131,0.6416,-1.1478{)}^{\top }$.

3. For unique finite time points, $t_1,\dots ,t_p$, let ${\hat{S}}_1,\dots ,{\hat{S}}_p$ be the empirical non-parametric estimates for the interval censored data (Lindsey & Ryan, 1998). Using Euclidean distance measure and some approximations, we obtain ${\mathit{x}}_j={\left(x_{1j},x_{2j},x_{3j}\right)}^{\top }$ for $j=1,\dots ,p$ from the smoking cessation data. By our assumption of proportional hazards structure, we get ${\hat{S}}_j={\left\{S_{0j}\right\}}^{\text{exp}\left\{{\mathit{x}}_j^{\top }\mathit{\gamma }\right\}}$, and hence, the transformed equation obtained is

   $$\text{ln}\left\{-\text{ln}\left\{{\hat{S}}_j\right\}\right\}={\mathit{x}}_j^{\top }\mathit{\gamma }+\text{ln}\left\{-\text{ln}\left\{S_{0j}\right\}\right\}$$

   (21)

   where $S_{0j}$, for $j=1,\dots ,p$, is the corresponding baseline survival function values. Applying multiple linear regression to the data

   $$\left\{\left({\hat{S}}_j,{\mathit{x}}_j\right);j=1,\dots ,p\right\},$$

   we obtain ${\mathit{\gamma }}_0=(-0.2227,-0.0042,0.0185{)}^{\top }$.

4. From eqn. (21), we obtain ${\hat{S}}_{0j}$, and hence, we obtain an estimate of the baseline cumulative hazard function value ${\hat{H}}_{0j}=-\text{ln}\left\{{\hat{S}}_{0j}\right\}$, for $j=1,\dots ,p$, by eliminating the effect of ${\mathit{x}}_j^{\top }{\mathit{\gamma }}_0$. From the set of paired numbers $\left\{\left(t_j,{\hat{H}}_{0j}\right);j=1,\dots ,p\right\}$, we may obtain another set of paired numbers $\left\{\left(t_j,{\hat{h}}_{0j}\right);j=1,\dots ,(p-1)\right\}$ where ${\hat{h}}_{0j}=\frac{{\hat{H}}_{0j\overline{+}1}-{\hat{H}}_{0j}}{t_{j+1}-t_j}$. Here, ${\hat{h}}_{0j}$ is an approximation to the baseline hazard function value at $t_j$.

5. We may plot $\left\{\left(t_j,{\hat{h}}_{0j}\right);j=1,\dots ,(p-1)\right\}$ to decide on our choice of the cut-points $\left\{{\tau }_0,\dots ,{\tau }_B\right\}$ and corresponding set of ${\hat{h}}_{0j}\mathrm{s}$ as initial parameter estimate ${\mathit{\psi }}_0={\left({\psi }_{00},\dots ,{\psi }_{B0}\right)}^{\top }$ for fixed $B.{\tau }_0=0$, and ${\tau }_B$ is taken as the maximum value of the observed finite left and right limits of the time intervals. The cut-points are decided based on the points of inflections of the plotted graph. Other strategy that one may take is by considering suitable quantiles of the observed time interval limit values.

It is observed that the maximized log-likelihood values increase with an increase in $B$ (see Table 8). This model behavior is expected since the cut-points are chosen in a way that the resulting model becomes a part of the family of nested models, i.e., BCTM with $B=b$ is nested within BCTM with $B=b^{\prime }$ for $b<b^{\prime }$. AIC, on the other hand, does not show monotonic behavior necessarily. Findings for the BCTM models with $B>6$ are not presented here since increasing $B$ further does not show significant decrease in AIC. Alongside BCTM models with PLA of the baseline hazard function corresponding to the latency part, BCTM models with Weibull and gamma distribution assumptions of the baseline hazard function are also considered for comparison. The BCTM model with gamma baseline hazard function ($\hat{l}=-194.3569$, AIC=408.7138) provides the minimum AIC value for the given smoking cessation data set. However, note that, the maximum log-likelihood is observed for the BCTM model with PLA $(B=6)$ of the baseline hazard function ($\hat{l}=-194.1805$, AIC=418.3610). As a justification in favor of estimating the baseline hazard function corresponding to the latency part of the BCTM with PLA even though the number of parameter involved is large, we would like to point out that our proposed model is more flexible since it avoids considering a more rigid assumption of any parametric distribution for the baseline hazard function.

The model parameter estimates and corresponding standard errors for the best model among the BCTM models with PLA of the baseline hazard function are presented in Table 9. Note that $x_1$ and $x_2$ show negative relationship while $x_3$ shows a positive relationship with the cure rate (or, in this case, non-relapse rate). Most of the parameter and standard error estimates seem reasonable. However, a higher variability is observed to be associated with ${\psi }_5$, that is for the baseline hazard function value corresponding to the last cut-point. This point has already been discussed in details in Section 4. The ML estimate of $\alpha$ is obtained as 1.0000 (SE = 0.0001). This indicates that for the smoking cessation interval censored data, the MCM with $B=6$ gives the best fit in terms of AIC. Further, it is observed that baseline hazards increase in the intermediate cut-points, however, shows decreasing trend in the later cut-points, except, for the last cut-point which provides the highest hazard value. Findings are also presented in Table 10 for the overall best model among the candidate models with respect to AIC, i.e., the BCTM with gamma baseline hazard function.

Figure 1 represents estimated population survival functions corresponding to our proposed model, and it shows that the group with smoking intervention is likely to have a greater cure rate than the group with usual care.

6. Concluding remarks

In this article, for interval censored data, a flexible family of transformation cure model is proposed that unifies mixture cure model and promotion time cure model. The transformation considered is a Box-Cox transformation taken on the population survival function, and is linked to incidence part and latency part of the cure model. The transformation parameter $\alpha$ provides flexible characteristics to the BCTM, thereby, giving rise to a flexible family of cure models. To introduce non-homogeneity with respect to time-to-events for susceptible subjects, a proportional hazards framework is considered where the baseline hazard function is captured in a distribution-free or non-parametric manner by incorporating a piecewise linear function similar to smoothing splines. Estimation of the model parameters is carried out by implementing the expectation-maximization algorithm to accommodate the missing cure status present in the data. It is to be noted that all parameters including $\alpha$ are considered simultaneously in the maximization step of the EM algorithm as opposed to applying the profile-likelihood technique considered by authors in previous literature. The performance of our proposed method is studied numerically, and established through a simulation study under various parameter settings (with respect to true $\alpha =0.0$, 0.5 and 1.0) and sample size ($n=200$ and 400) consideration, and an application to real-life interval censored data obtained from a smoking cessation study.

The estimation results obtained from the simulation study are overall satisfactory. The parameter estimates obtained are accurate, consistent with relatively lesser SEs and RMSEs with CPs being closer to the nominal level of 95% for most of the parameters. The parameter estimate of $\alpha$ shows lesser accuracy, higher SE and RMSE, with CP showing under-coverage for most of the settings. The flatness of the likelihood function surface with respect to $\alpha$ is considered as accountable for this issue. However, simultaneous-maximization approach is found to be 5 – 10 times faster than previously used profile-likelihood based approach to estimate $\alpha$ even for basic profile likelihood search range of $\alpha \in \{0.0,0.1,\dots ,1.0\}$. The SE corresponding to ${\psi }_B$ is found to be relatively higher when compared with other ${\psi }_b,b=0,\dots ,B-1$. This problem is attributed to lesser precision of estimation of ${\psi }_B$ owing to greater Euclidean distance between ${\tau }_{B-1}$ and ${\tau }_B$, and approximating the true exponential baseline hazard function in $[0,\infty )$ with PLA in the interval $\left[{\tau }_0,{\tau }_B\right]$. Due to nested nature of the BCTM incorporating PLA in its framework, increasing $B$ is associated with higher likelihood value (even though marginal). Analyzing the real-life example of smoking cessation data with our proposed method, the BCTM with PLA $(B=6)$ provides the highest maximized log-likelihood value. An estimated population survival function curve based on the same model reflects the presence of significant cure rates in the data, and that the smoking intervention group performs better than the usual care group of subjects. On comparing with parametric counterpart, the BCTM with gamma baseline hazard function corresponding to the latency of the model provides lowest AIC for the smoking cessation data.

As a potential future work, it is of great interest to extend the current framework to accommodate covariates with measurement error Barui & Yi (2020). In this context, it may be also be worth exploring the use of machine learning techniques in the incidence part to capture non-linearity in the data, which may result in improved predictive accuracy of cure Pal, Peng, Aselisewine, & Barui (2023). We are currently working on these problems and hope to report the findings in future manuscripts.

Supplementary Material