tests.qmd

# Statistical tests {#sec-tests}

```{r}
#| label: setup4
#| include: false
#| cache: false

source("before-each-chapter.R")
```

The early history of anomaly detection methods involved statistical tests and assumptions about the distribution of the data.
There are much better methods available now, but these old methods are still widely used.
So they are included here for historical interest, and to point out their flaws to those still using them.
Apart from the six examples introduced in @sec-examples, which are reused in @sec-boxplots, nothing in later chapters assumes knowledge of this chapter.
So if you want to skip on to more useful methods, read @sec-examples and then head to @sec-boxplots.

Only univariate methods are considered in this chapter.
There have been multivariate parametric tests for anomaly detection, but they are much less widely used and are also best avoided.

## Examples {#sec-examples}

The methods are tested on six examples: two involving real data, and four involving simulated data.

1. The cricket batting averages from the `cricket_batting` data set.
   As discussed in @sec-scatterplots, these appear to contain one genuine anomaly, the batter Don Bradman.
2. Old Faithful eruption durations between 2017 and 2023, also discussed in @sec-scatterplots.
   For these data, there were three extreme anomalies --- with durations up to 1 minute --- and several "inliers" in an area of low density between 140 and 180 seconds.
3. The first 18 rows of the second variable in the `n01` data, along with the values 4.0 and 4.5.
   The latter two are anomalies as they are unlikely to arise from the N(0,1) distribution.
   We save this data set as `n01b`.
   @fig-n01b shows the data for this example.

```{r}
#| label: fig-n01b
#| code-fold: false
#| fig-asp: 0.2
#| fig-cap: "Eighteen observations from a N(0,1) distribution along with two anomalies."

n01b <- tibble(y = c(n01$v2[1:18], 4, 4.5))
n01b |>
  ggplot(aes(x = y, y = 1)) +
  geom_jitter(width = 0, alpha = 0.5) +
  scale_y_discrete() +
  labs(y = "", x = "Synthetic data")
```

The other three simulated examples will use 1000 observations from each of the following distributions:

4. an N(0,1) distribution (the first variable in the `n01` data set);
5. a $\text{t}_3$ distribution, which has heavy tails (but finite variance);
6. a $\chi^2_4$ distribution, which is positively skewed.

The density functions for these three distributions can be seen in @fig-testdist.

```{r}
#| label: fig-testdist
#| fig-cap: "Three distributions used to test the anomaly detection methods in this chapter."
#| message: false

c(t3 = dist_student_t(3), X4 = dist_chisq(4), N01 = dist_normal()) |>
  gg_density() +
  labs(y = "Probability density") +
  scale_colour_discrete(
    breaks = c("t3", "X4", "N01"),
    labels = c(
      TeX("$t_{3}$"),
      TeX("$\\chi^2_4$"),
      "N(0,1)"
    )
  ) +
  theme(legend.text = element_text(hjust = 0)) +
  coord_cartesian(xlim = c(-5, 10))
```

Data for the last two examples are generated below.

```{r}
#| label: testdist-data
#| code-fold: false

set.seed(1)
t3 <- tibble(y = rt(1000, df = 3))
chisq4 <- tibble(y = rchisq(1000, df = 4))
```

All methods considered in this chapter assume the underlying data follow a Normal distribution.
In our examples, only #4 is from a Normal distribution, so the methods should work well in that case, but perhaps not in the other cases.

Good anomaly detection methods should pick up one anomaly in the cricket batting example, at least three anomalies in the Old Faithful example, and the two anomalies in the `n01b` example.
They should identify no anomalies in the remaining examples.

```{r}
#| label: check_numbers
#| include: false
#| dependson: ["fig-n01b", "testdist-data"]

nanomalies <- tibble(
  Example = c(
    "Cricket batting",
    "Old Faithful duration",
    "N(0,1) + 2 outliers",
    "N(0,1)",
    "t(3)",
    "Chi-squared(4)"
  ),
  `N observations` = c(
    NROW(filter(cricket_batting, Innings >= 20)),
    NROW(oldfaithful),
    NROW(n01b),
    NROW(n01),
    NROW(t3),
    NROW(chisq4)
  ),
  `Expected anomalies` = c(
    1,
    3,
    2,
    0,
    0,
    0
  )
)
```

## Z-scores {#sec-zscores}

Many parametric methods are based on $z$-scores (@sec-scaling) and assume that the data come from a Normal distribution.
Recall that if our data are given by $y_1,\dots,y_n$, then their $z$-scores are given by
$$
  z_i = (y_i - \bar{y})/s_y
$$
where $\bar{y}$ is the mean and $s_y$ is the standard deviation of the observations.

Some books recommend that observations be identified as anomalies when the absolute value of the corresponding $z$ is above some threshold (usually 3).
This is a bad idea for several reasons.

1. If the data really do come from a Normal distribution that has been contaminated with anomalies, then the estimated mean ($\bar{y}$) and standard deviation ($s_y$) will also be affected by those anomalies.
   Any anomaly detection method should be relatively robust to the anomalies in the data.
2. Having a fixed threshold regardless of sample size means that the probability of a spurious anomaly increases with the sample size.
3. The assumption of Normality is unlikely to be satisfied with most real data.

### Examples

We now examine what happens when we apply this test to the six examples.

```{r}
#| label: zscore1
#| code-fold: false
#| dependson: ["fig-n01b", "testdist-data"]

cricket_batting |>
  filter(Innings >= 20) |>
  mutate(z = (Average - mean(Average)) / sd(Average)) |>
  filter(abs(z) > 3) |>
  select(Player, Average, z)
oldfaithful |>
  mutate(z = (duration - mean(duration)) / sd(duration)) |>
  filter(abs(z) > 3)
n01b |>
  mutate(z = (y - mean(y)) / sd(y)) |>
  filter(abs(z) > 3)
n01 |>
  select(v1) |>
  mutate(z = (v1 - mean(v1)) / sd(v1)) |>
  filter(abs(z) > 3)
t3 |>
  mutate(z = (y - mean(y)) / sd(y)) |>
  filter(abs(z) > 3)
chisq4 |>
  mutate(z = (y - mean(y)) / sd(y)) |>
  filter(abs(z) > 3)
```

