---
title: "Vacuum Vignette"
author: "Ron Sielinski"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Vacuum Vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup, include = FALSE}
library(vacuum)
library(ggplot2)
library(magrittr)
library(tidyr)
library(dplyr)
```

## A tidy implementation of John Tukey's vacuum cleaner 

One of John Tukey’s landmark papers, “The Future of Data Analysis” [1], contains a set of analytical procedures that focus on the analysis of two-way (contingency) tables. **Vacuum** contains an implementation of the three procedures:

* FUNOP (**FU**ll **NO**rmal **P**lot): Identifies outliers within a vector by calculating the slope (*z*) of every element, relative to the vector's median.

* FUNOR-FUNOM (**FU**ll **NO**rmal **R**ejection-**FU**ll **NO**rmal **M**odification): Treats a two-way (contingency) table for outliers by
isolating residuals from the table's likely systemic effects, which
are calculated from the table's grand, row, and column means.

* Vacuum cleaner: Removes systemic effects from values in a two-way table, returning a set of residuals that can be used for further examination (e.g., analysis of variance). 

### FUNOP
The base prodecure is FUNOP, which identifies outliers in a vector. Here’s Tukey’s definition:

>(b1) Let *a~i|n~* be a typical value for the *i*th ordered observation in a sample of *n* from a unit normal distribution.

>(b2) Let *y*~1~ ≤ *y*~2~ ≤ … ≤ *y*~n~ be the ordered values to be examined. Let *y̍* be their median (or let *y̆*, read “y trimmed”, be the mean of the *y~i~* with \frac{1}{3}*n* < *i* ≤ \frac{1}{3}(2*n*).

>(b3) For *i* ≤ ⅓*n* or > ⅓(2*n*) only, let *z~i~* = (*y~i~* – *y̍*)/*a~i|n~* (or let *z~i~* = (*y~i~* – *y̆*) /*a~i|n~*).

>(b4) Let *z̍* be the median of the *z*’s thus obtained (about \frac{1}{3}(2*n*) in number).

>(b5) Give special attention to *z*’s for which both |*y~i~* – *y̍*| ≥ A · *z̍* and *z~i~* ≥ B · *z̍* where A and B are prechosen. 

>(b5*) Particularly for small *n*, *z~j~*’s with *j* more extreme than an *i* for which (b5) selects *z~i~* also deserve special attention…  (p 23).

The basic idea is very similar to a Q-Q plot. 

Tukey gives us a sample of 14 data points. On a normal Q-Q plot, if data are normally distributed, they form a straight line. But in the chart below, based upon data from Tukey’s example, we can clearly see that a couple of the points are relatively distant from the straight line. They’re outliers. 

```{r echo = FALSE, fig.width = 7, fig.height = 4}
# Normal Q-Q plot for Tukey's sample data
data.frame(y = table_1) %>% 
ggplot(aes(sample = y)) +
  stat_qq() +
  stat_qq_line() +
  #geom_abline(aes(intercept = 0, slope = 1)) +
  labs(title = 'Normal Q-Q plot of sample data')
```  

The goal of FUNOP, though, is to eliminate the need for visual inspection by automating interpretation.

The first variable in the FUNOP procedure (*a~i|n~*) simply gives us the theoretical distribution, where *i* is the ordinal position in the range 1..*n* and Gau^-1^ is the quantile function of the normal distribution (i.e., the “Q” in Q-Q plot): 

$$
a_{i|n}=\ {\rm Gau}^{-1}\left[\frac{\left(3i-1\right)}{\left(3n+1\right)}\right]
$$

The key innovation of FUNOP is to calculate the slope of each point, relative to the median. 

``` {r echo = FALSE, fig.width = 7, fig.height = 4}
x <- table_1
n <- length(x)
y_split <- median(x)

dat <- data.frame(i = 1:n, 
                  y_sorted = sort(x),
                  a = a_qnorm(1:n, n), 
                  z = (sort(x) - y_split) / a_qnorm(1:n, n)
)

# find the point with the greatest slope
key_pt <- dat[which.max(dat$z), ]

