Introduction

Standard statistical methods for time-to-event analysis—often referred to as survival models—typically assume that every subject in the study population is susceptible to the event of interest and will eventually experience the event if the follow-up is sufficiently long. However, in practice, a fraction of individuals may never experience the event. These subjects are considered as having infinite survival times and are referred to as cured (Peng and Yu 2021).

Over the last few decades, specific methods have been developed to deal with this type of data, collectively known as cure model analysis. Several reviews and methodological contributions summarize the theory and application of these models (Maller et al. 2024; Peng and Yu 2021; Patilea and Van Keilegom 2020; Mailis Amico and Van Keilegom 2018; Morbiducci, Nardi, and Rossi 2003).

Cure models explicitly account for a subpopulation that is not susceptible to the event, allowing direct modeling of the cure proportion and the effects of covariates on it. Additionally, these models provide inference on the survival distribution of the uncured (susceptible) subjects. Most existing cure models are modifications of standard survival models that incorporate the cure probability. Based on how this incorporation is done, cure models are typically categorized as either mixture or nonmixture cure models. A key advantage of mixture cure models is that they allow covariates to affect both the cured and uncured groups differently. This package focuses on mixture cure models.

Let \(Y\) denote the time to event. It is assumed that subjects are subject to right censoring, and the censoring time \(C\) and the event time \(Y\) are conditionally independent given covariates \(\boldsymbol{X}\). The conditional distribution and survival functions of \(Y\) are given by \[ F(t|\boldsymbol{x}) = \mathbb{P}(Y \leq t \mid \boldsymbol{X} = \boldsymbol{x}), \quad S(t|\boldsymbol{x}) = 1 - F(t|\boldsymbol{x}). \] Under right censoring, only the pair \((T, \delta)\) is observed, where \(T = \min(Y, C)\) and \(\delta = \mathbb{I}(Y \leq C)\) is the censoring indicator. The conditional distributions of \(C\) and \(T\) are denoted by \(G(t|\boldsymbol{x})\) and \(H(t|\boldsymbol{x})\), respectively.

Define the cure indicator \(\nu = \mathbb{I}(Y = \infty)\), where \(\nu = 0\) for susceptible individuals and \(\nu = 1\) for cured individuals. The probability of not being cured (incidence) is \(p(\boldsymbol{x}) = \mathbb{P}(\nu = 0 \mid \boldsymbol{X} = \boldsymbol{x})\), and the conditional survival function of the uncured group (latency) is \[ S_0(t|\boldsymbol{x}) = \mathbb{P}(Y > t \mid \nu = 0, \boldsymbol{X} = \boldsymbol{x}). \] Then, the mixture cure model can be written as: \[ S(t|\boldsymbol{x}) = 1 - p(\boldsymbol{x}) + p(\boldsymbol{x}) S_0(t|\boldsymbol{x}). \]

Nonparametric hypothesis tests for the cure rate

The proposed hypotheses are the following:

\[\begin{equation}\label{Eq:Test} \mathcal{H}_0 : \mathbb{E}(\nu | \boldsymbol{X}) \,\, \equiv \,\, 1 - p \quad \text{a.s.} \quad \text{vs} \quad \mathcal{H}_1 : \mathbb{E}(\nu | \boldsymbol{X}) \,\, \not\equiv \,\, 1 - p \quad \text{a.s.}, \end{equation}\]

which tests whether the covariate vector \(\boldsymbol{X}\) has an effect on the cure probability. The main problem with the hypotheses \(\eqref{Eq:Test}\) is that the response variable (the cure indicator, \(\nu\)) is only partially observed due to censoring.

In this package, this challenge is addressed by estimating the cure indicator using the methodology proposed by M. Amico, Van Keilegom, and Han (2021). The idea is as follows: define \(\tau\) as the upper limit of the support of the lifetime for a susceptible individual, where \(\tau = \sup_{x} \tau(x)\) and \(\tau(x) = \inf\{t:S_0(t|x) = 0\}\). We assume that \(\tau < \infty\) and that the follow-up time is long enough so that \(\tau < \tau_{G(x)}\) for all \(x\), where \(\tau_{G(x)} = \inf\{t:G(t|x) = 1\}\).

Therefore, it is reasonable to consider that all the individuals with a censored observed time greater than \(\tau\) can be categorized as cured \((\nu = 1)\). Then, the proxy of the cure rate \(\nu\) is defined as:

\[\begin{equation}\label{eq:nu_tilde_0} \tilde{\nu} = \mathbb{I}(T > \tau) + (1-\delta)\mathbb{I}(T \leq \tau) \, \frac{1 - p(\boldsymbol{X})}{1-p(\boldsymbol{X}) + p(\boldsymbol{X})S_0(T|\boldsymbol{X})}. \end{equation}\]

