Title:	Survival Control Charts Estimation Software
Version:	1.1.1
Maintainer:	Daniel Gomon <dgstatsoft@gmail.com>
Description:	Quality control charts for survival outcomes. Allows users to construct the Continuous Time Generalized Rapid Response CUSUM (CGR-CUSUM) <doi:10.1093/biostatistics/kxac041>, the Biswas & Kalbfleisch (2008) <doi:10.1002/sim.3216> CUSUM, the Bernoulli CUSUM and the risk-adjusted funnel plot for survival data <doi:10.1002/sim.1970>. These procedures can be used to monitor survival processes for a change in the failure rate.
License:	GPL (≥ 3)
URL:	https://github.com/d-gomon/success
Depends:	ggplot2, pbapply, plotly, Rfast, R (≥ 3.5.0)
Imports:	survival, ggrepel, Matrix, matrixcalc
Suggests:	gridExtra, knitr, rmarkdown, testthat (≥ 3.1.10), covr
Config/testthat/edition:	3
Encoding:	UTF-8
LazyData:	true
VignetteBuilder:	knitr
RoxygenNote:	7.2.3
NeedsCompilation:	no
Packaged:	2025-05-15 09:59:24 UTC; danie
Author:	Daniel Gomon [aut, cre], Mirko Signorelli [ctb]
Repository:	CRAN
Date/Publication:	2025-05-15 11:40:02 UTC

Estimate arrival rate of a Poisson points process

Description

In a Poisson point process, subjects arrive with exponentially distributed inter-arrival times with rate \psi. This function can be used to estimate the parameter \psi.

Usage

arrival_rate(data)

Arguments

Value

A (named) vector containing the estimated arrival rate in the data, or for each unit in the data.

Author(s)

Daniel Gomon

Examples

arrival_rate(surgerydat)

Average Run Length for Bernoulli CUSUM

Description

This function allows to estimate the Average Run Length (ARL) of the risk-adjusted Bernoulli CUSUM (see bernoulli_cusum()) through a Markov Chain Approach (Brook & Evans(1972) & Steiner et al. (2000)) or exploiting the relationship with the Sequential Probability Ratio Test (Kemp (1971)). The function requires the specification of one of the following combinations of parameters as arguments to the function:

• glmmod & theta

• p0 & theta

• p0 & p1

Average run length of lower-sided Bernoulli CUSUM charts can be determined by specifying theta < 0.

Usage

bernoulli_ARL(h, n_grid, glmmod, theta, theta_true, p0, p1, method = c("MC",
  "SPRT"), smooth_prob = FALSE)

Arguments

Details

The average run length of a CUSUM chart S_n is given by E[\tau_n], where \tau_n is defined as:

\tau_n = \inf\{n \geq 0: S_n \geq h\}.

When method = "MC", the average run length will be determined by the Markov Chain approach described in Brook & Evans (1972), using the risk-adjustment correction proposed in Steiner et al. (2000). The idea is to discretize the domain (0, h) into n_{grid} -1 state spaces, with E_0 of width w/2 and E_1, \ldots, E_{n_{grid}-1} of width w, such that E_{n_{grid}} is an absorbing state. This is done using the following steps:

• w is determined using the relationship \frac{2h}{2t-1}.

• Transition probabilities between the states are determined and 'transition matrix' R is constructed.

• The equation (\bm{I}-\bm{R}) \bm{ARL} = \bm{1} is solved to find the ARL starting from each of the states.

When method = "SPRT", the average run length will be determined by the relationship between the SPRT and CUSUM described in Kemp (1971), using the risk-adjustment correction proposed in Steiner et al. (2000). If N is the run length of a SPRT, P(0) the probability of a SPRT terminating on the lower boundary of zero and R the run length of a CUSUM, then:

E[R] = \frac{E[N]}{1 - P(0)}.

E[N] and P(0) are completely determined by

G_n(z) = \int_0^h F(z-w) dG_{n-1}(w)

with F(x) the cdf of the singletons W_n. The integral can be approximated using the generalized trapezoidal quadrature rule:

G_n(z) = \sum_{i=0}^{n_{grid}-1} \frac{F(z-x_{i+1}) + F(z-x_{i})}{2} \left(G_{n-1}(x_{i+1}) - G_{n-1}(x_{i}) \right)

Value

A list containing:

• ARL_0: A numeric value indicating the average run length in number of outcomes when starting from state E_0.