# plot the points and the slope of the point key point
ggplot(dat, aes(x = a, y = y_sorted)) +
  geom_point(color = 'darkgrey') + 
  geom_point(dat = key_pt, aes(x = a, y = y_sorted), color = 'black') + 
  geom_point(aes(x = 0, y = y_split), color = 'red') +
  # example line
  geom_segment(aes(x = 0, y = y_split, xend = key_pt$a, yend = key_pt$y_sorted)) +
  labs(title = 'Sorted values of y plotted against theoretical distribution a(y)', y = 'y')

```
  
If *y̍* is the median of the sample, and we presume that it’s located at the midpoint of the distribution (where *a*(*y*) = 0), then the slope of each point can be calculated as:

$$
z_i=\frac{y_i - \acute{y}}{a_{i|n}}
$$

The chart above illustrates how slope of one point (1.2, 454) is calculated, relative to the median (0, 33.5).

$$
z\ =\ \frac{\Delta y}{\Delta a}\ =\ \frac{\left(454.0-\ 33.5\right)}{\left(1.2-0\right)}=350.4
$$

Any point that has a slope significantly steeper than the rest of the population is necessarily farther from the straight line. To do this, FUNOP simply compares each slope (*z~i~*) to the median of all *calculated* slopes (*z̍*). 

Note, however, that FUNOP calculates slopes for the top and bottom thirds of the sorted population only, in part because *z~i~* won’t vary much over the inner third, but also because the value of *a~i|n~* for the inner third will be close to 0 and dividing by ≈0 when calculating *z~i~* might lead to instability. 

Significance---what Tukey calls “special attention”---is partially determined by *B*, one of two predetermined values (or hyperparameters). For his example, Tukey recommends a value between 1.5 and 2.0, which means that FUNOP simply checks whether the slope of any point, relative to the midpoint, is 1.5 or 2.0 times larger than the median. 

The other predetermined value is *A*, which is roughly equivalent to the number of standard deviations of *y~i~* from *y̍* and serves as a second criterion for significance. 

The following chart shows how FUNOP works. 

``` {r echo = FALSE, , fig.width = 7, fig.height = 4}
A <- 0
B <- 1.5
tmp <- funop(table_1, A, B)
# need z & special for every point
tmp$z <- (tmp$y - median(tmp$y)) / a_qnorm(tmp$i, n)
tmp$special[is.na(tmp$special)] <- FALSE
# need outer to get special and outer to add up as a factor
tmp$outer <- !tmp$middle
z_split <- attr(tmp, 'z_split')
y_split <- attr(tmp, 'y_split')
extreme_B <- B * z_split
extreme_A <- A * z_split

ggplot(tmp, aes(x = y, y = z)) +
  geom_point(aes(color = as.factor(outer + special))) +
  geom_hline(yintercept = z_split) +
  geom_hline(yintercept = extreme_B,
             color = 'blue',
             alpha = 0.5) +
  geom_vline(xintercept = y_split,
             color = 'green3',
             alpha = 0.75) +
  geom_rect(
    xmin = -Inf,
    xmax = Inf,
    ymin = extreme_B,
    ymax = Inf,
    alpha = 0.0075,
    fill = 'blue'
  ) +
  geom_rect(
    xmin = y_split - extreme_A,
    xmax = y_split + extreme_A,
    ymin = -Inf,
    ymax = Inf,
    alpha = 0.01,
    fill = 'green'
  ) +
  labs(title = 'FUNOP: slope (z) plotted against original values') +
  labs(x = 'original values') +
  scale_color_manual(
    name = 'values',
    labels = c(
      'excluded',
      'undistinguished',
      'special attention'
    ),
    values = c('green3', 'black', 'red')
  )


