Chapter 12 SPC Charts for Rare Events

When dealing with potentially serious or fatal events, the number of occurrences is often, fortunately, very low. This creates challenges for traditional SPC charts designed for count data.

The difficulty arises because charts such as the P and U charts assume a sufficiently frequent and reasonably stable occurrence of events. When counts are low, the data become sparse and highly variable, which makes these charts less useful.

Sparse data can produce lower control limits that are censored at zero, thereby invalidating the 3-sigma test in that direction. Run charts may also become problematic when more than half of the data points are zero, because the median then becomes zero as well. This undermines the runs analysis by producing one long run above the median.

One way of dealing with rare events is to plot the number of opportunities or the time between events rather than event proportions or rates. In other words, we turn the indicator around and focus on the gaps between occurrences rather than on the occurrences themselves.

Another approach is to examine the cumulative sums (CUSUM) of binary data.

In this chapter, we introduce the G chart for the number of opportunities between events, the T chart for time between events, and the Bernoulli CUSUM chart for binary data.

12.1 Introducing the birth dataset

The Robson group 1 births dataset is a data frame containing 2193 observations from Robson group 1 deliveries, that is, first pregnancy, single baby, head first, and gestational age of at least 37 weeks.

# import data
births <- read.csv('data/robson1_births.csv',
                   comment.char = '#',
                   colClasses   = c('POSIXct',
                                    'Date',
                                    'logical',
                                    'logical',
                                    'factor',
                                    'integer',
                                    'integer',
                                    'integer',
                                    'double',
                                    'logical'))

# make sure the rows are sorted in time order
births <- births[order(births$datetime), ]
# add dummy variable, case, for counting the denominator for the proportion charts
births$case <- 1
# show data structure
str(births)
## 'data.frame':    2193 obs. of  11 variables:
##  $ datetime: POSIXct, format: "2016-01-04 06:13:00" "2016-01-04 07:20:00" ...
##  $ biweek  : Date, format: "2016-01-04" "2016-01-04" ...
##  $ csect   : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ cup     : logi  FALSE FALSE FALSE FALSE TRUE TRUE ...
##  $ sex     : Factor w/ 2 levels "female","male": 2 1 2 2 2 2 1 2 2 2 ...
##  $ length  : int  55 52 50 47 55 50 50 47 50 53 ...
##  $ weight  : int  3872 3750 3458 3145 4224 2730 3130 2755 3270 3645 ...
##  $ apgar   : int  10 10 10 10 10 10 10 10 10 10 ...
##  $ ph      : num  7.1 7.26 7.36 7.31 7.16 7.14 7.23 7.21 7.17 7.2 ...
##  $ asphyxia: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ case    : num  1 1 1 1 1 1 1 1 1 1 ...

In this chapter, we are interested in the asphyxia variable, which is a logical vector that is TRUE when the baby had either a low Apgar score or a low umbilical cord pH, suggesting lack of oxygen during delivery. Asphyxia is a very serious condition that may result in permanent brain damage or death. In total, there are 16 cases corresponding to 0.7%.

Figure 12.1 is a run chart of the biweekly counts.

# run chart of counts
qic(biweek, asphyxia, 
    data    = births, 
    agg.fun = 'sum')  # use sum() function to aggregate data by subgroup
## Subgroup size > 1. Data have been aggregated using sum().
Run chart of number of deliveries with neonatal asphyxia

Figure 12.1: Run chart of number of deliveries with neonatal asphyxia

Because more than half of the data points are zero, the centre line, that is, the median, is also zero. This invalidates the runs analysis, because the chart then contains only one very long run. As a result, the chart may falsely suggest one or more shifts in the data.

Figure 12.3 plots proportions rather than counts and leads to the same conclusion.

# run chart of proportions
qic(biweek, asphyxia, case,  # use the case variable for denominators
    data = births)
Run chart of proportion deliveries with neonatal asphyxia

Figure 12.2: Run chart of proportion deliveries with neonatal asphyxia

Plotting the data on a P chart provides control limits and allows for a more meaningful runs analysis, because in this case the average is a better representation of the process centre (Figure 12.3).

# P chart
qic(biweek, asphyxia, case, 
    data  = births, 
    chart = 'p')
P chart of percent deliveries with neonatal asphyxia

Figure 12.3: P chart of percent deliveries with neonatal asphyxia

