Data

X <- readRDS(system.file("testdata", "three-pops.rds", package = "bigutilsr"))
pca <- prcomp(X, scale. = TRUE, rank. = 10)
U <- pca$x

library(ggplot2)
theme_set(bigstatsr::theme_bigstatsr(0.8))
qplot(U[, 1], U[, 2]) + coord_equal()

Measuring outlierness

The standard way to detect outliers in genetics is the criterion of being “more than 6 standard deviations away from the mean”.

apply(U, 2, function(x) which( abs(x - mean(x)) > (6 * sd(x)) ))

## integer(0)

Here, there is no outlier according to this criterion. Let us make some fake one.

U2 <- U
U2[1, 1] <- 30
qplot(U2[, 1], U2[, 2]) + coord_equal()

apply(U2, 2, function(x) which( abs(x - mean(x)) > (6 * sd(x)) ))

## integer(0)

Still not an outlier..

U3 <- U2
U3[1, 1] <- 80
qplot(U3[, 1], U3[, 2]) + coord_equal()

library(magrittr)
apply(U3, 2, function(x) which( abs(x - mean(x)) > (6 * sd(x)) )) %>%
  Reduce(union, .)

## [1] 1

Now, the first sample is considered as an outlier by this criterion.

A more robust variation

Note that you might want to use median() instead of mean() and mad() instead of sd() because they are more robust estimators. This becomes

ind.out <- apply(U3, 2, function(x) which( (abs(x - median(x)) / mad(x)) > 6 )) %>%
  Reduce(union, .) %>%
  print()

## [1]   1 516

We get a new outlier.

col <- rep("black", nrow(U3)); col[ind.out] <- "red"
qplot(U3[, 1], U3[, 3], color = I(col), size = I(2)) + coord_equal()

A continuous view of this criterion

This criterion flag an outlier if it is an outlier for at least one principal component (PC). This corresponds to using the max() (infinite) distance (in terms of number of standard deviations) from the mean.

dist <- apply(U3, 2, function(x) abs(x - median(x)) / mad(x)) %>%
  apply(1, max)
qplot(U3[, 1], U3[, 3], color = dist, size = I(3)) + coord_equal() + 
  scale_color_viridis_c(trans = "log", breaks = c(1, 3, 6))

qplot(y = sort(dist, decreasing = TRUE)) +
  geom_hline(yintercept = 6, color = "red")

Robust Mahalanobis distance

Instead of using the infinite distance, Mahalanobis distance is a multivariate distance based on all variables (PCs here) at once. We use a robust version of this distance, which is implemented in packages {robust} and {robustbase} (Gnanadesikan and Kettenring 1972, Yohai and Zamar (1988), Maronna and Zamar (2002), Todorov, Filzmoser, and others (2009)) and that is reexported in {bigutilsr}.

dist2 <- covRob(U3, estim = "pairwiseGK")$dist
qplot(dist, sqrt(dist2))

cowplot::plot_grid(
  qplot(U3[, 1], U3[, 2], color = dist2, size = I(2)) + coord_equal() + 
    scale_color_viridis_c(trans = "log", breaks = NULL),
  qplot(U3[, 3], U3[, 7], color = dist2, size = I(2)) + coord_equal() + 
    scale_color_viridis_c(trans = "log", breaks = NULL),
  rel_widths = c(0.7, 0.4), scale = 0.95
)

This new criterion provides similar results for this data. These robust Mahalanobis distances are approximately Chi-square distributed, which enables deriving p-values of outlierness.

pval <- pchisq(dist2, df = 10, lower.tail = FALSE)
hist(pval)

is.out <- (pval < (0.05 / length(dist2)))  # Bonferroni correction
sum(is.out)

## [1] 33

qplot(U3[, 3], U3[, 7], color = is.out, size = I(3)) + coord_equal()

Local Outlier Factor (LOF)

LOF statistic (Breunig et al. 2000) has been cited more than 4000 times. Instead of computing a distance from the center, it uses some local density of points. We make use of the fast K nearest neighbours implementation of R package {nabor} (Elseberg et al. 2012) to implement this statistic efficiently in {bigutilsr}.

llof <- LOF(U3)  # log(LOF) by default
qplot(dist2, llof)

The fake outlier that we introduced is now clearly an outlier. The other points, not so much.

