racusum_arl: ARL of RA-CUSUM charts

In vlad: Variable Life Adjusted Display and Other Risk-Adjusted Quality Control Charts

Description

Compute the ARL of risk-adjusted CUSUM charts.

Usage

racusum_arl_mc(h, pmix, RA, RQ, scaling = 600, rounding = "p", method = "Toep")

racusum_arl_sim(h, pmix, r, RA = 2, RQ = 1, yemp = FALSE)

Arguments

h	Double. h is the control limit (>0).
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.
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. Odds ratio of death under the null hypotheses is 1. RQ. Use RQ = 1 to compute the in-control ARL and other values to compute the out-of-control ARL.
RQ	Double. Defines the true performance of a surgeon with the odds ratio ratio of death RQ. Use RQ = 1 to compute the in-control ARL and other values to compute the out-of-control ARL.
scaling	Double. The scaling parameter controls the quality of the approximation, larger values achieve higher accuracy but increase the computation burden (larger transition probability matrix).
rounding	Character. If rounding = "p" a paired rounding implementation of Knoth et al. (2019) is used, if rounding = "s" a simple rounding method of Steiner et al. (2000) is used.
method	Character. If method = "Toep" a combination of Sequential Probability Ratio Test and Toeplitz matrix structure is used to calculate the ARL. "ToepInv" computes the inverted matrix using Toeplitz matrix structure. "BE" solves a linear equation system using the classical approach of Brook and Evans (1972) to calculate the ARL.
r	Integer. Number of runs.
yemp	Logical. If TRUE use observed outcome value, if FALSE use estimated binary logistc regression model.

Value

Returns a single value which is the Average Run Length for "racusum_arl_mc" and the Run Length for "racusum_arl_sim".

Author(s)

Philipp Wittenberg

References

Steiner SH, Cook RJ, Farewell VT and Treasure T (2000). Monitoring surgical performance using risk-adjusted cumulative sum charts. Biostatistics, 1(4), pp. 441–452.

Knoth S, Wittenberg P and Gan FF (2019). Risk-adjusted CUSUM charts under model error. Statistics in Medicine, 38(12), pp. 2206–2218.

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.

Brook D and Evans DA (1972) An approach to the probability distribution of CUSUM run length. Biometrika, 59(3), pp. 539–549

Webster RA and Pettitt AN (2007) Stability of approximations of average run length of risk-adjusted CUSUM schemes using the Markov approach: comparing two methods of calculating transition probabilities. Communications in Statistics - Simulation and Computation 36(3), pp. 471–482

Examples

## Not run: 
library(vlad)
library(dplyr)
data("cardiacsurgery", package = "spcadjust")

## Markov Chain
## preprocess data to 30 day mortality and subset phase I (In-control) of surgeons 2
SALLI <- cardiacsurgery %>% rename(s = Parsonnet) %>%
  mutate(y = ifelse(status == 1 & time <= 30, 1, 0),
        phase = factor(ifelse(date < 2*365, "I", "II"))) %>%
  filter(phase == "I") %>% select(s, y)

## estimate risk model, get relative frequences and probabilities
mod1 <- glm(y ~ s, data = SALLI, family = "binomial")
fi  <- as.numeric(table(SALLI$s) / length(SALLI$s))
usi <- sort(unique(SALLI$s))
pi1 <- predict(mod1, newdata = data.frame(s = usi), type = "response")
pi2 <- tapply(SALLI$y, SALLI$s, mean)

## set up patient mix (risk model)
pmix1  <- data.frame(fi, pi1, pi1)

## Average Run Length for detecting deterioration RA = 2:
racusum_arl_mc(pmix = pmix1, RA = 2, RQ = 1, h = 4.5)

## Average Run Length for detecting improvement RA = 1/2:
racusum_arl_mc(pmix = pmix1, RA = 1/2, RQ = 1, h = 4)

## set up patient mix (model free)
pmix2  <- data.frame(fi, pi1, pi2)

## Average Run Length for detecting deterioration RA = 2:
racusum_arl_mc(pmix = pmix2, RA = 2, RQ = 1, h = 4.5)

## Average Run Length for detecting improvement RA = 1/2:
racusum_arl_mc(pmix = pmix2, RA = 1/2, RQ = 1, h = 4)

## compare results with R-code function 'findarl()' from Steiner et al. (2000)
source("https://bit.ly/2KC0SYD")
all.equal(findarl(pmix = pmix1, R1 = 2, R = 1, CL = 4.5, scaling = 600),
         racusum_arl_mc(pmix = pmix1, RA = 2, RQ = 1, h = 4.5, scaling = 600, rounding = "s"))

## Monte Carlo simulation
set.seed(1234)
SALLI <- cardiacsurgery %>% mutate(s = Parsonnet) %>%
  mutate(y = ifelse(status == 1 & time <= 30, 1, 0),
         phase = factor(ifelse(date < 2*365, "I", "II"))) %>%
  filter(phase == "I") %>% select(s, y)

## estimate risk model, get relative frequences and probabilities
mod1 <- glm(y ~ s, data = SALLI, family = "binomial")
y <- SALLI$y
pi1 <- fitted.values(mod1)

## set up patient mix (risk model)
pmix <- data.frame(y, pi1, pi1)
h <- 2.75599

m <- 1e4
RLS <- sapply(1:m, racusum_arl_sim, h=h, pmix=pmix, RA=2)
data.frame(cbind(ARL=mean(RLS), ARLSE=sd(RLS)/sqrt(m), h, m))

## End(Not run)

vlad documentation built on Feb. 15, 2021, 5:12 p.m.