• ARL: A data.frame containing the average run length (#outcomes) depending on the state in which the process starts (E_0, E_1, \ldots, E_{n_{grid}-1})

  start_val:

  Starting value of the CUSUM, corresponding to the discretized state spaces E_{i};

  #outcomes:

  ARL for the CUSUM with initial value start_val;

• R: A transition probability matrix containing the transition probabilities between states E_0, \ldots, E_{t-1}. R_{i,j} is the transition probability from state i to state j.

• h: Value of the control limit.

The value of ARL_0 will be printed to the console.

Author(s)

Daniel Gomon

References

Brook, D., & Evans, D. A. (1972). An Approach to the Probability Distribution of CUSUM Run Length. Biometrika, 59(3), 539-549. doi:10.2307/2334805

Steiner, S. H., Cook, R. J., Farewell, V. T., & Treasure, T. (2000). Monitoring surgical performance using risk-adjusted cumulative sum charts. Biostatistics, 1(4), 441-452. doi:10.1093/biostatistics/1.4.441

Kemp, K. W. (1971). Formal Expressions which Can Be Applied to Cusum Charts. Journal of the Royal Statistical Society. Series B (Methodological), 33(3), 331-360. doi:10.1111/j.2517-6161.1971.tb01521.x

See Also

bernoulli_cusum, bernoulli_control_limit

Examples

#Determine a risk-adjustment model using a generalized linear model.
#Outcome (failure within 100 days) is regressed on the available covariates:
glmmodber <- glm((survtime <= 100) & (censorid == 1)~ age + sex + BMI,
                  data = surgerydat, family = binomial(link = "logit"))
#Determine the Average Run Length in number of outcomes for
#control limit h = 2.5 with (0, h) divided into n_grid = 200 segments
ARL <- bernoulli_ARL(h = 2.5, n_grid = 200, glmmod = glmmodber, theta = log(2))
#Calculate ARL, but now exploiting connection between SPRT and CUSUM:
#n_grid now decides the accuracy of the Trapezoidal rule for integral approximation
ARLSPRT <- bernoulli_ARL(h = 2.5, n_grid = 200, glmmod = glmmodber,
theta = log(2), method = "SPRT")

Average run length for Bernoulli CUSUM using Markov Chain methodology

Description

Internal function that discretizes grid and solves matrix equation involving transition matrix for Markov Chain methodology

Usage

bernoulli_ARL_MC(n_grid, R, h)

Arguments

Average run length for Bernoulli CUSUM using Integral Equation methodology

Description

Internal function that calculates the ARL using the connection between the ARL of a Wald SPRT and a CUSUM.

Usage

bernoulli_ARL_SPRT(h, n_grid, Wncdf, glmmod, theta, theta_true, p0,
  tol = 1e-06)

Arguments

Cumulative distribution function (cdf) of Run Length for Bernoulli CUSUM

Description

Calculate the cdf of the Run Length of the Bernoulli CUSUM, starting from initial value between 0 and h, using Markov Chain methodology.

Usage

bernoulli_RL_cdf(h, x, n_grid, glmmod, theta, theta_true, p0, p1,
  smooth_prob = FALSE, exact = TRUE)

Arguments

Details

Let K denote the run length of the Bernoulli CUSUM with control limit h, then this function can be used to evaluate P(K \leq x).

The formula on page 543 of Brook & Evans (1972) is used if exact = TRUE. When exact = FALSE, formula (3.9) on page 545 is used instead, approximating the transition matrix using its Jordan canonical form. This can save computation time considerably, but is not appropriate for small values of x.

Value

A list containing:

• Fr_0: A numeric value indicating the probability of the run length being smaller than x.

• Fr: A data.frame containing the cumulative distribution function of the run length depending on the state in which the process starts (E_0, E_1, \ldots, E_{n_{grid}-1})

  start_val:

  Starting value of the CUSUM, corresponding to the discretized state spaces E_{i};

  P(K <= x):

  Value of the cdf at x for the CUSUM with initial value start_val;

• R: A transition probability matrix containing the transition probabilities between states E_0, \ldots, E_{t-1}. R_{i,j} is the transition probability from state i to state j.

The value of ARL_0 will be printed to the console.

References

Brook, D., & Evans, D. A. (1972). An Approach to the Probability Distribution of CUSUM Run Length. Biometrika, 59(3), 539-549. doi:10.2307/2334805

Steiner, S. H., Cook, R. J., Farewell, V. T., & Treasure, T. (2000). Monitoring surgical performance using risk-adjusted cumulative sum charts. Biostatistics, 1(4), 441-452. doi:10.1093/biostatistics/1.4.441

Examples

#Determine a risk-adjustment model using a generalized linear model.
#Outcome (failure within 100 days) is regressed on the available covariates:
glmmodber <- glm((survtime <= 100) & (censorid == 1)~ age + sex + BMI,
                  data = surgerydat, family = binomial(link = "logit"))
#Determine probability of run length being less than 600
prob600 <- bernoulli_RL_cdf(h = 2.5, x = 600, n_grid = 200, glmmod = glmmodber, theta = log(2))

Average run length for Bernoulli CUSUM using Markov Chain methodology

Description

Internal function that discretizes grid and solves matrix equation involving transition matrix for Markov Chain methodology

Usage

bernoulli_cdf_MC(n_grid, R, r, h, exact = TRUE)

Arguments

Determine control limits for the Bernoulli CUSUM by simulation

Description

This function can be used to determine control limits for the Bernoulli CUSUM (bernoulli_cusum) procedure by restricting the type I error alpha of the procedure over time.

Usage

bernoulli_control_limit(time, alpha = 0.05, followup, psi, n_sim = 200,
  glmmod, baseline_data, theta, p0, p1, h_precision = 0.01, seed = 1041996,
  pb = FALSE, assist)

Arguments

Details

This function performs 3 steps to determine a suitable control limit.

• Step 1: Generates n_sim in-control units (failure rate as baseline). If data is provided, subject covariates are resampled from the data set.

• Step 2: Determines chart values for all simulated units.

• Step 3: Determines control limits such that at most a proportion alpha of all units cross the control limit.

The generated data as well as the charts are also returned in the output.

Value

A list containing three components:

• call: the call used to obtain output;

• charts: A list of length n_sim containing the constructed charts;

• data: A data.frame containing the in-control generated data.

• h: Determined value of the control limit.

Author(s)

Daniel Gomon

See Also

bernoulli_cusum

Other control limit simulation: bk_control_limit(), cgr_control_limit()

Examples

#We consider patient outcomes 100 days after their entry into the study.
followup <- 100

#Determine a risk-adjustment model using a generalized linear model.
#Outcome (failure within 100 days) is regressed on the available covariates:
exprfitber <- as.formula("(survtime <= followup) & (censorid == 1)~ age + sex + BMI")
glmmodber <- glm(exprfitber, data = surgerydat, family = binomial(link = "logit"))

#Determine control limit restricting type I error to 0.1 over 500 days
#using the risk-adjusted glm constructed on the baseline data.
a <- bernoulli_control_limit(time = 500, alpha = 0.1, followup = followup,
 psi = 0.5, n_sim = 10, theta = log(2), glmmod = glmmodber, baseline_data = surgerydat)

print(a$h)

Risk-adjusted Bernoulli CUSUM

Description

This function can be used to construct a risk-adjusted Bernoulli CUSUM chart for survival data. It requires the specification of one of the following combinations of parameters as arguments to the function:

• glmmod & theta

• p0 & theta

• p0 & p1

Usage

bernoulli_cusum(data, followup, glmmod, theta, p0, p1, h, stoptime, assist,
  twosided = FALSE)

Arguments

Details

The Bernoulli CUSUM chart is given by

S_n = \max(0, S_{n-1} + W_n),

where

W_n = X_n \ln \left( \frac{p_1 (1-p_0)}{p_0(1-p_1)} \right) + \ln \left( \frac{1-p_1}{1-p_0} \right)

and X_n is the outcome of the n-th (chronological) subject in the data. In terms of the Odds Ratio:

W_n = X_n \ln \left( e^\theta \right) + \ln \left( \frac{1}{1-p_0 + e^\theta p_0} \right)

For a risk-adjusted procedure (when glmmod is specified), a patient specific baseline failure probability p_{0i} is modelled using logistic regression first. Instead of the standard practice of displaying patient numbering on the x-axis, the time of outcome is displayed.

Value

An object of class bercusum containing:

• CUSUM: A data.frame containing the following named columns:

  time:

  times at which chart is constructed;

  value:

  value of the chart at corresponding times;

  numobs:

  number of observations at corresponding times.

• call: the call used to obtain output;

• glmmod: coefficients of the glm() used for risk-adjustment, if specified;

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

Author(s)

Daniel Gomon

See Also

plot.bercusum, runlength.bercusum

Examples

#We consider patient outcomes 100 days after their entry into the study.
followup <- 100
#Determine a risk-adjustment model using a generalized linear model.
#Outcome (failure within 100 days) is regressed on the available covariates:
exprfitber <- as.formula("(survtime <= followup) & (censorid == 1)~ age + sex + BMI")
glmmodber <- glm(exprfitber, data = surgerydat, family = binomial(link = "logit"))
#Construct the Bernoulli CUSUM on the 1st hospital in the data set.
bercus <- bernoulli_cusum(data = subset(surgerydat, unit == 1), glmmod = glmmodber,
 followup = followup, theta = log(2))
#Plot the Bernoulli CUSUM
plot(bercus)

Determine control limits for BK-CUSUM by simulation

Description

This function can be used to determine control limits for the BK-CUSUM (bk_cusum) procedure by restricting the type I error alpha of the procedure over time.

Usage

bk_control_limit(time, alpha = 0.05, psi, n_sim = 200, theta, coxphmod,
  baseline_data, cbaseh, inv_cbaseh, interval = c(0, 9e+12),
  h_precision = 0.01, seed = 1041996, pb = FALSE, chartpb = FALSE,
  assist)

Arguments

Details

This function performs 3 steps to determine a suitable control limit.

• Step 1: Generates n_sim in-control units (failure rate as baseline). If data is provided, subject covariates are resampled from the data set.

• Step 2: Determines chart values for all simulated units.

• Step 3: Determines control limits such that at most a proportion alpha of all units cross the control limit.

The generated data as well as the charts are also returned in the output.

Value

A list containing three components:

• call: the call used to obtain output;

• charts: A list of length n_sim containing the constructed charts;

• data: A data.frame containing the in-control generated data.

• h: Determined value of the control limit.

• achieved_alpha: Achieved type I error on the sample of n_sim simulated units.

Author(s)

Daniel Gomon

See Also

bk_cusum

Other control limit simulation: bernoulli_control_limit(), cgr_control_limit()

Examples

require(survival)

#Determine a cox proportional hazards model for risk-adjustment
exprfit <- as.formula("Surv(survtime, censorid) ~ age + sex + BMI")
tcoxmod <- coxph(exprfit, data= surgerydat)

#Determine a control limit restricting type I error to 0.1 over 500 days
#with specified cumulative hazard function without risk-adjustment
a <- bk_control_limit(time = 500, alpha = 0.1, theta = log(2),
cbaseh = function(t) chaz_exp(t, lambda = 0.02),
inv_cbaseh = function(t) inv_chaz_exp(t, lambda = 0.02), psi = 0.5,
n_sim = 10)

#Determine a control limit restricting type I error to 0.1 over 500 days
#using the risk-adjusted cumulative hazard determined using coxph()
b <- bk_control_limit(time = 500, alpha = 0.1, theta = log(2),
coxphmod = tcoxmod, psi = 0.5, n_sim = 10)

print(a$h)
print(b$h)

Continuous time BK-CUSUM

Description

This function implements the BK-CUSUM procedure based on the Biswas & Kalbfleisch (2008) CUSUM. To construct the Biswas & Kalbfleisch (2008) CUSUM, set C = 1 (years) or C = 365 (days). For detection purposes, it is sufficient to determine the value of the chart at the times of failure. This can be achieved by leaving ctimes unspecified. The function requires the specification of theta and 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

Usage

bk_cusum(data, theta, coxphmod, cbaseh, ctimes, h, stoptime, C,
  twosided = FALSE, pb = FALSE, assist)

Arguments

Details

The BK-CUSUM can be used to test the alternative hypothesis of an instant change of fixed size e^\theta in the subject specific hazard rate from h_i(t) to h_i(t) e^\theta. The parameter C can be used to ignore the contributions of subjects, C time units after their entry into the study. The BK-CUSUM is constructed as

G(t) = \max_{0 \leq k \leq t} \left( \theta N(k,t) - \left( e^\theta -1 \right) \Lambda(k,t) \right),

where \theta is the log expected hazard ratio,

N(k,t) = N(t) - N(k)

with N(t) the counting process of all failures at time t, and

\Lambda(k,t) = \Lambda(t) - \Lambda(k)

with \Lambda(t) the summed cumulative intensity of all subjects at time t.

Value

An object of class bkcusum containing:

• BK: a data.frame containing the following named columns:

  time:

  times at which chart is constructed;

  value:

  value of the chart at corresponding times.

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

• call: the call used to obtain output;

• h: Specified value for the control limit.

Author(s)

Daniel Gomon

References

Biswas P. and Kalbfleisch J.D. (2008), A risk-adjusted CUSUM in continuous time based on the Cox Model, Statistics in medicine 27, 3452-3452. doi:10.1002/sim.3216

See Also

plot.bkcusum, runlength.bkcusum

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

Examples

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

#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 <- Surv(survtime, censorid) ~ age + sex + BMI
tcoxmod <- coxph(exprfit, data= surgerydat)

#Determine the values of the chart
bk <- bk_cusum(data = tdat, theta = log(2), coxphmod = tcoxmod, cbaseh = tcbaseh, pb = TRUE)
#plot the BK-CUSUM (exact hazard)
plot(bk)

#Alternatively, cbaseh can be left empty when specifying coxphmod through coxph()
bk_cox <- bk_cusum(data = tdat, theta = log(2), coxphmod = tcoxmod, pb = TRUE)
#plot the BK-CUSUM (estimated hazard from coxph)
plot(bk_cox)

Survival after breast cancer surgery

Description

Data about patient survival after their breast cancer surgery procedure performed at one of the 15 units participating in a cancer treatment study. The data is based on a trial performed at the European Organisation for Research and Treatment of Cancer (EORTC).

Usage

breast

Format

A data.frame with 2663 rows and 11 variables:

entrytime

Chronological time of entry of patient into study/time of surgery (numeric)

survtime

Time from entry until failure of patient (numeric)

censorid

Censoring indicator (0 - right censored, 1 - observed) (integer)

unit

Unit number at which patient received treatment (integer)

var1-7

Covariates associated with patient (factor)

Source

Based on trial data from the European Organisation for Research and Treatment of Cancer, https://www.eortc.org/

Examples

#Determine the estimated arrival rate for all units in the data
arrival_rate(breast)

#Plot Quality Control charts for unit 11 in the study
library(survival)
phmodbreast <- coxph(Surv(survtime, censorid) ~ . - entrytime - unit ,
data = breast)
glmmodbreast <- glm((survtime <= 36) & (censorid == 1) ~ . - entrytime - unit,
 data = breast, family = binomial(link = "logit"))

par(mfrow = c(1, 3))
p1 <- plot(cgr_cusum(data = subset(breast, unit == 11),  coxphmod = phmodbreast)) +
ggtitle("CGR-CUSUM")
p2 <- plot(bk_cusum(data = subset(breast, unit == 11),  coxphmod = phmodbreast,
 theta = log(2))) + ggtitle("BK-CUSUM")
p3 <- plot(bernoulli_cusum(data = subset(breast, unit == 11), followup = 36,
 glmmod = glmmodbreast, theta = log(2))) + ggtitle("Bernoulli CUSUM")
p4 <- plot(funnel_plot(data = breast, glmmod = glmmodbreast, followup = 36 )) +
ggtitle("Funnel plot")

if(require("gridExtra")){
 grid.arrange(p1, p2, p3, p4, nrow = 2)
}

Transition probability matrix for Bernoulli CUSUM

Description

Calculates the transition probability matrix for the Bernoulli CUSUM described in Brook & Evans (1972).

Usage

calc_MC_trans_matrix(h, n_grid, Wncdf, glmmod, p0, theta, theta_true)

Arguments

Calculate cdf of singletons W_n for CUSUM

Description

Internal function to calculate cdf of singletons W_n of the Bernoulli CUSUM chart. The cdf is used to create the transition matrix when Markov Chain methodology is used or to determine the integral equation/probabilities of a Wald test when integral equation or Kemp's methodology is used.

Usage

calc_Wncdf(glmmod, theta, theta_true, p0, smooth_prob = FALSE)

Arguments

Calculate the Cox Proportional hazards relative risk associated with the covariates of subjects

Description

This function can be used to calculate the change in the relative risk of a subject pertaining to their covariates under a specified Cox proportional hazards model.

Usage

calc_risk(data, coxphmod = NULL)

Arguments

Details

The subject-specific relative risk is e^{\beta Z_i}, where \beta is a vector of regression coefficients and Z_i a vector of covariates for subject i.

Value

A vector of nrow(data) specifying the increased/decreased risk of failure for each subject.

Author(s)

Daniel Gomon

Examples

#Small example data
crdat <- data.frame(age = rnorm(10, 40, 5), BMI = rnorm(10, 24, 3))
#Example risk-adjustment list (can also specify coxphmod)
crlist <- list(formula = as.formula("~age + BMI"), coefficients = c("age"= 0.02, "BMI"= 0.009))
#Calculate the increase or decrease of the relative risk for the subjects
#in crdat.
calc_risk(crdat, crlist)

Determine control limits for CGR-CUSUM by simulation

Description

This function can be used to determine control limits for the CGR-CUSUM (cgr_cusum) procedure by restricting the type I error alpha of the procedure over time.

Usage

cgr_control_limit(time, alpha = 0.05, psi, n_sim = 20, coxphmod,
  baseline_data, cbaseh, inv_cbaseh, interval = c(0, 9e+12),
  h_precision = 0.01, ncores = 1, seed = 1041996, pb = FALSE,
  chartpb = FALSE, detection = "upper", maxtheta = log(6), assist)

Arguments

Details

This function performs 3 steps to determine a suitable control limit.

• Step 1: Generates n_sim in-control units (failure rate as baseline). If data is provided, subject covariates are resampled from the data set.

• Step 2: Determines chart values for all simulated units.

• Step 3: Determines control limits such that at most a proportion alpha of all units cross the control limit.

The generated data as well as the charts are also returned in the output.

Value

A list containing three components:

• call: the call used to obtain output;

• charts: A list of length n_sim containing the constructed charts;

• data: A data.frame containing the in-control generated data.

• h: Determined value of the control limit.

• achieved_alpha: Achieved type I error on the sample of n_sim simulated units.

Author(s)

Daniel Gomon

See Also

cgr_cusum

Other control limit simulation: bernoulli_control_limit(), bk_control_limit()

Examples

require(survival)

#Determine a cox proportional hazards model for risk-adjustment
exprfit <- as.formula("Surv(survtime, censorid) ~ age + sex + BMI")
tcoxmod <- coxph(exprfit, data= surgerydat)

#Determine a control limit restricting type I error to 0.1 over 500 days
#with specified cumulative hazard function without risk-adjustment
a <- cgr_control_limit(time = 500, alpha = 0.1, cbaseh = function(t) chaz_exp(t, lambda = 0.02),
inv_cbaseh = function(t) inv_chaz_exp(t, lambda = 0.02), psi = 0.5, n_sim = 10)


#Determine a control limit restricting type I error to 0.1 over 500 days
#using the risk-adjusted cumulative hazard found using coxph()
b <- cgr_control_limit(time = 500, alpha = 0.1, coxphmod = tcoxmod, psi = 0.5, n_sim = 10,
baseline_data = subset(surgerydat, unit == 1))

print(a$h)
print(b$h)

Continuous time Generalized Rapid response CUSUM (CGR-CUSUM)

Description

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

Usage

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

Arguments

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) e^\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),

