Monitoring time to event in registry data using CUSUMs based on relative survival models

Abstract

An aspect of interest in surveillance of diseases is whether the survival time distribution changes over time. By following data in health registries over time, this can be monitored, either in real time or retrospectively. With relevant risk factors registered, these can be taken into account in the monitoring as well. A challenge in monitoring survival times based on registry data is that the information related to cause of death might either be missing or uncertain. To quantify the burden of disease in such cases, relative survival methods can be used, where the total hazard is modelled as the population hazard plus the excess hazard due to the disease.

We propose a cumulative sum (CUSUM) procedure for monitoring for changes in the survival time distribution in cases where the use of excess hazard models is relevant. The CUSUM chart is based on a survival log-likelihood ratio and extends previously suggested methods for monitoring of time to event data to the excess hazard setting. The procedure takes into account changes in the population risk over time, as well as changes in the excess hazard which is explained by observed covariates. Properties, challenges and an application to cancer registry data will be presented.

Keywords

Excess hazard relative survival survival time monitoring risk-adjusted monitoring cancer registry data

1. Introduction

In health registries, such as cancer registries, patients are routinely registered with key disease related and demographic data at diagnosis of disease, and outcome data like survival times are usually added later depending on the characteristics and clinical course of the key disease. Based on such data, an aspect of interest could be to monitor whether the distribution of the time to an outcome of interest changes over time, for instance if the survival times of cancer patients due to the disease itself change over time while adjusting for known prognostic factors. Such monitoring could be of interest both for real-time monitoring of incoming data and for retrospective analyses to pinpoint when in time important changes took place. This could be of value both to monitor for unexpected changes and to monitor if changes in treatment policies can be associated with better outcomes, for example survival benefit for cancer patients.

A common challenge in monitoring survival times based on such registry data is that time to death, but not necessarily cause of death, is registered. In addition, even if this information is available, it is often unreliable at population level. Therefore, to quantify the burden of disease in such cases, relative survival methods can be used. Under this framework, the frequent standard assumption is that the total hazard is modelled as the population hazard plus the excess hazard due to the disease, see for example Dickman et al.,¹ Perme et al.² and Perme et al.³ for an overview of notions in relative survival and regression models regarding this topic. The population hazard is usually retrieved from national population life tables, stratified on demographic variables.

Methods from statistical process control have been essential in many different fields and applications, originally suggested for monitoring processes in various industries. One of the most widely used control charts is the CUSUM (cumulative sum) presented by Page.⁴ Compared to other charts, for example Shewhart charts, the CUSUM chart is known to perform better for detecting smaller shifts that are persistent over time, thus being suitable for a range of medical settings. Over the years, extensions and adaptions of CUSUM charts to a number of different scenarios have been developed. CUSUM procedures for monitoring of ordinary time to event data were first proposed by Biswas and Kalbfleisch,⁵ Sego et al.,⁶ Gandy et al.⁷ and Steiner and Jones.⁸ These have since been extended to contexts like frailty models,^9,10 cure models,¹¹ illness-death models¹² and queue models.¹³ Impact of estimation error in time to event data monitoring was studied by Zhang et al.¹⁴ and effectiveness versus periodic evaluations was studied by Massarweh et al.¹⁵ Applications have been demonstrated in, for instance, monitoring of perioperative mortality.^16,17

We propose a CUSUM procedure for monitoring for changes in the time to event distribution when the use of excess hazard models is relevant, for example for monitoring based on registry data with uncertain or missing cause of death information, a typical example being cancer registry data. The CUSUM chart is based on a survival log-likelihood ratio and extends the literature discussed above on monitoring of time to event to the relative survival setting. The procedure takes into account changes in the population risk over time, as well as changes in the excess hazard which are explained by observed covariates. Properties, challenges and an application to cancer registry data will be presented.

The structure of the paper is as follows: Section 2 introduces the set-up and notation and presents the proposed method. In Section 3, numerical simulations and experiments are carried out to demonstrate the use of the method and how it performs under different scenarios, including studying the impact of estimation error. An application of the method on a real data set obtained from the Norwegian Cancer Registry is illustrated in Section 4. Finally, Section 5 provides some concluding remarks. R-code implementing the proposed methods and for running the simulations in Section 3 can be found at https://github.com/jihut/cusum_relative_survival_simulations.

2. CUSUM chart based on excess hazard models in relative survival setting

2.1. Set-up and notation

The setting considered is monitoring time to event under the assumption of an additive hazard model $h (\cdot) = h_{P} (\cdot) + h_{E} (\cdot)$ . Here, $h_{E} (\cdot)$ denotes the excess hazard due to the disease of interest. In theory, $h_{P} (\cdot)$ represents the hazard related to other causes. However, it is almost always the case that this quantity is substituted by the population hazard, which is assumed to be known and usually extracted from population life tables. In real time, we would like to monitor for changes in $h_{E} (\cdot)$ over calendar time.

For ease of notation, consider a monitoring system starting from a given calendar date which is now defined as the origin time point $t = 0$ of the monitoring period, for instance the first day of a particular year. Let $0 \leq B_{1} \leq B_{2} \leq \dots$ denote the arrival times after the starting point at which individuals enter the monitoring system. This could for instance be when a patient is diagnosed and enters the cancer registry system. Let $Z_{1}, Z_{2}, \dots$ denote the corresponding vectors of demographic variables determining $h_{P} (\cdot)$ . Age, gender and calendar year are always included as a part of $Z_{i}$ , but other variables could exist as well depending on different population life tables. Similarly, $X_{1}, X_{2}, \dots$ represent the covariate vectors affecting $h_{E} (\cdot)$ . Usually, $X_{i}$ contains available and relevant information regarding a patient, for example cancer stage and treatment, where $Z_{i}$ may be a subset of $X_{i}$ .

Further, denote $T_{E i}$ as the time to event due to the disease of interest. Here, the event of interest will be death. Similarly, $T_{P i}$ represents the time to event due to other causes appearing in the general population. With the additive hazard model and the assumption of $T_{E i}$ and $T_{P i}$ being conditionally independent given the covariates,³ let $T_{i} = min (T_{E i}, T_{P i})$ denote the overall time to event and $C_{i}$ the censoring time for individual $i$ . We do not observe the cause of event, only $min (T_{i}, C_{i})$ and whether the observation is censored or not. We then define the time at risk up to time $t$ after the start calendar time of the monitoring system for individual $i$ as $A_{i} (t) = min (T_{i}, C_{i}, max (t - B_{i}, 0))$ and the event indicator $δ_{i} (t) = I (T_{i} = A_{i} (t))$ .

Let $h_{0 i} (\cdot) = h_{P} (\cdot, Z_{i}) + h_{E, 0} (\cdot, X_{i})$ be the in-control hazard rate, that is the reference (or baseline) hazard based on the current situation at the start of the monitoring. This is the situation we want to be able to quickly detect deviations from. The known population hazard $h_{P} (\cdot)$ might change over calendar time, but this is suppressed in the notation. The excess hazard rate in the in-control (baseline) situation is assumed to be constant over calendar time. At some time $η$ after the monitoring start, the hazard switches to some out-of-control hazard rate $h_{1 i} (\cdot) = h_{P} (\cdot, Z_{i}) + h_{E, 1} (\cdot, X_{i})$ , which is a hazard rate that deviates from the reference period, and we would like to quickly detect this change. Finally, let the cumulative hazard be $H_{j i} (t) = \int_{0}^{t} h_{j i} (s) d s$ for $j = 0, 1$ .

2.2. The proposed CUSUM procedure

Following Gandy et al.,⁷ we define a new continuous time CUSUM procedure for the setting of relative survival based on the survival likelihood ratio. Moustakides¹⁸ has shown optimality properties for CUSUMs based on likelihood ratios. The (partial) likelihood¹⁹ based on $h_{j i}$ up to time $t$ after the start of monitoring is

L_{j} (t) = \prod_{i : B_{i} \leq t} h_{j i} (T_{i})^{δ_{i} (t)} \exp [- H_{j i} (A_{i} (t))] .

The (partial) log-likelihood ratio between the in- and out-of-control likelihood up to time

t

thus becomes

R (t) = \sum_{i : B_{i} \leq t} δ_{i} (t) \log (\frac{h_{1 i} (T_{i})}{h_{0 i} (T_{i})}) - \sum_{i : B_{i} \leq t} [H_{1 i} (A_{i} (t)) - H_{0 i} (A_{i} (t))] .

