bernoulli_control_limit {success}

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

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

p_1 = \frac{p_0 e^\theta}{1-p_0 +p_0 e^\theta}.

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

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)