vlad (version 0.1.0)

racusum_arl_sim: Compute ARLs of RA-CUSUM control charts using simulation

Description

Compute ARLs of RA-CUSUM control charts using simulation.

Usage

racusum_arl_sim(r, coeff, h, df, R0 = 1, RA = 2, yemp = TRUE)

Arguments

r

Integer vector. Number of runs.

coeff

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

h

double. Control Chart limit for detecting deterioration/improvement.

df

DataFrame. First column are Parsonnet Score values within a range of zero to 100 representing the preoperative patient risk. The second column are binary (0/1) outcome values of each operation.

R0

double. Odds ratio of death under the null hypotheses.

RA

double. Odds ratio of death under the alternative hypotheses. Detecting deterioration in performance with increased mortality risk by doubling the odds Ratio RA=2. Detecting improvement in performance with decreased mortality risk by halving the odds ratio of death RA=1/2.

yemp

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

Value

Returns a single value which is the Run Length.

References

Steiner SH, Cook RJ, Farewell VT and Treasure T (2000). <U+201C>Monitoring surgical performance using risk-adjusted cumulative sum charts.<U+201D> Biostatistics, 1(4), pp. 441-452. doi: 10.1093/biostatistics/1.4.441.

Examples

# NOT RUN {
library("vlad")
library("spcadjust")
set.seed(1234)
data("cardiacsurgery")
df1 <- subset(cardiacsurgery, select=c(Parsonnet, status))
coeff1 <- round(coef(glm(status ~ Parsonnet, data=df1, family="binomial")), 3)

## Parallel Simulation 1: y = random (10^4 runs, RA=2)
m <- 10^4; h_vec <- 2.7; yemp <- FALSE
no_cores <- parallel::detectCores()
cl <- parallel::makeCluster(no_cores)
parallel::clusterExport(cl, c("h_vec", "racusum_arl_sim", "coeff1", "df1", "yemp"))
time <- system.time( {
  ARL <- array(NA, dim=c( length(h_vec), m))
  for (h in h_vec) {
    ARL[which(h_vec==h), ] <- parallel::parSapply(cl, 1:m, racusum_arl_sim, h=h, coeff=coeff1,
                                                 df=df1, yemp=yemp, USE.NAMES=FALSE) }
} )
simMean <- apply(ARL, c(1), mean)
simSE <- sqrt(apply(ARL, c(1), var)/m)
print(list(simMean, simSE, time))
parallel::stopCluster(cl)
df.sim1 <- data.frame("RA"=2, "h"=h, "ARL"=simMean, "ARLSE"=simSE, "nsim"=m)

## Parallel Simulation 2: y = empirical (10^4 runs, RA=2)
m <- 10^4; h_vec <- 2.7
no_cores <- parallel::detectCores()
cl <- parallel::makeCluster(no_cores)
parallel::clusterExport(cl, c("h_vec", "racusum_arl_sim", "coeff1", "df1"))
time <- system.time( {
  ARL <- array(NA, dim=c( length(h_vec), m))
  for (h in h_vec) {
    ARL[which(h_vec==h), ] <- parallel::parSapply(cl, 1:m, racusum_arl_sim, h=h, coeff=coeff1,
                                                 df=df1, USE.NAMES=FALSE) }
} )
simMean <- apply(ARL, c(1), mean)
simSE <- sqrt(apply(ARL, c(1), var)/m)
print(list(simMean, simSE, time))
parallel::stopCluster(cl)
df.sim2 <- data.frame("RA"=2, "h"=h, "ARL"=simMean, "ARLSE"=simSE, "nsim"=m)

rbind(df.sim1, df.sim2)
# }