cgr_cusum function

Continuous time Generalized Rapid response CUSUM (CGR-CUSUM)

This function performs the CGR-CUSUM procedure described in Gomon et al. (2022). For detection purposes, it suffices to determine the value of the chart at the times of failure. This can be achieved by leaving ctimes unspecified. The function has two vital parameters, at least one of which must be specified:

• coxphmod:: Cox proportional hazards model to be used for risk-adjustment. If cbaseh is not specified, it will be determined from coxphmod numerically.
• cbaseh:: The cumulative baseline hazard rate to use for chart construction. If specified with coxphmod missing, no risk-adjustment will be performed

cgr_cusum(data, coxphmod, cbaseh, ctimes, h, stoptime, C, pb = FALSE,
  ncores = 1, cmethod = "memory", dependencies, detection = "upper",
  assist, maxtheta = log(6))

Arguments

• data: A data.frame with rows representing subjects and the following named columns:

  • entrytime:: time of entry into study (numeric);
  • survtime:: time from entry until event (numeric);
  • censorid:: censoring indicator (0 = right censored, 1 = observed), (integer).

  and optionally additional covariates used for risk-adjustment.

• coxphmod: A Cox proportional hazards regression model as produced by the function coxph(). Suggested:

  coxph(Surv(survtime, censorid) ~ covariates, data = data).

  Alternatively, a list with the following elements:

  • formula:: a formula() in the form ~ covariates;
  • coefficients:: a named vector specifying risk adjustment coefficients for covariates. Names must be the same as in formula and colnames of data.

• cbaseh: A function that returns the unadjusted cumulative baseline hazard $H_0(t)$. If cbaseh is missing but coxphmod has been specified as a survival object, this baseline hazard rate will be determined using the provided coxphmod.

• ctimes: (optional): Vector of construction times at which the value of the chart should be determined. When not specified, the chart is constructed at all failure times.

• h: (optional): Value of the control limit. The chart will only be constructed until the value of the control limit has been reached or surpassed.

• stoptime: (optional): Time after which the value of the chart should no longer be determined. Default = max(failure time). Useful when ctimes has not been specified.

• C: (optional): A numeric value indicating how long after entering the study patients should no longer influence the value of the chart. This is equivalent to right-censoring every observation at time entrytime + C.

• pb: (optional): A boolean indicating whether a progress bar should be shown. Default is FALSE.

• ncores: number of cores to use to parallelize the computation of the CGR-CUSUM chart. If ncores = 1 (default), no parallelization is done. You can use detectCores() to check how many cores are available on your computer.

• cmethod: Method to calculate chart values. One of the following:

  • "memory" (default): matrix formulation of the problem (faster for high volume/long time construction)
  • "CPU" calculates the value of the CGR-CUSUM for every time point from scratch. Recommended for small data volume (lower initialization time).

• dependencies: (optional): When ncores > 1, specify a list of variables/functions/other dependencies to be exported to the core clusters for parallel computation.

• detection: Should an "upper" or "lower" CGR-CUSUM be constructed. Upper CUSUMs can be used to monitor for an increase in the failure rate, while lower CUSUMs can be used to monitor for a decrease in the failure rate.

• assist: (optional): Output of the function parameter_assist()

• maxtheta: (optional): Maximum value the maximum likelihood estimate for parameter $\theta$ can take. If detection = "lower", -abs(theta)

  will be the minimum value the maximum likelihood estimate for parameter $\theta$ can take. Default is log(6), meaning that at most a 6 times increase/decrease in the odds/hazard ratio is expected.

Returns

An object of class "cgrcusum" containing:

• CGR: a data.frame with named columns:

  • time:: times at which chart is constructed;
  • value:: value of the chart at corresponding times;
  • exp_theta_t:: value of MLE $e^(\theta_t)$;
  • S_nu: time from which patients are considered for constructing the chart.

• call: the call used to obtain output;

• stopind: indicator for whether the chart was stopped by the control limit;