R (t)

is the continuous time analogue to a CUSUM. By restarting this sum whenever it drops to 0, we get the continuous time CUSUM control chart

Ψ (t) = R (t) - min_{s \leq t} R (s)

(1)

that signals at a time

τ = inf {t : Ψ (t) > c}

for some threshold

c > 0

. The interpretation is that when a signal occurs, evidence has been accumulated that the out-of-control model is a better model for the incoming observations than the in-control model. Thus, one should reject the in-control alternative and conclude with the out-of-control alternative. After such a signal, the monitoring will often be restarted with a new in-control model.

What differs from the procedure in Gandy et al.⁷ is that we here work with an excess hazard model and are interested in the change of $h_{E} (\cdot)$ , and this in particular implies that it is more challenging to determine $c$ and that we need to keep track of the population hazard in the first term of $R (t)$ . Usually, the population hazard can be found in national life tables. Also, in most practical applications, one has to estimate $h_{E, 0} (\cdot, X_{i})$ . A popular excess hazard model in relative survival is the proportional excess hazard model, that is $h_{E, 0} (\cdot, X_{i}) = h_{0} (\cdot) \exp (β X_{i})$ , where $h_{0}$ corresponds to the baseline excess hazard. For instance, all the available models implemented in the popular R package for relative survival relsurv²⁰ inherit this assumption, as well as models in Stata.²¹ For our purposes, a piecewise constant baseline excess hazard model, estimated by a GLM-approach with Poisson error structure,¹ and a semiparametric excess hazard model with a smooth nonparametric estimate of $h_{0} (\cdot)$ , fitted by an EM-algorithm,² will be the relevant models in the upcoming sections. Nevertheless, the method can be extended to further more complicated and flexible models of the excess hazard as long as the upcoming log-likelihood ratio can be computed, for example the framework of flexible excess hazard modelling of Eletti et al.²²

2.3. Out-of-control alternatives

In this section, we consider different alternative models for the out-of-control hazard $h_{1 i} (\cdot)$ . As motivated in the previous sections, the main focus here is to monitor the change in the excess hazard part related to the disease of interest. Therefore, the alternative considered in Gandy et al.,⁷ which implies monitoring change in the total hazard, is not suitable for the current purpose. However, based on the same idea, a proportional alternative directly on the excess hazard, that is $h_{E, 1} (\cdot, X_{i}) = ρ h_{E, 0} (\cdot, X_{i})$ for some $ρ > 0$ , can be considered. The corresponding log-likelihood ratio after a monitoring time period of length $t$ becomes

R (t) = \sum_{i : B_{i} \leq t} δ_{i} (t) \log (\frac{h_{P} (T_{i}, Z_{i}) + ρ h_{E, 0} (T_{i}, X_{i})}{h_{P} (T_{i}, Z_{i}) + h_{E, 0} (T_{i}, X_{i})}) - (ρ - 1) \sum_{i : B_{i} \leq t} [H_{E, 0 i} (A_{i} (t))] .

(2)

Here,

ρ

is the proportional factor between the out-of-control and in-control hazard. If

ρ > 1

, one can interpret this as a situation where the out-of-control hazard has increased by

(ρ - 1) \cdot 100 %

compared to the reference hazard. Similarly, if

ρ < 1

, the burden of disease has been decreased by

(1 - ρ) \cdot 100 %

with respect to the in-control hazard.

The proportional alternative implies assuming a larger change in absolute value of the hazard for patients with a higher hazard. However, this might not always be reasonable, and an alternative model for the change could be an additive model where the hazard changes by the same absolute amount for all patients. This motivates an additive out-of-control alternative for the excess hazard. Assuming that we work with non-negative excess hazard, which is usually the case when dealing with cancer patients, the additive alternative is given as $h_{E, 1} (\cdot, X_{i}) = max (0, h_{E, 0} (\cdot, X_{i}) + γ)$ for some $- \infty < γ < \infty$ so that $H_{E, 1} (t, X_{i}) = \int_{0}^{t} max (0, h_{E, 0} (u, X_{i}) + γ) d u$ and

\begin{aligned} R (t) = & \sum_{i : B_{i} \leq t} δ_{i} (t) \log (\frac{h_{P} (T_{i}, Z_{i}) + max (0, h_{E, 0} (T_{i}, X_{i}) + γ)}{h_{P} (T_{i}, Z_{i}) + h_{E, 0} (T_{i}, X_{i})}) \\ - \sum_{i : B_{i} \leq t} ((\int_{0}^{A_{i} (t)} max (0, h_{E, 0} (u, X_{i}) + γ) d u) - H_{E, 0} (A_{i} (t), X_{i})) . \end{aligned}

(3)

The additive change parameter

γ

can therefore be interpreted as the change in the excess hazard in absolute measure. If

γ < 0

, the out-of-control hazard is smaller than the in-control, which implies that the disease mortality has been reduced. On the other hand,

γ > 0

represents an increase of the excess hazard. Here, the use of the maximum function is to ensure that we do not get a negative excess hazard in the out-of-control scenario for some individuals when

γ < 0

. In the case where

γ > 0

, the maximum function can be omitted and the log-likelihood ratio simplifies to

R (t) = \sum_{i : B_{i} \leq t} δ_{i} (t) \log (\frac{h_{P} (T_{i}, Z_{i}) + h_{E, 0} (T_{i}, X_{i}) + γ}{h_{P} (T_{i}, Z_{i}) + h_{E, 0} (T_{i}, X_{i})}) - γ \sum_{i : B_{i} \leq t} A_{i} (t) .

(4)

The set-up can in theory be extended to negative excess hazard as well by rather requiring that the total hazard is non-negative, that is

h_{1 i} (\cdot) = max (h_{P} (\cdot, Z_{i}) + h_{E, 0} (\cdot, X_{i}) + γ, 0)

. For situations where

h_{1 i} (T_{i}) = 0

, the corresponding observations will have no events. The hazard term can thus be omitted for these individuals. However, we will not pursue this somewhat odd situation any further.

Further, similar to what is considered in Gandy et al.⁷ for ordinary survival models, time-transformation alternatives can be used to model changes in the hazard. This is for instance done in order to incorporate more general changes that depend on the definition of the hazard function, for example non-proportional change between the in-control and out-of-control hazard. We consider the following linear accelerated time alternative specified as $H_{E, 1} (u, X_{i}) = H_{E, 0} (k u, X_{i})$ for some $k > 0$ . This leads to the following log-likelihood ratio for the linear accelerated time alternative:

\begin{aligned} R (t) = & \sum_{i : B_{i} \leq t} δ_{i} (t) \log (\frac{h_{P} (T_{i}, Z_{i}) + k h_{E, 0} (k T_{i}, X_{i})}{h_{P} (T_{i}, Z_{i}) + h_{E, 0} (T_{i}, X_{i})}) - \sum_{i : B_{i} \leq t} [H_{E, 0 i} (k A_{i} (t)) - H_{E, 0 i} (A_{i} (t))] . \end{aligned}

(5)

This parameterisation implies in most cases an easy interpretation of the alternative - a value of

k < 1

yields

H_{E, 1 i} (u) = H_{E, 0 i} (k u) < H_{E, 0 i} (u)

, which implies that the survivor function is larger at the same time point

u

in the out-of-control setting compared to the in-control scenario, and the other way around for

k > 1

. In other words,

k > 1

accelerates time, indicating a faster increase in the cumulative excess hazard so that the out-of-control excess survival is smaller than the in-control period. In contrast,

k < 1

slows time down and can be interpreted as an improvement in the burden of disease. One could also choose the parameterisation

h_{E, 1} (u, X_{i}) = h_{E, 0} (k u, X_{i})

, but the interpretation is less clear as it will in general depend on the shape of

h_{E, 0} (u, X_{i})

as a function of

u

and thus will not be considered further here.

A possible challenge with the linear accelerated time alternative occurs if the monitoring system is used to monitor survival up to a certain time point, for instance 10-year survival, and a $k > 1$ alternative is of interest. Then, estimation of the hazard function up to $10 k$ years will be required for calculating $h_{E, 0} (k u, X_{i})$ and $H_{E, 0} (k u, X_{i})$ . If nonparametric modelling is used, this requires that survival data up to $10 k$ years must be available for estimating the hazard function. If survival data beyond 10 years are not available, a parametric hazard model can be used, but this will require extrapolations beyond the range of the data used to estimate the model.