While this P chart is still useful, the lower control limit is zero, which makes it impossible to detect a signal of improvement using the 3-sigma rule. We are therefore left to rely on runs analysis alone when looking for improvement.

One obvious solution is to increase subgroup size by aggregating data over longer periods, such as months or quarters. However, this comes at the cost of slower detection of both improvement and deterioration.

As an alternative, we now introduce the G chart.

12.2 G chart for opportunities between cases

The G chart plots the number of opportunities between cases, for example the number of deliveries between cases of neonatal asphyxia. This kind of data often follows a geometric distribution with standard deviation

\[\sigma = \sqrt{\bar{x}(\bar{x}+1)}\]

and control limits

\[\bar{x}\pm3\sqrt{\bar{x}(\bar{x}+1)}\]

where \(\bar{x}\) is the average number of opportunities between cases.

To obtain the number-between variable, we calculate the differences between the indices of the cases. As always, the data must first be sorted in time order.

# get indices of asphyxia cases
asph_cases <- which(births$asphyxia)

# calculate the number of deliveries between cases
asph_g <- diff(asph_cases)

# plot G chart
qic(asph_g, chart = 'g')
G chart of deliveries between neonatal asphyxia

Figure 12.4: G chart of deliveries between neonatal asphyxia

Because the geometric distribution is highly skewed, the average is not ideal for runs analysis, which assumes that the data are distributed symmetrically around the centre. For this reason, qicharts2 uses the median by default in the runs analysis for G charts (Figure 12.4).

As with other count charts, negative lower control limits are rounded up to zero, since values below zero are not possible.

Notice that process improvement, that is, fewer cases, appears as the curve moving upwards. In this example, a 3-sigma signal would require at least 518 consecutive deliveries without asphyxia.

G charts are therefore useful alternatives to P charts when events are rare and the main interest is in detecting improvement. The two charts are complementary rather than competing. The P chart is often more familiar and more useful for signalling deterioration, whereas the G chart is more helpful when looking for improvement.

12.3 T chart for time between events

The T chart plots the time between events, which is a continuous variable. Although such data could in principle be plotted on an I chart, this may not be ideal. If events occur according to a Poisson process, as is often assumed, then the time between events follows an exponential distribution, which is highly skewed. One way to deal with this is to transform the data before plotting them on an I chart. A suitable transformation is

\[x = y^{(1 / 3.6)}\]

where \(x\) is the transformed variable and \(y\) is the original time-between-events variable.

Control limits can then be calculated for the transformed data in the usual way. Because transformed values are difficult to interpret, it is often preferable to plot the original values together with back-transformed centre and control lines. The back-transformation is

\[y=x^{3.6}\].

# get time between events
asph_t <- diff(c(births$datetime[asph_cases]))
# plot T chart
qic(asph_t, chart = 't')
T chart of time (days) between deliveries with neonatal asphyxia

Figure 12.5: T chart of time (days) between deliveries with neonatal asphyxia

T charts do not allow zero time between events. This may become a problem when time is recorded with insufficient precision. For example, if two events occur on the same day and time of day is not recorded, the variable is not truly continuous. In such situations, a G chart plotting days between events may be more appropriate.

In this example, the control limits on the T chart in Figure 12.5 appear rather wide. The upper control limit indicates that about 190 days, or 27 weeks, must pass without any asphyxia before a signal is triggered. This suggests that these asphyxia cases, being discrete cases rather than random events, may not be especially well suited to analysis with a T chart.

12.4 The Bernoulli CUSUM chart for binary data

The Bernoulli CUSUM chart, also known as the B chart, is a specialised control chart for detecting small shifts in the proportion of binary outcomes over time (Neuburger et al. 2017). It operates on a logical vector (TRUE/FALSE), where each observation contributes to a trace statistic calculated as the cumulative sum of deviations from a predefined target proportion.

The trace value \(s_i\) at time \(i\) is updated recursively. When the current observation \(x_i\) is TRUE, the trace increases; when \(x_i\) is FALSE, the trace decreases. The size of each increment depends on both the target proportion and the size of the shift the chart is designed to detect. In this way, the B chart accumulates evidence over time and becomes sensitive to sustained deviations from the target.

If \(x_i\) = TRUE:

\[s_i=s_{i-1}+logOR-log(1+p_0(OR-1))\]

If \(x_i\) = FALSE:

\[s_i=s_{i-1}-log(1-p_0(OR-1))\]

