cgr_control_limit {success} | R Documentation |
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
time |
A numeric value over which the type I error |
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
|
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 greatly increase computation times. Default is 20. |
coxphmod |
(optional): A cox proportional hazards regression model as
produced by
the function
|
baseline_data |
(optional): A
and optionally additional covariates used for risk-adjustment. Can only be specified
in combination with |
cbaseh |
(optional): A function that returns the unadjusted cumulative
baseline hazard |
inv_cbaseh |
(optional): A function that returns the unadjusted inverse cumulative
baseline
hazard |
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. |
ncores |
(optional): Number of cores to use to parallelize the computation of the
CGR-CUSUM charts. If ncores = 1 (default), no parallelization is done. You
can use |
seed |
(optional): A numeric seed for survival time generation. Default = my birthday. |
pb |
(optional): A boolean indicating whether a progress bar should
be shown. Default is |
chartpb |
(optional): A boolean indicating whether progress bars should
be displayed for the constructions of the charts. Default is |
detection |
Should the control limit be determined for an
|
maxtheta |
(optional): Maximum value of maximum likelihood estimate for
parameter |
assist |
(optional): Output of the function |
Details
This function performs 3 steps to determine a suitable control limit.
Step 1: Generates
n_sim
in-control units (failure rate as baseline). Ifdata
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 lengthn_sim
containing the constructed charts; -
data
: Adata.frame
containing the in-control generated data. -
h
: Determined value of the control limit. -
achieved_alpha
: Achieved type I error on the sample ofn_sim
simulated units.
Author(s)
Daniel Gomon
See Also
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)