Another type of out-of-control models that could in principle be used is models allowing for different changes depending on the value of some of the covariates. For proportional alternatives, this might for instance be an out-of-control model of the form $h_{E, 1} (\cdot, X_{i}) = ρ^{I (X_{i} \in D_{i})} h_{E, 0} (\cdot, X_{i})$ , where $D_{i}$ is a subset of the covariate space. As an example, this could be used to monitor for a change only expected to affect some individuals, for example a change in treatment that is only likely to alter the survival of patients in certain disease stages. Similar adaptions could be constructed for other types of alternatives. A further refinement could be to assume graded effects in different subsets.

A final note is that one needs to choose a value of $ρ$ , $γ$ or $k$ in order to calculate the charts. Choosing the values of these parameters can be seen as the equivalence of specifying the alternative hypothesis in a hypothesis testing setup. Thus, a specific value of the parameters corresponds to a specific alternative, that is a specific out-of-control scenario needs to be determined. In practical settings, if a clear specification of one single alternative of interest is not given, an option is to set up several charts to monitor varying degrees and directions of change. This is done as an illustration in the cancer data example in Section 4.

2.4. Criteria for determining the signal threshold

The signal threshold $c$ needs to be tuned to obtain the desired performance of the CUSUM charts. Usually, the threshold is chosen so that the charts achieve a certain performance under the in-control alternative. For example, one could relate $c$ with the in-control average run length (ARL) defined as $ARL = E (τ ∣ η = \infty)$ , which corresponds to the expected time until the charts cross the threshold $c$ when the hazard never changes to the out-of-control state. Then, $c$ is chosen such that the in-control ARL is equal to a desired value. Another common strategy is to choose $c$ such that the probability of a false signal until a given time point $t_{m}$ is at a desired level $α$ , that is $p_{hit} = P (τ \leq t_{m} ∣ η = \infty) = α$ . This is often called the in-control hitting probability.

Gandy et al.⁷ explored the setting of choosing $c$ with respect to the expected number of events until hitting or the probability that the number of events at hitting time point does not exceed a given amount when in-control, that is $E (N (τ) ∣ η = \infty)$ and $P (N (τ) \leq N_{max} ∣ η = \infty)$ . In the ordinary survival model setting studied in Gandy et al.,⁷ if a proportional change scenario is considered, $c$ can be analytically computed using a method based on a discrete time Markov chain. However, this methodology does not apply for the setting of monitoring the excess hazard as the jump of the chart is no longer constant across observations. Nevertheless, for the proportional change alternative, it is possible to approximate the current setting to a situation where the Markov chain method is applicable if the population hazard is much smaller than the excess hazard. This was examined in the master thesis of Tran,²³ where the thresholds obtained from the method yield approximately the desired performances in circumstances where the excess hazard dominates the population hazard. Otherwise, for general excess hazard scenarios studied in the present paper, it is seemingly not possible to compute $c$ analytically. In the remaining sections, the thresholds will thus be calculated via different simulation approaches depending on the setting.

2.5. Simulating the signal threshold

Motivated by the preceding section, we opt for the standard strategy of simulating the threshold so that a desired value of the false signal probability (in-control hitting probability) is achieved. In Algorithm 1, an overview of the general approach to simulate the threshold $c$ based on a predefined $p_{hit}$ is described.

We will consider two different scenarios for simulating the covariates $X_{i}$ and $Z_{i}$ . For the first simulations performed in Section 3, we assume that the true covariate distribution is known and therefore simulate $X_{i}$ and $Z_{i}$ directly from these distributions. In real-life applications, this is usually not the case. We then opt for a nonparametric bootstrap procedure that estimates the baseline/in-control distribution in such scenarios. Also, the interim censoring rate is estimated from baseline data.

2.6. Monitoring and updating schemes

The proposed method is best suited for detecting changes in the setting where all observations in the system will experience the out-of-control hazard whenever $t > η$ , for example if a new regime for post-operative treatment is introduced for all patients. However, in certain cases, only the individuals arriving after the change point $η$ will be affected by the change in the hazard, for instance if a new surgery procedure is introduced. Our method should still be able to detect the change in this latter scenario, but there will be more delay in detecting the change.

With the definition $A_{i} (t) = min (T_{i}, C_{i}, max (t - B_{i}, 0))$ , we assume that information on the status of a patient is continuously available from arrival and onward. This is realistic if cases are registered to the database more or less in real time and information about censoring or events is swiftly added. Thus, an important motivation for this setup is for real-time monitoring when information of the observations is registered immediately in the system. Immediate information availability will also be the case when the method is used for retrospective analysis, which we will illustrate in Section 4.

Although the current method works for real-time monitoring, it is often the case that observations are not registered immediately, and modifications related to the definition of the at-risk time are needed to alleviate this problem. Adjustments to the procedure for various information arrival scenarios are discussed in the following.

For instance, a scenario could be that information about each case becomes available only after censoring or an event. For this setting, Gandy et al.⁷ suggested to redefine $A_{i} (t)$ such that the observation contributes to the likelihood ratio only after an event or censoring has occurred:

A_{i} (t) = min (T_{i}, C_{i}) I (min (T_{i}, C_{i}) + B_{i} \leq t) .

This can also be applied to our set-up. It is clear that the chart will not have a continuous drift with this specification, but only jumps in value whenever an observation has experienced an event or censoring.

Another circumstance that could appear is when the information is only available at certain time points. An example is when a patient is diagnosed in the middle of a year, but the observation is only included in the database at the end of the year. A possible solution is to again redefine $A_{i} (t)$ so that the observation is included in the system periodically, for example at the start of every new year so that

A_{i} (t) = min (T_{i}, C_{i}, t - B_{i}) I (⌈ B_{i} ⌉ \leq t),

where

⌈ \cdot ⌉

is the ceiling function.

Similar to the above, we can extend the first situation where the information of an observation is only available periodically after the observation experiences an event or has been censored, for example when a patient either dies or drops out of a study in the middle of a year and the database is updated with the given observation at the start of the new year. Then, $A_{i} (t)$ can be defined as follows:

A_{i} (t) = min (T_{i}, C_{i}) I (⌈ min (T_{i}, C_{i}) + B_{i} ⌉ \leq t) .

Another scenario is when the entire database is updated in regular intervals, adding all information accumulating since the last update. Then, an option would be to run the CUSUM retrospectively by running the CUSUM for the last period once the data for that period become available. This would of course lead to a certain delay in signalling versus when running the chart in real time. Overall, the limitations are not related to the method, but rather the data registrations process. If observations cannot be included in real time, the method can be modified to accommodate this scenario, although that implies that one does not have the trait of real-time monitoring.

In some cases, it is of interest to just monitor for the excess hazard related to survival up to a certain time after diagnosis $t_{D}$ , for instance $t_{D} = 5$ years survival in cancer patients. A benefit with such monitoring is that we only need baseline data over a time period long enough to be able to estimate $h_{E, 0} (\cdot)$ up to $t_{D}$ . Furthermore, with this set-up, there is no restriction as to how long we can run the monitoring as we will only need to evaluate $h_{E, 0} (t)$ for $t \leq t_{D}$ . Thus, it is not required to have an estimate of $h_{E, 0} (t)$ for $t \geq t_{D}$ . Estimation of the in-control model is further discussed in the next subsection.

2.7. Estimating the in-control model

In practical applications, the true in-control excess hazard will be unknown and has to be estimated from baseline data. To run the monitoring procedure, we need to be able to evaluate $h_{E, 0} (t)$ and $H_{E, 0} (t)$ for any survival time $t$ that could be observed during the monitoring. In practice, this means that if nonparametric estimates of $h_{0} (t)$ and $H_{0} (t)$ are used, baseline data for at least as long period as the maximum monitoring period would be required to estimate these functions, if the monitoring is run without any upper limit on the time to event of interest. If a parametric estimate is used, one could in principle use the estimated model beyond the event horizon observed in the baseline data.

When estimating the in-control excess hazard, the estimation error will propagate to the achieved in-control performance. The impact of this, as well as the impact of model misspecifications, will be illustrated in Section 3.5.

3. Simulation study

In this section, properties of the suggested procedures are studied by simulations.

