vlad (version 0.1.0)

eocusum_adoc_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_adoc_sim(r, k, h, df, coeff, coeff2, QS = 1, side = "low",
  type = "cond", m = 50)

Arguments

r

int. Number of of simulation runs.

k

double. Reference value of the CUSUM control chart.

h

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

df

DataFrame. First column Parsonnet Score and second column outcome of each operation.

coeff

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

coeff2

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

QS

double. Defines the performance of a surgeon with the odds ratio ratio of death Q.

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

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.

Taylor HM (1968). <U+201C>The Economic Design of Cumulative Sum Control Charts.<U+201D> Technometrics, 10(3), pp. 479-488. doi: 10.1080/00401706.1968.10490595.

Crosier R (1986). <U+201C>A new two-sided cumulative quality control scheme.<U+201D> Technometrics, 28(3), pp. 187-194. doi: 10.2307/1269074.

Examples

# NOT RUN {
library("vlad")
library("spcadjust")
## Datasets
data("cardiacsurgery")
s5000 <- dplyr::sample_n(cardiacsurgery, size=5000, replace=TRUE)
df1 <- subset(cardiacsurgery, select=c(Parsonnet, status))
df2 <- subset(s5000, select=c(Parsonnet, status))
## estimate coefficients from logit model
coeff1 <- round(coef(glm(status~Parsonnet, data=df1, family="binomial")), 3)
coeff2 <- round(coef(glm(status~Parsonnet, 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, parsonnetscores=df1$Parsonnet, coeff=coeff1)
# 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_adoc_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()
# }