A New Cure Rate Model with Discrete and Multiple Exposures

1. Introduction

The study of models for survival or time-to-event data with a surviving fraction, known as cure rate models, have been extensively studied during the past few decades. Cure rate models have become more prevalent as new treatments for diseases cause a larger population of patients who are permanently cured of their disease. The first cure rate model was proposed by Boag (1949) and later revised by Berkson & Gage (1952). They proposed a two-component mixture model to analyze data where some subjects in the study were not susceptible to the event of interest. These subjects are called immunes, cured, or long-term survivors. Yakovlev et al. (1993) created a new formulation for mixture models to more effectively study cancer metastases. Tsodikov et al. (1998) further advanced the study of mixture models with the development of bounded cumulative hazard models. Other recent works of note include Davies et al. (2021), Pal (2021), Milienos (2022), Pal et al. (2022), Wang & Pal (2022), Pal et al. (2023), Pal & Aselisewine (2023) and Aselisewine & Pal (2023), amongst others. Interested readers may also look at Peng & Yu (2021) for a book-length account on mixture cure models. However, in the works previously mentioned, the models described primarily dealt with cancer metastasis and the spread of tumors in patients. Until recently, the application of these models in the context of infectious disease modeling had not been well researched. Some of the first work looking into the application of these models in the context of infectious diseases was done by Tournoud & Ecochard (2007), who studied a cure model with multiple discrete infection occasions. Specifically, the authors proposed a cure model where the number of pathogens affecting the $i^{th}$ patient, $i=\mathrm{1,2},\dots ,n$, at exposure time $t_k,k=\mathrm{0,1},2,\dots ,T$ ($k=0$ representing the time to initial exposure), was modeled using a Poisson random variable with mean ${\theta }_{t_k}$. This allowed them to study mother to child transmission of HIV-1 virus and the spread of nosocomial infections. Although the Poisson distribution was used to capture the latent number of pathogens at exposure time $t_k$, there was no justification provided for the suitability of the Poisson distribution, given that it cannot handle over- and under-dispersion that we commonly encounter with respect to the count in the number of pathogens. To advance their work, Tournoud & Ecochard (2008) again proposed a cure model with multiple infection occasions, but considered different distributions to model the number of pathogens affecting the $i^{th}$ patient at time $t_k$. They considered modeling the number of pathogens using a compound Poisson distribution to consider inter-subject heterogeneity. The issue with the compound Poisson distribution is the difficulty in interpreting a distribution in which a discrete random variable is attained using the sum of continuous random variables. In order to overcome these drawbacks, we propose the use of the Conway-Maxwell Poisson (COM-Poisson) distribution to model the number of pathogens at each infection occasion $t_k,k=\mathrm{0,1},2,\dots ,T$, assuming multiple infection occasions. This results in a flexible model to study infectious diseases with multiple discrete exposures and is more realistic compared to the model studied by Tournoud & Ecochard (2007) and Tournoud & Ecochard (2008). Furthermore, if we consider the special case of just one exposure, i.e. $k=0$, our model reduces to the COM-Poisson cure rate model of Rodrigues et al. (2009) and Balakrishnan & Pal (2016). To develop the likelihood inference, we do not rely on the numerical maximization of the log-likelihood function. Instead, motivated by the introduction of missing data arising from the set of censored observations, we develop the steps of the expectation maximization (EM) algorithm to calculate the maximum likelihood estimates (MLEs) of the model parameters.

The COM-Poisson distribution was first discussed in Conway & Maxwell (1962) to help solve problem related to queuing systems. Since its introduction, many authors have studied its properties including Shmueli et al. (2005), Boatwright et al. (2003), and more recently Li et al. (2020). The COM-Poisson is a generalized discrete distribution that allows for under- and over- dispersion of count data. As such, the COM-Poisson distribution allows for far more flexibility in modeling count data than other commonly used distributions. The Poisson, geometric, and Bernoulli distributions are all special cases of the COM-Poisson distribution. Rodrigues et al. (2009) were one of the first to explore the uses of the COM-Poisson distribution in the context of cure rate model. In the work of Rodrigues et al. (2009), the authors used the COM-Poisson distribution as a way to form a more flexible version of the unified cure rate model that was originally formulated by Yin & Ibrahim (2005). For other recent works on the COM-Poisson distribution in the context of cure rate models together with flexible modeling of lifetime data, interested readers may look at Balakrishnan & Pal (2013, 2015b), Pal & Balakrishnan (2017b), Pal et al. (2018), Wiangnak & Pal (2018), and Majakwara & Pal (2019), among others.

