In meyera/rcusum: Tools for the easy creation cusum analysis charts

source("functions.R")
# pandoc Manual.latex --output "Manual RCusum.pdf" --smart --standalone --number-sections

RCusum - An Implementation of Cusum Control Charts as an R package

Author: Alex Meyer
Date: r format(Sys.time(), "%d.%m.%Y")
Version: 1.00

This R Package implements the charts as described by Rogers et. al 2004. Cusum control charts are devices for monitoring arbitrary processes for a prespecified event. This package and manual uses the term failure instead of event, although any other type of event such as successes could be monitored as well. The processes of our interest are typically surgical procedures.

To create a basic control chart 4 parameters have to be specified:

$p0$: risk of failure when the process is in control (acceptable failure rate for the procedure)
$p1$, where $p1>p0$: failure rate considered unacceptable
$\alpha$: Type I error (the probability of concluding that the failure rate has increased, when in fact it has not)
$\beta$: Type II error (the probability of concluding that the failure rate has not increased, when in fact it has)

Anatomy of a Cusum Control Chart

A cusum chart shows the cumulative sum of failures for a given process. There are two boundaries shown on the graph.

If the cusum curve...

... lies within the boundaries, there is not enough evidence for a significant conclusion, thus there is still need for monitoring!
... crossed the upper boundary, we can conclude that the failure rate has increased to the unacceptable event rate $p1$ (at the specified significance level).
... crossed the lower boundary, we can conclude that the failure rate is less than or equal to the acceptable event rate $p0$ (at the specified significance level)

"The natural progression of the graph for an individual or institution with acceptable performance is toward the lower boundary for this method of constructing control limits." (Rogers et. al, 2004)

Unadjusted Cumulative Failure Charts

For the following examples a data set is simulated by binomial sampling with three specified probabilties. The simulated dataset contains data for three different surgeons.

Unadjusted Cumulative Failure Chart

The parameters are specified as:

$p0 = .1$
$p1 = .2$
$\alpha = .01$
$\beta = .01$

cusum_plot = cusum(
              failure_indicator=df$is_failure,
              p0=.10,
              p1=.20,
              alpha=0.01,
              beta=0.01,
              by=df$by
              )
print(cusum_plot)

Unadjusted Cumulative log-likelihood ratio test charts

"An alternative but equivalent presentation of the data involves graphing a modified CUSUM against the operation number. Interpretation of the graph in relation to the boundary lines is the same as for the cumulative failures chart. If performance is acceptable, the graph will tend downward toward the lower boundary; it will not follow the horizontal axis." (Rogers et. al, 2004)

For the example below, the dataset and parameter specification is the same as aboved.

cusum_plot = cusum(
              failure_indicator=df$is_failure,
              p0=.10,
              p1=.20,
              alpha=0.01,
              beta=0.01,
              by=df$by,
              loglike_chart=TRUE
              )
print(cusum_plot)

Cumulative observed minus expected failure chart

With the previously mentioned graphs, the expected failure rate will naturally progress towards the lower boundary, which is somewhat counter-intuitive.

By using the Cumulative observed minus expected failure graph, the graph will roughly oscillate around the horizontal axis if the performance does not deviate from the acceptable rate of failures. If any major trend is visible, most likely significant deviance is observed.

This graph does not show any inferential boundaries, so it should be interpreted only in conjunction with the already mentioned graph types. Be careful to not overinterpret this graph.

This graph requires only $p0$ specification. Below the same dataset as in the examples before is used, $p0$ was specified as before: $p0 = .1$

cusum_plot = cusum.obs_minus_exp(
              failure_indicator=df$is_failure,
              p0=.10,
              by=df$by
              )
print(cusum_plot)

Risk-adjusted Cumulative Failure Charts

So far we have discussed only unadjusted cumulative failure charts. So basically we assumed the acceptable and the unacceptbale failure rates are the same among all patients. In the following we adjust for the case-mix by supplying an individual acceptable failure rate for each patient calculated by using an empirically derived risk model or by using an accepted risk model such as the ESII or STS-Score.

When using risk adjusted calculation, we can control for confounders and risk determining factors to facilitate comparison.

However, Rogers et. al (2004) mention "Monitoring a health-care process is not the same as monitoring a manufacturing process; case mix represents a fundamental difference, and risk adjustment is imperfect and cannot remove all confounding.".

Comparison of unadjusted und risk adjusted Cumulative observed minus expected failure charts

