Mahalanobis Distance and its Limitations

library(dplyr)
library(ggplot2)
library(mvtnorm)
library(robust)

Here are presented some situations using artificially generated data where the Mahalanobis distance does a good and a bad job detecting outliers.

Spherical Data

We know that for spherical data (no correlation between variables) the euclidean distance works ok.

sample = rmvnorm(1000, sigma=diag(2))
center = colMeans(sample)
distances = euclidean_norm(sample-center)
cutoff = quantile(distances, 0.90)
custom_plot(sample, distances>cutoff, center)

Non-Spherical Data

For non-spherical data the euclidean norm does not work as well.

sigma = matrix(c(1, 0.8, 0.8, 1), ncol=2)
sample = rmvnorm(1000, sigma=sigma)
center = colMeans(sample)
distances = euclidean_norm(sample-center)
cutoff = quantile(distances, 0.90)
custom_plot(sample, distances>cutoff, center)

The Mahalanobis distance works better since it considers the correlation structure of the data.

sigma = matrix(c(1, 0.8, 0.8, 1), ncol=2)
sample = rmvnorm(1000, sigma=sigma)
center = colMeans(sample)
distances = mahalanobis(sample, center=center, cov=cov(sample))
cutoff = quantile(distances, 0.90)
custom_plot(sample, distances>cutoff, center)

Mahalanobis Distance Robustness

The Mahalanobis distance to the center of the data defined as:

\[ d(x) = (x-\mu)^t \Sigma^{-1} (x-\mu)\]

Where \(\mu\) is an estimation of the center of the cloud of points, generally the mean and \(\Sigma\) is an estimation of the covariance matrix. Therefore the Mahalanobis distance is as good as these estimations are.

We’ll now see some shortcomings of the Mahalanobis distance and how to deal with them.

Better location estimator

We already know that the mean is not very robust to outliers. We can improve our location estimator to a more robust alternative defined for a set of points \(x_1, x_2 \dots x_n\) as:

\[\arg \min_{1 \leq i \leq n} \sum_{j=1}^n \| x_i - x_j \|_2\]

Or in other words, the point from the set of points that minimizes the euclidian distance to all the other points.

geometric_median <- function (sample) {
  euclidian_distances = apply(sample, 1, function (x) {sum(euclidean_norm(sample - x))})
  sample[which.min(euclidian_distances),]
}

sigma1 = matrix(c(1, 0.8, 0.8, 1), ncol=2)
sigma2 = matrix(c(1, 0.3, 0.3, 1), ncol=2)
sample1 = rmvnorm(1000, sigma=sigma1)
sample2 = rmvnorm(100, mean=c(-3, 8), sigma=sigma2)
sample = rbind(sample1, sample2)
center = colMeans(sample)
custom_plot_no_outliers(sample, center)

sigma1 = matrix(c(1, 0.8, 0.8, 1), ncol=2)
sigma2 = matrix(c(1, 0.3, 0.3, 1), ncol=2)
sample1 = rmvnorm(1000, sigma=sigma1)
sample2 = rmvnorm(100, mean=c(-3, 8), sigma=sigma2)
sample = rbind(sample1, sample2)
center = geometric_median(sample)
custom_plot_no_outliers(sample, center)

We can see that this estimate of the center of the cloud of data is more robust to noise.

Better scale estimator

Minimum Covariance Determinant

A classical approach to computing a robust covariance matrix is the use of the Minimum Covariance Determinant (MCD). This method not only provides estimators for the covariance matrix but also location estimates.

The location estimate of MCD, \(\hat{\mu}\) is the mean of the \(h\) observations for which the determinant of the covariance matrix is minimal.

\(\hat\Sigma\) is the covariance matrix of the \(h\) observations multiplied by a consistency factor \(c_0\) (needed to achieve consistency on normal distributions).

The main parameter of this method is \(\alpha\), this value roughly represents the value \(h/n\), the proportion of observations used for the estimation. Since the value of \(h\) must lay between \((n + p + 1)/2\) and \(n\) (for \(n\) observations of \(p\) variables), \(\alpha\) must be between 0.5 and 1. A low \(\alpha\) translates in more robustness but low consistency at the normal distribution, whereas an \(\alpha\) of 1 is as bad as the classical sample covariance matrix.

A good in depth explanation of MCD can be found here

Problems with the Sample Covariance Matrix

So we’ve introduced some noise and we can see that the Mahalanobis distance does not properly identify the outliers. This is because the heavy presence of outliers inflates the covariance estimates allowing outliers to pass undetected.

sigma1 = matrix(c(1, 0.9, 0.9, 1), ncol=2)
sigma2 = matrix(c(1, 0.3, 0.3, 1), ncol=2)
sample1 = rmvnorm(1000, sigma=sigma1)
sample2 = rmvnorm(100, mean=c(-2, 3), sigma=sigma2)
sample = rbind(sample1, sample2)
center = colMeans(sample)
distances = mahalanobis(sample, center=center, cov=cov(sample))
cutoff = quantile(distances, 0.90)
custom_plot(sample, distances>cutoff, center)

sigma1 = matrix(c(1, 0.9, 0.9, 1), ncol=2)
sigma2 = matrix(c(1, 0.3, 0.3, 1), ncol=2)
sample1 = rmvnorm(1000, sigma=sigma1)
sample2 = rmvnorm(100, mean=c(-2, 3), sigma=sigma2)
sample = rbind(sample1, sample2)
cov.obj = covRob(sample, estim="mcd", alpha=0.7)
center = cov.obj$center
cov = cov.obj$cov
distances = mahalanobis(sample, center=center, cov=cov)
cutoff = quantile(distances, 0.90)
custom_plot(sample, distances>cutoff, center)

Other issues

This section exhibits some limitations of the Mahalanobis distance that I do not know how to overcome.

Non Linearity

The Mahalanobis distance is only able to capture linear relationships between the variables. This case presents a cuadratic relationship so the outlier cut here is not quite right.

x = rnorm(1000, mean=5, sd=2)
y = rnorm(1000, mean=x^2, sd=3)
sample = cbind(x, y)
cov.obj = covRob(sample, estim="mcd", alpha=0.7)
center = cov.obj$center
cov = cov.obj$cov
distances = mahalanobis(sample, center=center, cov=cov)
cutoff = quantile(distances, 0.90)
custom_plot(sample, distances>cutoff, center)

Heteroskedacity

We can find heteroskedacity or heterogeneity of the variance when the variance of a variable depends on the value of another. In this example we can see that when \(x\) increases not only \(y\) increases but also the variability of \(y\) is greater.

x = rgamma(1000, shape=1, scale=1)
y = rnorm(1000, mean=x, sd=0.2*x)
sample = cbind(x, y)
cov.obj = covRob(sample, estim="mcd", alpha=0.7)
center = cov.obj$center
cov = cov.obj$cov
distances = mahalanobis(sample, center=center, cov=cov)
cutoff = quantile(distances, 0.90)
custom_plot_heteroskedacity(sample, distances>cutoff, center)

A differs more significantly than B, because points with higher \(x\) should be expected to vary more than points near \(x=0\). But since the Mahalanobis distance is only able to define a tolerance ellipsoid, it does not do such a good job in this case.

Other properties of the data that I believe the Mahalanobis distance would not be able to address well but are not showcased here are: skewness and multimodality.