Since under \(\mathcal{H}_0\) the cure probability does not depend on \(\boldsymbol{X}\), the approximation of \(\nu\) becomes under \(\mathcal{H}_0\):

\[\begin{equation}\label{eq:nu_tilde} \tilde{\nu}_0 = \mathbb{I}(T > \tau) + (1-\delta)\mathbb{I}(T \leq \tau) \, \frac{1 - p}{1-p + pS_0(T|\boldsymbol{X})}. \end{equation}\]

One issue with the proposal of \(\tilde{\nu}_0\) is that it cannot be computed in practice since \(p\), \(\tau\), and the latency \(S_0(t|\boldsymbol{X})\) are unknown. An estimator that handles this problem is proposed as:

\[\begin{equation}\label{eq:nu_hat} \hat{\nu}_{h} = \mathbb{I}(T > \hat{\tau}) + (1-\delta)\mathbb{I}(T \leq \hat{\tau}) \, \frac{1 - \hat{p}}{1-\hat{p} + \hat{p}\hat{S}_{0,h}(T|\boldsymbol{X})} \end{equation}\]

where \(\hat{p}\) is the nonparametric estimator proposed by Laska and Meisner (1992), \(\hat{\tau}\) is the largest observed uncensored survival time, and \(\hat{S}_{0,h}(T|\boldsymbol{X})\) is the nonparametric estimator based on the Beran estimator proposed by López-Cheda, Jácome, and Cao (2017).

Two statistics for the covariate hypothesis test, based on the martingale difference correlation between the estimated cure rate \(\hat{\nu}_h\) and the covariate \(\boldsymbol{X}\), are proposed:

\[ \text{MDCV}_n(\hat{\nu}_h|\boldsymbol{X})^2 \quad \text{and} \quad \text{MDCU}_n(\hat{\nu}_h|\boldsymbol{X})^2, \]

where \(\text{MDCV}_n\) and \(\text{MDCU}_n\) correspond to the biased and bias-corrected estimators, respectively.

In practice, we do not have access to the true null distribution of these statistics. To address this, the null distribution is approximated using two methods: a permutation test (see Gretton et al. (2008), Pfister et al. (2018)) and a chi-square test.

The permutation test compares the statistic computed from the original data to the distribution of the statistic computed over multiple randomly permuted datasets. While this approach is nonparametric and robust, it can be computationally demanding, especially for large sample sizes or when many permutations are required.

To provide a faster alternative, a chi-square approximation is implemented for the \(\text{MDCU}_n^2\) test statistic. This approach extends the idea proposed by Shen, Panda, and Vogelstein (2022) and offers a computationally efficient solution for hypothesis testing.

Goodness-of-Fit Test

As discussed by Müller and Van Keilegom (2019), most studies in the cure model literature either assume that no cure fraction exists, or, if it does, that it follows a logistic model. However, both assumptions may be unrealistic in practice.

First, the existence of a cure fraction has been debated in the survival analysis literature. Second, even when a cure fraction is present, there is no theoretical reason to expect the cure probability to follow a monotonic function of the covariates—let alone a logistic form.

To address this issue, Müller and Van Keilegom (2019) propose goodness-of-fit tests for evaluating whether a parametric model for the cure probability is consistent with the observed data. These tests are flexible and well-suited to situations where the functional form of the cure model is uncertain or possibly misspecified. Notably, their test is the first to allow the cure rate to vary with a covariate.

The hypotheses tested are:

\[ \mathcal{H}_0 : p(x) = p_{\theta}(x) \quad \text{for some } \theta \in \Theta \quad \text{vs.} \quad \mathcal{H}_1 : p(x) \neq p_{\theta}(x) \quad \forall \theta \in \Theta, \]

where \(\Theta\) is a finite-dimensional parameter space, and \(p_{\theta}(x)\) is a known parametric form up to the parameter vector \(\theta\).

The proposed test statistic is a weighted \(L_2\)-distance between a nonparametric estimator \(\hat{p}_h(x)\) of the cure probability and its parametric counterpart \(p_{\hat{\theta}}(x)\) estimated under \(\mathcal{H}_0\). The test statistic is defined as:

\[ \mathcal{T}_n = nh^{1/2} \int \left( \hat{p}_h(x) - p_{\hat{\theta}}(x) \right)^2 \pi(x)\, dx, \]

where:

• \(\hat{p}_h(x)\) is a nonparametric estimator with bandwidth \(h\),
• \(\hat{\theta}\) is an estimator of the parameter \(\theta\),
• \(\pi(x)\) is a user-specified weight function (often a density on the covariate space).