The upper graph is plotted using a fixed $p0 = .10$ acceptable failure rate. The lower is the same plot, only using a simulated failure rate $p0_i \sim N(\mu=0.10, \sigma=0.03)$ for each individual procedure. The failure rate is simulated using a probability of failure $P(F)=0.1$

cp1 = cusum.obs_minus_exp(rbinom(200,1,0.10), 0.10) + ggtitle("Unadjusted Plot - p0=0.10")
cp2 = cusum.obs_minus_exp(rbinom(200,1,0.10), c(rnorm(100, 0.10, 0.03),rnorm(100, 0.10, 0.03))) + ggtitle("Risk-adjusted Plot - p0i ~ N(µ=0.10, sd=0.03)")

library(gridExtra)
grid.arrange(cp1,cp2)

The next set of plots show a simulated set of operations with the first 100 cases having a simulated failure rate of 10 %, the 100 consecutive cases however having a raising rate up to 20 %.

The upper graph is plotted again using a fixed $p0 = .10$ acceptable failure rate for all cases.
The lower plot however accounts for a raising case complexity using a simulated failure rate $p0_i \sim N(\mu=0.10, \sigma=0.03)$ for the first 100 patients, and then $p0_i \sim N(\mu=0.20, \sigma=0.03)$ for each individual procedure.

cp1 = cusum.obs_minus_exp(c(rbinom(100,1,0.10), rbinom(100,1,0.20)), 0.10, scale_ylim=30) + ggtitle("Unadjusted Plot")
cp2 = cusum.obs_minus_exp(c(rbinom(100,1,0.10), rbinom(100,1,0.20)), c(rnorm(100, 0.10, 0.03),rnorm(100, 0.20, 0.03)), scale_ylim=30) + ggtitle("Risk-adjusted Plot")

library(gridExtra)
grid.arrange(cp1,cp2)

Risk adjusted Cumulative log-likelihood ratio test charts (SPRT charts)

The risk adjusted counterparts of Cumulative log-likelihood ratio test charts are the risk-adjusted sequential probability ratio test (SPRT) charts. They also come with control limits as described before, which facilitate formal hypothesis testing.

For the unadjusted chart, increase in risk is defined in terms of reaching a constant unacceptable failure rate.
However, when risk for each patient varies, it does not make sense to have a common unacceptable rate applied across all operations. A variable unacceptable rate is achieved by defining the increase in terms of a relative risk (ie, odds ratio), rather than a specific rate: an odds ratio of 2, for example, would equate approximately to a doubling of patientspecific risk of failure, an odds ratio of 1.5 to a 50% increase in failure risk, and so on.

As before: "The natural progression of the risk-adjusted graph for an individual or institution with acceptable performance is toward the lower boundary." (Rogers et al., 2004)

So the following parameters need to be given to construct an risk-adjusted cumulative log-likelihood ratio test chart:

$p0_i$: the procedure specific risk-adjusted probability of failure
$OR$: the Odds Ratio --- a relative increase of the patient specific risk considered to be unacceptable
$\alpha$: Type I error (the probability of concluding that the failure rate has increased, when in fact it has not)
$\beta$: Type II error (the probability of concluding that the failure rate has not increased, when in fact it has)

Below two graphs are shown. The upper is the same unadjusted cumulative log-likelihood chart as prevousily shown:

Surgeon A crosses the lower boundary, so we can conclude that there is enough evidence that his failure rate is equal or below the acceptable failure rate.
Surgeon B almost crosses the lower boundary, so he is still in the monitoring area and formally no conclusion is yet possible.
Surgeon C crossed the upper limit. We conclude that evidence suggests that his failure rate reached the unacceptable rate.

However, with this chart, we assume all procedures performed had the same risk.

The lower plot shows the same data risk-adjusted. Procedure specific risk was simulated to distribute randomly around 10 % with a standard deviation of 3 %. For Surgeon C, who reached the unacceptable failure rate before, the procedure specific risk, however, was increasingly rising in a random fashion up to 20 %. This change is reflected in the risk adjusted chart.

cp1 = cusum(
              failure_indicator=df$is_failure,
              p0=.10,
              p1=.20,
              alpha=0.01,
              beta=0.01,
              by=df$by,
              loglike_chart=TRUE
              )
cp2 = cusum.sprt(df$is_failure, df$p0, 1.5, by=df$by)

library(gridExtra)
grid.arrange(cp1 + ggtitle("Unadjusted Cusum log-likelihood"), cp2 + ggtitle("Risk-adjusted Cusum log-likelihood"))