Here, \(s_i\) is the CUSUM statistic at time i, \(p_0\) is the target or baseline proportion, and \(OR\) is the odds ratio representing the smallest shift the chart is designed to detect. For example, setting OR = 2 configures the chart to detect a doubling of the target proportion.

If the trace remains close to zero, the process is operating near target.

bchart(births$asphyxia, target = 0.007, or = 2, limit = 3.5)
B chart of newborns with asphyxia

Figure 12.6: B chart of newborns with asphyxia

Figure 12.6 shows the asphyxia data plotted in a B chart using the bchart() function from qicharts2. The target corresponds to the known baseline proportion of asphyxia cases, and the default settings are used for odds ratio (= 2) and control limits (= ±3.5).

Several features are worth noting:

The chart displays two cumulative traces. The upper trace is configured to detect an increase from the target proportion, specifically a doubling in the odds, while the lower trace is configured to detect a decrease, corresponding to an odds ratio of 0.5.

The trace values themselves have no direct meaning beyond their direction of movement. They increase in response to TRUE values and decrease in response to FALSE values. Their magnitude reflects the cumulative departure from target, not the absolute severity or frequency of events.

The traces are reset to zero whenever they cross the horizontal axis or a control limit, so that the chart remains responsive to new patterns of deviation.

The control limits are user-defined rather than statistically derived. Their placement therefore reflects a compromise between sensitivity, that is, the ability to detect true shifts, and specificity, that is, resistance to false alarms. By default, the limits are set at ±3.5.

Figure 12.7 compares the data with a target proportion of 2%.

bchart(births$asphyxia, target = 0.02)
B chart of newborns with asphyxia, target = 0.2%

Figure 12.7: B chart of newborns with asphyxia, target = 0.2%

Here, the lower trace signals repeatedly, suggesting that the process is not centred around a target of 2%. More specifically, the chart suggests that the current process is performing better than that target.

Conversely, when the target is set to 0.1%, the chart suggests that the process is performing worse than the target (Figure 12.8).

bchart(births$asphyxia, target = 0.001)
B chart of newborns with asphyxia, target = 0.01%

Figure 12.8: B chart of newborns with asphyxia, target = 0.01%

Control limits for B charts are user-defined and reflect a trade-off between sensitivity and specificity. A common default is ±3.5, which is suitable for detecting a doubling or halving of the event rate. Tighter limits increase sensitivity but may lead to more false alarms; wider limits reduce false alarms but may miss smaller shifts. The most appropriate choice depends on the context, including the amount of available data, the size of shift of interest, and the consequences of missing a signal or raising a false one.

In practice, limits are often chosen on the basis of simulation, historical data, and domain knowledge, especially in relation to the consequences of missed detections and false alerts. See Neuburger et al. (2017) for details on the selection of control limits.

Looking at Figure 12.6, which shows a historically stable process based on more than 2000 observations, with most data points clustered around the centre line, one might judge it reasonable to tighten the control limits to ±3 or even ±2.5, as shown in Figure 12.9.

bchart(births$asphyxia, target = 0.007, limit = 2.5)
B chart of newborns with asphyxia, target = 0.01%

Figure 12.9: B chart of newborns with asphyxia, target = 0.01%

12.5 Selecting the rigth chart for rare events

Having reviewed several alternatives for rare events data – U and T charts for low event rates, and P, G, and B charts for low case proportions – we may ask which chart should be preferred. The answer is that it depends. No single chart is universally best, and in many situations two or more charts provide complementary insights.

In particular, we recommend pairing a T or G chart with its corresponding U or P chart. Together, these pairs provide a fuller picture by helping to detect both deterioration and improvement.

The B chart is a particularly powerful tool, but it requires careful attention and specialist knowledge to construct and interpret. Its diagnostic performance depends strongly on the chosen target, odds ratio, and control limits, and these choices should be made by people with a good understanding of the underlying process. That said, the default settings – an odds ratio of 2 and control limits of ±3.5 – are generally reasonable for processes with a baseline rate around 1%, where the aim is to detect a doubling or halving of the event rate relative to that baseline.

References

Neuburger, Jenny, Kate Walker, Chris Sherlaw-Johnson, Jan van der Meulen, and David A Cromwell. 2017. “Comparison of Control Charts for Monitoring Clinical Performance Using Binary Data.” BMJ Quality & Safety 26 (11): 919–28. https://doi.org/10.1136/bmjqs-2016-005526.