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