where

N(\geq \nu)(t) = \sum_{i \geq \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),

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.

Value

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.

Author(s)

Daniel Gomon

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 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, doi:10.32614/rj-2023-095

See Also

plot.cgrcusum, runlength.cgrcusum

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

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)

Exponential hazard, cumulative hazard and inverse cumulative hazard

Description

Functions which return the hazard, cumulative hazard and inverse cumulative hazard at time t for an exponential distribution with parameter \lambda and true hazard ratio \mu.

Usage

haz_exp(t, lambda, mu = log(1))

chaz_exp(t, lambda, mu = log(1))

inv_chaz_exp(t, lambda, mu = log(1))

Arguments

Details

The hazard function of an exponential distribution is given by:

h(t| \lambda, \mu) = \lambda e^\mu

The cumulative hazard (with true hazard ratio \mu) is given by:

H(t| \lambda, \mu) = \lambda t e^\mu

The inverse cumulative hazard (with true hazard ratio \mu) by:

H^{-1}(t| \lambda, \mu) = \frac{t}{\lambda e^\mu}

Value

Value of specified function at time t.

Extract (inverse) cumulative baseline hazard from Cox PH model

Description

Extracts a function which returns the (inverse) cumulative baseline hazard from a coxph() call.

Usage

extract_hazard(coxphmod)

