Biometrical Journal 58 (2016) 6, 1485–1505 DOI: 10.1002/bimj.201500094 1485 Risk-adjusted monitoring of time to event in the presence of long-term survivors Jocelânio W. Oliveira∗,1, Dione M. Valença2, Pledson G. Medeiros2, and Magaly Marçula3 1 Instituto de Matemática e Estatı́stica, USP, 05508-090 São Paulo, Brazil 2 Departamento de Estatı́stica, CCET, UFRN, 59078-970 Natal, Brazil 3 Instituto do Coração, FMUSP, 05403-900 São Paulo, Brazil Received 26 May 2015; revised 20 April 2016; accepted 5 May 2016 The use of control charts for monitoring schemes in medical context should consider adjustments to incorporate the specific risk for each individual. Some authors propose the use of a risk-adjusted survival time cumulative sum (RASTCUSUM) control chart tomonitor a time-to-event outcome, possibly right censored, using conventional survival models, which do not contemplate the possibility of cure of a patient. We propose to extend this approach considering a risk-adjusted CUSUM chart, based on a cure rate model. We consider a regression model in which the covariates affect the cure fraction. The CUSUM scores are obtained for Weibull and log-logistic promotion time model to monitor a scale parameter for nonimmune individuals. A simulation study was conducted to evaluate and compare the performance of the proposed chart (RACUF CUSUM) with RAST CUSUM, based on optimal control limits and average run length in different situations. As a result, we note that theRASTCUSUM chart is inappropriate when applied to data with a cure rate, while the proposed RACUF CUSUM chart seems to have similar performance if applied to data without a cure rate. The proposed chart is illustrated with simulated data and with a real data set of patients with heart failure treated at the Heart Institute (InCor), at the University of São Paulo, Brazil. Keywords: Cure rate; Risk-adjusted CUSUM; Statistical process control; Survival analysis.  Additional supporting information including source code to reproduce the resultsmay be found in the online version of this article at the publisher’s web-site 1 Introduction Statistical Process Control (SPC) is a set of statistical tools widely used in the industrial context to monitor and control a variety of processes. However, according to Woodall (2006), there are many applications of SPC for health-related monitoring, especially in the use of control charts. Examples include monitoring the performance of hospitals based on variables such as infection rates, rates of patient falls, or time until an event of interest happens. The CUSUM (Cumulative Sum) chart (Page, 1954) is one of the main tools to monitor such situations where there is an interest in detecting small and persistent changes in quality. Several authors argue that, unlike usual applications in industry, where individuals (equipment, products, goods, etc.) are considered as homogeneous, monitoring schemes in medical context should consider adjustments to incorporate the specific risk for each individual (risk adjustment). Thus, it is important to take into account, for example, sex, age, health conditions, or risk factors associ- ated with each patient before undergoing surgery. Grigg and Farewell (2004) provide an overview ∗Corresponding author: e-mail: jwesley@ime.usp.br ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 1486 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event on the subject and describe several works proposing the use of the risk-adjusted CUSUM (RA CUSUM) control chart, mainly for monitoring binary outcomes. Steiner et al. (2000, 2001) de- veloped CUSUM methods for monitoring quality based on a binary variable, which essentially indicates whether the surgical procedure was successful or not for a given individual. These au- thors propose a monitoring scheme of surgical performance of patients (the chance of each in- dividual in surviving more than 30 days), using a likelihood score based on a logistic regression model. Some proposals for risk-adjusted CUSUM charts based on time-to-event models for right censored data have recently been described. Biswas and Kalbfleisch (2008) consider a Cox model for a failure time outcome and provide a method to obtain approximations to the ARL (Average Run Length). Gandy et al. (2010) propose a more general CUSUM chart based on the partial likelihood ratio between an out-of-control hazard rate and an in-control hazard rate. Sego et al. (2009) developed an RA CUSUM control chart to monitor a time-to-event outcome using accelerated failure time models (Lawless, 2003). They call it RAST CUSUM (risk-adjusted survival time CUSUM) and express the CUSUM score for the Weibull and log-logistic models. All these mentioned works demonstrate in simulation studies that the time-to-event approach tends to provide a faster detection when compared to the logistic regression based methods. Somemore recent studies about CUSUMprocedures based on time-to-event data, can be cited. Sun and Kalbfleisch (2013) propose a risk-adjusted Observed−Expected CUSUM chart with monitoring bands to simultaneously monitor for failure time outcomes that are worse or better than expected. Phinikettos and Gandy (2014) studied a nonparametric method for monitoring time-to-event data. Sun, Kalbfleisch and Schaubel (2014) developed a weighted cumulative sum to monitor medical outcomes with dependent censoring and Zhang et al. (2016) investigate the effect of estimation error on the performance of risk-adjusted survival time CUSUM. In this context, it seems natural that monitoring schemes in medical procedures could consider the success of procedures as the cure of the patients. However, the conventional survival models do not contemplate this possibility. In fact, these models assume that the event (death due to the procedure) will actually occur for everybody, although the time until this occurrence is not observable for some individuals (censoring). A natural extension should be a survival model incorporating a cure rate. Cure rate models have been extensively studied in the literature for datasets where part of the population will never experience the event. The book by Maller and Zhou (1996) presents an updated review of the first models proposed, called standard mixture model. These authors argue that the occurrence of a high percentage of censoring at the end of the study in a sufficient follow-up time is an indication of a cure rate in the population. Another class of models with many applications to cancer clinical trials has been proposed by Yakovlev et al. (1993) and Chen et al. (1999). Rodrigues et al. (2009) refer to this model as promotion time cure model. In this article, a risk-adjusted CUSUM chart based on a cure rate model is proposed. We call this RACUF CUSUM (risk-adjusted with cure fraction CUSUM) and develop a sim- ulation study to evaluate the detection performance of the proposed chart when compared with the RAST CUSUM (Sego et al., 2009) based on the ARL performance, in different scenarios. This article is organized as follows: In Section 2 we describe the general form of the CUSUM chart, the RA CUSUM chart and its extension to deal with time until the occurrence of an event based on survival models (RAST CUSUM). In Section 3 we present our proposal for the use of CUSUM charts applied to survival data with cure fraction (RACUF CUSUM). In Section 4 a simulation study is presented. Section 5 presents two illustrations of the results with a simu- lated and a real dataset. We end this article with some discussions about the results in Section 6. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1487 2 Risk-adjusted CUSUM charts Ordinary control charts have been widely used as a tool to check the stability of industrial processes (Montgomery, 2007). The CUSUM chart is characterized by evaluating the process as a whole, accumulating information of all individuals, rather than using each observation alone. The CUSUMproposed by Page (1954) was based on theWald sequential test (Grigg et al., 2003) and involves a sample weight (score) assigned to each observation. To define this weight, let i = 1, 2, . . . the index observations (patients) whose outcomes will be monitored in chronological order of execution of the medical procedure. We assume a probability model for the outcome and consider L(ξ |Di) a contribution to likelihood function of the i-th observation, where ξ is a vector of parameters and Di represents the outcome and the characteristics of patient i. Patient risk is taken into account by using a regression model. Thus, in the control state, it is expected that ξ0 is the true value of ξ . The model fitted for historical data under control is used to estimate ξ0. A common assumption is that the model for in-control data is known and the estimation error is negligible. In this regard, some authors argue that the impact of the error should not be ignored because it can affect the performance of the chart (e.g., Jones et al., 2004; Gandy and Kvaløy, 2013; Zhang et al., 2016). Here we make this usual assumption but we briefly investigate the estimation error in Section 4. The weight for C[USUM de]fined as log-likelihood ratio and denoted by Wi is given by = L(ξ1|Di)Wi log = (ξL(ξ |D ) 1|Di)− (ξ0|Di), (1)0 i where ξ1 is the nominal value of the parameter for out-of-control state and (.) is the log-likelihood function. The change of ξ = ξ0 (under control) to ξ = ξ1 (out-of-control), must be interpretable with respect to process quality. The CUSUM procedure involves plotting Zi = max(0,Zi−1 +Wi), i = 1, 2, . . . , (2) where Z0 = 0. An alarm is raised when Zi > h, where h is a positive control limit. Moustakides (1986) shows that this procedure is asymptotically optimal for detecting changes. Note that, the higher the score Wi in Eq. (1), the more the CUSUM accumulates errors in the statistic Zi, which can exceed the upper limit, indicating that the observed process is out-of-control. Besides, note that the CUSUM procedure has a barrier at zero. We describe below the proposed RAST CUSUM chart by Sego et al. (2009), based on accelerated failure time (AFT) models. 2.1 RAST AFT models CUSUM chart Let Ti be the event time for individual i, and let xi = (xi1, ..., xip) be a fixed covariate vector. The AFT model can be represented by  logTi = μ+ γ xi + σVi, i = 1, . . . , n, (3) where Vi are independent and identically distributed random errors with a distribution independent of xi. The vector γ = (γ1, . . . , γ p) , μ and σ are unknown parameters. Survival times are possibly subject to right censoring. Here, the censoring times are represented by the independent random variables Ci, for i = 1, . . . , n, which are assumed to be independent of T1, . . . ,Tn. The censoring mechanism is assumed to be noninformative, that is, the distribution of the Ci’s does not depend on unknown parameters. Let δi = 1, if the observation for individual i is a failure time, and δi = 0, if it is a censoring time and letYi = min(Ti,Ci). The observations can be represented by the pairs of random variables (Yi, δi), and the covariate vectors x  i , for i = 1, . . . , n. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com 1488 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event The log-likelihood function for the unknown parameters vector ηwith respect to the i-th individual, is given by log[ f (ti|η)δi S(t |η)1−δi i ], (4) whereS(t|η) denotes the survival functions ofTi, given xi, and f (t|η) = − ∂t S(t|η) is the corresponding∂ probability density function.  To obtain a Weibull AFT model, it is assumed that Ti ∼ Weibull (α, λ exp(γ xi)), where α = 1/σ and λ = eμ are the form and scale parameters of Ti (when the influence of covariates is ignored), γ is a vector of regression coefficients, and Vi in Eq. (3) has a standard extreme value distribution.  Similarly, the log-logistic AFT model is obtained assuming that Ti ∼ log-logistic(α, λ exp(γ xi)), where Vi in Eq. (3) has a standard logistic distribution. Survival functions for the Weibull and log- logistic AFT models ar{e res(pectively give)n b}y [ ( ) t α t α ]−1 S(ti|η) = exp − i  and S(ti|η) = 1+ iexp x exp  , (5)λ (γ i) λ (γ xi) where η = (α, λ, γ ). Considering similar interpretations for these models, Sego et al. (2009) assume that when the process is in-control, the scale parameter takes the value λ0, and that the interest is to detect a change for λ1 =  ρ1λ0, as new individuals are observed. The vector η0 = (α, λ0, γ ) is estimated from historical data (training) from a process (theoretically) under control. The authors also assume that the model used is accurate, and the error in the parameter estimation is negligible. It is easy to verify that a reduction in λ (0 < ρ1 < 1) corresponds to a reduction in the average (and median) survival time for Weibull and log-logistic models, as well as an increase in λ (ρ1 > 1) leads to an increase in these measurements. The Weibull and log-logistic CUSUM scores on Sego et al. (2009), based on Eq. (1), are given respectively by ( ) t α W = (1− ρ−α ) ii 1  − δexp x iα log ρ1 (6)λ0 (γ i) and { [ ( )α] [ ( )α]} W = −αδ log ρ + 2δ ti log 1+ i tii i 1 − log 1+ . (7) λ0 exp   (γ xi) ρ1λ0 exp(γ xi) Finally, the monitoring of the parameter λ can then be performed from the statistical CUSUM described in Eq. (2). 3 The RACUF CUSUM chart In this section, the unified cure rate model given by Rodrigues et al. (2009) is presented, and the CUSUM chart is developed from this model. Suppose that for an individual i in the population, we have Mi causes or risks competing for the occurrence of a single event of interest, with probability function pθ (mi) = Pθ (Mi = mi). Given Mi = mi, the random variables Ri1,Ri2, . . . ,Rim , are independent and identically distributed (i.i.d .)i representing (unobserved) times to event occurrence for the i-th individual due to j-th cause ( j = 1, . . . ,Mi), with common distribution function F (z|η) and survival function S(z|η) = 1− F (z|η) where η is a vector of parameters. LetTi be a random variable representing the time until the occurrence of the event, defined as Ti = min{Ri0,Ri1, . . . ,RiM } where P(Ri0 = ∞) = 1. This assumption permitsi the occurrence of an infinite lifetime in immune individuals, since Mi = 0 means that there are no causes or risks for the occurrence of the event. To exemplify the latent variables Mi and Ri within a ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1489 biological context, consider that Mi denotes the number of metastatic-competent tumour cells for that individual i left active after the initial treatment. A metastatic-competent tumour cell is a tumour cell that has a potential for metastasizing. If Mi = 0 the individual is cured. If Mi = m (m > 0), then this individual has m tumour cells and Ri1,Ri2, . . . ,Rim are random times for each metastatic-competent cell to produce a detectable tumour. For different applications with this model see, for example, Chen et al. (1999), Tournoud and Ecochard (2007), and Fonseca et al. (2013). It can be shown that the common long-term survival function for Ti is given by ∑∞ Sp(t) = P(Ti > t) = pθ (0)+ pθ (m mi)S(t|η) i . (8) mi=1 The survival function S(t|η) is such that limt−→∞ S(t|η) = 0. Hence Sp(t) is an improper sur- vival function, that is, limt−→∞ Sp(t) > 0. The proportion of cured individuals (cure rate) is given by limt−→∞ Sp(t) = pθ (0). The probability density function of the random variables Ti can be obtained from Eq. (8) as ∑∞ fp(t) = f (t|η) mi pθ (mi) [S(t|η)]mi−1. mi=1 3.1 Likelihood for unified approach In addition, consider that, for the individual i (i = 1, . . . , n), Yi = min{Ti,Ci} is the observable life- time, where Ci is the right censoring time (random and uninformative), independent of Ti, and let δi be the failure/censoring indicator, with δi = 1 if Ti ≤ Ci and δi = 0 if Ti > Ci. Also consider xi = (xi1, xi2, . . . , xip)′ a vector of associated covariates. To simplify the notation, the n-dimensional vectors of observations y = (y1, y2, . . . , y ′n) , δ = (δ1, δ2, . . . , δn)′ and m = (m1,m ′′ 2, . . . ,mn) and thecovariate matrix X = (x1, x2, . . . ,xn) of dimension n × p is defined. Hence, the complete dataset is de- noted by Dc = (n, y, δ,m,X), and the data without the latent variables is denoted by D = (n, y, δ,X). The covariates can be included in the model through some relation: θ ≡ θ (x′iβ), where β = (β1, ..., βp)′ is the vector of regression coefficients. Thus, the vector of unknown parameters in the model is denoted by φ = (β′, η′)′, and, after some algebraic manipulations it can be shown that the likelihood for the complete data Dc is given by ( ∏n L φ;D ) = [ ] [ ]c m f y | δi S y | m −δi ( i η) ( i η) i i pθ (mi). (9) i=1 Note that the likelihood (9) is not observable, since it depends on the latent variables. The marginal likelihood for the observed data is obtained by summing overall possible values for the variables Mi, i = 1, . . . , n. Therefore, the logarithm of the marginal likelihood function is given by ∑n (φ;D) = δi log[ fp(yi|φ)]+ (1− δi) log[Sp(yi|φ)]. (10) i=1 3.2 Promotion time model When each random variable Mi follows a Poisson distribution with parameter θi, the unified model comes down to promotion time model (Yin and Ibrahim, 2005). An usual relation used between the ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com 1490 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event parameter θ and the covariates in this model is given by ( ) θi = ′ exp xiβ . (11) From Eq. (8), and using the relation (11), it is derived that the common survival function and the probability density function for the random variables Ti are, respectively, given by{ ′ } Sp(ti|φ) = exp −exiβ [1− S(ti|η)] (12) and | = ′ { } f (t φ) ex βp i i f (ti|η) exp − ′ exiβF (ti|η) . (13) In this case, the cure probability of the i-th observation is related to the covariates through the ′ expression pθ (0) = exp[− exp(xiβ)].i Finally, the logarithm of the marginal likelihood function for this model is obtained by substituting Eqs. (12) and (13) in the logarithm of the marginal likelihood function (10). 3.3 The Weibull and log-logistic RACUF CUSUM chart For Weibull and log-logistic promotion time models, it is assumed that failure times of susceptible individuals follow distributions given in Eq. (5). Here, the covariables are considered to affect only the cure rate, that is, in Eq. (3) we assume that γ = 0. When a process is in-control the random variables Ri j , with j = 1, . . . ,Mi, and i = 1, . . . , n are assumed to follow aWeibull or a log-logistic distribution, with shape parameter α and scale parameter λ = λ0. For out-of-control processes, we suppose there is a change only in the scale parameter for λ1 = ρ1λ0, with other parameters held fixed. We consider 0 < ρ1 < 1 (a reduction in λ) or ρ1 > 1 (an increase in λ) to mean, respectively, a reduction or an increasing in the median lifetime of individuals not immune. From Eq. (1), the CUSUM score for the models is given by Wi = (φ1;Di)− (φ0;Di), (14) where φl = (β′, η ′)′l for l = 0, 1, is a vector of parameters associated with the distribution of Rj , with η0 = (α, λ0) and η1 = (α, λ1) representing the value of the monitored parameter for in-control and out-of-control situations, respectively. Here, Di denotes the data for an individual i. Note that the marginal log-likelihood function for complete data is given by ∑n (φ;D) = (φ;Di). (15) i=1 3.3.1 Weibull promotion time cure model Consider the Weibull distribution to model the lifetime of individuals not immune, with survival and density functions given{, res(pec)tiv}ely, by | = − t α ( )α−1 { ( )α} S(t η) exp i α t and f (t |η) = ii i exp − ti . (16) λ λ λ λ Thus, the long-term survival function Sp(ti|φ), and the subdensity function associated fp(ti|φ) for the Weibull promotion time cure model are obtained by substituting Eq. (16) to Eqs. (12) and (13). The log-likelihood function for a single observation is given by (φ;Di) = ′ δi[xiβ − α log(λ)+ log(αyα−1i )− (yi/λ)α]− ′ exiβ [1− exp(−(yi/λ)α)]. (17) ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1491 Finally, the weight to obtain the RACUF CUSUM chart is [ ( )α ] { [ ( )α] [ ( )α]} W = δ − yα log ρ + i (1− ′ρ−α ) − ex yiβ i yii i 1 λ 1 exp − − exp − . (18)0 λ0 ρ1λ0 3.3.2 Log-logistic promotion time cure model Consider the following survival and density functions to model susceptible individuals [ (t )α]−1 α (t )α−1 [ (t )α]−2 S(ti|η) = 1+ i and f (t |η) = i 1+ ii . (19)λ λ λ λ The long-term survival function and density for the log-logistic promotion time as given in Eqs. (12) and (13) is obtained by the use of the functions given above (19). Hence, the log-likelihood function for a single individual is { ′ ( ) [ (y )α]} (φ;Di) = δi xiβ + lo α−1{g αyi − α log λ−}2 log 1+ i + [ ( ) λα]−1 − ′exp(x β) 1− y1+ ii , (20)λ and the log-likelihood function for complete data under log-logistic is obtained as in Eq. (15). The weight to obtain the log-logistic RACUF CUSUM chart is ⎨⎧ ⎡ ( )y α⎤⎫1+ {i ⎬ [ ( )α]−1 [ ( )α]− }1 Wi = ρ λ δi ⎩−α log(ρ1)− 2 log⎣ ( 1 0) ⎦α ⎭− ′ex β + yi − yi 1 1+ i .(21)1+ yi λ0 ρ1λ0 λ0 Themaximum-likelihood estimates for the parameters inWeibull or log-logistic models are obtained by maximizing the full log-likelihood (15), using the routine optim of R software, through the BFGS method (Broyden-Fletcher-Goldfarb-Shanno). For computational reasons, in the Weibull case (17) the reparametrizations α = exp{α∗} and ξ = −α log(λ) are considered, and, in log-logistic case (20) we use α = exp{−σ ∗} and μ = log(λ). Thus, we obtain parameters that take values in the whole real line. 4 Simulation study The average run length is a common criterion used in literature to assess a CUSUM chart’s perfor- mance.When the process is in-control, the average run length before a false alarm (ARL0) is considered. The out-of-control ARL, called ARL1, considers the average run length to detect a real change in process, and represents the detection speed of the chart. Generally, the ARL0 is fixed (and h obtained according to this), and then the ARL1 is measured for a chart with the same control limit. An extensive simulation study was conducted to perform comparisons of RACUF CUSUM and RAST CUSUM charts. The charts were calibrated to have a fixed ARL0 and the ARL1 of the charts was compared, where early detection of a change in the process is essential for chart performance. We consider samples generated with and without cure fraction and we investigate the effect on the detection speed by using RAST CUSUM (which does not incorporate the cure rate) in data presenting cured individuals, as well as the effect of using RACUF CUSUM, which allows the existence of the cure rate when, in fact, this is not present in the data. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com 1492 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event 4.1 Control limit, ARL0 and ARL1 According to Sego et al. (2009), the control limit h for CUSUM charts can be determined by simu- lations, by setting up an in-control average run length. In this case, an in-control sample of n = 100 survival times was generated from a distribution (conditioned to covariates), involving a censoring proportion and, possibly, a cure rate. Such observations are used as the training sample, to obtain estimates of parameters to in-control process (η0). For a given (true) model, the estimate is used as the true value of the parameter vector, to generate new in-control data. Then, the values Zi (i = 1, . . . , n) of the CUSUM charts (RAST and RACUF) are calculated based on estimates of parameters associated to the two models (AFT and promotion time, respectively). For each chart, we record the first position in which Zi > h, that is, the number of observations for a false alarm on that sample. The idea is to obtain the value of the control limit h, such that the average for K replicas of these observed values (estimated ARL0) is approximately equal to the nominal ARL0 desired. In this work, the nominal ARL0 = 1000 is used, and, for the two models and different parameter combinations, an optimal control limit h is sought, which leads to this value of ARL0. To do so, 1000 in-control samples are generated, each one with 10,000 observations, and based on this data structure, the function f (h) = |1000−ARL0(h)| is minimized, where ARL0(h) is the observed average of the number of observations until a false alarm (Zi > h) in each of the 1000 samples. This is done for both RACUF and RAST CUSUM charts. The golden section search method (Pan and Jarrett, 2013) was used to obtain the optimal control limit h, which minimizes f (h), for both charts in different scenarios. According to Pan and Jarrett (2013) this method reduces the computational effort of minimizing f (h), in comparison with some usual optimization methods. In addition, to check the optimal control limit achieved, the ARL0 of each chart is estimated, considering 100,000 in-control samples of length 10,000, as the average number of observations until a false alarm (Zi > h) in each sample. Samples that do not have false alarm were discarded (Koshti, 2011). Although this discard can cause bias and a more appropriate procedure would be to generate more samples until the signal is reached, we verified in our simulations that this occurrence was quite rare (around 12 occurrences per 100,000 samples). It is hoped that the estimated ARL0 obtained here are close to the desired value 1000. Finally, to evaluate the detection power of the chart, the ARL1 is estimated as the average number of observations until a true alarm, based on 100,000 out-of-control samples with 5000 observations. A detailed description of each step of the simulation procedure is given in Appendix. 4.2 Simulation scenarios The survival times (Ti or Ri j) were generated under aWeibull and a log-logistic model (given in Section 2 or 3), with the covariate xi considered as a scalar and drawn from a Bernoulli random variable with success probability of 0.5 for the i = 1, . . . , n. For data without a cure rate, the survival times Ti were generated as Weibull(α, λe γxi ) or log-logistic(α, λeγxi ), with parameters α = 4, γ = −0.5 and λ = λ0 = 40 for the in-control process and λ = ρ140 for the out-of-control process. For the cure rate data, the latent variable Mi was generated using a Poisson distribution with mean θ = eβxi i , representing the number of risk possibilities for the i-th individual, where different cure rates in the sample are obtained by changing the value of β. For the i-th noncured individual (Mi > 0), random samples from Ri j ∼ Weibull (α, λ) and Ri j ∼ log-logistic(α, λ) of size Mi were generated with parameters α = 4 and λ = λ0 = 40 for the in-control process, and with parameters α = 4 and λ = ρ140 for the out-of-control process. Hence, the failure time is denoted by ti = min{Ri j, j = 1, . . . ,Mi}, i = 1, . . . , n. Six possibilities are considered for the intensity change ρ1, with {0.5; 0.7; 0.9} indicating reduction in survival time, and {1.1; 1.3; 1.5} indicating an increase in survival time. Censoring times Ci j were generated as independent uniform random variables on [0, κ]. Two different censoring proportions ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1493 50% and 70% were taken, which resulted by adequately choosing the value of κ. It is considered here the censoring proportion with respect to all units at risk. Hence, the censoring proportion used in the simulation study can be represented as number of censored or cured units . number of population units Then, none cured rate was considered for AFT model and a cure rate of 50% (β = −0, 77) for the promotion time (PT) model, with censoring proportions of 50% (κ = 58 for AFTmodel and κ = 1000 for PT model) and 70% (κ = 48 for AFT model and κ = 80 for PT model). For data with cure rate, in the case of censoring proportion 50% all censored individuals are actually cured, and, in the case of censoring proportion 70%, around 20% of the data are in fact censoring of individuals susceptible to the occurrence of the event. The observed times and status are yi = min{ti, ci} and δi (δi = 1 if ti ≤ ci and δi = 0 if ti > ci), with i = 1, . . . , n. Maximum-likelihood estimates for the parameter (in-control data) are obtained by maximizing the likelihood functions given in Eq. (4) or (15), which requires a numerical solution. The algorithm was implemented using the R software (R Development Core Team, 2012). Maximization for the promotion time model was obtained considering the intrinsic function optim. For AFT models we use the survreg function in the survival package in R (R Development Core Team, 2012), which uses the Newton-Raphson algorithm (see Therneau and Lumley, 2014). 4.3 Results for samples with cure rate Table 1 presents simulation results obtained for both models (log-logistic and Weibull), when the use of RACUF and RAST CUSUM charts for monitoring data with cure rate of 50% are compared. It is noted there is not a fixed control limit to be used in all cases. The variation of optimal control limit h depends on the intensity of change in the process (ρ1), the censoring proportion in the data and the model considered. For the RACUF CUSUM, most of the optimal limits obtained lie between 4 and 5, with some of them less than 3. For the RAST CUSUM, the variability of control limits is greater, with one of them greater than 6 (in case ρ1 = 0.3, for Weibull data with censoring proportions of 50%) and some of them less than 2 (most of the cases with ρ1 = 0.9). With respect to the average run length performance, when the process is in-control, the ARL0 stays close to 1000 in all situations, as expected in the simulation. Here, a better performance for RACUF CUSUM chart to detect a real change in a process is observed, since this chart presents ARL1 with lower values than those obtained from the RAST CUSUM chart, especially for slight changes in the process (ρ1 = 0.9). So, it is understood that, if the data contain cured individuals, the use of RAST CUSUM chart is not efficient to detect a change in the survival quality, whereas the RACUFCUSUM chart is highly suggested. This result is obviously expected, since the RACUF CUSUM chart uses a survival model that considers the existence of a cure rate in the data, and so it is more suitable. In the next section, the performance of the charts in data without a cure rate (generated from an AFTmodel) is evaluated. 4.4 Results for samples without cure rate The results obtained in the simulation of data without cure rates are presented in Table 2. As noted before, the optimal control limits also change in the different scenarios for data without cure rate. Although most of the control limits obtained are close to 5, for ρ1 = 0.9, these values are smaller. The values of ARL0 are close to 1000, but, unlike the results observed for the lifetime data with cure rate, a similar performance for both charts is seen in the detection of a real change in survival quality. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com 1494 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event Table 1 Optimal control limits (h), estimated average run lengths (ARL0,ARL1) and their respective standard error (SE0, SE1) for RAST and RACUFCUSUM, based onWeibull and log-logistic models, considering 100,000 replications of the procedure with samples generated with a cure rate of 50%. Censoring/model CUSUM chart ρ1 h ARL0 SE0 ARL1 SE1 RACUF 0.5 4.72 1036 3.21 11 0.02 0.7 4.21 1042 3.25 27 0.05 0.9 2.64 982 2.94 134 0.28 1.1 2.33 988 2.91 174 0.37 1.3 3.85 1007 3.09 44 0.08 50% 1.5 4.03 1035 3.17 37 0.07 log-logistic RAST 0.5 3.59 1005 3.03 161 0.36 0.7 3.06 1038 2.86 315 0.69 0.9 1.38 1054 2.60 678 1.48 1.1 1.22 998 2.37 678 1.40 1.3 1.82 988 2.91 428 1.09 1.5 4.54 1035 2.75 290 0.51 RACUF 0.5 4.78 966 3.02 8 0.02 0.7 4.26 973 3.02 21 0.04 0.9 2.40 977 2.89 157 0.35 1.1 2.15 1010 2.91 226 0.47 1.3 3.71 1031 3.14 61 0.10 70% 1.5 4.41 996 3.10 28 0.04 log-logistic RAST 0.5 4.22 952 3.06 16 0.04 0.7 4.31 1003 3.15 59 0.13 0.9 1.94 1038 3.10 241 0.63 1.1 2.56 983 2.69 342 0.70 1.3 2.51 1012 3.10 166 0.40 1.5 3.51 976 3.00 82 0.17 RACUF 0.5 5.06 1011 3.15 6 0.01 0.7 4.86 1028 3.19 12 0.02 0.9 3.44 1017 3.07 79 0.13 1.1 2.66 1009 3.05 109 0.25 1.3 3.95 1008 3.16 15 0.04 50% 1.5 3.79 962 3.05 8 0.02 Weibull RAST 0.5 3.59 1026 3.03 319 0.91 0.7 1.95 989 3.03 525 1.54 0.9 1.53 1054 2.50 841 1.88 1.1 1.50 982 2.23 818 1.74 1.3 2.38 1022 2.90 628 1.68 1.5 4.39 1029 2.82 538 1.22 RACUF 0.5 4.78 965 3.02 8 0.01 0.7 4.55 1006 3.15 19 0.03 0.9 3.31 956 2.89 75 0.14 1.1 2.49 966 2.91 132 0.30 1.3 3.57 986 3.04 37 0.08 ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1495 Table 1 Continued Censoring/model CUSUM chart ρ1 h ARL0 SE0 ARL1 SE1 70% 1.5 4.12 1035 3.22 29 0.06 Weibull RAST 0.5 5.00 987 3.10 38 0.10 0.7 4.75 1031 3.19 72 0.19 0.9 3.90 1035 2.90 319 0.73 1.1 2.64 965 2.80 352 0.82 1.3 6.14 1050 3.08 154 0.25 1.5 4.41 978 3.00 69 0.14 In fact, the values of ARL1 for RAST and RACUF CUSUM charts are similar in every situation. As in the previous case, the main differences observed for ARL1 occur when ρ1 = 0.9. But, even in such case, the differences are far smaller than those obtained with the cure rate data. 4.5 Studying estimation error To evaluate the estimation error we made a brief simulation, where we ran the chart with the estimated parameters (the same obtained in Tables 1 and 2), but generated the data from the distribution as we generated the original training sample. The results considering only censoring proportion of 70% are shown in Table 3. It can be seen that, although the obtained values differ from the values in the previous simulation due to estimation error (mainly when ρ1 is near 1), the conclusions about the detection speed (ARL1) are similar. For example, for the log-logistic case with ρ1 = 0.7 and a 50% cure rate, the values of ARL1 obtained with RACUF and RAST are, respectively, 33 and 73 in Table 3, and 21 and 59 in Table 1, indicating a better performance for RACUF in both tables. In the similar case without a cure rate (log-logistic, ρ1 = 0.7), the values of ARL1 obtained with RACUF and RAST are, respectively, 23 for the two charts in Table 3, and 24 and 21 in Table 2, indicating in both cases a similar result for RACUF and RAST. 5 Applications To better illustrate the performance of the RACUFCUSUM chart and compare with RASTCUSUM two examples are presented. In the first, simulated data is used to show graphically one specific scenario presented in the simulation study. In the second case, a retrospective analysis based on real data of patients with heart failure who were followed for almost two decades is considered. 5.1 Illustration with simulated data For a simulated dataset, we took the first 200 observations for the in-control sample and the other 100 for out-of-control. The in-control sample was used as historical data to estimate the parameters. Considering initially simulated data from a promotion time Weibull model, a cure rate of 50% is simulated, with a censoring proportion of 71% and considering ρ1 = 0.7, which represents a reduction of 30% in the parameter λ. It is noted the Weibull RACUF CUSUM chart for simulated data (Fig. 1a) can signal the change quickly, while theWeibull RASTCUSUMperforms poorly (Fig. 1b). Thus, when the existence of cured individuals is ignored, and this method of monitoring is applied, it can lead to a lower performance in detecting changes in the quality of process. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com 1496 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event Table 2 Optimal control limits (h), estimated average run lengths (ARL0,ARL1) and their respective standard error (SE0, SE1) for RAST and RACUF CUSUM, based on the Weibull and log-logistic models, considering 100,000 replications of the procedure with samples generated without cure rate. Censoring/model CUSUM chart ρ1 h ARL0 SE0 ARL1 SE1 RACUF 0.5 5.14 964 3.06 7 0.01 0.7 4.71 1040 3.12 16 0.03 0.9 2.84 993 3.00 128 0.28 1.1 4.03 1001 2.87 154 0.27 1.3 5.09 1019 3.17 28 0.05 50% 1.5 4.33 1035 3.23 27 0.05 log-logistic RAST 0.5 5.22 1006 3.18 6 0.01 0.7 4.71 996 3.09 16 0.03 0.9 3.01 996 3.01 103 0.21 1.1 2.76 999 2.98 133 0.27 1.3 4.57 1057 3.29 24 0.04 1.5 4.54 1052 3.26 26 0.04 RACUF 0.5 4.85 1004 3.23 8 0.02 0.7 4.36 1000 3.12 27 0.05 0.9 2.65 1060 3.20 127 0.28 1.1 6.17 1051 2.61 234 0.33 1.3 3.66 963 2.93 70 0.14 70% 1.5 4.76 998 3.11 32 0.06 log-logistic RAST 0.5 5.08 1016 3.20 7 0.01 0.7 4.34 984 3.09 25 0.05 0.9 2.88 1023 3.08 116 0.25 1.1 2.47 1006 2.95 173 0.35 1.3 3.44 950 2.89 66 0.12 1.5 4.34 1004 3.11 32 0.06 RACUF 0.5 5.56 998 3.13 5 0.01 0.7 5.20 987 3.08 12 0.02 0.9 3.36 1009 3.20 83 0.19 1.1 2.09 963 2.77 221 0.47 1.3 4.27 1010 3.13 24 0.06 50% 1.5 3.86 1007 3.13 11 0.03 Weibull RAST 0.5 5.42 1028 3.24 5 0.01 0.7 5.16 961 3.00 11 0.01 0.9 3.72 1008 3.09 63 0.11 1.1 1.80 1012 2.99 297 0.74 1.3 4.22 1004 3.15 19 0.05 1.5 4.39 964 3.02 10 0.02 RACUF 0.5 5.15 975 3.04 6 0.01 0.7 4.65 938 2.99 18 0.03 0.9 3.17 971 2.95 88 0.17 1.1 2.86 995 3.03 89 0.19 1.3 3.47 997 3.02 29 0.06 ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1497 Table 2 Continued Censoring/model CUSUM chart ρ1 h ARL0 SE0 ARL1 SE1 70% 1.5 4.27 1011 3.16 16 0.03 Weibull RAST 0.5 5.29 986 3.09 5 0.01 0.7 4.83 975 3.07 16 0.03 0.9 3.31 1016 3.10 84 0.16 1.1 3.04 1027 3.11 86 0.18 1.3 4.00 987 3.03 28 0.06 1.5 4.34 1000 3.14 15 0.03 On the other hand, making the comparison between procedures on simulated data without a cure rate (based on an AFTmodel) and with a censoring proportion of 72%, theWeibull RACUFCUSUM chart (Fig. 2a) has similar performance to theWeibull RAST CUSUM chart (Fig. 2b), which indicates that, in this example, the proposed chart can be efficient, even without the presence of a cure rate in the data. 5.2 Monitoring the survival of patients with heart failure In a study conducted at the Heart Institute (InCor) of FMUSP (Faculdade de Medicina da Universi- dade de São Paulo), Brazil, two cohorts of patients with heart failure treated at the General Outpatient Clinic were analyzed (see Marçula et al., 2011). From March 1991 to June 2003, 1353 patients were evaluated (cohort 1) and 2124 consecutive patients were treated from July 2003 to June 2007 (cohort 2). Besides to find prognostic variables from demographic, clinical and laboratory data obtained in the initial evaluation of these patients, other objective was to compare the two cohorts. Considering the Kaplan–Meier curve of the data for the first 975 patients of cohort 1 (Fig. 3), the occurrence of about 20% of censoring at the end of the study was noted in a reasonable follow-up time, which is an indication of cure rate in the data (Maller and Zhou, 1996). When a promotion time Weibull model is fit to this data, 15.77% is found to be an estimative of the general cure rate. In order to investigate any change in survival of patients over time, and the association with the two cohorts, RACUF CUSUM and RAST CUSUM were used to monitor retrospectively the data. As historical in-control data, the lifetimes (in years) of the patients from March 1991 to December 1998 in cohort 1 were considered, and the lifetimes of the other patients in cohort 1 and cohort 2 were monitored. To incorporate the specific risk for each individual we use the covariates given in Marçula et al. (2011). For RACUF CUSUM, we assume the promotion time Weibull model (18) to model the lifetimes. For the RAST CUSUM chart risk adjustment is made through an accelerated failure time Weibull model (6). The reduction (ρ1 = 0.9) or increase (ρ1 = 1.1) in the average lifetime of the (not immune) individuals was evaluated. We use optimal control limits h obtained for each case, which lead to an ARL0 about 2000 (since the historical data have around 1000 observations, we can be more rigorous on the false alarm rate). To obtain the control limit in this application (as described in Section 4.1)we should simulate samples considering censoring, cure fraction and the 9 covariates used in the study (Marçula et al., 2011), what would be computationally hard. Instead, we chose to use here a simpler resampling procedure. We select at random and with replacement 1000 samples of size 10,000 from historical data, then we use the procedure given in Section 4.1 to obtain the control limit. To validate this approach we run some additional simulations (Table 4). For that we consider some covariates (drawn from Bernoulli random ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com 1498 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event Table 3 Optimal control limits (h), estimated average run lengths (ARL0,ARL1) and their respective standard error (SE0, SE1) for RAST and RACUF CUSUM based on the Weibull and log-logistic models to study the estimation error, considering the censoring proportion of 70%. Model CUSUM chart ρ1 h ARL0 SE0 ARL1 SE1 RACUF 0.5 4.82 994 3.14 12 0.02 0.7 5.19 1025 3.19 33 0.07 0.9 2.92 995 3.00 175 0.43 1.1 1.15 1011 3.10 256 0.74 1.3 5.16 1042 3.10 65 0.10 Log-logistic 1.5 5.41 1008 3.11 38 0.07 (cure rate of 50%) RAST 0.5 3.41 1027 3.23 20 0.05 0.7 3.59 1026 3.18 73 0.19 0.9 2.64 993 2.91 333 0.74 1.1 1.56 1016 3.05 413 1.16 1.3 2.38 1011 3.06 138 0.33 1.5 3.65 1017 3.14 79 0.17 RACUF 0.5 4.44 1016 3.19 10 0.02 0.7 3.63 986 3.08 23 0.05 0.9 1.73 962 2.94 157 0.45 1.1 2.36 1036 2.90 181 0.32 1.3 2.51 1020 3.14 71 0.19 Log-logistic 1.5 3.88 976 3.07 30 0.06 (without cure rate) RAST 0.5 4.82 986 3.09 7 0.01 0.7 4.25 973 3.01 23 0.04 0.9 1.79 1012 3.11 155 0.44 1.1 3.09 990 2.65 183 0.27 1.3 2.28 1019 3.14 66 0.17 1.5 3.64 1039 3.23 29 0.05 RACUF 0.5 5.09 1005 3.16 11 0.02 0.7 7.12 983 3.05 24 0.04 0.9 3.75 984 2.95 128 0.26 1.1 1.38 1039 3.22 166 0.47 1.3 3.95 977 3.06 39 0.09 Weibull 1.5 3.33 997 3.07 24 0.05 (cure rate of 50%) RAST 0.5 5.03 977 3.02 35 0.09 0.7 5.65 1021 3.11 105 0.27 0.9 2.31 1010 3.00 359 0.97 1.1 2.28 972 2.86 425 1.12 1.3 4.49 956 2.90 145 0.32 1.5 4.26 1007 3.06 73 0.16 RACUF 0.5 4.01 1020 3.20 7 0.01 0.7 4.81 984 3.05 16 0.02 0.9 6.35 1010 2.80 114 0.14 1.1 2.97 982 2.93 110 0.22 1.3 3.24 1020 3.18 35 0.06 ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1499 Table 3 Continued Model CUSUM chart ρ1 h ARL0 SE0 ARL1 SE1 Weibull 1.5 3.91 992 3.09 17 0.04 (without cure rate) RAST 0.5 4.14 1003 3.16 6 0.01 0.7 4.98 988 3.08 15 0.02 0.9 5.37 991 2.91 97 0.14 1.1 2.93 981 2.95 108 0.23 1.3 3.78 1017 3.16 33 0.06 1.5 4.20 1024 3.21 16 0.04 0 50 100 150 200 250 300 0 50 100 150 200 250 300 Observation Observation Figure 1 Monitoring by the charts compared in a Weibull process with cure rate (50%), taking λ = 40, α = 4, ρ1 = 0.7, β = −0.77 and censoring proportion of 71%. Control limits: h = 4.55 (RACUF CUSUM chart) and h = 4.75 (RAST CUSUM chart). variables with success probability of 0.5) and a scenario similar to what we observe in the data (20% cure fraction and samples of size 1000 as historical data). In Table 4 the optimal control limit was obtained as described in Section 4 (h) and using the resampling approach (h∗). We note that the results of h and h∗ were similar to each other, which justifies the resampling procedure. In Fig. 4, the charts show the CUSUM scores for an increase in average lifetime (points plotted above zero), and for a reduction (which are presented as negative values, plotted below zero, with the control limit h considered as negative). The horizontal dashed lines represent the control limits and the vertical dotted line indicates the last observation used as historical data. We hope that the chart signal only after this observation. The charts clearly signal only an increase in the average lifetime of patients (of 10%). For the RACUF CUSUM chart, the first detection occurs in the observation 1032 (August 6, 1999) and for the RAST CUSUM chart, the first detection occurs in the observation 1052 (October 1, 1999). The RASTCUSUMchart presents a false alarm (monitoring an increase, which signals for the observation 627). The RACUF CUSUM chart did not present a false alarm and detected the increase in survival 20 observations before the RAST CUSUM chart did, so it was more effective. In fact, similar plots (not shown here) were constructed to evaluate a reduction or increase of 50% (ρ1 = 0.5 or ρ1 = 1.5) and they exhibit the same type of behavior. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Zi 0 5 10 15 Zi 0 1 2 3 4 5 6 1500 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event 0 50 100 150 200 250 300 0 50 100 150 200 250 300 Observation Observation Figure 2 Monitoring by charts compared in aWeibull process without a cure rate, taking λ = 40, α = 4, ρ1 = 0.7, γ = −0.5) and censoring proportion of 72%. Control limits: h = 4.65 (RACUF CUSUM chart) and h = 4.83 (RAST CUSUM chart). 0 5 10 15 Time (years) Figure 3 Empirical survivor function (Kaplan–Meier) for the first 975 patients of cohort 1 of the data. 6 Discussion We propose a risk-adjusted CUSUM chart based on a cure rate model to monitor lifetime data in which part of the population studied will never experience the event of interest. Using the promotion time regression model, a CUSUM score is obtained when the promotion times follow a Weibull and a log-logistic distribution to monitor their scale parameters, and considering the covariables affect only the cure rate. A simulation study to make comparisons between the proposed ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Zi 0 5 10 15 20 25 Survival 0.0 0.2 0.4 0.6 0.8 1.0 Zi 0 5 10 15 20 25 30 Biometrical Journal 58 (2016) 6 1501 Table 4 Comparison of control limits obtained by simulation (h) and control limits obtained by the resampling procedure (h∗) for Weibull model, considering a censoring proportion of around 50%, and training data samples of size n0 = 1000. Cure rate CUSUM chart ρ1 One covariate Two covariates Three covariates h h∗ h h∗ h h∗ 20% RACUF 0.5 5.18 5.24 5.96 5.18 5.18 5.08 0.7 4.80 4.95 5.90 5.08 4.94 5.04 0.9 3.26 3.24 4.38 3.70 3.42 3.53 1.1 2.69 2.92 3.45 3.33 3.04 3.17 1.3 4.01 4.23 4.80 4.33 3.70 3.40 1.5 4.48 4.53 5.02 4.44 4.13 4.51 RAST 0.5 5.23 5.18 5.96 5.57 5.57 5.69 0.7 5.26 5.05 5.52 5.38 5.02 5.14 0.9 3.72 3.36 3.59 3.48 3.01 3.25 1.1 3.22 3.24 2.96 3.04 3.02 2.89 1.3 5.48 5.50 4.16 4.21 3.33 3.52 1.5 6.42 6.34 4.50 4.53 4.26 4.27 0% RACUF 0.5 5.75 5.35 4.55 5.01 5.11 4.75 0.7 5.51 5.07 4.71 4.48 4.72 4.66 0.9 4.00 3.51 3.12 3.10 3.06 3.27 1.1 2.60 3.11 2.62 2.81 3.19 2.82 1.3 3.90 4.26 3.54 4.16 4.04 3.64 1.5 3.59 3.83 3.54 4.76 4.01 3.59 RAST 0.5 5.58 5.26 4.79 5.04 5.47 4.98 0.7 5.53 5.17 4.82 4.58 5.13 4.77 0.9 4.25 3.71 3.32 3.14 3.47 3.30 1.1 2.82 3.22 2.61 2.86 3.00 2.84 1.3 4.25 4.74 3.57 4.04 3.98 3.69 1.5 4.18 4.26 3.73 4.49 4.11 3.70 h = 1.76 h = 2.03 h = −1.88 h = −2.21 0 500 1500 2500 3500 0 500 1500 2500 3500 Observation Observation Figure 4 Weibull RACUF and RAST CUSUM charts for monitoring of reduction or increase of 10% in the mean survival time. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Zi 0 10 20 30 Zi 0 10 20 30 1502 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event RACUF CUSUM and the RAST CUSUM charts is considered. As the first result, we note that, with respect of ARL1 of the charts, the RACUF CUSUM performs much better than RAST CUSUM to monitor data with cure rate, and seems to be equivalent in monitoring data without a cure rate. In fact, we also consider some simulations involving the mixture cure model (not shown here) and the conclusions were the same. We think this advantage is due to the added flexibility of cure rate models to fit long-term survivors. So we believe that if there is any doubt on the existence of cured in the data, it is worth using the CUSUMRACUF to monitor survival times. In addition, we emphasize the importance of obtaining optimal control limits for each situation, since they may vary according to specific data or model assumption. In the application with heart failure data, it was possible to show the increase in survival (median) of patients by monitoring data with risk-adjusted CUSUM charts. Although expected, this fact had not been identified in other studies. It is pointed out that, even with a small estimated cure rate in this application, the RACUF detected the change more quickly than the RAST. In fact, the simulation suggests that the higher the cure rate, the greater the advantage of the RACUF. In the simulation study, we generate data from an AFT model (where the covariates act on the scale parameter) to be monitored by a CUSUM chart based on a PT model (where the covariates act on the cure rate), and vice versa. Thus we are apparently violating the assumption that the correct model is used and that the maximum likelihood estimate of the parameter (using historical data) is very close to the true value (for an in-control process). However, in many practical situations it is not possible to assess whether the used model is correct and, consequently, if the assumption about the form of action of the covariates is correctly specified. By White (1982), it could be argued that, once valid the regularity conditions, the maximum likelihood estimators γ̂ and β̂ are, respectively, consistent estimators for the values γ ∗ and β∗ that minimize “our ignorance” about the real model. Furthermore, as the “wrong choice” favors a model in each case (see Tables 1 and 2), the RACUF CUSUM seems to be more robust for model misspecification. Although the existence of cured individuals is considered, in this approach only a parameter as- sociated to non-immune individuals is monitored, and the cure rate is considered unaffected by the change in process. In practice, some parameter associated to cure rate might need to be monitored. For example, when a medical procedure starts to have a bad response, a reduction of the cure rate in patients can be expected. Thus, an important topic for a future study is to monitor parameters associated to cure rate. Regarding the estimation error problem, we have made a brief simulation based on just one training sample to study it (see Table 3).However, an anonymous referee observes that to get a good understand- ing of the impact of estimation error we need to repeat the estimation from different training samples from the same underlying true model several times. This is an important topic to be investigated in future works. Also, a control chart might need to be used in a prospective context, which is not considered in this work. Thus, developing an updating control chart as in Biswas and Kalbfleisch (2008) and Steiner and Jones (2009) is a topic for future research. Besides, we intend to provide a friendly routine in R language’s C interface to obtain the optimal control limit for the RACUF CUSUM. Acknowledgments This research was supported in part by the Coordenação de Aperfeiçoamento de Pessoal de Nı́vel Superior (CAPES) and Conselho Nacional de Desenvolvimento Cientı́fico e Tecnológico (CNPq), Brazil. The authors thank the referees for their careful reading and important suggestions. Conflict of interest The authors have declared no conflict of interest. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1503 Appendix Description of simulation procedure for a fixed value of ρ1 1. Obtaining optimal control limits for models: (a) Generate a sample of size n0, according to the real model (the correctly specified model), denoted MR(φR0), where φR0 is a vector of parameters associated with this model when the process is under control. In practice, these data represent the historical data. (b) Get φ̂R0 the maximum-likelihood estimator of φR0, based on the historical data. (c) Replicate mh samples of size Nh from the model MR(φ̂R0), obtaining the data matrix Dh = Dh(φ̂R0), with dimension mh × Nh. (d) For each element Dh of Dhi j and considering the estimated φ̂R0, calculate the CUSUM score for the real model, denoted by Z hRi j = ZR(Di j, φ̂R0, ρ1), thus obtaining the mh × Nh matrix of scores, ZR = [ZRi j ]. (e) Based on ZR, obtain the optimal control limit for the real model: hR = hR(Dh,ZR) .1 (f) Given the false model (the improperly specified model), denoted by MF (φF ), where φF is a vector of parameters associated with this model, let φ̂F be the maximum-likelihood estimator of φF based on historical data (n0). (g) For each element Dh hi j of D and considering the estimated φ̂F , calculate the CUSUM score for the false model, denoted by Z hF i j = ZF (Di j, φ̂F , ρ1), obtaining a mh × Nh matrix of scores ZF = [ZF i j ]. (h) Based on ZF , obtain the optimal control limit for the false model: h h 1 F = hF (D ,ZF ) . 2. Checking the closeness of ARL0(h) with respect to the nominal value M: (a) Replicate m0 samples of size N0 from the model MR(φ̂R0), obtaining the data matrix D 0 = D0(φ̂R0) with dimension m0 × N0. (b) For each element D0 of D0i j calculate the CUSUM score for the real model, denoted by Z0Ri j = Z0R(D0i j, φ̂R0, ρ1), obtaining the m0 × N0 matrix of scores, Z0 0R = [ZRi j ]. (c) Compare each line of Z0R with hR, to calculate the ARL R 0 of the real model, as the average number of observations to a false alarm, based on m0 lines of this matrix. (d) For each element D0i j of D 0 calculate the CUSUM score for the false model, denoted by Z0 = Z0 (D0F i j F i j, φ̂F , ρ1), obtaining the m0 × N0 matrix of scores, Z0F = [Z0F i j ]. (e) Compare each line of Z0F with hF , to calculate the ARL F 0 of the false model, as the average number of observations to a false alarm, based on the m0 lines of this matrix. 3. Calculation of ARL1 for models: (a) Replicate m1 samples of size N1 from the model MR(φ̂R1) where φ̂R1 represents the value of the estimated parameter for an out-of-control process, that is φ̂R1 = ρ1 ∗ φ̂R0, obtaining the data matrix D1 = D1(φ̂R1) with dimension m1 × N1. (b) For each element D1 of D1i j calculate the CUSUM score for the real model, denoted by Z1Ri j = Z1R(D1 1 1i j, φ̂R0, ρ1), obtaining the m1 × N1 matrix of scores, ZR = [ZRi j ]. (c) Compare each line of Z1R with hR, to calculate the ARL R 1 of the real model, as the average run length to detect the real change in the process, based on m1 lines of this matrix. (d) For each element D1i j of D 1 calculate the CUSUM score for the false model, denoted by Z1 1 1F i j = ZF (Di j, φ̂F , ρ1), obtaining the m1 × N1 matrix of scores, Z1 = [Z1F F i j ]. 1Golden section search method to minimize f (h) = |M −ARL0(h)|, where M is the nominal ARL0. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com 1504 J. W. Oliveira et al.: Risk-adjusted monitoring of time to event (e) Compare each line of Z1F with hF , to calculate the ARL F 1 of the false model, as the average run length to detect the real change in the process, based on m1 lines of this matrix. References Biswas, P. and Kalbfleisch, J. D. (2008). A risk-adjusted CUSUM in continuous time based on the Cox model. Statistics in Medicine 27, 3382–3406. Chen, M. H., Ibrahim, J. G. and Sinha, D. (1999). A new Bayesian model for survival data with a surviving fraction. Journal of the American Statistical Association 94, 909–919. Fonseca, R. S., Valença, D. M. and Bolfarine, H. (2013). Cure rate survival models with missing covariates: a simulation study. Journal of Statistical Computation and Simulation 83, 97–113. Gandy, A. andKvaløy, J. T. (2013). Guaranteed conditional performance of control charts via bootstrapmethods. Scandinavian Journal of Statistics 40, 647–668. Gandy, A., Kvaløy, J. T., Bottle, A. and Zhou, F. (2010). Risk-adjusted monitoring of time to event. Biometrika 97, 375–388. Grigg, O. and Farewell, V. (2004). An overview of risk-adjusted charts. Journal of the Royal Statistical Society: Series A (Statistics in Society) 167, 523–539. Grigg, O. A., Farewell, V. T. and Spiegelhalter, D. J. (2003). Use of risk-adjusted CUSUM and RSPRTcharts for monitoring in medical contexts. Statistical Methods in Medical Research 12, 147–170. Jones, L. A., Champ, C. W. and Rigdon, S. E. (2004). The run length distribution of the CUSUM with estimated parameters. Journal of Quality Technology 36, 95–108. Koshti, V. V. (2011). Cumulative sum control chart. International Journal of Physics and Mathematical Sciences 1, 28–32. Lawless, J. F. (2003). Statistical Models and Methods for Lifetime Data 2nd edn. John Wiley, New York, NY. Maller, R. A. and Zhou, X. (1996). Survival Analysis with Long-term Survivors. Wiley, New York, NY. Marçula, M., Cuoco, M. A. R., Yamada, A. T., Freitas, H. F. G., DePaula, R. S. T., Uemura, R. T., Ferreira, P., Alencar, A. P., Barretto, A. C. P. and Mansur, A. J. (2011). Comparison of prognostic variables between two cohorts of outpatients with heart failure in the last two decades. Paris European Heart Journal (Abstract Supplement) 32, 459. Montgomery, D. C. (2007). Introduction to Statistical Quality Control. John Wiley & Sons, New York, NY. Moustakides, G. V. (1986). Optimal stopping times for detecting changes in distributions. The Annals of Statistics 14, 1379–1387. Page, E. S. (1954). Continuous inspection schemes. Biometrika 41, 100–115. Pan, X. and Jarrett, J. E. (2013). On quality control chart construction and simulation. College of Business Working Paper Series, University of Rhode Island. Phinikettos, I. and Gandy, A. (2014). An omnibus CUSUM chart for monitoring time to event data. Lifetime Data Analysis 20, 481–494. R Development Core Team (2012). R: A Language and Environment for Statistical Computing. Vienna, AT. Available at: . Rodrigues, J., Cancho, V. G., de Castro,M. and Louzada-Neto, F. (2009). On the unification of long-term survival models. Statistics & Probability Letters 79, 753–759. Sego, L. H., Reynolds, M. R. and Woodall, W. H. (2009). Risk-adjusted monitoring of survival times. Statistics in Medicine 28, 1386–1401. Steiner, S. H., Cook, R. J. and Farewell, V. T. (2001). Risk-adjusted monitoring of binary surgical outcomes. Medical Decision Making 21, 163–169. Steiner, S. H., Cook, R. J., Farewell, V. T. and Treasure, T. (2000). Monitoring surgical performance using risk-adjusted cumulative sum charts. Biostatistics 1, 441–452. Steiner, S. H. and Jones, M. (2009). Risk-adjusted survival time monitoring with an updating exponentially weighted moving average (EWMA) control chart. Statistics in Medicine 29, 444–454. Sun, R. J., and Kalbfleisch, J. D. (2013). A risk-adjusted OE CUSUM with monitoring bands for monitoring medical outcomes. Biometrics 69, 62–69. Sun, R. J., Kalbfleisch, J. D. and Schaubel, D. E. (2014). A weighted cumulative sum (WCUSUM) to monitor medical outcomes with dependent censoring. Statistics in Medicine 33, 3114–3129. Therneau, T. and Lumley, T. (2014). Survival: survival analysis, including penalised likelihood. S original by Terry Therneau and ported by Thomas Lumley. R package version 2.37-7. R Foundation for Statistical ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com Biometrical Journal 58 (2016) 6 1505 Computing, Vienna. CRAN. Available at: . Tournoud, M. and Ecochard, R. (2007). Application of the promotion time cure model with time-changing exposure to the study of HIV/AIDS and other infectious diseases. Statistics in Medicine 26, 1008–1021. White, H. (1982). Maximum likelihood estimation of misspecified models. Econometrica: Journal of the Econo- metric Society 50, 1–25. Woodall, W. H. (2006). The use of control charts in health-care and public-health surveillance. Journal of Quality Technology 38, 89–104. Yakovlev, A. Y., Asselain, B., Bardou, V. J., Fourquet, A., Hoang, T., Rochefediere, A. and Tsodikov, A. D. (1993). A simple stochastic model of tumor recurrence and its applications to data on premenopausal breast cancer. Biometrie et Analyse de Donnees Spatio-temporelles 12, 66–82. Yin,G. and Ibrahim, J.G. (2005). Cure ratemodels: a unified approach.Canadian Journal of Statistics 33, 559–570. Zhang, M., Xu, Y., He, Z. and Hou, X. (2016). The effect of estimation error on risk-adjusted survival time CUSUM chart performance. Quality and Reliability Engineering International, 32, 1445–1452. ©C 2016 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.biometrical-journal.com