vlad (version 0.2.0)

racusum_scores: Compute CUSUM scores based on the log-likelihood ratio statistic

Description

Compute CUSUM scores based on the log-likelihood ratio statistic.

Usage

racusum_scores(wt1, wt2, reset = FALSE, h1 = NULL, h2 = NULL)

Arguments

wt1

Double. Log-likelihood ratio scores from function llr_score for upper CUSUM.

wt2

Double. Log-likelihood ratio scores from function llr_score for lower CUSUM.

reset

Logical. If FALSE CUSUM statistic is not reset. If TRUE CUSUM statistic is reset to 0 after a signal is issued.

h1

Double. Upper control limit of the CUSUM chart.

h2

Double. Lower control limit of the CUSUM chart.

Value

Returns a list with two components for the CUSUM scores.

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.

Parsonnet V, Dean D, Bernstein AD (1989). A method of uniform stratification of risk for evaluating the results of surgery in acquired adult heart disease. Circulation, 79(6):I3-12.

Rigdon SE and Fricker RD (2015). Health Surveillance. In Chen DG and Wilson J (eds) Innovative Statistical Methods for Public Health Data, pp. 203--249. Springer, Cham.

Examples

# NOT RUN {
#' library(vlad)
# patient Cusum values with different odds ratios, see Rigdon and Fricker p. 225, 226
coeff <- c("(Intercept)" = -3.67, "Parsonnet" = 0.077)
wt1 <- round(llr_score(df = data.frame(19L, 0), coeff = coeff, R0 = 1, RA = 2), 4)
wt2 <- round(llr_score(df = data.frame(19L, 0), coeff = coeff, R0 = 1, RA = 1/2), 5)
all.equal(racusum_scores(wt1 = wt1, wt2 = wt2), list(s1 = 0, s1l = 0.05083))

library("dplyr")
library("tidyr")
library(ggplot2)
data("cardiacsurgery", package = "spcadjust")

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

## subset phase I (In-control)
SI <- filter(SALL, phase == "I") %>% select(s, y)

## estimate coefficients from logit model
coeff1 <- round(coef(glm(y ~ s, data = SI, family = "binomial")), 3)

## subset phase II of surgeons 2
S2II <- filter(SALL, phase == "II", surgeon == 2) %>% select(s, y)
n <- nrow(S2II)

## CUSUM statistic without reset
wt1 <- sapply(1:n, function(i) llr_score(S2II[i, c("s", "y")], coeff = coeff, RA = 2))
wt2 <- sapply(1:n, function(i) llr_score(S2II[i, c("s", "y")], coeff = coeff, RA = 1/2))
cv <- racusum_scores(wt1 = wt1, wt2 = wt2)
s1 <- cv$s1; s1l <- cv$s1l
dm1 <- data.frame(cbind("n" = 1:length(s1), "Cup" = s1, "Clow" = -s1l, "h1" = 2, "h2" = -2))

## CUSUM statistic reset after signal
cv <- racusum_scores(wt1 = wt1, wt2 = wt2, reset = TRUE, h1 = 2, h2 = 2)
s1 <- cv$s1; s1l <- cv$s1l
dm2 <- data.frame(cbind("n" = 1:length(s1), "Cup" = s1, "Clow" = -s1l, "h1" = 2, "h2" = -2))

## plot
dm3 <- bind_rows(dm1, dm2, .id = "type")
dm3$type <- recode_factor(dm3$type, `1`="No resetting", `2`="Resetting")
dm3 %>%
  gather("CUSUM", value, c(-n, - type)) %>%
  ggplot(aes(x = n, y = value, colour = CUSUM, group = CUSUM)) +
  geom_hline(yintercept = 0, colour = "darkgreen", linetype = "dashed") +
  geom_line(size = 0.5) +
  facet_wrap( ~ type, ncol = 1, scales = "free") +
  labs(x = "Patient number n", y = "CUSUM values") + theme_classic() +
  scale_y_continuous(sec.axis = dup_axis(name = NULL, labels = NULL)) +
  scale_x_continuous(sec.axis = dup_axis(name = NULL, labels = NULL)) +
  guides(colour = "none") +
  scale_color_manual(values = c("blue", "orange", "red", "red"))
# }