vlad (version 0.2.0)

eocusum_arl_sim: Compute ARLs of EO-CUSUM control charts using simulation

Description

Compute ARLs of EO-CUSUM control charts using simulation.

Usage

eocusum_arl_sim(r, k, h, df, coeff, yemp = TRUE, side = "low")

Arguments

r

Integer. Number of of simulation runs.

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.

df

Data Frame. First column Parsonnet Score and second column outcome of each operation.

coeff

Numeric Vector. Estimated coefficients \(\alpha\) and \(\beta\) from the binary logistic regression model.

yemp

Logical. If TRUE use observed outcome value, if FALSE use estimated binary logistc regression model.

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.

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.

Examples

# NOT RUN {
library("dplyr")
library("tidyr")
library(ggplot2)

## Datasets
data("cardiacsurgery", package = "spcadjust")
cardiacsurgery <- cardiacsurgery %>% rename(s = Parsonnet) %>%
  mutate(y = ifelse(status == 1 & time <= 30, 1, 0))
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)

## set up
RNGkind("L'Ecuyer-CMRG")
m <- 10^3
kopt <- optimal_k(QA = 2, df = S2I, coeff = coeff1, yemp = FALSE)
h <- eocusum_arloc_h_sim(L0 = 370, df = df1, k = kopt, m = m, side = "low", coeff = coeff1,
                         coeff2 = coeff2, nc = 4)

## Serial simulation
RLS <- do.call(c, lapply(1:m, eocusum_arloc_sim, h = h, k = kopt, df = df1, side = "low",
                         coeff = coeff1, coeff2 = coeff2))
data.frame(cbind(ARL = mean(RLS), ARLSE = sd(RLS)/sqrt(m)))

## Parallel simulation (FORK)
RLS <- simplify2array(parallel::mclapply(1:m, eocusum_arloc_sim, h = h, k = kopt, df = df1,
                                         side = "low", coeff = coeff1, coeff2 = coeff2,
                                         mc.cores = parallel::detectCores()))
data.frame(cbind(ARL = mean(RLS), ARLSE = sd(RLS)/sqrt(m)))

## Parallel simulation (PSOCK)
no_cores <- parallel::detectCores()
cl <- parallel::makeCluster(no_cores)
side <- "low"
h_vec <- h
QS_vec <- 1
k <- kopt
parallel::clusterExport(cl, c("h_vec", "eocusum_arloc_sim", "df1", "coeff1", "coeff2",
                              "QS_vec", "side", "k"))
time <- system.time( {
  RLS <- array(NA, dim = c( length(QS_vec), length(h_vec), m))
  for (h in h_vec) {
    for (QS in QS_vec) {
      cat(h, " ", QS, "\n")
      RLS[which(QS_vec==QS), which(h==h_vec), ] <- parallel::parSapply(cl, 1:m, eocusum_arloc_sim,
                                                                       side = side, QS = QS, h = h,
                                                                       k = k, df = df1,
                                                                       coeff = coeff1,
                                                                       coeff2 = coeff2,
                                                                       USE.NAMES = FALSE)
    }
  }
} )
ARL <- apply(RLS, c(1, 2), mean)
ARLSE <- sqrt(apply(RLS, c(1, 2), var)/m)
print(list(ARL, ARLSE, time))
parallel::stopCluster(cl)
# }