Arguments

Details

The baseline hazard is extracted from the coxph() call using the basehaz() function. The baseline hazard function is then smoothed using approxfun() to obtain the linear interpolant. If required, the inverse baseline hazard is determined using root linear interpolation. For this, a function written by Zheyuan Li (see references) is used.

Value

A list containing:

• cbaseh: A function which returns the cumulative baseline hazard at specified time;

• inv_cbaseh: A function which returns the inverse cumulative baseline hazard at specified time.

• max_time: maximal time at which cbaseh is known;

• max_haz: value of maximal hazard (at maximum time).

Author(s)

Daniel Gomon

References

Zheyuan Li: How to estimate x value from y value input after approxfun in R? (accessed: 09/10/2023)

See Also

coxph

Examples

require(survival)
exprfit <- as.formula("Surv(survtime, censorid) ~ age + sex + BMI")
tcoxmod <- coxph(exprfit, data= surgerydat)
tcox_hazard_fcts <- extract_hazard(tcoxmod)

Risk-adjusted funnel plot

Description

This function allows to construct a risk-adjusted funnel plot for comparing survival proportion between units, see Spiegelhalter (2005).

Usage

funnel_plot(data, ctime, p0, glmmod, followup, predlim = c(0.95, 0.99),
  assist)

Arguments