The rest of this paper is organized as follows. In Section 2, we introduce the COM-Poisson cure rate model in the context of infectious diseases with discrete multiple exposures to infection. In Section 3, we develop the estimation algorithm based on the EM algorithm to calculate the MLEs of the model parameters. In Section 4, we carry out an extensive Monte Carlo simulation study, where we present the results of both model fitting and model discrimination. To carry out the model discrimination, we use both likelihood ratio test and information-based criteria. In Section 5, we illustrate our proposed model using a recently collected data on COVID-19. Finally, in Section 6, we make some concluding remarks and discuss some potential future research work in the context of infectious diseases.

2. COM-Poisson cure rate model with discrete multiple exposures

Let $n$ represent the number of patients in the study. Let $t=\left\{t_0,\dots ,t_k,\dots ,t_T\right\}$ denote the multiple and successive moments of infections with $t_0$ denoting the initial moment of infection. At each moment of infection, let $M_{i,t_k}$, for $1\le i\le n,t_0\le t_k\le t_T$, denote the number of pathogens infecting the $i^{th}$ patient at time $t_k$. Note that $M_{i,t_k}$ are unobservable variables and we assume it to follow a discrete distribution with probability mass function (pmf) $p_{m_{i,t_k}}$. We assume the number of pathogens at each exposure time to be independent, which is along the lines of Tournoud & Ecochard (2008). It is important to note that we do not assume the distributions of $M_{i,t_0},M_{i,t_1},\dots ,M_{i,t_T}$ to be identical for the $i^{\mathit{\text{th}}}$ patient. Rather, at each exposure time $t_k$, we assume a suitable parameter of the distribution to depend on the exposure time. Let $Z_{i,j,t_k}$, for $1\le i\le n,t_0\le t_k\le t_T$ and $0\le j\le M_{i,t_k}$, be the $j^{\mathit{\text{th}}}$ pathogen promotion time of the $i^{th}$ patient at time $t_k$, in other words, the time taken by the $j^{\mathit{\text{th}}}$ pathogen to produce the event, which in the context of infectious diseases may represent the first biological sign of infection. For a given $M_{i,t_k}$, we assume $Z_{i,j,t_k}$ to be distributed with cumulative distribution function (cdf) $F\left(z\right)=F(z\mid \gamma )=1-S(z\mid \gamma )$, where $\gamma$ is the associated vector of parameters and $S(\cdot )$ is the corresponding survival function. Since this model follows data with multiple exposures, we need to re-parameterize the time scale of infection on the original scale. Let $y=z+t_k$, then, we can re-parameterize the cdf and survival function as follows:

and

Given that patients are exposed at different discrete time points and at each exposure time there are several competing pathogens, the time-to-event, which can be the first biological sign of infection or time to recovery from an infectious disease, can be defined as

where $Z_{i,0,t_k}$ is such that

A patient is termed “immune” of an infection if there are no competing pathogens at each exposure time point, and its probability, termed as cure rate, is defined as

We will now introduce the COM-Poisson distribution as a means of modeling the number of competing pathogens at each discrete exposure time point. For this purpose, we will first introduce the COM-Poisson distribution. If a random variable $M$ follow a COM-Poisson distribution, then, its pmf is given by

where $\theta >0,\nu >0$, and $Z(\theta ,\nu )$ is the normalizing constant calculated by

From the above equation, we can see the cured fraction, in other words the part of the population that is not affected by the pathogen, denoted by $p_0$, is given by