• h: Specified value for the control limit.

Details

The CGR-CUSUM can be used to test for a change of unknown positive fixed size $\theta$

in the subject-specific hazard rate from $h_i(t)$ to $h_i(t) exp(\theta)$

starting from some unknown subject $\nu$. The starting time of the first subject who had an increase in failure rate as well as the estimated increase in the hazard rate are shown in the output. The CGR-CUSUM is determined as

$$\max_{1 \leq \nu \leq n} \left( \hat{\theta}_{\geq \nu}(t) N_{\geq \nu}(t) - \left( \exp\left( \hat{\theta}_{\geq \nu}(t) \right) - 1 \right) \Lambda_{\geq \nu}(t)\right),max{1<=\nu<=n} (\theta_{>=\nu}(t)N_{>=\nu}(t)) - (exp(\theta_{>=\nu}(t))-1) \Lambda_{>=\nu}(t)),$$

where

$$N(\geq \nu)(t) = \sum_{i \geq \nu} N_i(t),N_{>=\nu}(t) = \sum_{i>=\nu} N_i(t),$$

with $N_i(t)$ the counting process for the failure at time $t$ of subject $i$

and

$$\Lambda_{\geq \nu}(t) = \sum_{i \geq \nu} \Lambda_i(t),\Lambda_{>=\nu}(t) = \sum_{i>=\nu}\Lambda_i(t),$$

where $\Lambda_i(t)$ is the cumulative intensity of subject $i$ at time $t$.

When maxtheta is specified, the maximum likelihood estimate of $\theta$

is restricted to either abs(maxtheta) (upper sided CGR-CUSUM) or -abs(maxtheta) (lower sided CGR-CUSUM). Choosing this value smaller leads to smaller control limits and therefore quicker detection times, but can cause delays in detection if the true increase in failure rate is larger than the cut-off. The default of expecting at most a 6 times increase in hazard/odds ratio can therefore be the wrong choice for your intended application area. If unsure, the most conservative choice is to take maxtheta = Inf.

Examples

require(survival)
#Select only the data of the first year of the first hospital in the surgerydat data set
tdat <- subset(surgerydat, unit == 1 & entrytime < 365)

#We know that the cumulative baseline hazard in the data set is
#Exponential(0.01). If you don't know the cumulative baseline, we suggest
#leaving the cbaseh argument empty and determining a coxphmod (see help)
tcbaseh <- function(t) chaz_exp(t, lambda = 0.01)

#Determine a risk-adjustment model using a Cox proportional hazards model.
#Outcome (survival) is regressed on the available covariates:
exprfit <- as.formula("Surv(survtime, censorid) ~ age + sex + BMI")
tcoxmod <- coxph(exprfit, data= surgerydat)

#Determine the values of the chart
cgr <- cgr_cusum(data = tdat, coxphmod = tcoxmod, cbaseh = tcbaseh, pb = TRUE)
#plot the CGR-CUSUM (exact hazard)
plot(cgr)

#Alternatively, cbaseh can be left empty when specifying coxphmod through coxph()
cgr_cox <- cgr_cusum(data = tdat, coxphmod = tcoxmod, pb = TRUE)
#plot the CGR-CUSUM (estimated hazard from coxph)
plot(cgr_cox)

References

Gomon, D. and Putter, H. and Nelissen, R. G. H. H. and van der Pas, S. (2022), CGR-CUSUM: A Continuous time Generalized Rapid Response Cumulative Sum chart, Biostatistics tools:::Rd_expr_doi("10.1093/biostatistics/kxac041")

Gomon D. and Fiocco M and Putter H and Signorelli M, SUrvival Control Chart EStimation Software in R: the success Package, The R Journal, vol. 15, no. 4, pp. 270–291, Dec. 2023, tools:::Rd_expr_doi("10.32614/rj-2023-095")

See Also

plot.cgrcusum, runlength.cgrcusum

Other quality control charts: bk_cusum(), funnel_plot()

Author(s)

Daniel Gomon