Description Usage Arguments Value Author(s) References Examples
Compute ARLs of EO-CUSUM control charts using simulation.
1 | eocusum_arl_sim(r, pmix, k, h, RQ = 1, yemp = FALSE, side = "low")
|
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 |
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
|
yemp |
Logical. If |
side |
Character. Default is |
Returns a single value which is the Run Length.
Philipp Wittenberg
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 | ## 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)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.