```

Our original values are plotted along the *x*-axis. The points in the green make up the inner third of our sample, and we use them to calculate *y̍*, the median of just those points, indicated by the green vertical line. 

The points not in green make up the outer thirds (i.e., the top and bottom thirds) of our sample, and we use them to calculate *z̍*, the median slope of just those points, indicated by the black horizontal line.

Our first selection criterion is *z~i~* ≥ *B* · *z̍*. In his example, Tukey sets *B* = 1.5, so our threshold of interest is 1.5*z̍*, indicated by the blue horizontal line. We’ll consider any point above that line (the shaded blue region) to “deserve special attention”. We have only one such point, colored red. 

Our second criterion is |*y~i~* – *y̍*| ≥ *A* · *z̍*. In his example, Tukey sets *A* = 0, so our threshold of interest is |*y~i~* – *y̍*| ≥ 0 or (more simply) *y~i~* ≠ *y̍*. Basically, any point not on the green line. Our one point in the shaded blue region isn’t on the green line, so we still have our one point. 

Our final criterion is any *z~j~*’s with *j* more extreme than any *i* selected so far. Basically, that’s any value more extreme than the ones already identified. In this case, we have one value that’s larger (further to the right on the *x*-axis) than our red dot. That point is also colored red, and we flag it as deserving special attention. 

The two points identified by FUNOP are the same ones that we identified visually in Chart 1. 

```{r}
# Table 1 contains example data from Tukey's paper
table_1

result <- funop(table_1)
result

# value of the outliers
result[which(result$special == TRUE), 'y']

```

### FUNOR-FUNOM

A common reason for identifing outliers is to do something about them, often by trimming or Winsorizing the dataset. The former simply removes an equal number of values from upper and lower ends of a sorted dataset. Winsorizing is similar but doesn’t remove values. Instead, it replaces them with the closest original value not affected by the process. 

Tukey’s FUNOR-FUNOM offers an alternate approach. The procedure’s name reflects its purpose: FUNOR-FUNOM uses FUNOP to identify outliers, and then uses separate *rejection* and *modification* procedures to treat them. 

The technique offers a number of innovations. First, unlike trimming and Winsorizing, which affect all the values at the top and bottom ends of a sorted dataset, FUNOR-FUNOM uses FUNOP to *identify* individual outliers to treat. Second, FUNOR-FUNOM leverages statistical properties of the dataset to determine *individual* modifications for those outliers. 

FUNOR-FUNOM is specifically designed to operate on two-way (or contingency) tables. Similar to other techniques that operate on contingency tables, it uses the table’s grand mean (*x~..~*) and the row and column means (*x~j.~* and *x~.k~*, respectively) to calculate expected values for entries in the table. 

The equation below shows how these effects are combined. Because it’s unlikely for expected values to match the table’s actual values exactly, the equation includes a residual term (*y~jk~*) to account for any deviation. 

$$

\begin{align}

(actual \space value) &= (general \space \mathit{effect}) + (row \space \mathit{effect}) + (column \space \mathit{effect}) + (deviation)\\
x_{jk} &= (x_{..}) + (x_{j.} - x_{..}) + (x_{.k} - x_{..}) + (y_{jk})\\
&= x_{..} + x_{j.} - x_{..} + x_{.k} - x_{..} + y_{jk}\\
&= x_{j.} + x_{.k} - x_{..} + y_{jk}

\end{align}

$$

FUNOR-FUNOM is primarily interested in the values that deviate most from their expected values, the ones with the largest residuals. So, to calculate residuals, simply swap the above equation around:

$$
y_{jk}=x_{jk}-x_{j.}-\ x_{.k}+\ x_{..}\ 
$$

FUNOR-FUNOM starts by repeatedly applying FUNOP, looking for outlier residuals. When it finds them, it modifies the outlier with the greatest deviation by applying the following modification:

$$
x_{jk}\rightarrow x_{jk}-{\Delta x}_{jk}
$$

where

$$
\Delta x_{jk}=\ z_{jk} \cdot a_{a(n)}\cdot\frac{r\cdot c}{\left(r-1\right)\left(c-1\right)}
$$



Recalling the definition of slope (from FUNOP) 

$$
z_i=\frac{y_i - \acute{y}}{a_{i|n}}
$$
the first portion of the Δ*x~jk~* equation reduces to just *y~jk~* – *y̍*, the difference of the residual from the median. The second portion of the equation is a factor, based solely upon table size, meant to compensate for the effect of an outlier on the table’s grand, row, and column means. 

When Δ*x~jk~* is applied to the original value, the *y~jk~* terms cancel out, effectively setting the outlier to its expected value (based upon the combined effects of the contingency table) plus a factor of the median residual (~ *x~j.~* + *x~.k~* + *x~..~* + *y̍*).

FUNOR-FUNOM repeats this same process until it no longer finds values that “deserve special attention.” 

In the final phase, the FUNOM phase, the procedure uses a lower threshold of interest---FUNOP with a lower *A*---to identify a final set of outliers for treatment. The adjustment becomes

$$
\Delta x_{jk} = (z_{jk} - B_M \cdot z̍)a_{i|n}
$$

There are a couple of changes here. First, the inclusion of (–*B~M~* · *z̍*) effectively sets the residual of outliers to FUNOP’s threshold of interest, much like the way that Winsorizing sets affected values to the same cut-off threshold. FUNOM, though, sets *only the residual* of affected values to that threshold: The greater part of the value is determined by the combined effects of the grand, row, and column means. 

Second, because we’ve already taken care of the largest outliers (whose adjustment would have a more significant effect on the table’s means), we no longer need the compensating factor. 

The chart below shows the result of applying FUNOR-FUNOM to the data in Table 2 of Tukey’s paper. 

```{r echo = FALSE, fig.width = 7, fig.height = 4}
# cleanse the table
cleansed_t2 <- funor_funom(table_2)