1. Don Bradman is correctly identified as the only anomaly in the cricket batting averages.
2. The three obvious anomalies in the Old Faithful data have been identified, but none of the unusual eruptions between 140 and 180 seconds in length.
3. Neither of the two anomalies added to the 18 N(0,1) observations has been correctly identified.
4. One spurious anomaly has been identified in the 1000 N(0,1) observations.
5. Eighteen spurious anomalies have been identified in the 1000 $\text{t}_3$ observations.
6. Fifteen spurious anomalies have been identified in the 1000 $\chi^2_4$ observations.

```{r}
#| label: check2
#| include: false
#| dependson: check_numbers

nanomalies <- nanomalies |>
  mutate(
    Zscore = c(
      cricket_batting |>
        filter(Innings >= 20) |>
        mutate(z = (Average - mean(Average)) / sd(Average)) |>
        filter(abs(z) > 3) |>
        NROW(),
      oldfaithful |>
        mutate(z = (duration - mean(duration)) / sd(duration)) |>
        filter(abs(z) > 3) |>
        NROW(),
      n01b |>
        mutate(z = (y - mean(y)) / sd(y)) |>
        filter(abs(z) > 3) |>
        NROW(),
      n01 |>
        mutate(z = (v1 - mean(v1)) / sd(v1)) |>
        filter(abs(z) > 3) |>
        NROW(),
      t3 |>
        mutate(z = (y - mean(y)) / sd(y)) |>
        filter(abs(z) > 3) |>
        NROW(),
      chisq4 |>
        mutate(z = (y - mean(y)) / sd(y)) |>
        filter(abs(z) > 3) |>
        NROW()
    )
  )
stopifnot(nanomalies$Zscore == c(1, 3, 0, 1, 18, 15))
```

### Probability of spurious anomalies

Even if the data did come from a Normal distribution, the probability of an observation being more than 3 standard deviations from the mean is $`r round(2 * pnorm(-3),4)`$, so we would expect to see approximately 1 in every $1 / `r round(2 * pnorm(-3),4)` = `r round(.5/pnorm(-3))`$ regular observations being identified as an anomaly using this approach.

In fact, we can be more accurate than this.
If our $n$ observations come from a Normal distribution, then the $z$-scores approximately follow a N(0,1) distribution.
(The approximation arises because we replace the mean and standard deviation by their sample estimates.)
When we identify anomalies as those points with $|z| > c$, then the probability of finding at least one spurious anomaly is the probability of the maximum being above $c$ or the minimum being below $-c$.
Since the distribution is symmetric, these probabilities are the same.
Let $M_n$ be the maximum of $n$ $z$-scores computed from the data; then $\Pr(M_n \le c) \approx (\Phi(c))^n$ where $\Phi$ is the standard Normal cdf.
So the probability of at least one spurious anomaly is approximately
$$
  1 - (\Phi(c))^{2n}.