3.1. Simulation set-up

For the following, we will use a proportional excess hazard model, that is $h_{E, 0} (\cdot, X_{i}) = h_{0} (\cdot) \exp (β X_{i})$ . Inspired by data from the Norwegian Cancer Registry, Table A.1 in Appendix A presents the covariates that we will consider in the simulations, with corresponding parameter values. These parameter values are chosen based on the results obtained from fitting a piecewise constant baseline excess hazard model¹ using a given subset of the real data from the Norwegian Cancer Registry.

The monitoring is run for 10 years (in one case 5 years) and is using the population hazard of Norway from the beginning of 2010 to the end of 2019. The baseline excess hazard is chosen to be a piecewise constant function of the form $h_{0} (t) = \exp [\sum_{k} χ_{k} I_{k} (t)]$ , where $I_{k} (t) = 1$ whenever $t$ lies in the $k$ -th band of the 10-year follow-up interval. Observations that have not experienced an event 10 years after the diagnosis will be administratively censored at 10 years. The interim censoring time, that is the time until censoring due to other causes like dropout etc., is set to be exponentially distributed with the rate parameter equal to 0.000275. Here, we partition the follow-up interval into yearly bands during the first 5 years before defining one single band for $t \in [5, 10]$ . Depending on the cancer type and the ratio of events happening during the short period right after diagnosis, a finer partition at the start of the follow-up might be reasonable to capture rapid changes over time in the hazard function if this is indeed the case. However, we opt for the general recommendation of yearly bands¹ during the first few years in this illustration. We will use $χ = (- 1.4, - 1.6, - 1.8, - 2.0, - 2.1, - 3.0)$ based on the same model fitted to the subset of the real data mentioned previously.

In most simulations, we assume that both the distributions of the covariates and the parameter vectors $β$ and $χ$ are known. Thus, the charts are calculated using the true parameters given in Table A.1 for $β$ and the preceding paragraph for $χ$ , and hence the true functional form of $h_{0}$ . However, the impact of estimation error will be studied in Section 3.5.

The arrival of patients is simulated as a homogeneous Poisson process with an arrival rate $λ_{a}$ . To determine the threshold $c$ , we use the false signal probability during the 10 years of monitoring as criteria. The threshold is then computed by simulating 1000 or $10, 000$ CUSUMs from the in-control model, and we tune the threshold to achieve the desired false signal probability, which we specify to either 0.05 or 0.01. After determining the signal threshold $c$ , the signal probabilities for various out-of-control situations are simulated by simulating 1000 or $10, 000$ CUSUMs, the former being used for larger values of $λ_{a}$ to obtain the results in a reasonable time. Using 10 cores in parallel on a Macbook Pro M2 Pro, computing the thresholds for the setup considered in Section 3.2, more specifically Table 1, takes around 6 to 7 minutes. Running the $10, 000$ CUSUMs for the three out-of-control scenarios requires a similar amount of time in each case. For the setup in Table 2 where the arrival rate is larger, 12 to 13 minutes are spent to run the 1000 simulations needed to compute the thresholds, while the 1000 simulations for the out-of-control settings take 3 to 5 minutes per scenario. Next, a similar running time is observed for the situation considered in Table 3. Finally, the estimated time to run both the 1000 in-control simulations for the thresholds and each of the out-of-control scenarios is around 16 to 17 minutes.

Table 1.
The table reports proportions of cases when the CUSUM signals in a setting where individuals arrive over calendar time according to a Poisson process with rate $λ_{a} = 250$ , and the out-of-control model is $h_{1 i} (\cdot) = h_{P} (\cdot, Z_{i}) + ρ h_{E, 0} (\cdot, X_{i})$ with the latter having a piecewise constant baseline excess hazard.

Signal ratio

Scenario $η = 0$ $η = 5$ $η^{*} = 5$

$ρ$ $P (τ \leq 10 ∣ η = \infty)$

0.80 0.01 ( $c = 6.41$ ) 95.34% 78.39% 42.97%

0.05 ( $c = 5.02$ ) 98.50% 89.85% 64.17%

0.90 0.01 ( $c = 4.37$ ) 43.57% 21.49% 9.21%

0.05 ( $c = 3.41$ ) 65.28% 42.61% 24.44%

1.10 0.01 ( $c = 4.32$ ) 36.09% 17.39% 8.02%

0.05 ( $c = 3.36$ ) 59.78% 37.78% 21.40%

1.20 0.01 ( $c = 6.12$ ) 91.20% 67.78% 32.45%

0.05 ( $c = 4.80$ ) 96.75% 83.95% 54.09%

		Signal ratio
0.80	0.01 ( $c = 6.41$ )	95.34%	78.39%	42.97%
	0.05 ( $c = 5.02$ )	98.50%	89.85%	64.17%
0.90	0.01 ( $c = 4.37$ )	43.57%	21.49%	9.21%
	0.05 ( $c = 3.41$ )	65.28%	42.61%	24.44%
1.10	0.01 ( $c = 4.32$ )	36.09%	17.39%	8.02%
	0.05 ( $c = 3.36$ )	59.78%	37.78%	21.40%
1.20	0.01 ( $c = 6.12$ )	91.20%	67.78%	32.45%
	0.05 ( $c = 4.80$ )	96.75%	83.95%	54.09%

Different shifts in terms of $ρ$ and different shift points $η$ are considered. For the shift after 5 years of monitoring, both the setting that the shift affects all individuals ( $η = 5$ ) and the setting when the shift only affects new individuals ( $η^{*} = 5$ ) are considered. $10, 000$ simulations are performed for each combination of $ρ$ and scenario of shift, and for finding thresholds $c$ .

Table 2.

The table reports proportions of cases when the CUSUM signals in a setting where individuals arrive over calendar time according to a Poisson process with rate $λ_{a} = 3700$ , and the out-of-control model is $h_{1 i} (\cdot) = h_{P} (\cdot, Z_{i}) + ρ h_{E, 0} (\cdot, X_{i})$ with the latter having a piecewise constant baseline excess hazard.

		Signal ratio
	Scenario	$η = 0$	$η = 2.5$	$η^{*} = 2.5$
$ρ$	$P (τ \leq 10 ∣ η = \infty)$
0.90	0.01 ( $c = 8.00$ )	97.3%	84.9%	26.0%
	0.05 ( $c = 6.23$ )	99.2%	94.0%	46.0%
0.95	0.01 ( $c = 6.06$ )	37.2%	16.4%	3.0%
	0.05 ( $c = 4.96$ )	54.8%	33.3%	8.5%
1.05	0.01 ( $c = 6.08$ )	33.3%	13.8%	3.3%
	0.05 ( $c = 4.69$ )	57.6%	37.1%	10.1%
1.10	0.01 ( $c = 7.63$ )	96.0%	84.6%	22.2%
	0.05 ( $c = 6.06$ )	98.5%	92.5%	42.0%

The CUSUM is run for 5 years. Different shifts in terms of $ρ$ and different shift points $η$ are considered. For the shift after 2.5 years of monitoring, both the setting that the shift affects all individuals ( $η = 2.5$ ) and the setting when the shift only affects new individuals ( $η^{*} = 2.5$ ) are considered. 1000 simulations are performed for each combination of $ρ$ and scenario of shift, and for finding thresholds $c$ .

Table 3.

The table reports proportions of cases when the CUSUM signals in a setting where individuals arrive over calendar time according to a Poisson process with rate $λ_{a} = 3700$ , and the out-of-control model is $h_{1 i} (\cdot) = h_{P} (\cdot, Z_{i}) + max (0, h_{E, 0} (\cdot, X_{i}) + γ)$ with the latter having a piecewise constant baseline excess hazard.

		Signal ratio
	Scenario	$η = 0$	$η = 5$	$η^{*} = 5$
$γ$	$P (τ \leq 10 ∣ η = \infty)$
-0.002	0.01 ( $c = 4.93$ )	61.8%	49.3%	5.0%
	0.05 ( $c = 3.78$ )	79.3%	70.7%	17.6%
0.002	0.01 ( $c = 4.78$ )	60.6%	46.9%	5.6%
	0.05 ( $c = 3.83$ )	76.6%	64.7%	13.5%
0.005	0.01 ( $c = 6.81$ )	100.0%	100.0%	22.5%
	0.05 ( $c = 5.50$ )	100.0%	100.0%	40.0%

