vlad (version 0.2.2)

eocusum_ad_sim: Compute steady-state ARLs of EO-CUSUM control charts using simulation

Description

Compute steady-state ARLs of EO-CUSUM control charts using simulation.

Usage

eocusum_ad_sim(r, pmix, k, h, RQ = 1, side = "low", type = "cond", m = 50)

Arguments

r

Integer. Number of of simulation runs.

pmix

Data Frame. A three column data frame. First column is the operation outcome. Second column are the predicted probabilities from the risk model. Third column can be either the predicted probabilities from the risk model or average outcome.

k

Double. Reference value of the CUSUM control chart. Either 0 or a positive value. Can be determined with function optimal_k.

h

Double. Decision interval (alarm limit, threshold) of the CUSUM control chart.

RQ

Double. Defines the true performance of a surgeon with the odds ratio ratio of death RQ. Use RQ = 1 to compute the in-control ARL and other values to compute the out-of-control ARL.

side

Character. Default is "low" to calculate ARL for the upper arm of the V-mask. If side = "up", calculate the lower arm of the V-mask.

type

Character. Default argument is "cond" for computation of conditional steady-state. Other option is the cyclical steady-state "cycl".

m

Integer. Simulated in-control observations.

Value

Returns a single value which is the Run Length.

References

Wittenberg P, Gan FF, Knoth S (2018). A simple signaling rule for variable life-adjusted display derived from an equivalent risk-adjusted CUSUM chart. Statistics in Medicine, 37(16), pp 2455--2473.

Taylor HM (1968). The Economic Design of Cumulative Sum Control Charts. Technometrics, 10(3), pp. 479--488.

Crosier R (1986). A new two-sided cumulative quality control scheme. Technometrics, 28(3), pp. 187--194.

Examples

# NOT RUN {
data("cardiacsurgery", package = "spcadjust")
library("dplyr")

## preprocess data to 30 day mortality and subset phase I/II
cardiacsurgery <- cardiacsurgery %>% rename(s = Parsonnet) %>%
  mutate(y = ifelse(status == 1 & time <= 30, 1, 0),
         phase = factor(ifelse(date < 2*365, "I", "II")))

s5000 <- sample_n(cardiacsurgery, size = 5000, replace = TRUE)
df1 <- select(cardiacsurgery, s, y)
df2 <- select(s5000, s, y)

## estimate coefficients from logit model
coeff1 <- round(coef(glm(y ~ s, data = df1, family = "binomial")), 3)
coeff2 <- round(coef(glm(y ~ s, data = df2, family = "binomial")), 3)

## Number of simulation runs
m <- 10^3
## Number of cores
nc <- parallel::detectCores()
# steady state
RNGkind("L'Ecuyer-CMRG")
m <- 10^3
tau <- 50
kopt <- optimal_k(QA = 2, df = S2I, coeff = coeff1, yemp = FALSE)
# eocusum_arloc_h_sim(L0 = 370, df = df1, k = kopt, m = m, side = "low", coeff = coeff1,
 coeff2 = coeff2, nc = nc)
res <- sapply(0:(tau-1), function(i){
  RLS <- do.call(c, parallel::mclapply( 1:m, eocusum_ad_sim, k = kopt, QS = 2, h = 2.637854,
  df = df1, m = i, coeff = coeff1, coeff2 = coeff2, side = "low", mc.cores = nc))
  list(data.frame(cbind(ARL = mean(RLS), ARLSE = sd(RLS)/sqrt(m))))
} )
RES <- data.frame(cbind(M = 0:(tau-1), do.call(rbind, res)))
ggplot2::qplot(x = M, y = ARL, data = RES, geom = c("line", "point")) +
ggplot2::theme_classic()
# }