Value

An object of class "funnelplot" containing:

• data: A data.frame containing:

  unit:

  unit number/name;

  observed:

  observed number of failures at unit;

  expected:

  expected (risk-adjusted) number of failures at unit;

  numtotal

  total number of individuals considered at this unit;

  p:

  (risk-adjusted) proportion of failure at unit;

  predlimels:

  worse/in-control/better performance than expected at specified confidence levels.

• call: the call used to obtain output

• plotdata: data used for plotting confidence intervals

• predlim: specified confidence level(s)

• p0: (Estimated) baseline failure probability

Author(s)

Daniel Gomon

References

Spiegelhalter D. J. (2005). Funnel plots for comparing institutional performance. Statistics in medicine, 24(8), 1185-1202. doi:10.1002/sim.1970

See Also

plot.funnelplot, summary.funnelplot

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

Examples

#Determine a risk-adjustment model using a generalized linear model.
#Outcome (survival in first 100 days) is regressed on the available covariates:
exprfitfunnel <- as.formula("(survtime <= 100) & (censorid == 1)~ age + sex + BMI")
glmmodfun <- glm(exprfitfunnel, data = surgerydat, family = binomial(link = "logit"))
#Determine the necessary values to produce a funnel plot
funnel <- funnel_plot(data = surgerydat, ctime = 3*365, glmmod = glmmodfun, followup = 100)
#Produce a funnel plot!
plot(funnel)
## Not run: 
require(plotly)
#Create an interactive plot!
ggplotly(plot(funnel))