Different shifts in terms of $γ$ and different shift points $η$ are considered. For the shift after 5 years of monitoring, both the setting that the shift affects all individuals ( $η = 5$ ) and the setting when the shift only affects new individuals ( $η^{*} = 5$ ) are considered. 1000 simulations are performed for each combination of $γ$ and scenario of shift, and for finding thresholds $c$ .

We study three different scenarios for the shift from in-control to out-of-control: (i) all patients are in the out-of-control situation from the beginning ( $η = 0$ ), (ii) all patients change to the out-of-control hazard at some time point $η$ years after monitoring start (either $η = 2.5$ or $η = 5$ ), and (iii) only patients arriving after some time point $η^{*}$ have the out-of-control hazard (either $η^{*} = 2.5$ or $η^{*} = 5$ ). In all cases, the data are simulated from the correct specification of the CUSUM chart, that is the CUSUM charts are correctly specified for the incoming data. The reason for choosing two different change-point times and monitoring lengths is due to the differences in sample sizes simulated, resulting in differences in how fast a signal is obtained.

3.2. Proportional alternative

This subsection considers different scenarios with the proportional alternative represented by (2). The charts are computed using the correct value of $ρ$ that determines the various out-of-control situations. The first simulation set-up examines the scenario of four different values of $ρ$ in a situation with a moderate yearly arrival rate of $λ_{a} = 250$ . The true change-point occurs 5 years after the monitoring start.

The results from the described simulation set-up are given in Table 1. One can observe that in the situation when the excess hazard is in the out-of-control state from the start, almost all simulations yield a signal during the 10-year monitoring period for the largest shifts, that is $ρ = 0.80$ and $ρ = 1.20$ . When the change is less pronounced like $ρ = 0.90$ and $ρ = 1.10$ , the signal ratio is reduced, and of course, lower for the 0.01 false signal probability versus 0.05. If the change occurs later in the monitoring period, it is natural that the signal ratio is further decreased as there is less time for the charts to signal, hence the lower signal ratios for situations with $η = 5$ . Finally, the situation where only the individuals arriving after 5 years will experience the out-of-control state yields the smallest signal ratio. As expected, the proposed CUSUM charts have less power to detect these shifts.

Next, a situation with a much larger arrival rate, similar to what is observed in the cancer registry data used in Section 4, is considered by letting $λ_{a} = 3700$ . If one considers the same change-point time of 5 years after monitoring start and a monitoring period of 10 years, this situation will lead to very large detection ratios close to 100% in many cases. This is naturally intuitive as the power will increase with larger sample sizes. To balance out the increased power, we just run the CUSUM for 5 years and set an earlier change-point of 2.5 years after the start of monitoring and examine less prominent changes in the excess hazard by investigating values of $ρ$ that differ less from unity compared to the previous situation. From Table 2, the power of the monitoring system is relatively high if the excess hazard changes by 10% for both scenarios of $η$ . On the other hand, even with the much larger amount of arrivals, just above half of the charts will signal when the false signal probability is $0.05$ and the excess hazard changes by 5%.

3.3. Additive alternative

In this subsection, we perform a similar experiment as in the preceding for the additive alternative given by (3). By letting $λ_{a} = 3700$ , we examine three values of $γ$ : -0.002, 0.002 and 0.005. It turns out that combining the general simulation set-up described earlier with $γ = - 0.002$ yields individual out-of-control excess hazards that are non-negative so that (4) can also be used for faster calculations. For smaller shifts like $0.002$ in absolute value, we can see from Table 3 that the power of the charts is moderate when using the threshold giving a 5% probability of false signal. On the other hand, increasing the effect of change to $0.005$ leads to signals for all simulated charts for the cases where all individuals are affected by the shift. We also notice that the $γ = - 0.002$ and $γ = 0.002$ cases achieve roughly the same power. This is similar to the case with proportional alternative, if the CUSUM chart is set to monitor a large change and this is indeed the case in reality, the chart will detect and signal much faster compared to the situation with less pronounced change that is harder to identify. As before, the situation where only new arrivals are affected by the change has the smallest signal ratios for all values of $γ$ .

3.4. Linear accelerated time alternative

For the last set-up in this section, we investigate the performance of the linear accelerated time alternative (5) with four different values of $k$ . When $k > 1$ , we extend the upper limit of the last band of $h_{0} (t) = \exp [\sum_{k} χ_{k} I_{k} (t)]$ covering the time interval between 5 and 10 years to $10 k$ years. Table 4 shows that almost all charts signal for the situations when $η = 0$ and $η = 5$ across all four values of $k$ . For $η^{*} = 5$ , we can still observe that the procedure struggles more to capture these changes.

Table 4.
The table reports proportions of cases when the CUSUM signals in a setting where individuals arrive over calendar time according to a Poisson process with rate $λ_{a} = 3700$ , and the out-of-control model is $h_{1 i} (\cdot) = h_{P} (\cdot, Z_{i}) + k h_{E, 0} (k \cdot, X_{i})$ with the latter having a piecewise constant baseline excess hazard.

Signal ratio

Scenario $η = 0$ $η = 5$ $η^{*} = 5$

$k$ $P (τ \leq 10 ∣ η = \infty)$

0.90 0.01 ( $c = 7.83$ ) 100.0% 100.0% 44.7%

0.05 ( $c = 6.50$ ) 100.0% 100.0% 62.8%

0.95 0.01 ( $c = 7.29$ ) 99.9% 99.7% 10.6%

0.05 ( $c = 5.89$ ) 99.9% 100.0% 25.4%

1.05 0.01 ( $c = 7.18$ ) 99.9% 99.6% 14.3%

0.05 ( $c = 5.97$ ) 100.0% 99.9% 26.2%

1.10 0.01 ( $c = 8.49$ ) 100.0% 100.0% 35.3%

0.05 ( $c = 6.47$ ) 100.0% 100.0% 64.2%

		Signal ratio
0.90	0.01 ( $c = 7.83$ )	100.0%	100.0%	44.7%
	0.05 ( $c = 6.50$ )	100.0%	100.0%	62.8%
0.95	0.01 ( $c = 7.29$ )	99.9%	99.7%	10.6%
	0.05 ( $c = 5.89$ )	99.9%	100.0%	25.4%
1.05	0.01 ( $c = 7.18$ )	99.9%	99.6%	14.3%
	0.05 ( $c = 5.97$ )	100.0%	99.9%	26.2%
1.10	0.01 ( $c = 8.49$ )	100.0%	100.0%	35.3%
	0.05 ( $c = 6.47$ )	100.0%	100.0%	64.2%

Different shifts in terms of $k$ and different shift points $η$ are considered. For the shift after 5 years of monitoring, both the setting that the shift affects all individuals ( $η = 5$ ) and the setting when the shift only affects new individuals ( $η^{*} = 5$ ) are considered. 1000 simulations are performed for each combination of $k$ and scenario of shift, and for finding thresholds $c$ .

3.5. Estimation error and model misspecification

In the preceding simulation studies, the true covariate parameters and functional form of the baseline excess hazard are used when computing the CUSUM charts. In practice, these are all unknown quantities and need to be estimated as discussed in Section 2.7. To illustrate the effect of estimation error in achieving, for example the specified in-control false probability or average run length, we perform the simulation set-up presented in Algorithm 2.

Furthermore, we will also consider some model misspecification aspects. One would presume that the observed in-control false signal probability will be further off from the nominal value if the specified form of the baseline is not correct, for example when using a piecewise constant baseline excess hazard model when the true excess hazard is a continuous function. Consequently, we will consider two different types of baseline excess hazard in this subsection: a piecewise constant baseline and a flexible smooth nonparametric baseline hazard. Both methods can be found in the R package relsurv via the rsadd-function, see for example Perme and Pavlic.²⁰ A note regarding the smoothing procedure to obtain continuous nonparametric estimates of the baseline excess hazard is given in Appendix B.

In addition, we also investigate for each combination of estimation procedure and true excess hazard if the results obtained from bootstrapping the covariate distributions from the data set obtained in step 2 of Algorithm 2 differ from simulating directly from the true covariate distributions. We will only perform the simulation using the proportional alternative, but the idea is the same for the remaining out-of-control situations.

3.5.1. Piecewise constant baseline

