09_AdditionalResources.Rmd

---
title: "09_AdditionalResources"
output: pdf_document
---

# Additional Resources {#additional_resources}

## Time series EDA {#time_series}

The chapter [Exploratory Data Analysis](#eda) contains the basic exploratory data analysis you should do before using SPC tools. But there are many other time series-oriented analytic tools available that can help you understand the data more completely.  

There is usually far more information in a time series than is typically explored with basic SPC methods. You can create a variety of exploratory and diagnostic plots that help you understand the data more thoroughly. 

Because Rachel's data used in previous chapters has no time series patterns, we'll use the `beer` dataset provided in the `fpp2` package. It has clear time-related patterns to explore with EDA tools.

<br>
\vspace{12pt}


```{r beer_data}
# Use Australian beer data, trimmed to a 15 year subset
data(ausbeer, package = "fpp2")
beer <- window(ausbeer, start = 1990.00, end = 2005.75)
```

### Trend

The first thing to look for is whether there is a trend. The simplest way to let the data speak for this is by using a [loess smoother](https://en.wikipedia.org/wiki/Local_regression). 

The `autoplot` function in the `forecast` package provides several out-of-the-box plots for time series data, and since it's built over `ggplot2`, it can use those functions as well. 

There does seem to be an initial overall declining trend in the `beer` data that seems to flatten out.  


```{r loess_trend1, fig.height=3}
autoplot(beer) + 
  geom_smooth()
```

### Seasonplot

The seasonplot places each year as its own line over an x-axis of the sequential frequency, which defaults to the frequency of the time series. When there's no seasonal pattern across or within that frequency, the plot looks like spaghetti as the result of being driven by natural variation.  

When there is a pattern in the time series, patterns emerge. In this case, the fourth quarter increase above the other quarters is quite evident.   

```{r seasonplot2, fig.height=3}
ggseasonplot(beer)
```

### Monthplot

A monthplot puts all years into seasonal groups, where each line is a group (e.g., month) and each point in that line is an individual year. When there is a lengthy trend in the series, you can see it in a consistent up or down pattern in each seasonal group. You can also compare central tendencies across those groups with a mean or median line. 

Data with no inherent pattern shows up as noise. Whereas in a time series with temporal patterns, you can see both the higher levels in Q4 as compared with the other quarters, but you can also see that this quarter's values is declining over the years, a pattern echoed to lesser extent in the early years' values for the other quarters.  

```{r monthplot2, fig.height=3}
ggmonthplot(beer)
```

### Autocorrelation 

We've touched on autocorrelation in other portions of this book.

The `acf` function provides a graphical summary of the autocorrelation function, with each data point correlated with a value at increasing lagged distances from itself. Each correlation is plotted as a spike; spikes that go above or below the dashed line suggest that significant positive or negative autocorrelation, respectively, occurs at that lag (at the 95% confidence level). If all spikes occur inside those limits, it's safe to assume that there is no autocorrelation. If only one or perhaps two spikes exceed the limits slightly, it could be due simply to chance. Clear patterns seen in the acf plot can indicate autocorrelation even when the values do not exceed the limits. 

With the `beer` data, the patterning is obvious, especially at lags 2 (6 months apart) and 4 (1 year apart), and the correlation values are quite large.  

```{r acf2, fig.height=3}
# acf plot using the autoplot function instead of base for the ggplot look
autoplot(acf(beer, plot = FALSE))
```

<br>
\vspace{12pt}
  
The autocorrelation function is most concisely plotted with the approach above, but you can also plot the increasing lags against an initial value in individual scatterplots. If the points look like a shotgun target, there's no autocorrelation. Patterns in the points indicate autocorrelation in the data. Patterns strung along or perpendicular to the 1:1 dashed line suggest strong positive and negative correlation, respectively, though any sort of pattern is cause for concern.  

The lagplot for the `df_ts` data shows the shotgun target "pattern" that suggests that only random variation is present.   

Clear patterns emerge---especially at lag 4 (1 year apart)---in the lagplot for the `beer` data.  

```{r lagplot2}
# Scatterplot of beer data autocorrelation through first 8 lags
lag.plot(beer, lags = 8, do.lines = FALSE)
```


The `pacf` function gives you a partial autocorrelation plot, which is the correlation between the first value and each individual lag. It's the same information provided by the lag plot, only more compact as it only displays the correlation value itself. This can be quite useful in identifying cycles in data. 

Using the `beer` data shows the partial autocorrelation pattern. The spike at the second line indicates that there is a moderate negative relationship in values 6 months (2 quarters) apart, and the spike at the fourth line shows there's a strong positive relationship in values 1 year (4 quarters) apart.  

```{r fig.height=3}
autoplot(pacf(beer))
```

### Cycles

Periodograms allow you to explore a time series for cycles that may or may not be regular in timing (which makes it slightly distinct from seasonality). Sunspot cycles are a classic example at ~11 years, a time span that obviously doesn't correspond to calendar seasons and frequencies.  

Spikes in the periodogram designate possible cycle timing lengths, where the x-axis is based on frequency. The reciprocal of the frequency is the time period, so a spike in a periodogram for an annual series at a frequency of 0.09 suggests a cycle time of about 11 years.  

<br>
\vspace{12pt}

A clear spike occurs in the `beer` data at a frequency of 0.26, a time period of about 4. Since this is quarterly data, it confirms the annual pattern seen in several plots above.   

```{r periodicity2, fig.height=3}
TSA::periodogram(beer)
```

### Decomposition

The `decompose` function extracts the major pieces of a time series, while the `autoplot` function presents the results using `ggplot2` for a cleaner look. 


```{r decomp2}
autoplot(decompose(beer))
```

### Seasonal adjustment

The `seasonal` package uses the U.S. Census Bureau's X-13ARIMA-SEATS method to calculate seasonal adjustment. The `seas` function can be used to view or save the results into another object.

<br>
\vspace{12pt}


```{r seas}
# Convert ts to data frame
beer_df = tsdf(beer)

# Get seasonally-adjusted values and put into data frame
beer_season = seasonal::seas(beer)
beer_df$y_seasonal = beer_season$data[,3]

# Show top 6 lines of data frame
knitr::kable(head(beer_df))
```

<br>
\vspace{12pt}

If you just want to plot it on the fly, `ggseas` provides the `stat_seas` function for use with `ggplot2`. As with all ggplots, you need a data frame first, which the `tsdf` function provides.    

```{r seasonal2, fig.height=3, eval=FALSE}
# Plot original and seasonally adjusted data
ggplot(beer_df, aes(x, y)) + 
  geom_line(color="gray70") +
  stat_seas(color="blue")
```

### Residuals 

Residuals---the random component of the time series---can also be explored for potential patterns. Ideally, you don't want to see patterns in the residuals, but they're worth exploring in the name of thoroughness. 

```{r residuals_beer}
# Convert ts residuals to data frame
beer_df_rand = tsdf(decompose(beer)$random)

# Add quarter as a factor
beer_df_rand$qtr = factor(quarter(date_decimal(beer_df_rand$x)))

# Plot residuals, with custom colors
ggplot(beer_df_rand, aes(x, y)) + 
  geom_hline(yintercept=0, linetype="dotted") +
  geom_smooth(color = "gray70", alpha = 0.2) +
  geom_point(aes(color = qtr)) +
  scale_color_manual(values=c("#E69F00", "#56B4E9", "#009E73", "#000000"))
```

<br>
\vspace{12pt}

We can take that same information and facet by quarter for a different perspective. 

<br>
\vspace{12pt}

```{r residuals_faceted_beer}
# Residuals faceted by quarter
ggplot(beer_df_rand, aes(x, y)) + 
  geom_hline(yintercept=0, linetype="dotted") +
  geom_smooth(color = "gray70", alpha = 0.2) +
  facet_wrap(~ qtr) +
  geom_point(aes(color = qtr)) +
  scale_color_manual(values=c("#E69F00", "#56B4E9", "#009E73", "#000000"))
```

### Accumuluation plots 

You can use the EDA tools above on rates, numerators, and denominators alike to explore patterns. When you do have a numerator and a denominator that create your metric, you can also plot them against each other, looking at the accumulation of each over the course of a relevant time frame (e.g., a year).  

To illustrate, we'll create a new time series for monthly central line associated infections, set up so that the last two years of a 10 year series are based on a different process.  

```{r accumplot_data}
# Generate sample data
set.seed(54)
bsi_8yr = data.frame(Linedays = sample(1000:2000, 96),
                     Infections = rpois(96, 4))
bsi_2yr = data.frame(Linedays = sample(1200:2200, 24),
                     Infections = rpois(24, 3))
bsi_10yr = rbind(bsi_8yr, bsi_2yr)
bsi_10yr$Month = seq(as.Date("2007/1/1"), by = "month", length.out = 120)
bsi_10yr$Year = year(bsi_10yr$Month)
bsi_10yr$Rate = round((bsi_10yr$Infections / bsi_10yr$Linedays * 1000), 2)
```

<br>
\vspace{12pt}

First, calculate the cumulative sums for the numerator and denominator for the time period of interest. Here, we use years.  

<br>
\vspace{12pt}

```{r accumplot_calcs}
# Calculate cumulative sums by year
accum_bsi_df = bsi_10yr %>% 
  group_by(Year) %>% 
  arrange(Month) %>% 
  mutate(cuml_linedays = cumsum(Linedays),
         cuml_infections = cumsum(Infections))
```


Then, plot them against each other. Much like a seasonplot, a spaghetti "pattern" indicates that only random, common cause variation is acting on the variables. Strands (individual years) that separate from that mess of lines suggest that a different process is in place for those strands.  


```{r accumplot}
# Accumulation plot
ggplot(accum_bsi_df, aes(x = cuml_linedays, y = cuml_infections,
                         group = as.factor(Year))) +
  geom_path(aes(color = as.factor(Year)), size = 1) +
  geom_point(aes(color = as.factor(Year))) +
  scale_y_continuous(name = "Cummulative Infections",
                     breaks = seq(0,120,10)) +
  scale_x_continuous(name = "Cumulative Central Line Days",
                     breaks = seq(0,40000,5000)) +
  scale_colour_brewer(type = "div", palette = "Spectral") +
  guides(color = guide_legend(title = "Year")) +
  ggtitle("Infections vesus Central Line Days by Year")
```


## Custom SPC Function {#code_customSPC_function}

We've created a function that highlights points that may indicate special cause variation. Using this function does require some thought about setting up the variables, which was done on purpose---you should put as much care into the construction of run and control charts as your nurses put into patient care. To do any less is a disservice to the decision-makers who would rely on your work and the patients that rely on these decision-makers to provide the conditions that support the best care possible. 

```{r eval = F}
plotSPC <- function(subgroup, point, mean, sigma, k = 3,
                    ucl.show = TRUE, lcl.show = TRUE,
                    band.show = TRUE, rule.show = TRUE,
                    ucl.max = Inf, lcl.min = -Inf,
                    label.x = "Subgroup", label.y = "Value") {
  # Plots control chart with ggplot
  
  # Args:
  # subgroup: Subgroup definition (for x-axis)
  # point: Subgroup sample values (for y-axis)
  # mean: Process mean value (for center line)
  # sigma: Process variation value (for control limits)
  # k: Specification for k-sigma limits above and below center line,
  #    default is 3
  # ucl.show: Visible upper control limit? Default is true
  # lcl.show: Visible lower control limit? Default is true
  # band.show: Visible bands between 1-2 sigma limits? Default is true
  # rule.show: Highlight run rule indicators in orange? Default is true
  # ucl.max: Maximum feasible value for upper control limit
  # lcl.min: Minimum feasible value for lower control limit
  # label.x: Specify x-axis label
  # label.y: Specify y-axis label

  df <- data.frame(subgroup, point)
  df$ucl <- pmin(ucl.max, mean + k*sigma)
  df$lcl <- pmax(lcl.min, mean - k*sigma)
  warn.points <- function(rule, num, den) {
    sets <- mapply(seq, 1:(length(subgroup) - (den - 1)),
                  den:length(subgroup))
    hits <- apply(sets, 2, function(x) sum(rule[x])) >= num
    intersect(c(sets[,hits]), which(rule))
  }
  orange.sigma <- numeric()

  p <- ggplot(data = df, aes(x = subgroup)) +
    geom_hline(yintercept = mean, col = "gray", size = 1)
  if (ucl.show) {
    p <- p + geom_line(aes(y = ucl), col = "gray", size = 1)
  }
  if (lcl.show) {
    p <- p + geom_line(aes(y = lcl), col = "gray", size = 1)
  }
  if (band.show) {
    p <- p +
      geom_ribbon(aes(ymin = mean + sigma,
                      ymax = mean + 2*sigma), alpha = 0.1) +
      geom_ribbon(aes(ymin = pmax(lcl.min, mean - 2*sigma),
                      ymax = mean - sigma), alpha = 0.1)
    orange.sigma <- unique(c(
      warn.points(point > mean + sigma, 4, 5),
      warn.points(point < mean - sigma, 4, 5),
      warn.points(point > mean + 2*sigma, 2, 3),
      warn.points(point < mean - 2*sigma, 2, 3)
    ))
  }
  df$warn <- "blue"
  if (rule.show) {
    shift.n <- round(log(sum(point!=mean), 2) + 3)
    orange <- unique(c(orange.sigma,
                      warn.points(point > mean - sigma &
                                  point < mean + sigma, 15, 15),
                      warn.points(point > mean, shift.n, shift.n),
                      warn.points(point < mean, shift.n, shift.n)))
    df$warn[orange] <- "orange"
  }
  df$warn[point > df$ucl | point < df$lcl] <- "red"


 p +
    geom_line(aes(y = point), col = "royalblue3") +
    geom_point(data = df, aes(x = subgroup, y = point, col = warn)) +
    scale_color_manual(values = c("blue" = "royalblue3",
                                  "orange" = "orangered",
                                  "red" = "red3"),
                       guide = FALSE) +
    labs(x = label.x, y = label.y) +
    theme_bw()
}

```

<br>
\vspace{12pt}

The following code sets up the variables for replicating the control charts in the chapter [A Guide to Control Charts](#guide_controlCharts) using the simulated data and the custom SPC function instead of qiccharts2.

<br>
\vspace{12pt}

```{r eval=F}
# Calculate u chart inputs
subgroup.u <- unique(months)
point.u <- infections.agg / linedays.agg * 1000
central.u <- sum(infections.agg) / sum(linedays.agg) * 1000
sigma.u <- sqrt(central.u / linedays.agg * 1000)

# Plot u chart
plotSPC(subgroup.u, point.u, central.u, sigma.u, k = 3, lcl.min  = 0,
         label.x = "Month", label.y = "Infections per 1000 line days")



# Calculate p chart inputs
subgroup.p <- dates
point.p <- readmits / discharges
central.p <- sum(readmits) / sum(discharges)
sigma.p <- sqrt(central.p*(1 - central.p) / discharges)

# Plot p chart
plotSPC(subgroup.p, point.p, central.p, sigma.p,
         label.x = "Month", label.y = "Proportion readmitted")



# Calculate g chart inputs
subgroup.g <- seq(2, length(infections.index))
point.g <- linedays.btwn
central.g <- mean(point.g)
sigma.g <- rep(sqrt(central.g*(central.g+1)), length(point.g))

# Plot g chart
plotSPC(subgroup.g, point.g, central.g, sigma.g, lcl.show = FALSE,
         band.show = FALSE, rule.show = FALSE,
         lcl.min = 0, k = 3, label.x = "Infection number",
         label.y = "Line days between infections")




# Calculate IMR control chart inputs
subgroup.i <- seq(1, length(exit))
subgroup.mr <- seq(1, length(exit) - 1)

point.i <- exit - arrival
point.mr <- matrix(nrow = length(point.i) - 1)
for (i in 1:length(point.i) - 1) {
    point.mr[i] <- abs(point.i[i + 1] - point.i[i])
}

mean.i <- mean(point.i)
mean.mr0 <- mean(point.mr)
mean.mr <- mean(point.mr[point.mr<=3.27*mean.mr0])
sigma.i <- rep(mean.mr, length(subgroup.i))
sigma.mr <- rep(mean.mr, length(subgroup.mr))

# Plot MR chart
plotSPC(subgroup.mr, point.mr, mean.mr, sigma.mr, k = 3.27,
        lcl.show = FALSE, band.show = FALSE,
        label.x = "Test number",
        label.y = "Turnaround time (moving range)")

# Plot I chart
plotSPC(subgroup.i, point.i, mean.i, sigma.i, k = 2.66,
         lcl.min = 0, band.show = FALSE,
         label.x = "Test number", label.y = "Turnaround time")




# Calculate XbarS control chart inputs
subgroup.x <- as.Date(unique(months))
subgroup.s <- subgroup.x
point.x <- aggregate(dfw$waits, by = list(months), FUN = mean,
                     na.rm = TRUE)$x
point.s <- aggregate(dfw$waits, by = list(months), FUN = sd, na.rm = TRUE)$x
mean.x <- mean(waits)
mean.s <- sqrt(sum((sample.n - 1) * point.s ^ 2) /
        (sum(sample.n) - length(sample.n)))
sigma.x <- mean.s / sqrt(sample.n)
c4 <- sqrt(2 / (sample.n - 1)) * gamma(sample.n / 2) /
  gamma((sample.n - 1) / 2)
sigma.s <- mean.s * sqrt(1 - c4 ^ 2)

# Plot s chart
plotSPC(subgroup.s, point.s, mean.s, sigma.s, k = 3,
         label.x = "Month", label.y = "Wait times standard deviation (s)")

# Plot xbar chart
plotSPC(subgroup.x, point.x, mean.x, sigma.x, k = 3,
         label.x = "Month", label.y = "Wait times average (x)")



# Calculate t chart inputs
subgroup.t <- subgroup.g
point.t <- y
central.t <- mean(y)
sigma.t <- rep(mr_prime, length(point.t))

# Plot t chart
plotSPC(subgroup.t, point.t, central.t, sigma.t, lcl.show = FALSE,
         band.show = FALSE, rule.show = FALSE,
         lcl.min = 0, k = 2.66, label.x = "Infection number",
         label.y = "Line days between infections (transformed)")
```


## Code used to generate examples {#code_examples}
```{r eval = F}
### Rachel's Data
# Set seed for reproducibility
set.seed(2019)

# Generate fake infections data
dates <- strftime(seq(as.Date("2013/10/1"), by = "day", length.out = 730),
                  "%Y-%m-01")
linedays <- sample(30:60,length(dates), replace = TRUE)
infections <- rpois(length(dates), 2/1000*linedays)


# Aggregate the data by month
infections <- aggregate(infections, by = list(dates), FUN = sum,
                        na.rm = TRUE)$x
linedays <- aggregate(linedays, by = list(dates), FUN = sum, na.rm = TRUE)$x
months <- unique(dates)

# Create a tibble
rachel_data = tibble(months, infections, linedays)



### example control charts

# Set seed for reproducibility
set.seed(72)

#### u chart 
# Generate fake infections data
dates <- seq(as.Date("2013/10/1"), by = "day", length.out = 730)
linedays <- sample(30:60,length(dates), replace = TRUE)
infections <- rpois(length(dates), 2/1000*linedays)

# Aggregate the data by month
infections <- aggregate(infections, by = list(dates), FUN = sum,
                        na.rm = TRUE)$x
linedays <- aggregate(linedays, by = list(dates), FUN = sum, na.rm = TRUE)$x
months <- unique(dates)

# Create a tibble
uchart_data <- tibble(months, infections, linedays)


#### p chart
# Generate sample data
discharges <- sample(300:500, 24)
readmits <- rbinom(24, discharges, .2)
dates <- seq(as.Date("2013/10/1"), by = "month", length.out = 24)

# Create a tibble
pchart_data <- tibble(dates, readmits, discharges)


#### g chart
# Generate fake data using u-chart example data
infections.index <- replace_na(which(infections > 0)[1:30], 0)
dfind <- data.frame(start = head(infections.index,
                                 length(infections.index) - 1) + 1,
                   end = tail(infections.index,
                              length(infections.index) - 1))

linedays.btwn <- matrix(nrow = length(dfind$start))

for (i in 1:length(linedays.btwn)) {
  sumover <- seq(dfind$start[i], dfind$end[i])
  linedays.btwn[i] <- sum(linedays[sumover])
}

gchart_data <- tibble(inf_index = 1:length(linedays.btwn),
                      days_between = linedays.btwn)


#### IMR chart
# Generate fake data
arrival <- cumsum(rexp(24, 1/10))
process <- rnorm(24, 5)
exit <- matrix(nrow = length(arrival))
exit[1] <- arrival[1] + process[1]

for (i in 1:length(arrival)) {
  exit[i] <- max(arrival[i], exit[i - 1]) + process[i]
}

imrchart_data <- tibble(turnaround_time = exit - arrival,
                        test_num = 1:length(exit))



#### XbarS chart
# Generate fake patient wait times data
set.seed(777)
waits <- c(rnorm(1700, 30, 5), rnorm(650, 29.5, 5))
months <- strftime(sort(as.Date('2013-10-01') +
                        sample(0:729, length(waits), TRUE)), "%Y-%m-01")
sample.n <- as.numeric(table(months))
xbarschart_data <- tibble(months, waits)



##### t chart
# Generate sample data using g-chart example data
y <- linedays.btwn ^ (1/3.6)
mr <- matrix(nrow = length(y) - 1)
for (i in 1:length(y) - 1) {
  mr[i] <- abs(y[i + 1] - y[i])
}
mr <-  mr[mr <= 3.27 * mean(mr)]
tchart_data <- tibble(inf_index = 1:length(y), days_between = y)


##### change in process example
# Create fake data with change in process at 28 months
intervention = data.frame(date = seq(as.Date("2006-01-01"), by = 'month',
                                     length.out = 48),
                          y = c(rpois(28, 6), rpois(20, 3)),
                          n = round(rnorm(48, 450, 50)))
```



## Useful References {#useful_resources}
- For more information, a good overview of run charts can be found in Perla et al. 2011, [*The run chart: a simple analytical tool for learning from variation in healthcare processes*](http://www.med.unc.edu/cce/files/education-training/The%20run%20chart%20a%20simple%20analytical%20tool.pdf), BMJ Quality & Safety 20:46-51.  

<br>
\vspace{12pt}

- A straight-to-the-point reference/tool for doing run charts in R is Anhøj 2016, [*Run charts with R*](https://cran.r-project.org/web/packages/qicharts/vignettes/runcharts.html).

<br>
\vspace{12pt}

- Some good overview papers on control charts include Benneyan et al. 2003, [*Statistical process control as a tool for research and healthcare improvement*](http://qualitysafety.bmj.com/content/12/6/458.full.pdf), BMJ Quality & Safety 12:458-464; Mohammed et al. 2008, [*Plotting basic control charts: tutorial notes for healthcare practitioners*](https://www.researchgate.net/profile/William_Woodall/publication/5468089_Plotting_control_charts_Tutorial_notes_for_healthcare_practitioners/links/00b49521d1165f1f49000000.pdf), BMJ Quality & Safety 17:137-145; and Limaye et al. 2008, [*A Case Study in Monitoring Hospital--Associated Infections with Count Control Charts*](https://www.researchgate.net/profile/Christina_Mastrangelo/publication/233015368_A_Case_Study_in_Monitoring_Hospital-Associated_Infections_with_Count_Control_Charts/links/552c5b530cf29b22c9c44787/A-Case-Study-in-Monitoring-Hospital-Associated-Infections-with-Count-Control-Charts.pdf), Quality Engineering 20:404-413. [Wheeler 2010](http://www.qualitydigest.com/inside/quality-insider-column/individual-charts-done-right-and-wrong.html) covers why you shouldn't use 3$\sigma$ for control limits in *I* charts.    

<br>
\vspace{12pt}

- A straight-to-the-point reference/tool for doing control charts in R is Anhøj 2016, [*Control Charts with qicharts for R*](https://cran.r-project.org/web/packages/qicharts/vignettes/controlcharts.html).

<br>
\vspace{12pt}

- A good basic overview book is Carey and Lloyd 2001, [*Measuring Quality Improvement in Healthcare*](https://www.amazon.com/Measuring-Quality-Improvement-Healthcare-Applications/dp/0527762938/), American Society for Quality.  

<br>
\vspace{12pt}

- A good book that covers both basic and advanced topics is Provost and Murray 2011, [*The Health Care Data Guide*](https://www.amazon.com/Health-Care-Data-Guide-Improvement/dp/0470902582/), Jossey-Bass.  

<br>
\vspace{12pt}

- The papers that discuss the uselessness of the trend test in run and control charts include Davis & Woodall 1988, [*Performance of the control chart trend rule under linear shift*](http://asq.org/qic/display-item/index.pl?item=5597), Journal of Quality Technology 20:260-262, and Anhøj & Olesen 2014, [*Run charts revisited: A simulation study of run chart rules for detection of non-random variation in health care processes*](http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0113825), PLOS One 9(11): e113825.

<br>
\vspace{12pt}

- Finally, some important warnings about when control charts fail (and a useful alternative, GAMs) can be found in Morton et al. 2009, [*Hospital adverse events and control charts: the need for a new paradigm*](http://www.journalofhospitalinfection.com/article/S0195-6701(09)00340-5/abstract), Journal of Hospital Infection 73(3):225–231, as well as in Morton et al. 2007, [*New control chart methods for monitoring MROs in Hospitals*](https://www.researchgate.net/profile/Edward_Tong2/publication/43477704_New_control_chart_methods_for_monitoring_MROs_in_hospitals/links/5600d22e08aec948c4fa93cd.pdf), Australian Infection Control 12(1):14-18.    

<br>
\vspace{12pt}

- Wikipedia is a good place to start learning about probability distributions and their mean-variance relationships, e.g., (click the name to go to the link):   
    - [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution)
    - [binomial](https://en.wikipedia.org/wiki/Binomial_distribution)
    - [normal](https://en.wikipedia.org/wiki/Normal_distribution)
    - [geometric](https://en.wikipedia.org/wiki/Geometric_distribution)
    - [Weibull](https://en.wikipedia.org/wiki/Weibull_distribution)
    - [A gallery of distributions (NIST)](http://www.itl.nist.gov/div898/handbook/eda/section3/eda366.htm)
    - [Common probability distributions: the data scientist’s crib sheet (Cloudera)](http://blog.cloudera.com/blog/2015/12/common-probability-distributions-the-data-scientists-crib-sheet/)  
  
<br>
\vspace{12pt}

- The template used for the book can be found at https://github.com/sydneykpaul/blue-bookdown-latex-theme, and was derived from The Legrand Orange Book, which can be found here: https://www.latextemplates.com/template/the-legrand-orange-book. The photos used in this template are courtesy of https://www.pexels.com/.