data.frame(
  i = 1:length(table_2),
  original = as.vector(table_2),
  cleansed = as.vector(cleansed_t2)
) %>%
  ggplot(aes(x = i)) +
  geom_point(aes(y = original), color = 'black') +
  geom_segment(aes(
    x = i,
    y = original,
    xend = i,
    yend = cleansed
  ),
  color = 'black',
  linetype = 'dashed') +
  geom_point(aes(y = cleansed), color = 'grey') +
  labs(title = 'Changes to Table 2', y = 'value', x = '')

```

The grey dots represent the cleansed dataset. The black dots represent the values of affected points *prior* to treatment, and the dashed lines connect those points to their treated values, illustrating the extent to which those points changed. 

FUNOR handles the largest adjustments, which Tukey accomplishes by setting *A~R~* = 10 and *B~R~* = 1.5 for that portion of the process, and FUNOP handles the finer adjustments by setting *A~M~* = 0 and *B~M~* = 1.5.

Again, because the procedure leverages the statistical properties of the data, each of the resulting adjustments is unique. 

Here's an example on a smaller dataset.

```{r}
# create a random dataset, representing a two-way table 
set.seed(42)
dat <- matrix(rnorm(16), nrow = 4)
# add some noise
dat[2, 1] <- rnorm(1, mean = 10)

# cleanse the table of outliers
cleansed_dat <- funor_funom(dat)

# look at the two tables, before and after
dat
cleansed_dat 

# look at the difference between the two tables
# only one non-zero difference
dat - cleansed_dat