## End(Not run)

Generate arrival times according to a Poisson point process

Description

This function can be used to generate arrival times for a Poisson point process with rate psi up until time t.

Usage

gen_arriv_times(psi, t)

Arguments

Details

Exponential(\psi) interarrival times.

Value

A vector of arrival times up until time t.

Author(s)

Daniel Gomon

Examples

set.seed(123)
gen_arriv_times(psi = 0.3, t = 5)
gen_arriv_times(psi = 0.3, t = 20)

Generate survival times

Description

Generate survival times according to hazard rate h(t) \exp(\mu) with h(t) the hazard rate associated with the specified inverse cumulative hazard rate invchaz and \mu the specified true hazard ratio mu. See Bender et al. (2005).

Usage

gen_surv_times(invchaz, mu = log(1), data, coxphmod = NULL)

Arguments

Details

Sometimes it is desirable to generate survival times from an increased hazard rate

h(t, \mu) = h_0(t) e^\mu

with h_0 the baseline hazard rate. We call e^\mu the true hazard ratio.

Value

A vector of survival times from subject entry time.

Author(s)

Daniel Gomon

References

Bender, R., Augustin, T., & Blettner, M. (2005). Generating survival times to simulate Cox proportional hazards models. Statistics in Medicine, 24(11), 1713-1723. doi:10.1002/sim.2059

Examples

set.seed(123)
gen_surv_times(invchaz = function(t) inv_chaz_exp(t, lambda = 0.01), data = 5)

Generate units with specified failure rate

Description

Generate n_sim units with subjects arriving according to a Poisson process with rate psi until time. Failure rate is determined either from cbaseh or inv_cbaseh, or from specified coxphmod. Covariates will be resampled from baseline_data if specified.

Usage

generate_units(time, psi, n_sim = 20, cbaseh, inv_cbaseh, coxphmod = NULL,
  baseline_data, interval = c(0, 9e+12), mu = 0)

Arguments

Details

In a Poisson arrival process, inter-arrival times are exponentially distributed with parameter psi. If cbaseh is specified, the inverse baseline hazard will be determined using uniroot(). The times of failure are then determined using gen_surv_times().