In this part of the simulation, we consider the same piecewise constant baseline excess hazard set-up presented previously in this section. The arrival rate is set to $λ_{a} = 3750$ , the proportionality constant relating the out-of-control and in-control scenario is $ρ = 0.90$ . We let the desired in-control false probability during the first 5 years of monitoring be equal to $P (τ \leq 5 ∣ η = \infty) = 0.05$ . Steps 2 to 5 are performed 1000 times to obtain the distribution of in-control false probability under a given combination of estimation method and covariate simulation. For each iteration, 1000 simulations are also done in order to obtain the threshold value $c$ . Figure 1 shows the results when the model estimation is performed using the piecewise baseline excess hazard model. This resembles the situation where the estimated baseline excess hazard is correctly specified. We can observe that both the median and mean in-control signal ratio, regardless if the covariate distributions are simulated using bootstrap or the true distributions, only differ slightly from the desired nominal value of 5% due to the effect of the estimation error of $β$ and $χ$ .

Figure 1.

Histograms of the achieved signal probabilities under the in-control scenario with the true baseline excess hazard following a piecewise constant function. The model estimation is done using the model with piecewise constant baseline in step 3 of Algorithm 2. Here, the blue vertical line represents the median and the red line corresponds to the mean signal ratio. (a) True covariate distribution used in step 4 of Algorithm 2. Mean signal ratio is 5.73%, median signal ratio is 4.30% and (b) bootstrapping used to estimate the covariate distribution. Mean signal ratio is 5.65%, median signal ratio is 4.30%.

Figure 2.

Histograms of the achieved signal probabilities under the in-control scenario with the true baseline excess hazard following a piecewise constant function. The model estimation is done using a smooth nonparametric baseline excess hazard in step 3 of Algorithm 2. Here, the blue vertical line represents the median and the red line corresponds to the mean signal ratio. (a) True covariate distribution used in step 4 of Algorithm 2. Mean signal ratio is 10.23%, median signal ratio is 8.20% and (b) bootstrapping used to estimate the covariate distribution. Mean signal ratio is 10.10%, median signal ratio is 8.10%.

Figure 3.

Histograms of the achieved signal probabilities under the in-control scenario with the true baseline excess hazard following the hazard function of the Weibull distribution. The model estimation is done using the model with piecewise constant baseline in step 3 of Algorithm 2. Here, the blue vertical line represents the median and the red line corresponds to the mean signal ratio. (a) True covariate distribution used in step 4 of Algorithm 2. Mean signal ratio is 4.14%, median signal ratio is 3.00% and (b) bootstrapping used to estimate the covariate distribution. Mean signal ratio is 4.09%, median signal ratio is 2.90%.

Figure 4.

Histograms of the achieved signal probabilities under the in-control scenario with the true baseline excess hazard following the hazard function of the Weibull distribution. The model estimation is done using a smooth nonparametric baseline excess hazard in step 3 of Algorithm 2. Here, the blue vertical line represents the median and the red line corresponds to the mean signal ratio. (a) True covariate distribution used in step 4 of Algorithm 2. Mean signal ratio is 8.55%, median signal ratio is 6.80% and (b) bootstrapping used to estimate the covariate distribution. Mean signal ratio is 8.24%, median signal ratio is 6.50%.

On the other hand, when examining the results of the same procedure using the smooth baseline excess hazard estimate from Figure 2, we see that the median and mean in-control signal ratio, independent of the scenario considered in step 4 of Algorithm 2, are further away from the nominal value. Not only do we have the estimation error in the covariate parameters that affects the results, the incorrect specification of the form of the baseline excess hazard when using the smooth baseline estimate leads to larger deviation from the desired false signal ratio. Therefore, if the true baseline excess hazard is indeed a piecewise constant function, approximating it with a continuous and smooth function might cause undesired behaviour in the monitoring procedure.

3.5.2. Weibull baseline

In theory, the previous true form of baseline excess hazard favours the model with piecewise constant baseline as it is difficult for the estimated smoothed nonparametric baseline excess hazard to capture the stepwise behaviour. We now illustrate the opposite scenario by considering a Weibull baseline as the true form, where

h_{0} (t) = a b t^{a - 1} .

Thus, in these simulations, we simulate data from a continuous baseline excess hazard rather than a piecewise constant. Inspired by the data from the Norwegian Cancer Registry, we let

a = 0.65

and

b = 0.25

. Otherwise, the other quantities remain the same as in the piecewise constant baseline scenario. Results are reported in Figure 3. Compared to the situation with a piecewise constant baseline as the true form, the differences between the 5% nominal value and the mean and median signal ratio are now increased as expected due to the misspecification of the baseline excess hazard in the model estimation.

Figure 4 reports the results obtained when using the smooth nonparametric baseline excess hazard estimate. The in-control signal ratio obtained is slightly closer to the desired value of 5%, but the difference in absolute value is still larger compared to the piecewise constant model from above. We suspect that this is a consequence of the choice of the smoothing parameters in the experiment. This example shows a potential issue with the nonparametric baseline excess hazard estimate² - the smoothing parameters will have a large impact on how the charts will perform compared to the required in-control performance. An ideal situation is to obtain the full continuous estimate of $h_{0}$ at the last E-step of the EM-algorithm when fitting the Cox-type excess hazard model to avoid another smoothing operation on top of the already smoothed estimates (see for example Appendix B). However, this is not possible, and one could argue that this combination might give rise to oversmoothing. Nonetheless, this simulation set-up illustrates the challenges of a nonparametric baseline in the model estimation when using the proposed monitoring scheme.

4. Application to colorectal cancer data monitoring

In this section, we apply the proposed methods to a data set of patients diagnosed with colorectal cancer in Norway between the start of 1970 and the end of 2020, followed up until the end of 2020. We will focus on patients with specified and known position of the tumour according to the International Classification of Diseases (ICD), who received standard surgical resections for cure of stage I–III colorectal cancer. With reference to Table A.1, this means that the relevant subset is the patients in surgery group 0, with SEER stage equal to localised or regional and with ICD indicator equal to 0, 1, or 2. Furthermore, while colorectal cancer is mostly a disease of the elderly, there has been a significant rise in incidence of this disease among the younger population, defined as early-onset colorectal cancer, that is $\leq 50$ years of age. Currently, there is an increased focus on this subgroup of colorectal cancer patients (see for example REACCT Collaborative et al.²⁴ for an overview). Inspired by this fact, we divide the data into patients up to 50 years and patients above 50 years of age. This leaves us with $5321$ patients in the younger age group and $86 360$ patients in the older age group.

In this example, we explore monitoring in bands of 10 years, that is we choose a baseline period of 10 years and use the model fitted for this period to monitor for change during the next 10 years. This may provide important insight into the effects of new treatments and policy changes on outcomes for large patient populations on a macro health care level. More specifically, we define the first baseline period to be between 1970 and 1980. The patient cohort diagnosed in this period is used to fit an excess hazard model that will represent $h_{E, 0} (\cdot, X_{i})$ in the CUSUM chart. The covariate vector $X_{i}$ here contains gender, ICD indicator, morphology type and SEER stage. Because we are interested in examining each age group separately due to the interest in early-onset cancer, age at diagnosis is not directly included as covariate as this is partly taken care of by stratifying for the two age groups.

In addition, only a piecewise constant baseline excess hazard model is considered here. The partition of the follow-up interval is in line with the simulation setups in Section 3, that is yearly bands for the first 5 years and then a last band from 5 to 10 years of follow-up. With these results as the in-control period, we run the CUSUM charts for the next 10 years, that is the time period between 1980 and 1990. Subsequently, this period is then used as the baseline period for the new monitoring period of 1990 to 2000, and this is done until we monitor the period between 2010 and 2020.

For each period, we will for illustration purposes run several charts corresponding to monitoring for different degrees of improvement or worsening from the in-control baseline excess hazard. In terms of the interpretation of the method, the CUSUM charts for the period between 1980 and 1990 will signal if there is accumulation of evidence in the data of improvement or worsening, depending on the specified alternative, with respect to the baseline time period between 1970 and 1980. Similarly, the charts for the period between 1990 and 2000 will signal if there is accumulated evidence of a change in the excess hazard from the baseline period of 1980 and 1990, and so on.

4.1. Proportional alternative