In practice, the empirical version of the test statistic is used:

\[ \tilde{\mathcal{T}}_n = nh^{1/2} \frac{1}{n} \sum_{i=1}^n \left( \hat{p}_h(X_i) - p_{\hat{\theta}}(X_i) \right)^2. \] The critical values of the test are approximated using a bootstrap procedure.

MDCcure

library(MDCcure)

set.seed(123)
X <- matrix(rnorm(10 * 2), 10, 2)
Y <- matrix(rnorm(10 * 2), 10, 2)

mdc(X, Y, center = "U")
## [1] 0.2724937
mdc(X, Y, center = "D")
## [1] 0.3998999

set.seed(123)
X <- matrix(rnorm(10 * 2), 10, 2)
Y <- matrix(rnorm(10 * 2), 10, 2)
Z <- matrix(rnorm(10 * 2), 10, 2)

pmdd(X, Y, Z)
## [1] 0.1247584
pmdc(X, Y, Z)
## [1] 0.2555037

set.seed(1234)

theta0 = c(1,1)
gamma0 = 0.5
gamma1 = 0.5

cT = 0

n <- 200
maxT = 0.02*n

X = runif(n,-1,1)
Z = runif(n,-1,1)
phi = exp(theta0[1]+theta0[2]*X)
phi = phi/(1+phi)   
B = (runif(n) <= phi)

aT = (X+1)^cT
lambdaX = exp(gamma0+gamma1*X)
bT = lambdaX^(-1/aT)
tau = qweibull(0.9,shape=mean(aT),scale=mean(bT))
Y = rep(100000,n)
count = 0
for (j in 1:n)
{if (B[j]==1)
{stop = 0
while (stop==0)  
{Y[j] = rweibull(1,shape=aT[j],scale=bT[j])
if ((Y[j] > tau)*(count <= maxT)) 
{Y[j] = tau
count = count + 1}
if (Y[j] <= tau) stop = 1}}}

aC = 1
bC = 1.5
C = rweibull(n,shape=aC,scale=bC)
C = replace(C,C>tau,tau+0.001)  

T = apply(cbind(Y,C),1,min)
Delta = as.numeric(Y <= C)

#Covariate hypothesis test for the cure rate with one covariate
testcov(X, T, Delta, method = "All", P = 499)
## ============== Covariate hypothesis test for the cure rate ============
## Method: All 
## 
## MDCU: Martingale difference correlation unbiased with permutations
## Test statistic: 0.03986 
## p-value: 0.002 
## Number of permutations: 499 
## 
## MDCV: Martingale difference correlation biased with permutations
## Test statistic: 0.04746 
## p-value: 0.008 
## Number of permutations: 499 
## 
## FMDCU: Fast Chi-square test based on MDC unbiased
## Test statistic: 0.03986 
## p-value: 0.0027426 
## 
## GOFT 
## Test statistic: 0.769 
## p-value: 0.03006 
## Number of bootstrap replicates: 499 
## 
## ========================================================================

#Covariate hypothesis test for the cure rate with two covariates
testcov2(X, T, Z, Delta, P = 499)
## ---------------------------------------------------------------------
## Covariate hypothesis test for the cure rate with two covariates 
## ---------------------------------------------------------------------
## Hypotheses:
##   H0: The conditional mean of the cure status given X adjusting on Z does not depend on X,
##       i.e., E[nu|X,Z] = E[nu|Z].
##   H1: The conditional mean of the cure status depends on X adjusting on Z
## 
## Test Statistic: 0.02869 
## p-value: 0.004 
## Number of permutations: 499
testcov2(Z, T, X, Delta, P = 499)
## ---------------------------------------------------------------------
## Covariate hypothesis test for the cure rate with two covariates 
## ---------------------------------------------------------------------
## Hypotheses:
##   H0: The conditional mean of the cure status given Z adjusting on X does not depend on Z,
##       i.e., E[nu|Z,X] = E[nu|X].
##   H1: The conditional mean of the cure status depends on Z adjusting on X
## 
## Test Statistic: -0.002226 
## p-value: 0.744 
## Number of permutations: 499

#Goodness-of-fit test for the cure rate 
goft(X, T, Delta, model = "logit")
## ---------------------------------------------------------------------
## Goodness-of-fit test for the cure rate in a mixture cure model
## -------------------------------------------------------------------
## 
## Model: logit 
## 
## Hypotheses:
##   H0: The model fits the data, i.e., p(x) = p_theta(x)
##   H1: The model does not fit the data, i.e., p(x) != p_theta(x)
## 
## Test Statistic: 0.1192 
## p-value: 0.4168 
## Number of bootstrap replications: 499
# plotCure(X, T, Delta, density = FALSE)