Value

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

entrytime:

time of subject entry into the study;

survtime:

survival time of subject;

censorid:

censoring indicator, 0 = censored, 1 = observed;

unit:

unit number;

expmu:

exponent of the log hazard ratio used to generate survival times;

psival:

arrival rate at unit;

covariates:

covariates resampled from baseline_data.

Author(s)

Daniel Gomon

Examples

require(survival)
#Fit a Cox model
exprfit <- as.formula("Surv(survtime, censorid) ~ age + sex + BMI")
tcoxmod <- coxph(exprfit, data= surgerydat)

#Generate 30 hospitals with on average 2 patients per day arriving
#according to the Cox model determined above, with resampling from the
#original data set. The hazard rate at the hospitals is twice the baseline
#hazard.
generate_units(time = 50, psi = 2, n_sim = 30, coxphmod = tcoxmod,
baseline_data = surgerydat, mu = log(2))

Plot a list of CUSUM charts (interactive)

Description

Create an interactive plot visualizing a combination of control charts which can be created using this package.

Usage

interactive_plot(x, unit_names, scale = FALSE, group_by = c("none", "unit",
  "type"), highlight = FALSE, manual_colors = c(), ...)

Arguments

Value

An interactive plot will be produced in the current graphics device. For more information on the possibilities for interaction, see https://plotly.com/r/.

See Also

cgr_cusum, bk_cusum, bernoulli_cusum, funnel_plot

Examples

require(survival)
#Extract data to construct CUSUM charts on
tdat <- subset(surgerydat, unit == 1 & entrytime < 365)
tdat2 <- subset(surgerydat, unit == 2 & entrytime < 365)

#Determine model parameters
followup <- 100
tcbaseh <- function(t) chaz_exp(t, lambda = 0.01)
exprfit <- as.formula("Surv(survtime, censorid) ~ age + sex + BMI")
tcoxmod <- coxph(exprfit, data= surgerydat)
exprfitber <- as.formula("(survtime <= followup) & (censorid == 1)~ age + sex + BMI")
glmmodber <- glm(exprfitber, data = surgerydat, family = binomial(link = "logit"))


#Construct the charts
cgr <- cgr_cusum(data = tdat, coxphmod = tcoxmod, cbaseh = tcbaseh, pb = TRUE)
cgr$h <- 8.29
bk <- bk_cusum(data = tdat, theta = log(2), coxphmod = tcoxmod, cbaseh = tcbaseh, pb = TRUE)
bk$h <- 6.23
bercus <- bernoulli_cusum(data = subset(surgerydat, unit == 1), glmmod = glmmodber,
followup = followup, theta = log(2))
bercus$h <- 3.36
bk2 <- bk_cusum(data = tdat2, theta = log(2), coxphmod = tcoxmod, cbaseh = tcbaseh, pb = TRUE)
bk2$h <- 6.23

#Create the plot
interactive_plot(list(cgr, bk, bercus, bk2), unit_names =
c("hosp1", "hosp1", "hosp1", "hosp2"))

Assist users in parameter selection

Description

This function can be used to determine some of the vital parameters used to construct control charts in this package.

Usage

parameter_assist(baseline_data, data, formula, followup, theta = log(2),
  time, alpha = 0.05, maxtheta = log(6))

Arguments

Details

Depending on the specified arguments, the function will return parameters. If covariate_names is not specified, the returned risk-adjustment models will be trivial. If formula is not specified but covariate_names are, the function assumes the simplest form for the regression model (cov1 + cov2 + ...). If followup is not specified, no glmmod will be determined

Value

A list of parameters to feed to quality control charts in this package:

• call: The call used to obtain output.

• data: The data used in the call to the function.

• baseline_data: The baseline_data used in the call to the function

• glmmod: A glm() model which can be fed to the funnel_plot() and bernoulli_cusum() functions.

• coxphmod: A coxph() model which can be fed to the cgr_cusum() and cgr_cusum() functions.

• theta: Expected increase in the odds/hazard ratio.

• psi: Estimated Poisson arrival rate in data.

• time: Time frame over which to restrict type I error.

• alpha: Desired level of type I error for control limit determination.

• maxtheta: Maximum expected increase/decrease in the odds/hazard ratio.

Author(s)

Daniel Gomon

Examples

require(survival)

#Minimal example - no risk-adjustment
pars_min <- parameter_assist(baseline_data = surgerydat,
data = subset(surgerydat, unit == 1))

#Specifying all parameters
pars <- parameter_assist(baseline_data = surgerydat,
data = subset(surgerydat, unit == 1),
formula = formula("survtime ~ age + sex + BMI"), followup = 100)

Plot a quality control chart

Description

Plot a cgrcusum, bkcusum, bercusum or funnelplot chart, or a list containing a combination of 'bercusum', 'bkcusum' and 'cgrcusum' charts.

Usage

## S3 method for class 'cgrcusum'
plot(x, h, ...)

## S3 method for class 'bkcusum'
plot(x, h, ...)

## S3 method for class 'funnelplot'
plot(x, percentage = TRUE, unit_label = TRUE,
  label_size = 3, col_fill = "blue", ...)

