bernoulli_control_limit function

Determine control limits for the Bernoulli CUSUM by simulation

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.

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

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

• followup: The value of the follow-up time to be used to determine event time. Event time will be equal to entrytime + followup for each subject.

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

• glmmod: Generalized linear regression model used for risk-adjustment as produced by the function glm(). Suggested:

  glm(as.formula("(survtime <= followup) & (censorid == 1) ~covariates"), data = data).

  Alternatively, a list containing 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.

• theta: The $\theta$ value used to specify the odds ratio $e^\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). Note that

• p0: The baseline failure probability at entrytime + followup for individuals.
• p1: The alternative hypothesis failure probability at entrytime + followup for individuals.
• 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.
• 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.

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

#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)

See Also

bernoulli_cusum

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

Author(s)

Daniel Gomon