First, we consider the proportional alternative. For each combination of age group and monitoring period, the chart is computed using four different values of $ρ$ : 0.80, 0.90, 0.95, and 1.05, that is three scenarios corresponding to decreased and one to increased excess hazard. Figure 5 shows the calculated charts for different combinations of monitoring time period and $ρ$ in the two age groups. The dashed lines are the thresholds obtained by simulations using the bootstrap procedure mentioned in Section 3.5. The plots of the charts have different scales for distinct values of $ρ$ as the corresponding CUSUMs are not numerically comparable.

Figure 5.

Standardised CUSUM charts using the proportional alternative for each age group monitoring 10-year survival in four different time periods with the estimated output from the piecewise constant baseline excess hazard model fitted on the preceding 10-year time period as the in-control. Four different values of $ρ$ are explored here. $Ψ (t)$ corresponds to the value of the CUSUM chart at a time point $t$ given by (1) with $R (t)$ of the form (2). The charts are then standardised by dividing with the associated thresholds to get $Ψ (t) / c$ and a common threshold equal to 1. (a) Age $> 50$ and (b) Age $\leq 50$ .

It is clear from both plots that the charts fluctuate and do not signal for any combinations of age group and time period when considering the increased burden of disease scenario of $ρ = 1.05$ . In contrast, all charts signal for all the $ρ < 1$ scenarios, indicating that improvement versus the previous decade is observed for all the periods considered. For the elders, a striking observation is that in the 1990s, no improvement was seen during the first 5 years.

An explanation for this could be that no advances in treatment compared to the 1980s were observed during the first half of the 1990s. On the other hand, the chart signals fastest in the period 2010 to 2020. This could indicate early advances in this period, but the faster signal could also possibly be partly explained by more patients arriving during this time interval. For the younger, there seems to be an indication of a better treatment advance in the period 2000 to 2010 versus the previous decade than in the other periods. The generally later signals for the younger could be explained by much fewer patients in this group.

4.2. Additive alternative

Next, we investigate the setting of an additive alternative. As in the previous setting, we also here use three shift parameters that correspond to decreasing and one corresponding to increasing excess hazard, that is the values of $γ$ considered are $- 0.020$ , $- 0.010$ , $- 0.005$ and 0.005. The results are depicted in Figure 6. We observe that none of the charts yield a signal for the $γ = 0.005$ case that implies a worse outcome over time in this additive setting, although for the younger patients, the charts are actually close to a signal in the periods 1980 to 1990 and 1990 to 2000.

Figure 6.

Standardised CUSUM charts using the additive alternative for each age group monitoring 10-year survival in four different time periods with the estimated output from the piecewise constant baseline excess hazard model fitted on the preceding 10-year time period as the in-control. $Ψ (t)$ corresponds to the value of the CUSUM chart at a time point $t$ given by (1) with $R (t)$ of the form (3). Four different values of $γ$ are explored here. The charts are then standardised by dividing with the associated thresholds to get $Ψ (t) / c$ and a common threshold equal to 1. (a) Age $> 50$ and (b) Age $\leq 50$ .

Figure 7.

Standardised CUSUM charts using the linear accelerated time alternative for each age group monitoring 10-year survival in four different time periods with the estimated output from the piecewise constant baseline excess hazard model fitted on the preceding 10-year time period as the in-control. $Ψ (t)$ corresponds to the value of the CUSUM chart at a time point $t$ given by (1) with $R (t)$ of the form (5). Four different values of $k$ are explored here. The charts are then standardised by dividing with the associated thresholds to get $Ψ (t) / c$ and a common threshold equal to 1. (a) Age $> 50$ and (b) Age $\leq 50$ .

For the remaining values of $γ$ , all the charts related to the old age group signal with an increasing trend. The behaviours of the charts across all the values of $γ$ are very similar to the proportional alternative, with the latest two monitoring periods appearing to have the largest improvement.

Looking at the charts obtained from the younger age group, no signals are now obtained in the first two monitoring periods for $γ < 1$ unlike the conclusions from the proportional alternative. The values of the charts for the periods 2000 to 2010 and 2010 to 2020 when $γ = - 0.010$ and $γ = - 0.005$ are also much closer. In addition, for the largest shift $γ = - 0.020$ , the chart for the period 2010 to 2020 starts decreasing again after signalling at around 5 years. Therefore, the evidence of a larger reduction of $- 0.020$ in the excess hazard is much weaker for this period compared to the evidence for smaller changes like $- 0.010$ and $- 0.005$ . This shows the effect of setting a value of the change parameter to a larger change than the observed, and that monitoring for different values of the alternative might give useful insight.

4.3. Linear accelerated time alternative

Finally, we try a linear accelerated time alternative with the following four values of the acceleration parameter $k$ : 0.80, 0.90, 0.95 and 1.05. Therefore, this corresponds again to three improvement alternatives and one worsening situation. From Figure 7, for the older age group, we see a similar behaviour as in the previous models, although with a later signal in the case with the smallest improvement of $k = 0.95$ .

The most interesting observation from the results presented in Figure 7 regarding the younger age group is the charts corresponding to the periods 1990 to 2000 and 2010 to 2020. With the additive alternative, the charts monitoring towards improvement for 1990 to 2000 did not signal at all. However, with the linear accelerated time alternative, all charts for this period signal for all values of $k < 1$ , as was also the case for the corresponding charts with the proportional model. Further, the charts for 2010 to 2020 have a much later response than what was observed for the additive and proportional alternative. Also, when $k = 0.95$ , the chart does not even signal. This is another example of how different alternatives can potentially lead to different conclusions.

5. Summary

We have presented a CUSUM-based method for monitoring changes in excess hazard for relative survival settings where the population hazard is known. This complements the literature on monitoring based on time-to-event models. The proposed method can be used for real-time monitoring, semi real-time monitoring or for retrospective analyses. The method can be adapted to various data updating schemes.

We have considered proportional, additive and linear accelerated time change models and studied properties by simulations and in an application to cancer registry data. In particular, we have also considered the impact of estimation error and some forms of model misspecifications. Simulations indicate that model misspecifications might be a somewhat severe issue, while with a decent amount of baseline data the impact of estimation error is not so critical. To our knowledge, there is no formal testing procedure available to test the functional form of the baseline excess hazard. Fortunately, the deviations observed in the simulations are smallest in the situation that will most likely appear in practice where the true baseline hazard function is continuous and the fitted model is piecewise constant.

Furthermore, we have mainly focused on monitoring using proportional excess hazard models. For some cancer types, it is known that the excess hazard is not proportional. In contrast to the scenario related to misspecification of the baseline hazard, checking and performing formal hypothesis test for violation of this assumption can be done using Schoenfeld-like residuals proposed by Stare et al.²⁵ or martingale residuals proposed by Danieli et al.²⁶ The latter residuals can also be used to check the validity of the functional form of the covariates in the model.²⁶

For small sample sizes, the bootstrap approach for handling estimation error presented in Gandy and Kvaløy²⁷ could in principle be adapted, although with a substantial computational burden in the current setting. Applications to data illustrate that careful considerations need to be made when specifying the type of changes in the excess hazard to monitor against.

For the scenarios considered in this paper, year of diagnosis or variants of time with respect to the start of monitoring is not included as a predictor. However, one could fit a model where year of diagnosis is included as predictor to model a general secular trend, assume that this trend continues into the monitoring period, and use the proposed method to monitor drift from the secular trend.

Another extension that could be considered is to try to extend the approach to cure models, for instance by trying to adapt the CUSUM for cure models presented by Oliveira et al.¹¹ to the excess hazard setting. Here, this would correspond to some individuals having zero excess hazard. A different variant would be cure models where all individuals have zero excess hazard after some cure time point, this should work automatically with the current set-up if a proportional alternative is used.

Finally, it could be of interest to avoid the requirement of specifying a specific alternative in order to capture more general out-of-control situations. Phinikettos and Gandy²⁸ have proposed a method that avoids this necessity for the setting of non-risk-adjusted total hazard monitoring. Extending this methodology to the relative survival situation of interest here, but also allowing for the incorporation of covariates via regression models, could be an interesting direction in the future.

Footnotes

Acknowledgments

The authors would like to thank the three anonymous referees for detailed and helpful feedback that led to an improved manuscript. The authors would also like to express gratitude to Professor Axel Gandy for helpful comments regarding the manuscript and interesting suggestions that have been partly mentioned in the discussion.

ORCID iDs

Jimmy Huy Tran

Jan Terje Kvaløy

Hartwig Kørner

Ethical approval and informed consent statements

