---
title: "Run charts with R"
author: "Jacob Anhoej"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    fig_height: 3.5
    fig_width: 6
    fig_caption: true
    toc: true
vignette: >
  %\VignetteIndexEntry{runcharts}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteDepends{qicharts}
  %\usepackage[utf8]{inputenc}
---

## Introduction

Plotting data over time is a simple method to learn from trends, patterns, and variation in data and to study the effect of improvement efforts.

A run chart is a simple line graph of a measure over time with the median shown as a horizontal line dividing the data points so that half of the points are above the median and half are below.

```{r fig1, fig.cap='Figure 1'}
library(qicharts)
set.seed(9)         # Lock random number generator
y <- rpois(24, 16)  # Random values from to plot
qic(y)              # Plot run chart of y
```

The main purpose of the run chart is to detect process improvement or process degradation, which will turn up as non-random patterns in the distribution of data points around the median.

## Testing for non-random variation in run charts

If the process of interest shows only random variation, the data points will be randomly distributed around the median. _Random_ meaning that we cannot know if the next data point will fall above or below the median, but that the probability of each event is 50%, and that the data points are independent. _Independence_ means that the position of one data point does not influence the position of the next data point, that is, data are not auto-correlated. 

If the process shifts, these conditions are no longer true and patterns of non-random variation may be detected by statistical tests.

Non-random variation may present itself in several ways. If the process centre is shifting due to improvement or degradation we may observe _unusually_ long runs of consecutive data points on the same side of the median or that the graph crosses the median _unusually_ few times. The length of the longest run and the number of crossings in a random process are predictable within limits and depend on the total number of data points in the run chart (Anhoej 2014, Anhoej 2015).

- A **shift signal** is present if any run of consecutive data points on the same side of the median is longer than the prediction limit, `round(log2(n) + 3)`. Data points that fall on the median do not count, they do neither break nor contribute to the run (Schilling 2012).

- A **crossings signal** is present if the number of times the graph crosses the median is smaller than the prediction limit, `qbinom(0.05, n - 1, 0.5)` (Chen 2010).

`n` is the number of useful data points, that is, data points that do not fall on the median.

The shift and the crossings signals are based on a false positive signal rate around 5% and have proven useful in practice.

```{r fig2, fig.cap='Figure 2'}
y[13:24] <- rpois(12, 24)  # Introduce a shift in process mean
qic(y)                     # Plot run chart of y
```

Figure 2 show a run chart with 24 data points of which 22 are not on the median. The longest run of consecutive data points on the same side of the median is 10 (not counting the two data points that fall on the median); and the graph crosses the median 3 times. Since the longest run is longer that predicted (7) and the number of crossings is smaller than predicted (7), we may conclude that the process exhibits non-random variation.

The shift and crossings signals are two sides of the same coin and will often signal together. However, any one of them is diagnostic of non-random variation.

Signal limits may be tabulated like this:

```{r limits table}
n <- 10:30
data.frame(
  n.useful      = n,
  longest.run   = round(log2(n) + 3),
  min.crossings = qbinom(0.05, n - 1, 0.5))
```

## Analysis of before-and-after data

If data have been collected before and after a change, it may be useful to calculate the median only from the data points that belong to the before period.

```{r fig3, fig.cap='Figure 3'}
qic(y, freeze = 12)
```

If a significant change in process performance has occurred, it may be useful to split the graph in two.

```{r fig4, fig.cap='Figure 4'}
qic(y, breaks = 12)
```

## Plotting proportion and rates

If one needs to plot proportions or rates, the denominator may be provided as the second argument, `n`.

```{r fig5, fig.cap='Figure 5'}
y <- rbinom(24, 20, 0.5)                # Numerator
n <- sample(16:20, 24, replace = TRUE)  # Denominator
qic(y, n)                               # Plot run chart of y/n
```

## Using title, labels, and annotations

Tick mark labels for the x axis may be provided with the `x` argument. Chart title and labels for the x and y axis are provided the usual way using the `main`, `xlab`, and `ylab` arguments. Annotations can be added by the `notes` argument, which takes a character vector containing text to be added to individual data points.

```{r fig6, fig.cap='Figure 6'}
startdate <- as.Date('2014-1-6')
date      <- seq.Date(startdate,         # Dates for x axis labels
                      by = 'day',
                      length.out = 24)
notes     <- NA
notes[18] <- 'This is a note'            # Character vector of annotations
qic(y, n,
    x     = date,
    main  = 'Run Chart', 
    ylab  = 'Proportion',
    xlab  = 'Date',
    notes = notes)
```

## Automatic data aggregation by subgroups

Besides providing tick mark labels for the x axis, the `x` argument serves as a subgrouping vector. If, for instance, one collects data daily but wishes to aggregate data by week, this can be achieved by using weeks as the subgrouping, `x`, vector.

This example uses a data frame containing the numerator, the denominator, and the subgroups to plot.

```{r data frame}
date      <- seq.Date(startdate, by = 'day',       # 20 week long day sequence
                      length.out = 7 * 20)
n         <- sample(3:5, 7 * 20, replace = TRUE)   # Denominator vector
y         <- rbinom(7 * 20, n, 0.5)                # Numerator vector
week      <- as.Date(cut(date, 'week'))            # Subgrouping vector
d         <- data.frame(date, y, n, week)          # Data frame
head(d, 10)
```

By using the `data` argument, we may avoid the clumsy `$` notation. And by using the `week` column as the subgrouping vector, the `qic` function takes care of aggregating and plotting data by subgroups.

```{r fig7, fig.cap='Figure 7'}
qic(y, n, x = week, data = d)
```

## References

1. Jacob Anhoej, Anne Vingaard Olesen (2014). [Run Charts Revisited: A Simulation Study of Run Chart Rules for Detection of Non-Random Variation in Health Care Processes](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0113825). PLoS ONE 9(11): e113825.

2. Jacob Anhoej (2015). [Diagnostic Value of Run Chart Analysis: Using Likelihood Ratios to Compare Run Chart Rules on Simulated Data Series](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0121349). PLoS ONE 10(3): e0121349.

3. Mark F. Schilling (2012). The Surprising Predictability of Long Runs. Math. Mag. 85, 141-149

4. Zhenmin Chen (2010). A note on the runs test. Model Assisted Statistics and Applications 5, 73-77
