vlad (version 0.1.0)

racusum_arloc_h_sim: Compute alarm threshold (Out of Control ARL) of RA-CUSUM control charts using simulation

Description

Compute alarm threshold (Out of Control ARL) of RA-CUSUM control charts using simulation.

Usage

racusum_arloc_h_sim(L0, df, coeff, coeff2, R0 = 1, RA = 2, RQ = 1,
  m = 100, nc = 1, verbose = FALSE)

Arguments

L0

double. Prespecified in-control Average Run Length.

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.

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.

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.

RQ

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

m

integer. Number of simulation runs.

nc

integer. Number of cores.

verbose

boolean. If TRUE verbose output is included, if FALSE a quiet calculation of h is done.

Value

Returns a single value which is the control limit h for a given in-control ARL.

Details

The function racusum_arloc_h_sim determines the control limit for given in-control ARL (L0) by applying a multi-stage search procedure which includes secant rule and the parallel version of racusum_arloc_sim using mclapply.

Examples

# NOT RUN {
require("vlad")
# Set seed for reproducibility
RNGkind("L'Ecuyer-CMRG")
set.seed(1234)
parallel::mc.reset.stream()
# 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

# Deterioration:
# 1. Determine critical value for given ARL
racusum_arloc_h_sim(L0=370, df=df1, coeff=coeff1, coeff2=coeff2, m=m, RA=2, nc=6)
# h=2.030933

# 2. Determine ARL and Standard Error
RLS <- do.call(c, parallel::mclapply(1:m, racusum_arloc_sim, h=2.035, df=df1, RA=2, coeff=coeff1,
                                     coeff2=coeff2, mc.cores=6))
data.frame(cbind("ARL"=mean(RLS), "ARLSE"=sd(RLS)/sqrt(m)))
# ARL=371.125; ARLSE=11.36053

# Improvement:
# 1. Determine critical value for given ARL
racusum_arloc_h_sim(L0=370, df=df1, coeff=coeff1, coeff2=coeff2, m=m, RA=1/2, nc=6)
# h=1.710999
#
# 2. Determine ARL and Standard Error
RLS <- do.call(c, parallel::mclapply(1:m, racusum_arloc_sim, h=1.760, df=df1, RA=1/2, coeff=coeff1,
                                     coeff2=coeff2, mc.cores=6))
data.frame(cbind("ARL"=mean(RLS), "ARLSE"=sd(RLS)/sqrt(m)))
# ARL=399.613; ARLSE=10.7601
# }