Ethical approval was provided by the Regional Committee for Medical and Health Research Ethics, Western Norway (Ref. no. 343976). The data set was anonymised according to the guidelines from the Norwegian Data Protection Authority. Exemption from confidentiality was provided by the National Institute of Public Health.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Jimmy Huy Tran acknowledges finanacial support by Finansmarkedsfondet (Grant Number: 337601).

Declaration of conflicting interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Disclaimer

The study has used data from the Cancer Registry of Norway. The interpretation and reporting of these data are the sole responsibility of the authors, and no endorsement by the Cancer Registry of Norway is intended nor should be inferred.

A. Summary table of the covariates in the simulation set-up of Section 3

Table 5.

Overview of the covariates used in the simulations, with the corresponding parameter value used and proportion of observations simulated to be in each category for the categorical covariates. Age is simulated from a normal distribution with mean 75 and standard deviation 10, truncated at 50 and 105. See Tran²³ for details about the medical interpretations of the covariates.

Variable	$β$	Proportion
Gender $=$ Male	$\cdot$	50.00%
Gender $=$ Female	$0.005$	50.00%
ICD indicator $=$ 0 (Rectum)	$\cdot$	27.00%
ICD indicator $=$ 1 (Right colon or appendix)	$0.500$	44.00%
ICD indicator $=$ 2 (Left colon)	$0.200$	28.00%
ICD indicator $=$ 3 (Unspecified location or colon polyps)	$0.300$	1.00%
Morphology type $=$ Adenocarcinoma	$\cdot$	90.00%
Morphology type $=$ Mucinous carcinoma	$- 0.050$	10.00%
SEER stage $=$ Distant	$\cdot$	20.00%
SEER stage $=$ Localised	$- 3.000$	20.00%
SEER stage $=$ Regional	$- 1.750$	55.00%
SEER stage $=$ Unknown	$- 1.000$	5.00%
Surgery group $=$ 0 (Major surgery)	$\cdot$	82.75%
Surgery group $=$ 1 (Minor surgery)	$1.500$	17.15%
Surgery group $=$ 2 (No surgery)	$2.500$	0.10%

B. Smoothing procedure for the semiparametric excess hazard model in the CUSUM method

An important note regarding the semiparametric excess hazard model implemented in the R package relsurv²⁰ is that the user can specify a smoothing parameter bwin. Originally, the estimation of the model is done using the EM-algorithm by Dempster et al.²⁹ Perme et al.² proposed to apply kernel smoothing to the nonparametric baseline excess hazard estimates at each E-step of the algorithm. The argument bwin therefore controls the bandwidth in the smoothing process, with larger values implying larger bandwidth, and therefore smoother estimates. In the end, only the estimates of $h_{0}$ at the observed times to event are returned. However, in order to apply the proposed CUSUM procedures, $h_{0} (t)$ is required at all possible values of $t \in [0, 10]$ for continuous monitoring if the main interest is the 10-year excess survival. To get a continuous estimate of $h_{0} (t)$ , we must again apply some sort of manual smoothing on top of the returned estimates at the times to event observed in the data. For this purpose, smoothing splines (see for example James et al.³⁰) are used with the baseline excess hazard estimates as the response variable and the times to event as the sole predictor. Therefore, there are two extra sources that can impact the upcoming results: the smoothing parameter bwin at each E-step of the EM-algorithm when fitting the model and the effective degrees of freedom from the smoothing splines procedure. To get somewhat sensible results, we have decided to set bwin and the effective degrees of freedom to 25. It is therefore important to note that a different combination of these two parameters will lead to a different result than what we have obtained in Section 3.5.

References

Dickman

Sloggett

Hills

, et al. Regression models for relative survival. Stat Med 2004; 23: 51–64.

Perme

Henderson

Stare

. An approach to estimation in relative survival regression. Biostatistics 2009; 10: 136–146.

Perme

Stare

Estève

. On estimation in relative survival. Biometrics 2012; 68: 113–120.

Page

. Continuous inspection schemes. Biometrika 1954; 41: 100–115.

Biswas

Kalbfleisch

. A risk-adjusted CUSUM in continuous time based on the Cox model. Stat Med 2008; 27: 3382–3406.

Sego

Reynolds Jr

Woodall

. Risk-adjusted monitoring of survival times. Stat Med 2009; 28: 1386–1401.

Gandy

Kvaløy

Bottle

, et al. Risk-adjusted monitoring of time to event. Biometrika 2010; 97: 375–388.

Steiner

Jones

. Risk-adjusted survival time monitoring with an updating exponentially weighted moving average (EWMA) control chart. Stat Med 2010; 29: 444–454.

Begun

Kulinskaya

MacGregor

. Risk-adjusted CUSUM control charts for shared frailty survival models with application to hip replacement outcomes: a study using the NJR dataset. BMC Med Res Methodol 2019; 19: 217

10.

Keshavarz

Asadzadeh

Niaki

STA

. Risk-adjusted frailty-based CUSUM control chart for phase I monitoring of patients’ lifetime. J Stat Comput Simul 2021; 91: 334–352.

11.

Oliveira

Valença

Medeiros

, et al. Risk-adjusted monitoring of time to event in the presence of long-term survivors. Biomet J 2016; 58: 1485–1505.

12.

Liu

Lai

Wang

, et al. SCR-CUSUM: an illness-death semi-Markov model-based risk-adjusted CUSUM for semi-competing risk data monitoring. Comput Indust Eng 2023; 184: 109530.

13.

Kuang

Das

Sir

, et al. Likelihood ratio-based CUSUM charts for real-time monitoring the quality of service in a network of queues. IISE Trans Health Syst Eng 2023; 13: 344–354.

14.

Zhang

, et al. The effect of estimation error on risk-adjusted survival time CUSUM chart performance. Quality Reliab Eng Int 2016; 32: 1445–1452.

15.

Massarweh

Chen

Rosen

, et al. Comparative effectiveness of risk-adjusted cumulative sum and periodic evaluation for monitoring hospital perioperative mortality. Med Care 2021; 59: 639–645.

16.

Lai

Liu

, et al. A risk-adjusted approach to monitoring surgery for survival outcomes based on a weighted score test. Comput Indust Eng 2021; 160: 107568.

17.

Chen

Chidi

Dong

, et al. Risk-adjusted cumulative sum for early detection of hospitals with excess perioperative mortality. JAMA Surg 2023; 158: 1176–1183.

18.

Moustakides

. Optimal stopping times for detecting changes in distributions. Ann Stat 1986; 14: 1379–1387.

19.

Aalen

Borgan

Gjessing

. Survival and event history analysis: a process point of view. New York: Springer, 2008.

20.

Perme

Pavlic

. Nonparametric relative survival analysis with the R package relsurv. J Stat Softw 2018; 87: 1–27.

21.

Dickman

Coviello

. Estimating and modeling relative survival. Stata J 2015; 15: 186–215.

22.

Eletti

Marra

Quaresma

, et al. A unifying framework for flexible excess hazard modelling with applications in cancer epidemiology. J R Stat Soc Ser C: Appl Stat 2022; 71: 1044–1062.

23.

Tran

. Relative survival methods. Theory, applications and extensions to monitoring. Master’s Thesis, Department of Mathematics and Physics, University of Stavanger, 2022.

24.

Collaborative

Zaborowski

Abdile

, et al. Characteristics of early-onset vs late-onset colorectal cancer: a review. JAMA Surg 2021; 156: 865–874.

25.

Stare

Pohar

Henderson

. Goodness of fit of relative survival models. Stat Med 2005; 24: 3911–3925.

26.

Danieli

Bossard

Roche

, et al. Performance of two formal tests based on martingales residuals to check the proportional hazard assumption and the functional form of the prognostic factors in flexible parametric excess hazard models. Biostatistics 2017; 18: 505–520.

27.

Gandy

Kvaløy

. Guaranteed conditional performance of control charts via bootstrap methods. Scand J Stat 2013; 40: 647–668.

28.

Phinikettos

Gandy

. An omnibus CUSUM chart for monitoring time to event data. Lifetime Data Anal 2014; 20: 481–494.

29.

Dempster

Laird

Rubin

. Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc: Ser B (Methodol) 1977; 39: 1–22.

30.

James

Witten

Hastie

, et al. An introduction to statistical learning. Vol. 112. New York, NY: Springer, 2013.