ProfMohammed commented 1 month ago

Hi Jacob -

Below is one way to provide an overview of the family of SPC charts

• run charts • Shewhart control charts
• cumulative sum (CUSUM) charts

I don't think qicharts2 produces CUSUM plots? Is that correct? In any case my question is how should the book address the issue of CUSUM plots - they are widely used to monitor binary outcomes in surgery

anhoej commented 1 month ago

Check bchart().

ProfMohammed commented 1 month ago

That's great - binary cusum are likely to be very important. I was thinking about cusum for continious data -

here is the basic code with bp data -

CUSUM plots - v mask (no longer used) and tabular (popular)

https://www.itl.nist.gov/div898/handbook/pmc/section3/pmc323.htm

library(tidyverse)

sdf <- data.frame( systolic = c(169, 172, 175, 174, 161, 142, 174, 171, 168, 174, 180, 194, 161, 181, 175, 176, 186, 166, 157, 183, 177, 171, 185, 176, 181, 174), day = 1:26 )

use these names as functions below use y and id

sdf$y <- sdf$systolic sdf$id <- sdf$day

use df as functions below use df

df <- sdf

df$target <- mean(df$y)

sensible choices suggested by literature (eg HCDG Provost)

k <- 0.5 # 1/2 sigma slack, so shifts above +/- 0.5 sigma (or 1 sigma in total) h <- 4 # constant to make CUSUM chart equivalent to 3 sigma when k=0.5

N <- dim(df)[1]

df$deviation <- df$y-df$target df$cusum <- cumsum(df$deviation)

df <- df %>% mutate(mr = abs(y - lag(y)))

mean_mr <- mean(df$mr[!is.na(df$mr)]) sigma_est <- mean_mr/1.128

df$Vupper <- NA df$Vlower <- NA

df$Vupper[N] <- df$cusum[N] + hsigma_est df$Vlower[N] <- df$cusum[N] - hsigma_est

for (i in (N-1):1) { print(i) df$Vupper[i] <- df$Vupper[i+1] + ksigma_est df$Vlower[i] <- df$Vlower[i+1] - ksigma_est }

df

V mask

par(mfcol=c(1,2)) plot(cusum~day, df, type='b', ylim=c(min(df$Vlower), max(df$Vupper))) abline(h=0) points(Vupper~id, df, type='l') points(Vlower~id, df, type='l')

tabular CUSUM

target <- mean(df$y) deviation <- df$y - target cusum <- cumsum(deviation) mr <- abs(df$y - lag(df$y))

mean_mr <- mean(mr, na.rm=T) sigma_est <- mean_mr/1.128

C_plus <- NULL C_minus <- NULL for (i in 1:N) { if(i==1) { C_plus[1] <- 0 C_minus[1] <- 0 } else { C_plus[i] <- max(0, df$y[i]-(target+ksigma_est) + C_plus[i-1]) C_minus[i] <- min(0, (ksigma_est-target) + df$y[i] + C_minus[i-1]) } }

Ucl <- +hsigma_est Lcl <- -hsigma_est

plot(df$id, C_plus, type='b', pch=1,ylim=c(min(df$Vlower), max(df$Vupper))) points(df$id, C_minus, type='b', pch=2) abline(h=c(0,Lcl,Ucl), col='blue')

would it be possible to add the tabular cusum to qihcarts2?