## S3 method for class 'bercusum'
plot(x, h = x$h, ...)

Arguments

Value

A plot of the associated chart is displayed in the current graphics device.

Methods (by class)

• plot(cgrcusum): Plot a CGR-CUSUM

• plot(bkcusum): Plot a BK-CUSUM

• plot(funnelplot): Display a funnel plot

• plot(bercusum): Plot a Bernoulli CUSUM

Author(s)

Daniel Gomon

See Also

cgr_cusum, bk_cusum, bernoulli_cusum, funnel_plot

Determine run length of a CUSUM chart

Description

This function can be used to calculate the run length of a 'cgrcusum', 'bkcusum' or 'bercusum' chart when using control limit h

Usage

runlength(chart, h)

## S3 method for class 'cgrcusum'
runlength(chart, h, ...)

## S3 method for class 'bkcusum'
runlength(chart, h, ...)

## S3 method for class 'bercusum'
runlength(chart, h, ...)

Arguments

Value

The run length of the chart with the given control limit.

Methods (by class)

• runlength(cgrcusum): determines runlength of cgrcusum object

• runlength(bkcusum): determines runlength of bkcusum object

• runlength(bercusum): determines runlength of bercusum object

Author(s)

Daniel Gomon

Examples

exprfitber <- as.formula("(survtime <= 100) & (censorid == 1) ~ age + sex + BMI")
glmmodber <- glm(exprfitber, data = surgerydat, family = binomial(link = "logit"))
bercus <- bernoulli_cusum(data = subset(surgerydat, unit == 14), glmmod = glmmodber,
                   followup = 100, theta = log(2))
#Determine the run length of the above Bernoulli CUSUM when using a control limit
#of h = 1.
runlength(bercus, h = 1)

Summarizes S3 objects in this package.

Description

Prints a summary of the funnelplot object.

Usage

## S3 method for class 'funnelplot'
summary(object, ...)

Arguments

Value

A data.frame with:

• unit: unit number/identifier;

• observed: the observed amount of failures at respective unit;

• expected: the expected amount of failures at respective unit, given that the unit is performing at target;

• numtotal: total number of subjects at respective unit;

• p: estimated probability of failure at unit;

• '0.xx': better/normal/worse proportion of failure at specified confidence levels.

Methods (by class)

• summary(funnelplot): summarize instances detected by the funnelplot object

See Also

funnel_plot

Simulated data set with data of surgery procedures performed at multiple hospitals.

Description

Data about patients and their surgery procedure from 45 simulated hospitals with patient arrivals in the first 400 days after the start of the study.
Patient survival times were determined using a risk-adjusted Cox proportional hazards model with coefficients age = 0.003, BMI = 0.02 and sexmale = 0.2 and exponential baseline hazard rate h_0(t, \lambda = 0.01) e^\mu. The increase in hazard rate is sampled from a normal distribution for all hospitals:

• \theta \sim N(log(1), sd = 0.4)

This means that the average failure rate of hospitals in the data set should be baseline (\theta = 0), with some hospitals experiencing higher and lower failure rates. True failure rate can be found in the column exptheta.
The arrival rate \psi of patients at a hospital differs. The arrival rates are:

• Hospitals 1-5 & 16-20: 0.5 patients per day (small hospitals)

• Hospitals 6-10 & 21-25: 1 patient per day (medium sized hospitals)

• Hospitals 11-15 & 26-30: 1.5 patients per day (large hospitals)

These are then respectively small, medium and large hospitals.

Usage

surgerydat

Format

A data.frame with 12010 rows and 9 variables:

entrytime

Time of entry of patient into study (numeric)

survtime

Time from entry until failure of patient (numeric)

censorid

Censoring indicator (0 - right censored, 1 - observed) (integer)

unit

Hospital number at which patient received treatment (integer)

exptheta

True excess hazard used for generating patient survival (numeric)

psival

Poisson arrival rate at hospital which the patient was at (numeric)

age

Age of the patient (numeric)

sex

Sex of the patient (factor)

BMI

Body mass index of the patient (numeric)

Weibull hazard, cumulative hazard and inverse cumulative hazard

Description

Functions which return the hazard, cumulative hazard and inverse cumulative hazard at time t for a Weibull distribution with shape parameter \lambda, scale parameter \theta and true hazard ratio \mu.

Usage

haz_weib(t, lambda, theta, mu = log(1))

chaz_weib(t, lambda, theta, mu = log(1))

inv_chaz_weib(t, lambda, theta, mu = log(1))

Arguments

Details

The hazard function of a Weibull distribution is given by:

h(t| \lambda, \theta, \mu) = \frac{\lambda}{\theta} \left(\frac{t}{\theta} \right)^{\lambda -1} e^\mu

The cumulative hazard (with true hazard ratio \mu) is given by:

H(t| \lambda, \theta, \mu) = \left( \frac{t}{\theta} \right)^{\lambda} e^\mu

The inverse cumulative hazard (with true hazard ratio \mu) by:

H^{-1}(t| \lambda, \theta, \mu) = \theta \left( \frac{t}{e^\mu} \right)^{1/\lambda}

Value

Value of specified function at time t.