```

### Vacuum cleaner

FUNOR-FUNOM treats the outliers of a contingency table by identifying and minimizing outsized residuals, based upon the grand, row, and column means.

Tukey takes these concepts further with his vacuum cleaner, whose output is a set of residuals, which can be used to better understand sources of variance in the data and enable more informed analysis. 

To isolate residuals, Tukey’s vacuum cleaner uses regression to break down the values from the contingency table into their constituent components (p 51): 

$$
\mathit{
(original \space values) = (dual \space regression) \\
+ (deviations \space \mathit{of} \space row \space regression \space \mathit{from} \space dual \space regression) \\
+ (deviations \space \mathit{of} \space column \space regression \space \mathit{from} \space dual \space regression) \\
+ (residuals) \\
}
$$


The idea is very similar to the one based upon the grand, row, and column means. In fact, the first stage of the vacuum cleaner produces the same result as subtracting the combined effect of the means from the original values.  

To do this, the vacuum cleaner needs to calculate regression coefficients for each row and column based upon the values in our table (*y~rc~*) and a carrier---or regressor---for both rows (*a~r~*) and columns (*b~c~*). 

Below is the equation used to calculate regression coefficients for columns. 

$$
[y/a]_c= \frac{\sum_r{a_ry_{rc}}}{\sum_r a_r^2}
$$

Conveniently, the equation will give us the mean of a column when we set ar ≡ 1: 


$$
[y/a]_c= \frac{\sum_r{1 \cdot y_{rc}}}{\sum_r 1} = \frac{\sum_r{y_{rc}}}{n_r}
$$

where *n~r~* is the number of rows. Effectively, the equation iterates through every row (Σ*~r~*), summing up the individual values in the same column (*c*) and dividing by the number of rows, the same as calculating the mean (*y~.c~*).

Note, however, that *a~r~* is a vector. So to set *a~r~* ≡ 1, we need our vector to satisfy this equation:

$$
\sum_r a_r^2=1
$$

For a vector of length *n~r~* we can simply assign every member the same value:

$$
\sqrt{\frac{1}{n_r}}
$$
Our initial regressors end up being two sets of vectors, one for rows and one for columns, containing either √(1/*n~c~*) for rows or √(1/*n~r~*) for columns.

Finally, in the same way that the mean of all row means or the mean of all column means can be used to calculate the grand mean, either the row coefficients or column coefficients can be used to calculate a dual-regression (or “grand”) coefficient: 

$$
\frac{\sum_{c}{\left[y/a\right]_cb_c}}{\sum_{c} b_c^2}\ \equiv\ \frac{\sum\sum{a_ry_{rc}b_c}}{\sum\sum{a_r^2b_c^2}}\equiv\frac{\sum_{r}{a_r\left[y/b\right]_r}}{\sum_{r} a_r^2}\ =\ \left[y/ab\right]\
$$

The reason for calculating all of these coefficients, rather than simply subtracting the grand, row, and column means from our table’s original values, is that Tukey’s vacuum cleaner reuses the coefficients from this stage as regressors in the next. (To ensure *a~r~* ≡ 1 and *a~c~* ≡ 1 for the next stage, we normalize both sets of new regressors.) 

The second phase is the real innovation. Tukey takes one of his earlier ideas, one degree of freedom for non-additivity, and applies it separately to each row and column. This, Tukey tells us, “…extracts row-by-row regression upon ‘column mean minus grand mean’ and column-by-column regression on ‘row mean minus grand mean’” (p 53). 

The result is a set of residuals, vacuum cleaned of systemic effects.

``` {r}
# convert the output of vacuum_cleaner into a long data frame
dat <- vacuum_cleaner(table_8) %>% 
  as.data.frame() %>% 
  dplyr::mutate('row' = dplyr::row_number()) %>%
  tidyr::pivot_longer(cols = dplyr::starts_with('V'),
               names_to = 'col',
               values_to = 'value')

# since the table_8 was a two-way table, the rows are factors
dat$row <- as.factor(dat$row)

# conduct analysis of variance
aov_results <- aov(value ~ row + col, dat) 
summary(aov_results)
```


### Data sets and errata

There are three data sets in this package, which Tukey provided, one for each procedure. 

* **Table 1** contains example data for FUNOP (p 24)
* **Table 2** is a 36 x 15 two-way table and contains example data for FUNOR-FUNOM (p 26)
* **Table 8** is a 36 x 15 two-way table and contains example data for the vacuum cleaner (p 54). Note that Tukey's published table contains a typo, which has been corrected in this package. (See note below)

 "The Future of Data Analysis" contains a handful of errors that can make a detailed investigation of the paper difficult. Below, I've documented what I found in case anyone wants to dig more deeply themselves:

*	One of Tukey's Gaussian calculations contains an arithmetic error: *a*~4|14~ should be Gau^-1^(11/43) (p 23)
*	As mentioned, one of the values in Table 8 has an inverted sign: the value in row 23, column 10 should be 0.100, not -0.100 (p 54)
* The numerator for [y/ab], when based up the row coefficent, should be Σ*~r~* *a~r~*[*y/**b***]*~r~* instead of Σ*~r~* *a~r~*[*y/**a***]*~r~* (p 52)
*	Tukey swaps two values in when he plugs them into a key regression equation: 0.258 and 0.167 (p 55). The equation should read:

$$
0.126 - (0.014)(0.167) - (-0.459)(0.258) - (0.921)(0.258)(0.167) \space \doteq \space 0.20
$$

### References
Tukey, John W. "The Future of Data Analysis." *The Annals of Mathematical* Statistics, *33*(1), 1962, pp 1-67. *JSTOR*, <https://www.jstor.org/stable/2237638>.