$$ {#eq-spuriousz}
This approximation is accurate for large $n$, but not for small sample sizes.
Because each observation contributes to both $\bar{y}$ and $s_y$, an extreme $y_i$ pulls the mean towards itself and inflates the standard deviation, so its $z$-score cannot be very large; in fact, $|z_i| \le (n - 1)/\sqrt{n}$.
For example, with $n \le 10$ no $z$-score can exceed 3, whatever the data.
So, instead, the probabilities are computed by simulation; these are plotted in @fig-zspurious for several values of $c$ and $n$.

```{r}
#| label: pspurious
#| echo: false

# Function to compute probability of spurious anomaly for sample size n,
# distribution dist with parameters in ...
# c is zscore threshold, rep is number of replicates
pspurious <- function(n, dist, c = 3, rep = 1e5, ...) {
  if (length(n) == 1) {
    anomaly <- numeric(rep)
    for (j in seq(rep)) {
      anomaly[j] <- max(abs(scale(dist(n, ...)))) > c
    }
    return(mean(anomaly))
  } else {
    pv <- numeric(length(n))
    for (i in seq_along(n)) {
      pv[i] <- pspurious(n = n[i], dist = dist, c = c[i], rep = rep, ...)
    }
    return(pv)
  }
}
n <- c(2:10, round(exp(seq(log(15), log(3e4), l = 99))))
breaks <- sort(c(2:6, 8, 10^(1:4), 2.5 * 10^(1:3), 5 * 10^(1:3)))
```

```{r}
#| label: fig-zspurious
#| echo: false
#| dependson: pspurious
#| fig-cap: The probability of at least one spurious anomaly using $z$-scores and with a threshold $c$ for Normally distributed data.

zspuriousdf <- tidyr::expand_grid(n = n, c = seq(2, 3.5, by = 0.5)) |>
  mutate(pmax = pspurious(n, rnorm, c = c)) |>
  cache("zspuriousdf")
labels <- zspuriousdf |>
  group_by(c) |>
  filter(abs(pmax - 0.75) == min(abs(pmax - 0.75))) |>
  slice(1) |>
  mutate(label = TeX(sprintf("$c = %.1f$", c), output = "character"))
zspuriousdf |>
  ggplot(aes(x = n, y = pmax, group = c, col = factor(c))) +
  geom_line() +
  scale_x_log10(
    breaks = breaks,
    labels = sprintf("%d", breaks),
    minor_breaks = NULL
  ) +
  labs(
    y = "Probability of at least one spurious anomaly",
    x = TeX("Sample size ($n$)")
  ) +
  geom_label(
    aes(x = n * 0.9, y = 0.78, label = label, hjust = "right"),
    data = labels,
    parse = TRUE
  ) +
  guides(col = "none") +
  theme(
    axis.text.x = element_text(size = 8),
    plot.margin = margin_part(r = 12)
  ) +
  coord_cartesian(xlim = c(6, 7500))
```

Let's also compute the probability of spurious anomalies for different data distributions --- the $\text{t}_3$ and $\chi^2_4$ distributions shown in @fig-testdist.
The results are shown in @fig-zspurious2.

```{r}
#| label: fig-zspurious2
#| echo: false
#| dependson: pspurious
#| fig-cap: The probability of at least one spurious anomaly using $z$-scores for different data distributions and with a threshold of $c=3$.

zspurious2df <- tibble(n = n, c = 3) |>
  mutate(
    ApproxNormal = 1 - pt(c, n - 1)^(2 * n),
    Normal = pspurious(n, rnorm, c = c),
    `t(3)` = pspurious(n, rt, df = 3, c = c),
    `Chi-squared(4)` = pspurious(n, rchisq, df = 4, c = c)
  ) |>
  select(-c) |>
  tidyr::pivot_longer(
    -n,
    values_to = "Probability",
    names_to = "Distribution"
  ) |>
  cache("zspurious2df")
zspurious2df |>
  filter(Distribution != "ApproxNormal") |>
  ggplot(aes(
    x = n,
    y = Probability,
    group = Distribution,
    col = Distribution
  )) +
  geom_line() +
  labs(
    x = TeX("$n$ (sample size)"),
    y = "Probability of at least one spurious anomaly"
  ) +
  scale_x_log10(
    breaks = breaks,
    labels = sprintf("%d", breaks),
    minor_breaks = NULL
  ) +
  labs(
    y = "Probability of at least one spurious anomaly",
    x = TeX("Sample size ($n$)")
  ) +
  scale_colour_discrete(
    breaks = c("t(3)", "Chi-squared(4)", "Normal"),
    labels = c(
      TeX("$t_{3}$"),
      TeX("$\\chi^2_4$"),
      "N(0,1)"
    )
  ) +
  theme(legend.text = element_text(hjust = 0)) +
  theme(axis.text.x = element_text(size = 8)) +
  coord_cartesian(xlim = c(6, 7500))
```

Notice how the departures from Normality, either via skewness (for the $\chi^2$ distribution) or with heavier tails (for the t distribution), lead to much greater probabilities for spurious anomalies.
Any tests for anomaly detection that assume an underlying data distribution will be sensitive to the shape of that distribution.

Even if we were prepared to believe that the data come from a Normal distribution, we need to adjust the threshold to allow for the sample size.
This is the idea behind the next methods we will consider.

## Peirce's and Chauvenet's criteria {#sec-peirce-chauvenet}

### Peirce's criterion {-}

Benjamin Peirce was a Harvard mathematician in the mid-1800s who worked with astronomical data that were prone to anomalous observations.
He proposed the first known test based on what are now called $z$-scores [@Peirce1852].
His criterion was that any observations with $|z| > c$ should be "rejected", where $c$ is a complicated function depending on the sample size $n$ and the number of suspected anomalies.

@fig-zscorecriteria shows the value of $c$ as a function of sample size, when there is only one suspected anomaly.

```{r}
#| label: fig-zscorecriteria
#| warning: false
#| echo: false
#| fig-cap: "Thresholds $c$ for $|z|$ under Peirce's criterion (with one suspected anomaly) and Chauvenet's criterion, as functions of the sample size."
#| dependson: pspurious

threshold <- numeric(length(n))
for (i in seq_along(n)) {
  threshold[i] <- weird:::peirce_threshold(n[i])
}
tibble(n = n, c = threshold) |>
  ggplot(aes(x = n, y = c)) +
  geom_line() +
  geom_line(aes(y = qnorm(1 - 0.25 / n)), col = "red") +
  scale_x_log10(
    limits = c(3, 20000),
    breaks = breaks,
    labels = sprintf("%d", breaks),
    minor_breaks = NULL
  ) +
  labs(x = TeX("$n$ (sample size)"), y = "c (threshold)") +
  geom_label(aes(x = 1500, y = 3.95, label = "Peirce criterion")) +
  geom_label(
    aes(x = 2600, y = 3.4, label = "Chauvenet criterion"),
    col = "red"
  ) +
  coord_cartesian(xlim = c(3, 1.2e4), expand = FALSE)
```

The threshold increases with sample size $n$ to allow for the increasing likelihood of observations falling in the extreme tails of the distribution.

The `peirce_anomalies()` function returns a logical vector indicating which observations are anomalous under this criterion.
Let's apply it to the six examples.

```{r}
#| label: peirce_examples
#| code-fold: false
#| dependson: ["fig-n01b", "testdist-data"]

cricket_batting |>
  filter(Innings >= 20) |>
  filter(peirce_anomalies(Average))
oldfaithful |> filter(peirce_anomalies(duration))
n01 |>
  select(v1) |>
  filter(peirce_anomalies(v1))
n01b |> filter(peirce_anomalies(y))
t3 |> filter(peirce_anomalies(y))
chisq4 |> filter(peirce_anomalies(y))
```

```{r}
#| label: checkpeirce
#| include: false
#| dependson: check2

nanomalies <- nanomalies |>
  mutate(
    Peirce = c(
      cricket_batting |>
        filter(Innings >= 20) |>
        filter(peirce_anomalies(Average)) |>
        NROW(),
      oldfaithful |> filter(peirce_anomalies(duration)) |> NROW(),
      n01b |> filter(peirce_anomalies(y)) |> NROW(),
      n01 |> filter(peirce_anomalies(v1)) |> NROW(),
      t3 |> filter(peirce_anomalies(y)) |> NROW(),
      chisq4 |> filter(peirce_anomalies(y)) |> NROW()
    )
  )

stopifnot(nanomalies$Peirce == c(1, 2, 2, 1, 10, 6))
```

1. Don Bradman is correctly identified as the only anomaly in the cricket batting averages.
2. When applied to the Old Faithful eruption durations, it finds only the two most extreme durations.
3. Both anomalies have been correctly spotted amongst the 18 N(0,1) observations.
4. One spurious anomaly is identified in the 1000 N(0,1) observations.
5. Ten spurious anomalies have been identified in the 1000 $\text{t}_3$ observations.
6. Six spurious anomalies have been identified in the 1000 $\chi^2_4$ observations.

### Chauvenet's criterion {-}

Peirce's proposal was largely superseded by an alternative proposed by the astronomer William Chauvenet, which was much simpler to describe.
He suggested [@Chauvenet1863] replacing the threshold $c$ by the $(1 - 0.25/n)$ quantile from the standard Normal distribution.
This threshold is also shown in @fig-zscorecriteria.
A consequence of this choice is that Chauvenet's criterion will reject, on average, half an observation of genuine data from a Normal distribution regardless of the value of $n$.
However, for non-Normal data, there is no such guarantee, and many more genuine observations may be rejected.
Despite its flaws, the method is still widely used in some disciplines, especially engineering.

The `chauvenet_anomalies()` function can be used to implement this test.
For our six examples, it gives similar results to those above for `peirce_anomalies()` (with two additional spurious anomalies for the $\text{t}_3$ example, and one additional spurious anomaly for the $\chi^2_4$ example).

```{r}
#| label: checkchauvenet
#| include: false
#| dependson: checkpeirce

nanomalies <- nanomalies |>
  mutate(
    Chauvenet = c(
      cricket_batting |>
        filter(Innings >= 20) |>
        filter(chauvenet_anomalies(Average)) |>
        NROW(),
      oldfaithful |> filter(chauvenet_anomalies(duration)) |> NROW(),
      n01b |> filter(chauvenet_anomalies(y)) |> NROW(),
      n01 |> filter(chauvenet_anomalies(v1)) |> NROW(),
      t3 |> filter(chauvenet_anomalies(y)) |> NROW(),
      chisq4 |> filter(chauvenet_anomalies(y)) |> NROW()
    )
  )
stopifnot(nanomalies$Chauvenet == nanomalies$Peirce + c(0, 0, 0, 0, 2, 1))
```

## Grubbs' test {#sec-grubbs}

By the 1930s, hypothesis testing had been developed and the t distribution discovered, and both were applied to the identification of anomalies using $z$-scores.
Many tests were developed under different assumptions about the underlying distribution and what was assumed to be known [@Hawkins1980].
Only the two most widely used are described here.

Egon Pearson and C Chandra Sekar proposed [@Pearson1936] that an observation be considered an anomaly if its $z$-score is larger than a critical value based on the t distribution.
Later, this was modified by Frank Grubbs [@Grubbs1950] who proposed that anomalies be identified if $|z_i| > c_\alpha$ where
$$
  c_\alpha
  = \frac{(n - 1)t_{\alpha/(2n), n-2}}{\sqrt{n(n - 2 + t^2_{\alpha/(2n), n-2})}},
$$
and $t_{p,k}$ is the $(1 - p)$ quantile of the t distribution with $k$ degrees of freedom.

@fig-zscoretests shows the critical values of Grubbs' test at $\alpha = 0.05$, along with the corresponding threshold from @Chauvenet1863 for comparison.

```{r}
#| label: fig-zscoretests
#| echo: false
#| dependson: pspurious
#| fig-cap: Critical values for maximum $z$-score tests with $\alpha=0.05$.

alpha <- 0.05
df <- tibble(n = n) |>
  filter(n >= 3, n <= 10130) |>
  mutate(
    t1 = qt(1 - alpha / n, n - 2),
    t3 = qt(1 - alpha / (2 * n), n - 2),
    chauvenet = qnorm(1 - 0.25 / n),
    grubbs = (n - 1) / sqrt(n) * sqrt(t3^2 / (n - 2 + t3^2)),
  ) |>
  select(-t1, -t3) |>
  tidyr::pivot_longer(chauvenet:grubbs, values_to = "c", names_to = "test")
df |>
  ggplot(aes(x = n, y = c, col = test)) +
  geom_line() +
  scale_x_log10(
    breaks = breaks,
    labels = sprintf("%d", breaks),
    minor_breaks = NULL
  ) +
  labs(x = TeX("$n$ (sample size)"), y = TeX("$c_{\\alpha}$ (threshold)")) +
  guides(col = "none") +
  geom_label(
    aes(x = x, y = y, col = col, label = label, hjust = "left"),
    data = data.frame(
      x = rep(1.1e4, 2),
      y = c(4.6, 4.0),
      col = c("grubbs", "chauvenet"),
      label = c("Grubbs test", "Chauvenet criterion")
    )
  ) +
  coord_cartesian(xlim = c(4, 1e5))
```

From these, we can compute the probability of a spurious anomaly in a Normal distribution, using simulation rather than the approximation in @eq-spuriousz.
For Grubbs' test, this probability should be equal to $\alpha = 0.05$.

```{r}
#| label: fig-pzscoretests
#| echo: false
#| dependson: ["pspurious", "fig-zscoretests"]
#| fig-cap: True size of the tests (i.e., the probability of a spurious anomaly) for Normally distributed data with $\alpha=0.05$.

pzscoredf <- df |>
  mutate(pmax = pspurious(n, rnorm, c = c, rep = 1e6)) |>
  cache("pzscoredf")
pzscoredf |>
  ggplot(aes(x = n, y = pmax, group = test, col = test)) +
  geom_hline(aes(yintercept = 0.05), col = "grey") +
  geom_line() +
  labs(
    x = TeX("$n$ (sample size)"),
    y = "Probability of at least one spurious anomaly"
  ) +
  scale_y_continuous(
    breaks = c(0.05, seq(0.1, 0.8, by = 0.1)),
    minor_breaks = NULL
  ) +
  scale_x_log10(
    breaks = breaks,
    labels = sprintf("%d", breaks),
    minor_breaks = NULL
  ) +
  labs(
    y = "Probability of at least one spurious anomaly",
    x = TeX("Sample size ($n$)")
  ) +
  guides(col = "none") +
  theme(axis.text.x = element_text(size = 8)) +
  geom_label(
    aes(x = x, y = y, col = test, label = label, hjust = "left"),
    data = data.frame(
      x = rep(1.1e4, 2),
      y = c(0.05, 0.40),
      test = c("grubbs", "chauvenet"),
      label = c("Grubbs test", "Chauvenet criterion")
    )
  ) +
  coord_cartesian(xlim = c(4, 1e5))
```

@fig-pzscoretests shows that Grubbs' test has the correct size.
The Chauvenet criterion is much worse, finding far too many anomalies except in very small samples.

However, when we allow for different data distributions, Grubbs' test also gives poor results (@fig-pgrubbs), showing it is sensitive to the assumed data distribution.

```{r}
#| label: fig-pgrubbs
#| echo: false
#| dependson: ["pspurious", "fig-pzscoretests"]
#| fig-cap: True size of Grubbs' test (i.e., the probability of a spurious anomaly) for different data distributions with $\alpha=0.05$.

pgrubbsdf <- pzscoredf |>
  filter(test == "grubbs") |>
  mutate(
    Normal = pmax,
    ApproxNormal = 1 - pt(c, n - 1)^(2 * n),
    `t(3)` = pspurious(n, rt, df = 3, c = c),
    `Chi-squared(4)` = pspurious(n, rchisq, df = 4, c = c)
  ) |>
  select(-c, -test, -pmax) |>
  tidyr::pivot_longer(
    -n,
    values_to = "Probability",
    names_to = "Distribution"
  ) |>
  cache("pgrubbsdf")
pgrubbsdf |>
  filter(Distribution != "ApproxNormal") |>
  ggplot(aes(
    x = n,
    y = Probability,
    group = Distribution,
    col = Distribution
  )) +
  geom_hline(aes(yintercept = 0.05), col = "grey") +
  geom_line() +
  labs(
    x = TeX("$n$ (sample size)"),
    y = "Probability of at least one spurious anomaly"
  ) +
  scale_x_log10(
    limits = c(3, 2e4),
    breaks = breaks,
    labels = sprintf("%d", breaks),
    minor_breaks = NULL
  ) +
  labs(
    y = "Probability of at least one spurious anomaly",
    x = TeX("Sample size ($n$)")
  ) +
  theme(axis.text.x = element_text(size = 8)) +
  scale_colour_discrete(
    breaks = c("t(3)", "Chi-squared(4)", "Normal"),
    labels = c(
      TeX("$t_{3}$"),
      TeX("$\\chi^2_4$"),
      "N(0,1)"
    )
  ) +
  theme(legend.text = element_text(hjust = 0)) +
  theme(axis.text.x = element_text(size = 8)) +
  coord_cartesian(xlim = c(4, 7500))
```

### Examples

We can apply the test using the `grubbs_anomalies()` function to our six examples.

```{r}
#| label: grubbs-anomalies
#| code-fold: false
#| dependson: ["fig-n01b", "testdist-data"]

cricket_batting |>
  filter(Innings >= 20) |>
  filter(grubbs_anomalies(Average)) |>
  select(Player, Country, Average)
oldfaithful |>
  filter(grubbs_anomalies(duration))
n01b |>
  filter(grubbs_anomalies(y))
n01 |>
  filter(grubbs_anomalies(v1))
t3 |>
  filter(grubbs_anomalies(y))
chisq4 |>
  filter(grubbs_anomalies(y))
```

```{r}
#| label: checkgrubbs
#| include: false
#| dependson: checkchauvenet

nanomalies <- nanomalies |>
  mutate(
    Grubbs = c(
      cricket_batting |>
        filter(Innings >= 20) |>
        filter(grubbs_anomalies(Average)) |>
        NROW(),
      oldfaithful |> filter(grubbs_anomalies(duration)) |> NROW(),
      n01b |> filter(grubbs_anomalies(y)) |> NROW(),
      n01 |> filter(grubbs_anomalies(v1)) |> NROW(),
      t3 |> filter(grubbs_anomalies(y)) |> NROW(),
      chisq4 |> filter(grubbs_anomalies(y)) |> NROW()
    )
  )

stopifnot(nanomalies$Grubbs == c(1, 1, 0, 0, 6, 1))
```

The clear failure is with the $\text{t}_3$ example, which has no real anomalies.
A spurious anomaly is also detected in the $\chi^2_4$ example, and real anomalies are missed in the Old Faithful data and in `n01b`.

## Dixon's Q test {#sec-dixon}

The $z$-score tests above all measure the distance of an observation from the centre of the distribution.
A different approach looks instead at the gaps between consecutive order statistics (see @sec-gaps): if a value is anomalous, we expect a noticeably large gap between it and its nearest neighbour in the sorted data.
@Dixon1950 introduced the first widely used test of this type.

If $y_{(1)},\dots,y_{(n)}$ denote the ordered values of our sample, then Dixon's Q statistic is given by
$$
  Q = \frac{y_{(n)} - y_{(n-1)}}{y_{(n)} - y_{(1)}},
$$
the ratio of the difference between the two largest values to the range of the data.
If the largest value is the only anomaly, then $Q$ will take a larger value than expected.

The corresponding test for the minimum to be an anomaly uses $y_{(2)} - y_{(1)}$ in the numerator instead, the difference between the second smallest and minimum observations.
Both minimum and maximum values can be tested simultaneously using the two-sided test, where the numerator is the maximum of $y_{(n)} - y_{(n-1)}$ and $y_{(2)} - y_{(1)}$.

The test is flawed for two reasons.
First, if the two largest values are both anomalies of similar size, then $Q$ will be small and these anomalies will be missed.
Also, if both maximum and minimum values are anomalies, the denominator will be larger than expected, thereby reducing the size of $Q$.

As with the other tests considered here, the test assumes that the underlying data distribution is Normal, and anomalies are identified which appear inconsistent with that assumption.

Simulation can be used to compute the critical values for this test, assuming that the data come from a Normal distribution.
We can also use simulation to compute the probability of false positives, giving the results shown in @fig-spuriousdixon.
As with the other tests we have considered, Dixon's test is sensitive to the assumed data distribution, making it largely useless for real data analysis.

```{r}
#| label: fig-spuriousdixon
#| fig-cap: "True size of Dixon's test (i.e., the probability of a spurious anomaly) for different data distributions with $\alpha=0.05$."
#| dependson: "pspurious"
#| echo: false

# Function to compute probability of spurious anomaly for sample size n,
# distribution dist with parameters in ...
# rep is number of replicates
pspuriousdixon <- function(n, dist, rep = 1e6, ...) {
  if (length(n) == 1) {
    pv <- numeric(rep)
    for (j in seq(rep)) {
      pv[j] <- any(dixon_anomalies(dist(n, ...)))
    }
    return(mean(pv))
  } else {
    pv <- numeric(length(n))
    for (i in seq_along(n)) {
      # print(n[i])
      pv[i] <- pspuriousdixon(n = n[i], dist = dist, rep = rep, ...)
    }
    return(pv)
  }
}
dixondf <- tibble(n = c(3:10, seq(12, 26, by = 2), seq(30, 100, by = 5))) |>
  mutate(
    Normal = pspuriousdixon(n, rnorm),
    `t(3)` = pspuriousdixon(n, rt, df = 3),
    `Chi-squared(4)` = pspuriousdixon(n, rchisq, df = 4)
  ) |>
  cache("dixondf")
dixondf |>
  tidyr::pivot_longer(
    -n,
    values_to = "Probability",
    names_to = "Distribution"
  ) |>
  ggplot(aes(
    x = n,
    y = Probability,
    group = Distribution,
    col = Distribution
  )) +
  geom_hline(aes(yintercept = 0.05), col = "grey") +
  geom_line() +
  scale_y_continuous(
    breaks = c(0.05, seq(0.1, 0.8, by = 0.1)),
    minor_breaks = NULL
  ) +
  labs(
    y = "Probability of at least one spurious anomaly",
    x = TeX("Sample size ($n$)")
  ) +
  scale_colour_discrete(
    breaks = c("t(3)", "Chi-squared(4)", "Normal"),
    labels = c(
      TeX("$t_{3}$"),
      TeX("$\\chi^2_4$"),
      "N(0,1)"
    )
  ) +
  theme(legend.text = element_text(hjust = 0)) +
  theme(axis.text.x = element_text(size = 8)) +
  coord_cartesian(xlim = c(4, 96))
```

### Examples

We can apply the test using the `dixon_anomalies()` function, which does a two-sided test of both minimum and maximum observations.

```{r}
#| label: dixon-anomalies
#| code-fold: false
#| dependson: ["fig-n01b", "testdist-data"]

cricket_batting |>
  filter(Innings >= 20) |>
  filter(dixon_anomalies(Average)) |>
  select(Player, Country, Average)
oldfaithful |>
  filter(dixon_anomalies(duration))
n01b |>
  filter(dixon_anomalies(y))
n01 |>
  filter(dixon_anomalies(v1))
t3 |>
  filter(dixon_anomalies(y))
chisq4 |>
  filter(dixon_anomalies(y))
```

```{r}
#| label: checkdixon
#| include: false
#| dependson: checkgrubbs

nanomalies <- nanomalies |>
  mutate(
    Dixon = c(
      cricket_batting |>
        filter(Innings >= 20) |>
        filter(dixon_anomalies(Average)) |>
        NROW(),
      oldfaithful |> filter(dixon_anomalies(duration)) |> NROW(),
      n01b |> filter(dixon_anomalies(y)) |> NROW(),
      n01 |> filter(dixon_anomalies(v1)) |> NROW(),
      t3 |> filter(dixon_anomalies(y)) |> NROW(),
      chisq4 |> filter(dixon_anomalies(y)) |> NROW()
    )
  )

stopifnot(nanomalies$Dixon == c(1, 0, 0, 0, 1, 0))
```

For these examples, the test has failed to identify any of the anomalies in `oldfaithful` and `n01b`, and it has identified a spurious anomaly in `t3`.

## Gap-based methods {#sec-gap-tests}

Dixon's test only examines the gap adjacent to the minimum or maximum observation.
@Tukey1977-xd proposed looking for large gaps anywhere in the sorted data, while @Wainer1978 formalised this into a practical method for determining whether any such spacing is unexpectedly large for its position in the distribution.
The method is best thought of as detecting unusually large **separations** or **holes** in a one-dimensional data set.
A large tail gap may suggest an outlier; a large interior gap usually suggests a separation between two groups rather than a single anomalous observation.

### Tukey--Wainer--Schacht gap scores {-}

Consider the spacings defined (@sec-gaps) by
$$
  d_i = y_{(i+1)} - y_{(i)}, \qquad i = 1, \dots, n - 1,
$$
where $y_{(1)} \le y_{(2)} \le \cdots \le y_{(n)}$ are the sorted observations.

Raw spacings are not directly comparable across the sorted sample.
Even for data from a single unimodal distribution, spacings tend to be small near the centre of the distribution, where observations are dense, and larger in the tails, where observations are sparse.

Tukey's **gapping** scheme addresses this by weighting each gap according to its position in the ordered data.
Specifically, the gap between $y_{(i)}$ and $y_{(i+1)}$ is assigned weight
$$
  a_i = i(n - i), \qquad i = 1, \ldots, n - 1,
$$
which is largest near the middle of the sample and smallest at the extremes.
This emphasises unusually large gaps in the centre, where observations are usually close together, and downweights gaps in the tails, where larger spacings are expected.

@Wainer1978 developed this idea into a practical diagnostic by transforming and standardising the weighted gaps.
They proposed using the square root of the weighted gaps,
$$
  h_i = \sqrt{a_i d_i},
$$
which makes their distribution more symmetric, and then comparing each $h_i$ with a robust estimate of typical gap size, namely the 25% trimmed mean of the $h_i$ values (the mean of the middle half of the values):
$$
  z^*_i = \frac{h_i}{\bar h_{T25}}.
$$
Large values of $z^*_i$ identify gaps that are large relative to other gaps after allowing for their position in the sorted data.
Tukey suggested $z^*_i > 2.25$ as an exploratory cutoff; @Wainer1978 showed by simulation that, under a Normal distribution, this corresponds to a per-gap false positive rate of about 3%.
Because there are $n - 1$ gaps, this 3% rate does not apply to the data set as a whole: if many gaps are inspected, the chance of at least one crossing the threshold can be much larger than the per-gap rate.

The code below constructs a data frame of gaps between batter averages, and plots their scores in @fig-ws-gaps.

```{r}
#| label: ws-gaps
#| code-fold: false

batters <- cricket_batting |>
  filter(Innings >= 20) |>
  arrange(Average)
batting_gaps <- batters |>
  transmute(
    lower_player = Player,
    lower_average = Average,
    upper_player = lead(Player),
    upper_average = lead(Average),
    spacing = upper_average - lower_average,
    weight = row_number() * (n() - row_number())
  ) |>
  mutate(
    weighted_gap = sqrt(weight * spacing),
    ws = weighted_gap / mean(weighted_gap, trim = 0.25, na.rm = TRUE)
  ) |>
  filter(!is.na(ws))
```

```{r}
#| label: fig-ws-gaps
#| dependson: ws-gaps
#| fig-cap: >-
#|   Tukey--Wainer--Schacht gap scores for cricket batting averages (batters with at least 20 innings).
#|   Each point is plotted at the upper end of a gap, so the score shown at a batting average is for the gap between that batter and the next lower batting average.
#|   The dashed line marks Tukey's exploratory threshold of 2.25.

batting_gaps |>
  ggplot(aes(x = upper_average, y = ws)) +
  geom_hline(yintercept = 2.25, linetype = "dashed", colour = "gray50") +
  geom_point(alpha = 0.5) +
  labs(x = "Career batting average", y = "Tukey–Wainer–Schacht gap score")
```

Restricting attention to scores above the exploratory threshold of 2.25, we show the five largest gaps below.

```{r}
#| label: ws-top-gaps
#| code-fold: false
#| dependson: ws-gaps

batting_gaps |>
  filter(ws > 2.25) |>
  slice_max(spacing, n = 5) |>
  select(lower_player, upper_player, spacing, weight, ws)
```

The largest gap of `r sprintf("%.1f", batting_gaps |> filter(ws > 2.25) |> slice_max(spacing, n = 1) |> pull(spacing))` should be interpreted as the separation between Kamindu Mendis (who, at the time the data were collected, had the second highest batting average among batters with at least 20 innings) and Don Bradman, not as a score belonging to Bradman alone.
This illustrates that the method detects separations, not individual anomalous observations.

```{r}
#| label: ws-noise
#| include: false
#| dependson: ws-gaps

ws_flagged <- batting_gaps |> filter(ws > 2.25)
bradman_gap <- batting_gaps |> filter(upper_player == "DG Bradman")
bradman_rank <- sum(batting_gaps$ws > bradman_gap$ws) + 1
stopifnot(
  NROW(ws_flagged) / NROW(batting_gaps) > 0.03,
  median(ws_flagged$spacing) < 0.5,
  all(batting_gaps |> slice_max(ws, n = 5) |> pull(spacing) < 1),
  bradman_rank > 5
)
```

The example also shows the main weakness of the method.
In total, `r NROW(ws_flagged)` of the `r NROW(batting_gaps)` gaps (`r sprintf("%.1f", 100 * NROW(ws_flagged) / NROW(batting_gaps))`%) have scores above 2.25, more than the 3% expected for Normal data, and almost all of them are gaps of a fraction of a run between batters with unremarkable averages, in the middle of the sample where the weights $a_i$ are largest.
The gap below Bradman, which is by far the largest, has a score of only `r sprintf("%.2f", bradman_gap$ws)` (the `ws` column in the first row of the table above), ranking `r english::ordinal(bradman_rank)` among all gaps, because gaps at the extremes are given so little weight.
We found it above only because we sorted the flagged gaps by their raw size; the scores themselves are dominated by small gaps in the middle of the sample.

```{r}
#| label: check-second-batter
#| include: false
#| dependson: ws-gaps

# The prose names the batter below Bradman; fail if a data update changes this
stopifnot(
  batting_gaps |>
    filter(ws > 2.25) |>
    slice_max(spacing, n = 1) |>
    pull(lower_player) ==
    "PHKD Mendis"
)
```

### Significance testing {-}

The Tukey--Wainer--Schacht score does not yield a direct $p$-value.
A related way to construct a formal gap test is to calculate the uniform spacings (see @sec-gaps) of the data after transforming them using the fitted distribution function $F$:
$$
  \delta_i = F(y_{(i+1)}) - F(y_{(i)}), \quad i = 1, \dots, n - 1.
$$
Since the transformed spacings are approximately $\text{Exp}(n + 1)$ (@sec-gaps), and the Exponential distribution is in the Gumbel domain of attraction (@sec-evt), $\delta_{\max} = \max(\delta_1,\dots,\delta_{n-1})$ satisfies
$$
  t = (n + 1)\delta_{\max} - \log(n - 1) \approx \text{Gumbel}(0,1).
$$
This gives the approximate upper-tail probability
$$
  p \approx 1 - \exp(-e^{-t}).
$$

In practice, $F$ is unknown.
The example below uses a robustly fitted Normal distribution, with the median and MAD as estimates of location and scale.
This makes the fit less sensitive to the possible outlier, but the resulting test is still only an approximate parametric test.
The MAD is used rather than $Q_n$ (@sec-robust-univariate) because the test is very sensitive to the scale estimate: with the slightly smaller $Q_n$, the gap below Bradman is no longer significant.

```{r}
#| label: gap-bt-test
#| code-fold: false
#| dependson: ws-gaps

mu_hat <- median(batters$Average)
sigma_hat <- robustbase::s_mad(batters$Average)

batting_gaps |>
  mutate(
    transformed_gap = pnorm(upper_average, mu_hat, sigma_hat) -
      pnorm(lower_average, mu_hat, sigma_hat),
    t = (nrow(batters) + 1) * transformed_gap - log(nrow(batting_gaps)),
    p_value = 1 - exp(-exp(-t))
  ) |>
  filter(p_value < 0.05) |>
  select(
    lower_player,
    lower_average,
    upper_player,
    upper_average,
    spacing,
    p_value
  )
```

```{r}
#| label: gap-qn-check
#| include: false
#| dependson: gap-bt-test

# The prose says Qn is smaller than the MAD here, and that with Qn the gap below Bradman is not significant
sigma_qn <- robustbase::s_Qn(batters$Average)
qn_gaps <- batting_gaps |>
  mutate(
    t = (nrow(batters) + 1) *
      (pnorm(upper_average, mu_hat, sigma_qn) -
        pnorm(lower_average, mu_hat, sigma_qn)) -
      log(nrow(batting_gaps)),
    p_value = 1 - exp(-exp(-t))
  )
stopifnot(sigma_qn < sigma_hat, all(qn_gaps$p_value >= 0.05))
```

Only the gap below Bradman's average has a $p$-value below 0.05, suggesting that this gap is too large to be consistent with all batters coming from a single fitted Normal distribution.

Finally, we apply the same approach to all six examples considered in this chapter.
Here we associate the $p$-value with the larger observation in each gap pair.

```{r}
#| label: gap-tests
#| code-fold: false
#| dependson: ["fig-n01b", "testdist-data"]

gap_test_pvalue <- function(x) {
  mu_hat <- median(x)
  sigma_hat <- robustbase::s_mad(x)
  n <- length(x)
  ix <- order(x)
  y <- x[ix]
  lower <- y[-n]
  upper <- y[-1]
  transformed_gap <- pnorm(upper, mu_hat, sigma_hat) -
    pnorm(lower, mu_hat, sigma_hat)
  t <- (n + 1) * transformed_gap - log(n - 1)
  p_value <- 1 - exp(-exp(-t))

  p_value_original_order <- rep(NA_real_, n)
  p_value_original_order[ix[-1]] <- p_value
  p_value_original_order
}
cricket_batting |>
  filter(Innings >= 20) |>
  filter(gap_test_pvalue(Average) < 0.05) |>
  select(Player, Country, Average)
oldfaithful |>
  filter(gap_test_pvalue(duration) < 0.05)
n01b |>
  filter(gap_test_pvalue(y) < 0.05)
n01 |>
  filter(gap_test_pvalue(v1) < 0.05)
t3 |>
  filter(gap_test_pvalue(y) < 0.05)
chisq4 |>
  filter(gap_test_pvalue(y) < 0.05)
```

```{r}
#| label: checkgap
#| include: false
#| dependson: ["checkdixon", "gap-tests"]

nanomalies$Gap <- c(
  cricket_batting |>
    filter(Innings >= 20) |>
    filter(gap_test_pvalue(Average) < 0.05) |>
    NROW(),
  oldfaithful |> filter(gap_test_pvalue(duration) < 0.05) |> NROW(),
  n01b |> filter(gap_test_pvalue(y) < 0.05) |> NROW(),
  n01 |> filter(gap_test_pvalue(v1) < 0.05) |> NROW(),
  t3 |> filter(gap_test_pvalue(y) < 0.05) |> NROW(),
  chisq4 |> filter(gap_test_pvalue(y) < 0.05) |> NROW()
)
stopifnot(nanomalies$Gap == c(1, 68, 0, 0, 0, 0))
nof <- nanomalies$Gap[2]
# The prose says the Old Faithful flags are every whole second in a central range
of_flagged <- oldfaithful |>
  filter(gap_test_pvalue(duration) < 0.05) |>
  pull(duration) |>
  sort()
stopifnot(
  all(of_flagged == round(of_flagged)),
  all(diff(of_flagged) == 1),
  min(of_flagged) == 207,
  max(of_flagged) == 274
)
```

While the method worked well for the cricket batting averages, it identified `r nof` "anomalies" in the Old Faithful data, and it failed to identify either of the two clear anomalies in `n01b`.
The Old Faithful "anomalies" are every whole-second duration from `r min(of_flagged)` to `r max(of_flagged)` seconds; because the durations are rounded to the nearest second, the spacings near the centre of the data are much larger than the test expects for continuous data.

Like all methods in this chapter, both the Tukey--Wainer--Schacht gap score and the transformed-gap test depend on an assumed null distribution, and so the results depend on the adequacy of the model.
Related extreme-spacing ideas underlie the HDoutliers and stray algorithms discussed in @sec-hdoutliers and @sec-stray, where gaps in nearest-neighbour distances are used rather than gaps in the raw univariate data.

## Summary {#sec-tests-summary}

@tbl-nanomalies summarises the results of the tests on the six examples.

```{r}
#| label: tbl-nanomalies
#| echo: false
#| tbl-cap: >-
#|   Number of anomalies detected by each test.
#|   Cells in red indicate that the test has failed to identify all known anomalies (i.e., it has some false negatives), or it has identified some spurious anomalies (i.e., it has some false positives).
#| dependson: checkgap
#| cache: false

library(kableExtra)
nanomalies_table <- nanomalies |>
  mutate(
    Example = case_when(
      row_number() == 5 ~ "<span class='math inline'>\\text{t}_3</span>",
      row_number() == 6 ~ "<span class='math inline'>\\chi^2_4</span>",
      .default = Example
    ),
    Example = paste0(row_number(), ". ", Example)
  )
method_cols <- seq.int(4, ncol(nanomalies_table))
method_colours <- nanomalies[method_cols] |>
  purrr::map(\(x) {
    if_else(x != nanomalies$`Expected anomalies`, "red", "black")
  })
table <- nanomalies_table |>
  knitr::kable(escape = FALSE, format = "html") %>%
  add_header_above(c(" " = 3, "Anomaly detection method" = length(method_cols)))
purrr::reduce2(
  .x = method_cols,
  .y = method_colours,
  .init = table,
  .f = \(table, col, colour) column_spec(table, col, color = colour)
)
```

Every test failed on at least one of the examples in @tbl-nanomalies.
Several anomalies can hide one another: they inflate the standard deviation used to compute $z$-scores, and two anomalies of similar size leave no large gap for Dixon's test to find.
The $z$-score-based tests and Dixon's test also flag spurious anomalies in large samples that contain none, especially when the data are heavy-tailed or skewed.
The gap test is misled by data rounded to a fixed precision, because the spacings between distinct values are then larger than it expects for continuous data.

As noted at the start of this chapter, anomaly detection methods based on assumed data distributions are unreliable.
Until about 1975, they were the only viable methods given the lack of computing facilities available, because they could be implemented using tables and hand calculations.
However, there is little justification for continuing to use such methods.

They are particularly sensitive to the assumed data distribution, and the probability of detecting spurious anomalies is usually much higher in reality than under the ideal conditions in which the tests were conceived.

These tests are sometimes applied iteratively, where observations are removed from the data if determined to be anomalies, and the test re-applied to the remaining data.
This process continues until no more anomalies are found.
However, this procedure will clearly change the size of the test due to the problem of multiple comparisons.

Next, we will turn our attention to methods, developed in the latter part of the 20th century, that are based on quantiles and data depth rather than on any underlying parametric data distribution.