The COM-Poisson distribution has three special cases. When $\nu =1$, then, $Z(\theta ,\nu )=e^{\theta }$, which results in the Poisson distribution with mean $\theta$. As $\nu \to \mathrm{\infty },Z(\theta ,\nu )\to 1+\theta$, which means the COM-Poisson distribution converges in distribution to the Bernoulli distribution with $P[M=1;\theta ,v]=\frac{\theta }{1+\theta }$. When $\nu =0$ and $\theta <1,Z(\theta ,\nu )=\frac{1}{1-\theta }$, and the COM-Poisson reduces to a geometric distribution with parameter $1-\theta$. However, if $\nu =0$ and $\theta \ge 1$, then $Z(\theta ,\nu )$ does not converge. Therefore, the COM-Poisson distribution is undefined in this particular case. One advantage of the COM-Poisson distribution is that it allows for both under-dispersion and over-dispersion of data relative to the Poisson distribution. When $\theta >1$, the data is under-dispersed relative to the Poisson distribution, whereas when $\theta <1$, the data is over-dispersed relative to the Poisson distribution.

In this research, we will assume the number of competing pathogens to follow a COM-Poisson distribution with parameters ${\theta }_{t_k}$ and $\nu$ at exposure time $t_k$, for $k=\mathrm{0,1},2,\dots ,T$. Note that the parameter ${\theta }_{t_k}$ represents the infection intensity at exposure time $t_k$ and hence carries biological interpretation. In a practical scenario, the infection intensity at a given exposure time will differ across patients. To capture this heterogeneity in patient population, we link ${\theta }_{t_k}$ to a set of covariates ${\mathit{x}}_{t_k}$ at each exposure time $t_k$, for $k=\mathrm{0,1},2,\dots ,T$, using the log-linear link function ${\theta }_{t_k}=\mathrm{e}\mathrm{x}\mathrm{p}\left({\mathit{x}}_{t_k}^{\prime }{\mathit{\beta }}_{t_k}\right)$, where ${\mathit{\beta }}_{t_k}$ is the corresponding vector of regression coefficients. Now, let us consider two exposure times, $t_0$ and $t_1$, and derive the survival function of the random variable $Y$ in (3), known as the population survival function or the long-term survival function.

Theorem 2.1 Given two discrete exposure times $t_0$ and $t_1$, let $M_{t_0}$ and $M_{t_1}$ denote the number of pathogens at times $t_0$ and $t_1$. Let $M_{t_0}$ and $M_{t_1}$ both follow COM-Poisson distribution with parameters $\left({\theta }_{t_0},\nu \right)$ and $\left({\theta }_{t_1},\nu \right)$, respectively. Furthermore, let $S_{t_0}\left(y\right)$ and $S_{t_1}\left(y\right)$ denote the pathogen promotion time distribution at exposure times $t_0$ and $t_1$, respectively, as described in (2). Then, the overall survival function of the variable $Y$ in (3), denoted by $S_{\mathit{\text{pop}}}\left(y\right)$, is given by

A proof of Theorem 2.1 is provided in the Appendix.

If we generalize the above survival function with multiple exposure times $\left\{t_0,\dots ,t_k,\dots ,t_T\right\}$, then, we have

Note that in (5), if we just consider one exposure, then, $S_{pop}\left(y\right)=\frac{Z\left({\theta }_{t_0}{S}_{t_0}\left(y\right),\nu \right)}{Z\left({\theta }_{t_0},\nu \right)}$, which reduces to the model proposed by Rodrigues et al. (2009) and Balakrishnan & Pal (2016). The density function corresponding to (5), known as the long-term density function, is given by

where $f_{t_k}\left(y\right)$ is the density function associated with $S_{t_k}\left(y\right)$. Hence, the cure rate or cured proportion is given by

In Table 1, we present the long-term survival function and the cured proportion for the COM-Poisson cure rate model and its special cases with discrete multiple exposures. In Table 2, we present the corresponding long-term density functions.

3. Maximum likelihood estimation

We develop the EM algorithm to carry out the maximum likelihood estimation of model parameters; see Pal & Balakrishnan (2016, 2017a,c, 2018) and Treszoks & Pal (2022). First, let $F_{\mathit{\text{pop}}}$ and $F_1$ denote the cdf of the entire population and the susceptible population, respectively. Furthermore, let $S_1$ denote the survival function of the susceptible population. Define $J$ to be an indicator variable which takes on the value of 0 if the patient is immuned and 1 if the patient is susceptible. As such, we have $P[J=0]=p_0$ and $P[J=1]=1-p_0$. Furthermore, note that

