Content Disclaimer
Copyright @2020.
All Rights Reserved.

StatsToDo: CUSUM for Proportions with Binomial Distribution

Links : Home Index (Subjects) Contact StatsToDo


Explanations R Code
The example is a made up one to demonstrate the numerical process, and the data is generated by the computer. It purports to be from a quality control exercise in an obstetric unit in a large public hospital which delivers about 100 babies a week, and concerns the Caesarean Section Rate.
  • From records in the past, we established the benchmark Caesarean Section Rate to be 20% (0.2), and this can be capped if the junior staff and midwives are well trained and closely supervised.
  • With time however, experienced staff leave and replaced by the less experienced and trained. The standard of supervision would gradually deteriorate, resulting in an increase in the Caesarean Section rate.
  • We would like to trigger an alarm and reorganize the working and supervision framework when the Caesarean Section Rate increases to 30% (0.3) or more.
  • Given that Caesarean Section are performed by schedule and in emergency, it is more frequent during working hours, and less frequent at night and on weekends. It was decided therefore that the Caesarean Section Rate should be evaluated in groups of 100 births (ssiz=100), approximately the number of births per week.
  • As re-organizing working framework is time consuming and disruptive, we would like any false alarm to be no more frequent than once every 100 weeks, so the average run length ARL = 100 (or 10,000 births)

Step 1: Parameters for calculating Reference Value (k) and Decision Interval (h)

# Step 1: parameters and data
inControlProp = 0.2      # 20%
outOfControlProp = 0.3   # 30%
ssiz = 100               # sample size 
arl = 100 
theModel = "F"   #F for FIR,   Z for zero, S for steady state
Step 1 contains the parameters.

The first 4 lines sets the parameters required for the analysis.

The 5th line, the model has 3 options, which sets the first value of the CUSUM

  • F means Fast Initial Response, where the initial CUSUM value is set at half of the Decision Interval h. The rationale is that, if the situation is in control then CUSUM will gradually drift towards zero, but if the situation is already out of control, an alarm would be triggered early. The down side is that a false alarm is slightly more likely early on in the monitoring. As FIR is recommended by Hawkins, it is set as the default option
  • Z is for zero, and CUSUM starts at the baseline value of 0. This will lower the risk of false alarm in the early stages of monitoring, but will detect the out of control situation slower if it already exists at the begining.
  • S is for steady state, intended for when monitoring is already ongoing, and a new plot is being constructed. The CUSUM starts at the value when the previous chart ends.
  • Each model will make minor changes to the value of the decision interval h. The setting of the initial values is mostly intended to determine how quickly an alarm can be triggered if the out of control situation exists from the beginning.

Step 2: Calculate Reference Value (k) and Decision Interval (h)

# Step 2: Estimating k and h
#install.packages("CUSUMdesign")   # if not already installed
library(CUSUMdesign)
result <- getH(distr=4, ICprob=inControlProp, OOCprob=outOfControlProp, samp.size=ssiz, ARL=arl, type=theModel)
k = result$ref
h = result$DI
if(outOfControlProp<inControlProp)
{
  h = -h
}
cat("Reference Value k=",k,"\tDecision Interval h=", h, "\n")
Step 2 performs that statistical calculations using the parameters entered. The package CUSUMdesign needs to be alrady installed, and the library activated each time the program is used.

result is the object that contains the results of the analysis. The result required for this program are the reference value (k) and decision interval h. Please note that h is calculated as a positive value. If the CUSUM is designed to detect a decrease from in control value, then h needs to be changed to a negative value.

The last line displays the results we need 

Reference Value k= 24.75 	Decision Interval h= 5.5
Please note that, although the parameters in and out of control are entered as proportions, k and h are related to the number of positives.

Step 3: CUSUM Plot

Step 3 is divided into 2 parts. Step 3a calculates the cusum vector, and 3b plots the vector and h in a graph.

Step 3a: create the CUSUM vector 

# Step 3a: Calculating CUSUM
dat=c(20,26,18,6,25,24,18,26,22,21,22,26,25,27,25,28,26,24,23,25)    # number of positives in each sample
cusum <- vector()
cusumValue = 0
if(theModel=="F")
{
  cusumValue = h / 2
}
for(i in 1 : length(dat))
{
  cusumValue = cusumValue + dat[i] - k
  if(outOfControlProp>inControlProp) # Up
  {
    if(cusumValue<0)
    {
      cusumValue = 0
    }
  }
  else                               # down
  {
    if(cusumValue>0)
    {
      cusumValue = 0
    }
  }
  cusum[i] = cusumValue
}
cusum
The vector dat contains the number of positives in each sampling (in this case the number of Caesarean Sections in each 100 consecutive births).

The next 6 lines of code in step 3a creates the empty cusum vector and sets the initial cusum value. The remaining codes calculates the cusum value for each measurement, and places it in the cusum vector

The result CUSUM vector represents the number of excessive positives above 20 (20% in 100 births), and is as follows 

> cusum
 [1] 0.00 1.25 0.00 0.00 0.25 0.00 0.00 1.25 0.00 0.00 0.00 1.25 1.50 3.75 4.00 7.25 8.50
[18] 7.75 6.00 6.25
Step 3b: Plotting the CUSUM vector 
# Step 3b: Plot the cusum vector and h
plot(cusum,type="l")
abline(h=h)
In step 3b, the first line plots the cusum vector, and the second line the decision interval h. The result plot is shown to the right.

Step 4: Optional export of results

# Step 4: Optional export of results
#myDataFrame <- data.frame(dat,cusum)       #combine dat and cusum to dataframe
#myDataFrame                                #display dataframe
#write.csv(myDataFrame, "CusumBinom.csv")   # write dataframe to .csv file
Step 4 is optional, and in fact commented out, and included as a template only. Each line can be activated by removing the #

The first line places the two vectors, dat and cusum together into a dataframe

The second line displays the data, along with row numbers, in the console, which can then be copied and pasted into other applications for further processing

The third line saves the dataframe as a comma delimited .csv file. This is needed if the data is too large to handle by copy and paste from the console.

Javascript Plot