Most readers will have first come across anomaly detection using boxplots, proposed by John Tukey as a quick summary of medium-sized data sets (Kafadar, 2026; Tukey, 1975, 1977; Wickham & Stryjewski, 2011). Tukey never proposed that boxplots be a generic anomaly detection method, but they have become widely used for this purpose, and are often the first method that people think of when they want to identify anomalies in univariate data.
This chapter describes the original boxplot method, along with some variations that have been developed to address some of the limitations of the original approach.
For univariate data, the depth of an observation is a measure of how deeply buried it is when all observations are ordered. That is, how far would you need to count from either the smallest or largest observation until you encountered the observation of interest. So the minimum and maximum both have depth 1, while the sample median has the largest depth of (n + 1)/2.
The other sample quantiles as defined in Section 2.6 do not correspond exactly to specific depths, but Tukey (1977) introduced variations for some quantiles that are based on depths. He called these “letter values”.
Letter values (Hoaglin, 1983; Tukey, 1977) are order statistics with specific depths, defined recursively starting with the median. The depth of the median is d_1 = (1 + n)/2. The depths of successive letter values are defined recursively as d_i = (1 + \lfloor d_{i-1}\rfloor)/2, i = 2,3,\dots. The corresponding letter values are defined as L_i = y_{(\lfloor d_i\rfloor)} \qquad\text{and} \qquad U_i = y_{(\lfloor n-d_i+1\rfloor)}, when the depth is an integer. Otherwise, the depth is an integer plus 1/2, and the letter values are given by L_i = (y_{(\lfloor d_i\rfloor)} + y_{(\lfloor d_i\rfloor+1)})/2 \qquad\text{and} \qquad U_i = (y_{(\lfloor n-d_i+1\rfloor)} + y_{(\lfloor n-d_i+1\rfloor+1)})/2 . Rather than label these using integers (L_2,L_3,\dots), Tukey proposed using letters (L_F,L_E,L_D,\dots) where F = fourths, E = eighths, D = sixteenths, and so on.
Because each depth is roughly half the previous depth, the lower letter values provide estimates of the quantiles with probabilities p = \frac{1}{2},\frac14,\frac18,\dots, while the upper letter values provide estimates of the quantiles with probabilities p = \frac{1}{2},\frac34,\frac78,\dots. Hoaglin (1983, p44) showed that \hat{Q}(p) with type = 8 gives approximately the same result as the corresponding letter value.
Consider the batting averages from Section 2.6. In this example, n = 1254, so the depths of the first four letter values are given by
d_1 = 627.5,\quad d_2 = 314,
\quad d_3 = 157.5,\quad\text{and}\quad d_4 = 79,
and the corresponding letter values are given by \begin{align*}
L_1 & = U_1 = (y_{(627)} + y_{(628)})/2 = 26.43 &\\
L_F & = y_{(314)} = 16.19 \qquad
& U_F & = y_{(941)} = 36.68 \\
L_E & = (y_{(157)} + y_{(158)})/2 = 11.26 \qquad
& U_E & = (y_{(1097)} + y_{(1098)})/2 = 43.08 \\
L_D & = y_{(79)} = 7.79 \qquad
& U_D & = y_{(1176)} = 46.86
\end{align*} We can compute the letter values using the lvtable() function from the lvplot package.
batave <- cricket_batting |>
filter(Innings >= 20) |>
pull(Average)
batave |> lvplot::lvtable(k = 4) depth LV 2.5% 97.5%
Dl 79.0 7.792 7.176 8.286
El 157.5 11.257 10.404 11.682
Fl 314.0 16.192 15.435 17.130
M 627.5 26.427 25.569 27.414
Fu 314.0 36.682 35.696 37.473
Eu 157.5 43.085 42.218 44.075
Du 79.0 46.857 46.379 48.227The output also provides 95% confidence intervals for each letter value. These are formed from the order statistics either side of each depth, using a Normal approximation for the depth of each letter value, so (unlike Equation 2.1) they do not require an estimate of the density. The estimates are similar, but not identical, to the quantiles calculated using the quantile() function.
batave |> quantile(prob = c(0.5^(4:1), 1 - 0.5^(2:4)), type = 8) 6.25% 12.5% 25% 50% 75% 87.5% 93.75%
7.74 11.23 16.19 26.43 36.68 43.10 46.94 Figure 5.1 shows a boxplot of the cricket batting average data.
cricket_batting |>
filter(Innings >= 20) |>
ggplot(aes(x = Average)) +
geom_boxplot() +
scale_y_discrete() +
labs(y = "", x = "Career batting average")
In Tukey’s original version, the ends of the box were determined by the letter values corresponding to the fourths, and the middle line within the box showed the median. This is what is produced by the base R function boxplot(). In the ggplot2 version of the boxplot shown here, the ends of the box are the quartiles computed using type = 7.
Whichever variation is used, roughly half of all observations lie within the box. The length of the box is an estimate of the interquartile range (IQR). The fences lie 1.5 IQR beyond each end of the box, and any points outside the fences are shown as separate points in the plot. The “whiskers” that extend out each side of the box show the range of the remaining points.
Originally, Tukey proposed two levels of outliers — those more than 1.5 IQR beyond the box were labelled “outside” values, while those more than 3 IQR beyond the box were labelled “far out” values. Most software implementations of boxplots do not distinguish between these groups. When a boxplot is used for anomaly detection, all points outside the whiskers are called “anomalies” or “outliers”, even though they may not be genuine anomalies.
In this cricket batting example, the boxplot works well because the data set is not too large or small, and the distribution of points other than the anomaly is unimodal.
However, boxplots can be misleading, and they are limited in at least two respects.
To better understand the first problem, imagine if the data comprised n observations from a standard Normal distribution N(0,1). Figure 5.2 shows an example with 10000 points.
tibble(x = rnorm(10000)) |>
ggplot(aes(x = x)) +
geom_boxplot() +
scale_y_discrete() +
labs(y = "")
Many “anomalies” are shown, but since all observations come from a simple distribution, none of them are actually anomalies. For this distribution, Q(0.25) = -0.674, Q(0.75) = 0.674, \text{IQR} = 1.349, and so for a large sample size, the boxplot whiskers would be approximately -0.674 - 1.5\times 1.349 = -2.698 and 0.674 + 1.5\times 1.349 = 2.698. Any points outside the whiskers would be identified as “anomalous” by a boxplot and plotted as separate points. The probability of a standard Normal observation being at least 2.698 in absolute value (and therefore lying outside the boxplot whiskers) is 0.00698. So with 10000 observations, we would have about 70 anomalies identified, none of which would be a genuine anomaly.
The second problem is demonstrated using the Old Faithful eruption duration data (see Figure 1.3), with a boxplot shown in Figure 5.3.
oldfaithful |>
ggplot(aes(x = duration)) +
geom_boxplot() +
scale_y_discrete() +
labs(y = "", x = "Duration (seconds)")
Quite a few “anomalies” are shown, including the three abnormally short eruptions we identified earlier. But the remaining “anomalies” are not particularly unusual observations. All the points below 183 seconds are identified as anomalies, even though we know that observations in the region between 100 and 140 seconds are not unusual for this geyser. Because the boxplot does not allow for more than one mode, all the points in the second smaller cluster are identified as anomalies. The points around 300 seconds are also not really anomalies — these are just values in the upper tail of the distribution for eruptions.
Under Tukey’s boxplot approach to identifying outliers, a regular outlier (an “outside” value) is more than 1.5 IQR beyond the quartiles, while an extreme outlier (a “far out” value) is more than 3 IQR beyond the quartiles. For a Normal distribution, the probability of genuine observations lying beyond these thresholds is 0.0070 and 0.0000023 respectively, so for large sample sizes, many spurious anomalies will be identified. Even with 1000 observations, Tukey’s approach will find at least one spurious anomaly in a Normal distribution with probability 0.9991.
Barbato et al. (2011) proposed a modification to the boxplot approach to identifying outliers, where the IQR in these thresholds is replaced with IQR[1 + 0.1\log(n/10)]. This allows the limits to increase with the sample size, in a way that controls the probability of spurious anomalies. Figure 5.4 shows the probability of identifying at least one spurious anomaly in a Normal distribution, using Tukey’s boxplot approach compared to those obtained using the modified IQR approach of Barbato et al. (2011).
tibble(n = exp(seq(log(3), log(1e6), l = 100))) |>
mutate(
Tukey1 = 2 * (1 - pnorm(q3 + 1.5 * iqr)),
Tukey2 = 2 * (1 - pnorm(q3 + 3 * iqr)),
Barbato1 = 2 * (1 - pnorm(q3 + 1.5 * iqr * (1 + 0.1 * log(n / 10)))),
Barbato2 = 2 * (1 - pnorm(q3 + 3 * iqr * (1 + 0.1 * log(n / 10))))
) |>
tidyr::pivot_longer(
Tukey1:Barbato2,
names_to = "method",
values_to = "probability"
) |>
mutate(
level = stringr::str_extract(method, "\\d"),
level = if_else(level == "1", "Regular outlier", "Extreme outlier"),
method = stringr::str_extract(method, "[A-Za-z]*"),
level = factor(level, levels = c("Regular outlier", "Extreme outlier")),
method = factor(method, levels = c("Tukey", "Barbato")),
probability = 1 - (1 - probability)^n
) |>
ggplot(aes(x = n, y = probability)) +
geom_line() +
facet_grid(level ~ method) +
labs(y = "Probability of at least one spurious anomaly", x = "Sample size") +
scale_x_log10(
limits = c(3, 2e6),
breaks = 10^(1:6),
minor_breaks = NULL,
labels = format(10^(1:6), scientific = FALSE, trim = TRUE)
)
Even with a huge sample size, the probability of identifying a spurious regular anomaly using this modified approach is less than 1/2, and it is almost impossible to identify a spurious extreme anomaly.
Let’s apply this approach to the six examples introduced in Section 4.1. First we will write a short function to implement the idea.
barbato_anomaly <- function(y, extreme = FALSE) {
n <- length(y)
q1 <- quantile(y, 0.25, na.rm = TRUE, type = 8)
q3 <- quantile(y, 0.75, na.rm = TRUE, type = 8)
threshold <- (1.5 + 1.5 * extreme) * (q3 - q1) * (1 + 0.1 * log(n / 10))
return(y > q3 + threshold | y < q1 - threshold)
}
cricket_batting |>
filter(Innings >= 20) |>
filter(barbato_anomaly(Average))# A tibble: 1 × 15
Player Country Start End Matches Innings NotOuts Runs HighScore HighScoreNotOut
<chr> <chr> <int> <int> <int> <int> <int> <int> <int> <lgl>
1 DG Bradman Austra… 1928 1948 52 80 10 6996 334 FALSE
# ℹ 5 more variables: Average <dbl>, Hundreds <int>, Fifties <int>, Ducks <int>,
# Gender <chr>oldfaithful |> filter(barbato_anomaly(duration))# A tibble: 273 × 4
time recorded_duration duration waiting
<dttm> <chr> <dbl> <dbl>
1 2017-01-27 23:57:00 2m 1s 121 3900
2 2017-02-05 20:41:00 2m 3s 123 4080
3 2017-02-07 16:36:00 1m 56s 116 4260
4 2017-02-10 14:17:00 2m 5s 125 4140
5 2017-02-14 22:36:00 <2m 120 4140
6 2017-02-19 17:40:00 2m 5s 125 3840
7 2017-03-07 15:20:00 ~1.75m 105 3540
8 2017-03-22 11:49:00 ~1.75m 105 3540
9 2017-05-03 06:20:00 <2m 120 4680
10 2017-05-14 11:39:00 ~110s 110 3540
# ℹ 263 more rowsn01b <- tibble(y = c(n01$v2[1:18], 4, 4.5))
n01b |> filter(barbato_anomaly(y))# A tibble: 1 × 1
y
<dbl>
1 4.5n01 |> filter(barbato_anomaly(v1))# A tibble: 1 × 10
v1 v2 v3 v4 v5 v6 v7 v8 v9 v10
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 3.81 1.34 1.22 1.14 0.386 -0.158 1.88 -0.639 -0.636 1.96set.seed(1)
t3 <- tibble(y = rt(1000, df = 3))
t3 |> filter(barbato_anomaly(y))# A tibble: 23 × 1
y
<dbl>
1 -7.40
2 5.53
3 10.5
4 -5.82
5 -5.23
6 4.32
7 6.20
8 4.88
9 -4.11
10 -4.35
# ℹ 13 more rowschisq4 <- tibble(y = rchisq(1000, df = 4))
chisq4 |> filter(barbato_anomaly(y))# A tibble: 13 × 1
y
<dbl>
1 14.4
2 14.1
3 13.0
4 13.5
5 13.7
6 13.3
7 17.0
8 13.7
9 12.9
10 14.0
11 15.2
12 14.1
13 13.9This is an improvement on the original boxplot approach, and it finds the true anomaly in the cricket_batting data and one of the two anomalies in the n01b data. But it also identifies false anomalies in the other data sets: only one in the n01 data, but many in those that are not Normally distributed. The assumption of Normality means that the method is not a good choice for most data sets.
The problem that boxplots have with large data sets was also addressed by Hofmann et al. (2017) who introduced letter-value plots, a variation of boxplots that replaces the whiskers with a variable number of letter values. In these plots, each pair of letter values marks the boundaries of a box. The box bounded by the fourths is the same as the box of a boxplot; the additional boxes extend to successive letter values until the quantiles corresponding to the letter values can no longer be estimated sufficiently accurately from the available data.
These can be produced using the lvplot package, as shown in Figure 5.5.
library(lvplot)
cricket_batting |>
filter(Innings >= 20) |>
ggplot(aes(x = 1, y = Average)) +
geom_lv(aes(fill = after_stat(LV))) +
scale_x_discrete() +
coord_flip() +
labs(x = "", y = "Career batting average") +
theme(legend.key.height = unit(0.2, "cm"))
Here the median is given by M, the fourths by F, and so on. The middle box (F) is bounded by the fourths and contains all but 2/4 of the data; the next box (E) is bounded by the eighths and contains all but 2/8 of the data; then D is bounded by the sixteenths and contains all but 2/16 of the data; and so on.
In this example, the most extreme box (labelled Z) is bounded by the 1/256 letter values. So it contains all but 2/256 of the data, and shows 2 / 256 \times 1254 \approx 10 points individually.
The stopping rule used in the letter-value plot is to show the boxes up to letter value k, where k is the smallest value satisfying 0.5\sqrt{2d_k} z_{1-\alpha/2} > d_{k+1} and z_{1-\alpha/2} is the (1 - \alpha/2) quantile of a standard Normal distribution. This choice is based on the idea that the edges of the boxes are quantile estimates, and the confidence interval for each quantile estimate that is displayed should not overlap the subsequent quantile estimate. The Normal distribution arises because, for large n, the depth (rank) at which a sample quantile falls is approximately Normally distributed; this is the rank-based counterpart of Equation 2.1. This stopping rule means that, on average, there should be no more than about 2z^2_{1-\alpha/2} legitimate observations in each tail. By default, \alpha = 0.05, so that, on average, there should be no more than about 2 \times (1.96)^2 = 7.7 legitimate observations in each tail, regardless of the size of the data set.
Letter-value plots were not designed to detect anomalies, but to be a useful data visualisation tool for univariate distributions with large numbers of observations. So the display of legitimate observations in the tails of the distribution is by design, not a flaw.
In this cricketing example, it looks like there is one true anomaly (Don Bradman) and the remaining 9 observations displayed directly are simply in the tails of the distribution of the remaining data.
When applied to three of the other examples (Figure 5.6, Figure 5.7 and Figure 5.8), we see between 8 and 16 observations shown as individual points in each case.
oldfaithful |>
ggplot(aes(x = 1, y = duration)) +
geom_lv(aes(fill = after_stat(LV))) +
scale_x_discrete() +
coord_flip() +
labs(x = "", y = "Eruption durations (seconds)") +
theme(legend.key.height = unit(0.2, "cm"))
n01 |>
ggplot(aes(x = 1, y = v1)) +
geom_lv(aes(fill = after_stat(LV))) +
scale_x_discrete() +
coord_flip() +
labs(x = "") +
theme(legend.key.height = unit(0.2, "cm"))
n01b |>
ggplot(aes(x = 1, y = y)) +
geom_lv(aes(fill = after_stat(LV))) +
scale_x_discrete() +
coord_flip() +
labs(x = "")
In this last example, because there are only 20 observations, there are not enough data to estimate the quantiles beyond the fourths. So only the middle box is shown, and the 10 observations outside it are all plotted individually.
For multivariate data, there is no unique natural ordering of observations by size, so the concept of depth has to be thought about differently. It is still a measure of how centrally a point is located in a data set, but we need to think about what “central” means when there are multiple variables.
There are numerous ways to define depth for multivariate data, nicely summarised in Liu et al. (1999), and implemented in the depthR package for R. Only the simplest and oldest of these is considered here: Tukey depth, which is the basis of the bagplot (Section 5.6).
The first approach to this problem was (again) due to John Tukey (Tukey, 1975) who defined the depth of a point \bm{z} in a data set \{\bm{y}_1,\dots,\bm{y}_n\} as the smallest number of \bm{y}_i contained in any halfspace that contains \bm{z}. This is often called the Tukey depth of the point. (The point \bm{z} does not have to be one of the observations.)
For bivariate data, a halfspace is either of the two parts formed by splitting the plane with a straight line. A diagram will help illustrate the idea. Suppose we have the observations shown in Figure 5.9. These are the first 10 observations from the first two variables of n01. The red line divides the plane into two sections which are called halfspaces. We are interested in the depth of the orange point (which is not one of the observations).
n01 |>
head(10) |>
ggplot(aes(x = v1, y = v2)) +
geom_point() +
geom_abline(aes(intercept = -0.8, slope = -1.8), col = "red") +
geom_point(data = data.frame(v1 = 0, v2 = -1), col = "#D55E00") +
coord_fixed(xlim = c(-0.85, 1.6), ylim = c(-1.85, 1.1))
There are an infinite number of ways of dividing the plane into halfspaces, and the number of points in each halfspace will vary depending on where the dividing line falls.
To find the depth of the orange point, we need to find the line which divides the plane into two sections where the half containing the orange point includes as few observations as possible. In fact, the red line is one such solution which has the orange point in a halfspace containing only two observations. There are no dividing lines that would put the orange point in a halfspace on its own, or in a halfspace containing only one observation. So it has a Tukey depth of 2.
The depth region D_k is the set of all points \bm{z} with Tukey depth at least k. These form a series of nested convex polygons where D_{k+1} \subseteq D_k. The depth regions for our example are shown in Figure 5.10. The blue region is depth 1, the orange region is depth 2, the pink region is depth 3, and the green region is depth 4.
bagplot_fit <- aplpack::compute.bagplot(n01[1:10, 1:2])
z <- tidyr::expand_grid(
v1 = seq(-0.9, 1.6, l = 200),
v2 = seq(-1.95, 1.15, l = 200)
) |>
bind_rows(n01[1:10, 1:2])
z$depth <- round(
as.numeric(DepthProc::depthTukey(as.matrix(z), as.matrix(n01[1:10, 1:2]))) *
10
)
regions <- list()
depths <- sort(unique(z$depth))
depths <- depths[depths > 0]
for (i in depths) {
tmp <- z |> filter(depth == depths[i])
hull <- chull(tmp)
regions[[i]] <- tmp[c(hull, hull[1]), ]
}
cols <- c("#0072B2", "#D55E00", "#CC79A7", "#009E73")
p <- n01 |>
head(10) |>
ggplot(aes(x = v1, y = v2))
for (i in depths) {
p <- p +
geom_polygon(
aes(x = v1, y = v2),
data = regions[[i]],
fill = cols[i],
alpha = 0.8
)
}
p +
geom_point() +
coord_fixed(xlim = c(-0.85, 1.6), ylim = c(-1.85, 1.1)) +
geom_point(
data = data.frame(v1 = bagplot_fit$center[1], v2 = bagplot_fit$center[2]),
col = "yellow",
size = 2
)
A simple way to define a multivariate median is the “centre of gravity” of all points of maximum depth (Rousseeuw & Ruts, 1998). The centre of gravity (or “centroid”) of a shape is the average of all points within the shape. If the shape is convex, it is the point where the shape could be perfectly balanced on the tip of a pin if it were made of a uniform material. For example, the centre of gravity of a rectangle is the point in the middle of the rectangle, and the centre of gravity of a circle is the centre of the circle.
In Figure 5.10, the points of maximum depth are shown in the central green region. The centre of gravity of this region is the yellow point in the middle.
The bagplot was proposed by Rousseeuw et al. (1999) as a bivariate version of a boxplot, constructed using similar principles. Like a univariate boxplot, the bivariate bagplot has a central point (the depth median), an inner region (the bag), and an outer region (the loop), beyond which outliers are shown as individual points.
To define the bag, we first find the smallest depth region D_k containing at least \lfloor n/2 \rfloor of the observations. Then the bag is linearly interpolated between D_k and D_{k+1}, with the linear interpolation depending on the number of observations in each depth region. Because the bag is interpolated between depth regions, it is also a convex polygon. The procedure is slightly more complicated for small data sets where only D_1 might contain more than half of the observations.
To find the loop, we inflate the bag relative to the median by a factor of 3. This forms the “fence”. Then the loop is the convex hull of the points contained within the fence.
Figure 5.11 shows the bagplot from the same 10 observations as were used in the illustration of depth in the previous section, along with the observations themselves shown in black.
n01 |>
head(10) |>
gg_bagplot(v1, v2) +
geom_point(aes(x = v1, y = v2), data = n01[1:10, ])
A more interesting example is obtained in Figure 5.12, showing a bagplot of the durations and waiting times of Old Faithful eruptions.
oldfaithful |>
gg_bagplot(duration, waiting) +
labs(x = "Duration (seconds)", y = "Waiting time (seconds)")
This demonstrates the drawback of using a bagplot to identify anomalies. While the plot has identified all the genuine anomalies, it has also highlighted points on the edge of the data cloud which are probably not anomalies.
A useful variation of the bagplot, especially for larger data sets, displays a scatterplot of the observations, but coloured using the same colours as the bagplot, with the deepest observation in the strongest blue, the points within the bag in a lighter blue, and the points within the loop in the lightest colour. Outliers are shown in black if any exist.
oldfaithful |>
gg_bagplot(duration, waiting, show_points = TRUE) +
labs(x = "Duration (seconds)", y = "Waiting time (seconds)")
Figure 5.13 shows the value of this version of a bagplot, as you can see a lot of detail in the bag and loop that would not be visible otherwise.
The idea of depth regions and the depth median can be easily generalised to data with more than two dimensions (Rousseeuw & Struyf, 1998). Consequently, it would be possible to define a higher-dimensional version of the bagplot, although it would be difficult to plot it in four or more dimensions. In any case, approaches for high-dimensional data are discussed in Chapter 9.
We can compute HDRs (Section 2.7) from a kernel density estimate to find regions of the sample space where observations are unlikely to occur.
Let’s illustrate the idea with univariate data. In Section 2.9, we estimated the density of the duration of Old Faithful eruptions using a kernel density estimate, and plotted it in Figure 2.16. We can find the 50% and 99% HDRs of this density, which can then be used to form an HDR boxplot (Hyndman, 1996) as shown in Figure 5.14.
dist_kde(oldfaithful$duration) |>
gg_density(
hdr = "fill",
prob = c(0.5, 0.99),
show_points = TRUE,
jitter = TRUE,
show_mode = TRUE
) +
scale_y_continuous(breaks = NULL) +
labs(y = "") +
xlim(-10, 320)
The HDR boxplot shown here can also be produced using the gg_hdrboxplot() function, as shown in Figure 5.15.
oldfaithful |>
gg_hdrboxplot(duration, show_anomalies = FALSE) +
labs(x = "Duration (seconds)") +
xlim(-10, 320)
Points outside the 99% region are shown separately and jittered vertically to reduce overplotting. This type of boxplot has the advantage that it allows for multimodal distributions, and can identify “inliers” that occur in regions of low density between regions of high density.
By definition, about 1% of points will lie outside the 99% region, so the points shown are not necessarily anomalies, just observations that occur in the lower probability regions of the space. By setting show_anomalies = TRUE, some of the points are highlighted as potential anomalies, based on how unlikely each observation is under a density estimate computed without it (Section 7.5).
The idea naturally extends to higher dimensions. Figure 5.16 shows the bivariate HDR boxplot of the duration and waiting times for the Old Faithful data.
oldfaithful |>
gg_hdrboxplot(duration, waiting, show_anomalies = FALSE) +
labs(x = "Duration (seconds)", y = "Waiting time (seconds)")
There are 2097 observations used to compute this density estimate, and even with that many observations, computing the 99% HDR contour is difficult, as shown by the rather jagged boundary of the HDR region. This is simply a consequence of the sparsity of points in multidimensional space and in the tails of a distribution.
As with bagplots, there is a variation that shows the individual points coloured according to the HDR region in which they fall, as shown in Figure 5.17.
oldfaithful |>
gg_hdrboxplot(duration, waiting, show_anomalies = FALSE, show_points = TRUE) +
labs(x = "Duration (seconds)", y = "Waiting time (seconds)")
The plots in this chapter summarise a distribution by a series of nested central regions, and show the observations outside the outermost region individually. They differ in how “central” is defined. Boxplots and letter-value plots use univariate depth: a boxplot shows the median, a box bounded by the fourths (or quartiles), and whiskers extending to the most extreme observations within fences 1.5 IQR beyond the box, while a letter-value plot adds further letter values (eighths, sixteenths, and so on) until the next one can no longer be estimated reliably, which makes it better suited to large data sets. Bagplots extend the boxplot to bivariate data using Tukey depth, with the depth median at the centre, a bag containing the deepest half of the observations, and a loop enclosing the observations that lie within the bag inflated by a factor of 3. HDR boxplots define the regions by probability density instead, using highest density regions of a kernel density estimate.
The points shown individually are not necessarily anomalies. With fixed fences, a boxplot shows more spurious anomalies as the sample size grows; the modified IQR rule widens the fences with the sample size to control this, but only under Normality. Letter-value plots show a few legitimate observations in each tail by design, and an HDR boxplot shows about 1% of observations outside its 99% region.
Boxplots, letter-value plots, bagplots, and HDR boxplots are all extremely useful tools in exploratory data analysis, and deserve to be widely applied. They are particularly useful in summarising univariate and bivariate data distributions. However, as we have seen, none of them are effective at anomaly detection, and only HDR boxplots allow for data with multiple modes.
Nevertheless, we will find these tools useful when summarising results from anomaly detection algorithms, and in understanding why some points have been labelled as anomalies.