bk_control_limit function

Determine control limits for BK-CUSUM by simulation

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.

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

• time: A numeric value over which the type I error alpha must be restricted.

• alpha: A proportion between 0 and 1 indicating the required maximal type I error. Default is 0.05.

• psi: A numeric value indicating the estimated Poisson arrival rate of subjects at their respective units. Can be determined using parameter_assist().

• n_sim: An integer value indicating the amount of units to generate for the determination of the control limit. Larger values yield more precise control limits, but increase computation times. Default is 200.

• theta: The expected log-hazard ratio $\theta$ under the alternative hypothesis. If $\theta >= 0$, the chart will try to detect an increase in hazard ratio (upper one-sided). If $\theta < 0$, the chart will look for a decrease in hazard ratio (lower one-sided).

• 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.

• baseline_data: (optional): A data.frame used for covariate resampling with rows representing subjects and at least 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. Can only be specified in combination with coxphmod.

• 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.

• inv_cbaseh: (optional): A function that returns the unadjusted inverse cumulative baseline hazard $H_0^-1(t)$. If inv_cbaseh is missing, it will be determined from cbaseh numerically.

• interval: (optional): Interval in which survival times should be solved for numerically.

• h_precision: (optional): A numerical value indicating how precisely the control limit should be determined. By default, control limits will be determined up to 2 significant digits.

• seed: (optional): A numeric seed for survival time generation. Default is 01041996 (my birthday).

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

• chartpb: (optional): A boolean indicating whether progress bars should be displayed for the constructions of the charts. Default is FALSE.

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

Returns

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.

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.

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)

See Also

bk_cusum

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

Author(s)

Daniel Gomon