cowplot::plot_grid(
  qplot(U3[, 1], U3[, 2], color = llof, size = I(3)) + coord_equal() + 
    scale_color_viridis_c(breaks = NULL),
  qplot(U3[, 3], U3[, 7], color = llof, size = I(3)) + coord_equal() + 
    scale_color_viridis_c(breaks = NULL),
  rel_widths = c(0.7, 0.4), scale = 0.95
)

Tukey’s rule

Tukey’s rule (Tukey 1977) is a standard rule for detecting outliers. Here, we will apply it on the previously computed statistics. Note that we could use it directly on PCs, which is not much different from the robust version of the first criterion we introduced.

x <- rnorm(10000)
(tukey_up  <- quantile(x, 0.75) + 1.5 * IQR(x))

##     75% 
## 2.70692

(tukey_low <- quantile(x, 0.25) - 1.5 * IQR(x))

##       25% 
## -2.725665

hist(x); abline(v = c(tukey_low, tukey_up), col = "red")

mean(x < tukey_low | x > tukey_up)

## [1] 0.0057

where IQR(x) is equal to quantile(x, 0.75) - quantile(x, 0.25) (the InterQuartile Range).

However, there are two pitfalls when using Tukey’s rule:

1. Tukey’s rule assumes a normally distributed sample. When the data is skewed, it does not work that well.

   x <- rchisq(10000, df = 5)
   (tukey_up  <- quantile(x, 0.75) + 1.5 * IQR(x))

   ##      75% 
   ## 12.42084

   (tukey_low <- quantile(x, 0.25) - 1.5 * IQR(x))

   ##       25% 
   ## -3.232256

   hist(x, "FD"); abline(v = c(tukey_low, tukey_up), col = "red")

   

   mean(x < tukey_low | x > tukey_up)

   ## [1] 0.0294

   To solve the problem of skewness, the medcouple (mc) has been introduced (Hubert and Vandervieren 2008) and is implemented in robustbase::adjboxStats().

2. Tukey’s rule uses a fixed coefficient (1.5) that does not account for multiple testing, which means that for large samples, you will almost always get some outliers if using 1.5.

To solve these two issues, we implemented tukey_mc_up() that accounts both for skewness and multiple testing by default.

x <- rchisq(10000, df = 5)
(tukey_up  <- quantile(x, 0.75) + 1.5 * IQR(x))

##      75% 
## 12.48751

hist(x, "FD"); abline(v = tukey_up, col = "red")
abline(v = print(tukey_mc_up(x, coef = 1.5)), col = "blue")

## [1] 16.74215

abline(v = print(tukey_mc_up(x)), col = "green")  # accounts for multiple testing

## [1] 25.93299

Applying this corrected Tukey’s rule to our statistics:

tukey_mc_up(dist)

## [1] 6.406337

qplot(dist2, llof) +
  geom_vline(xintercept = tukey_mc_up(dist2), color = "red") +
  geom_hline(yintercept = tukey_mc_up(llof),  color = "red")

Histogram’s gap

This rule I come up with assumes that the “normal” data is somewhat grouped and the outliers have some gap (in the histogram, there is a bin with no value in it) with the rest of the data.

For example, for dist, there is a gap just before 6, and we can derive an algorithm to detect this:

hist(dist, breaks = nclass.scottRob)
str(hist_out(dist))

## List of 2
##  $ x  : num [1:515] 2.08 2.06 1.74 1.86 2.04 ...
##  $ lim: num [1:2] -Inf 5.75

abline(v = hist_out(dist)$lim[2], col = "red")

hist(dist2, breaks = nclass.scottRob)
abline(v = hist_out(dist2)$lim[2], col = "red")

hist(llof, breaks = nclass.scottRob)
abline(v = hist_out(llof)$lim[2], col = "red")

This criterion is convenient because it does not assume any distribution of the data, just that it is compact and that the outliers are not in the pack.

It could be used in other contexts, e.g. choosing the number of outlier principal components:

eigval <- pca$sdev^2
hist(eigval, breaks = "FD")  # "FD" gives a bit more bins than scottRob
abline(v = hist_out(eigval, breaks = "FD")$lim[2], col = "red")

sum(eigval > hist_out(eigval, breaks = "FD")$lim[2])

## [1] 3

pca_nspike(eigval)  # directly implemented in {bigutilsr}

## [1] 3

Note the possible use of bootstrap to make hist_out() and pca_nspike() more robust.