and hence

Using the form of $S_{pop}\left(y\right)$ as in (5), we can get an expression for $S_1$ as

where $p_0$ is defined as in (7). Let $\delta$ denote the censoring indicator such that when $\delta =1$, the actual time-to-event is observed and $\delta =0$ when the time-to-event is right censored. Let us define two sets, say, $I_1$ and $I_0$ as $I_1=\left\{i:{\delta }_i=1\right\}$ and $I_0=\left\{i:{\delta }_i=0\right\}$. Note that for the set $I_0$, the value of $J$ is unknown (since a right censored patient may be susceptible to the event of interest), and this introduces the missing data. Let $\mathit{\theta }$ denote the unknown parameter vector of our model. Then, the complete data likelihood function is given by

Given the above likelihood function, we can write the complete data log-likelihood function as

Now that we have an expression for the complete data log-likelihood function, we can begin discussion on the development of the EM algorithm; see McLachlan & Krishnan (2008). In the expectation step (E-step), we calculate the expectation of the complete log-likelihood function with respect to the distribution of the unobserved $J_i^{\prime }$s, given the model parameters and the observed data. We should note here that $J_i^{\prime }\mathrm{s}$ are Bernoulli random variables and are linear in the complete data $\mathrm{l}\mathrm{o}\mathrm{g}$-likelihood function. As such, at the $r^{\mathit{\text{th}}}$ iteration step, we simply need to calculate ${\pi }_i^{\left(r\right)}=E\left(J_i\mid {\mathit{\theta }}^{\left(r\right)}\right.,\mathrm{ }\mathrm{ }\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a})$, for $i\in I_0$, where ${\mathit{\theta }}^{\left(r\right)}$ denotes the current parameter value at the $r^{\mathit{\text{th}}}$ iteration of the EM algorithm. Hence, for the $i^{\mathit{\text{th}}}$ censored observation, we can calculate ${\pi }_i^{\left(r\right)}$ as

Therefore, in the E-step, we only replace $J_i$ in (8) with ${\pi }_i^{\left(r\right)}$ if the $i^{\mathit{\text{th}}}$ observation is censored. As done in McLachlan & Krishnan (2008), we will denote the conditional expectation of the complete data $\mathrm{l}\mathrm{o}\mathrm{g}$-likelihood function at the $r^{\mathit{\text{th}}}$ iteration as $Q\left(\mathit{\theta },{\mathit{\pi }}^{\left(r\right)}\right)$, where ${\mathit{\pi }}^{\left(r\right)}$ is the vector of ${\pi }_i^{\left(r\right)}$ values.

The next step in the EM algorithm is the maximization step (M-step). In this step, we maximize the $Q\left(\mathit{\theta },{\mathit{\pi }}^{\left(r\right)}\right)$ function with respect to $\mathit{\theta }$ over the parameter space $\mathbf{\Theta }$. In other words, we choose ${\mathit{\theta }}^{(r+1)}$ to be a value of $\mathit{\theta }\in \mathbf{\Theta }$ such that

see Pal & Roy (2021, 2022, 2023). We continue to repeat this algorithmic process until some convergence criterion is satisfied, for example, $\left|\frac{{\mathit{\theta }}^{r+1}-{\mathit{\theta }}^r}{{\mathit{\theta }}^r}\right|<ϵ$, where $ϵ$ is some pre-defined tolerance such as $ϵ=0.001$.

4. Simulation study

We shall now assume the pathogen promotion times to follow a Weibull distribution with probability density function given by

Thus, we now have $\mathit{\gamma }={\left({\gamma }_1,{\gamma }_2\right)}^{\prime }$. Note that one can also assume other parametric lifetime distributions here; see Balakrishnan & Pal (2012, 2015a,b). One can also model these promotion times assuming a semi-parametric framework or a completely parametric framework with respect to the pathogen promotion times; see Balakrishnan et al. (2016) and Balakrishnan et al. (2017).