# What is the best way to identify outliers in multivariate data?

Suppose I have a large set of multivariate data with at least three variables. How can I find the outliers? Pairwise scatterplots won't work as it is possible for an outlier to exist in 3 dimensions that is not an outlier in any of the 2 dimensional subspaces.

I am not thinking of a regression problem, but of true multivariate data. So answers involving robust regression or computing leverage are not helpful.

One possibility would be to compute the principal component scores and look for an outlier in the bivariate scatterplot of the first two scores. Would that be guaranteed to work? Are there better approaches?

• If a scatterplot matrix won't catch it, you could try a 3D scatterplot. That won't work out to 4D, of course, but then you could create a 4th dimension as time and make a movie. :) Jul 20, 2010 at 13:06
• What about (hierarchical) cluster analysis? Jul 29, 2010 at 18:48
• 5th D in color hue, 6th D in color intensity, 7th D in point size and we can go ;-)
– user88
Jul 31, 2010 at 12:36
• There is a nice collection of anomality detection resources at: https://github.com/yzhao062/anomaly-detection-resources. For multivariate data, see Section 3.1.
– np8
Dec 6, 2020 at 18:26

Have a look at the mvoutlier package which relies on ordered robust mahalanobis distances, as suggested by @drknexus.

• I personally don't believe in this package! I think this is just noise of science. This types of theories are very good for simulated data while just pffff for real data. I did use this package to for example detect outliers :-) sadly half of my data detected to be outliers :-D it is funny how people just publish some (I don't know what to call this ) . if you have any comment, I will be happy to see . We have many types of outliers but should we simply remove them from our data? I believe NOT. Feb 27, 2015 at 19:03
• Could you explain when and why this works? Jan 25, 2019 at 21:25
• This is an important topic and I would also like more information about this package, or for it to be removed as a reference, please. If you don't know what "robust Mahalanobis distances" are, it's hard to choose this over the classic calculations.
– P.P.
Mar 31, 2021 at 11:04
• @Learner welcome to the world of R packages :)
– bers
May 20, 2021 at 12:26

I think Robin Girard's answer would work pretty well for 3 and possibly 4 dimensions, but the curse of dimensionality would prevent it working beyond that. However, his suggestion led me to a related approach which is to apply the cross-validated kernel density estimate to the first three principal component scores. Then a very high-dimensional data set can still be handled ok.

In summary, for i=1 to n

1. Compute a density estimate of the first three principal component scores obtained from the data set without Xi.
2. Calculate the likelihood of Xi for the density estimated in step 1. call it Li.

end for

Sort the Li (for i=1,..,n) and the outliers are those with likelihood below some threshold. I'm not sure what would be a good threshold -- I'll leave that for whoever writes the paper on this! One possibility is to do a boxplot of the log(Li) values and see what outliers are detected at the negative end.

• Thanks for this answer (should be a comment ? ) As I already mentionned, as a comment to Rich answer's High dimension are not a problem (even 1000 could work) if you make parametric structural assumption. Jul 23, 2010 at 6:43
• @rob "I'm not sure what would be a good threshold " this would be the purpose of the multiple testing procedure I mentionned .... but I fully agree that things have to be filled in and I really like the outlier detection in the outlier detection ! who wants to writte a paper :) ? Jul 23, 2010 at 6:47
• (-1) It is a bit upsetting that the accepted answer to this question is neither useful (e.g. too vague and not implemented) nor scientifically correct (as in not even wrong: as far as I know, it is not even mentioned in the relevant literature). Dec 27, 2012 at 21:20
• (-1) this seems as an incorrect answer - this method will not detect outliers! it is easy to imagine a 3D pancake-like cloud of points projected by PCA on a 2D-plane and 1 outlier point that is far away from the plane but its projection is right in the middle of projected cloud, so "kernel density" of the point in the 2D-plane will be high although it is an outlier. Dec 13, 2013 at 5:26
• I would think PCA would be problematic here. PCA is an approximation which considers the main contributions to the variance. But here the very aim is to find outliers - wouldn't they "fall between the cracks"? Jul 23, 2018 at 5:28

You can find a pedagogical summary of the various methods available in (1)

For some --recent-- numerical comparisons of the various methods listed there, you can check (2) and (3).

there are many older (and less exhaustive) numerical comparisons, typically found in books. You will find one on pages 142-143 of (4), for example.

Note that all the methods discussed here have an open source R implementation, mainly through the rrcov package.

• (1) P. Rousseeuw and M. Hubert (2013) High-Breakdown Estimators of Multivariate Location and Scatter.
• (2) M. Hubert, P. Rousseeuw, K. Vakili (2013). Shape bias of robust covariance estimators: an empirical study. Statistical Papers.
• (3) K. Vakili and E. Schmitt (2014). Finding multivariate outliers with FastPCS. Computational Statistics & Data Analysis.
• (4) Maronna R. A., Martin R. D. and Yohai V. J. (2006). Robust Statistics: Theory and Methods. Wiley, New York.
• I would not recommend papers which are based on simulated data, I believe the first two papers are based on simulated data. I am working on real data and RARELY I those simulated strong outliers can be found in a data. I personally never found a data set with very strong outliers. If you know any real data, please don't hesitate to share Feb 27, 2015 at 19:15
• Hi Nemo: thanks for your comment. Actually, it turns out that for many (saddly not all) of those models the worst possible configuration of outliers are actually known and these are the configurations one would typically simulate from and evaluate those models on. Perhaps you can formulate this as a question? I would be more than happy to point to some references! Feb 27, 2015 at 21:22
• I asked a general question already. You can find my post here stats.stackexchange.com/questions/139645/… Feb 27, 2015 at 21:52

I would do some sort of "leave one out testing algorithm" (n is the number of data):

for i=1 to n

1. compute a density estimation of the data set obtained by throwing $X_i$ away. (This density estimate should be done with some assumption if the dimension is high, for example, a gaussian assumption for which the density estimate is easy: mean and covariance)
2. Calculate the likelihood of $X_i$ for the density estimated in step 1. call it $L_i$.

end for

sort the $L_i$ (for i=1,..,n) and use a multiple hypothesis testing procedure to say which are not good ...

This will work if n is sufficiently large... you can also use "leave k out strategy" which can be more relevent when you have "groups" of outliers ...

• How would this deal with cases where you don't know how many outliers you have, i.e. when the N-1 points still have a bias since they include outliers? Jul 21, 2010 at 11:58
• if n is sufficiently large and the number of outlier is small then this bias is negligible. If there are a large number of outliers then, maibe it is not outliers and anyway, as I mentionned you can use leave k out strategy ... (in this case, you have to find out a strategy to avoid tracking all configurations which may be NP hard ... ) and if you don't know k, you can try many values for k and keep the most relevent. Jul 21, 2010 at 12:07
• This is a nice idea, but if I have understood it correctly, it seems to enlarge the idea of "outlier" to include any value in a dataset that is distant from the others. For example, in the batch {-110[1]-90, 0, 90[1]110} of 43 integers, wouldn't your procedure identify the 0 (which is the median of these numbers!) as the unique "outlier"?
– whuber
Aug 23, 2010 at 22:36
• @whuber good question! without structural assumption the zero might be removed. However you can assume the distribution is unimodal and compute the density accordingly. In this case the structural assumption (unimodality) will help the procedure to see zero as a "normal" value. At the end you think it is normal because you have that type of "structural assumption" in your head :) ? Jun 14, 2011 at 19:05
• I am uneasy about some aspects of this procedure; my example just points out one of the reasons why. The $L_i$ are strongly interdependent, so it's unclear which "multiple hypothesis testing procedure" would be applicable. Moreover, the $L_i$ will depend on the density estimator--as you point out--so that opens up a whole battery of questions about how best to estimate the densities for this purpose. These are problems even with parametric univariate outlier detection procedures, so I don't expect them to be any easier to solve in the multivariate case.
– whuber
Jun 14, 2011 at 20:15

I didn't see anybody mention influence functions. I first saw this idea in Gnanadesikan's multivariate book.

In one dimension an outlier is either an extremely large or an extremely small value. In multivariate analysis it is an observation removed from the bulk of the data. But what metric should we use to define extreme for the outlier? There are many choices. The Mahalanobis distance is just one. I think that looking for every type of outlier is futile and counterproductive. I would ask why do you care about the outlier? In estimating a mean they can have a great deal of influence on that estimate. Robust estimators downweight and accommodate outliers but they do not formally test for them. Now in regression, the outliers--like leverage points--could have large effects on the slope parameters in the model. With bivariate data they can unduly influence the estimated correlation coefficient and in three or more dimensions the multiple correlation coefficient.

Influence functions were introduced by Hampel as a tool in robust estimation and Mallows wrote a nice unpublished paper advocating their use. The influence function is a function of the point you are at in n-dimensional space and the parameter. It essentially measures the difference between the parameter estimate with the point in the calculation and with the point left out. Rather than go to the trouble of doing the calculation of the two estimates and taking the difference, often you can derive a formula for it. Then the contours of constant influence tell you the direction that is extreme with respect to the estimate of this parameter and hence tell you where in the n-dimensional space to look for the outlier.

For more you can look at my 1983 paper in the American Journal of Mathematical and Management Sciences titled "The influence function and its application to data validation." In data validation we wanted to look for outliers that affected the intended use of the data. My feeling is that you should direct your attention to outliers that greatly affect the parameters you are interested in estimating and not care so much about others that don't.

• +1 Very nice idea, especially for the consistent focus on the purpose of the analysis.
– whuber
May 5, 2012 at 0:50
• A similar idea is Random Sample Consensus (RANSAC). Aug 20, 2016 at 15:41

You can find candidates for "outliers" among the support points of the minimum volume bounding ellipsoid. (Efficient algorithms to find these points in fairly high dimensions, both exactly and approximately, were invented in a spate of papers in the 1970's because this problem is intimately connected with a question in experimental design.)

I novel approach I saw was by IT Jolliffe Principal Components Analysis. You run a PCA on your data (Note: PCA can be quite a useful data exploration tool in its own right), but instead of looking at the first few Principal Components (PCs), you plot the last few PCs. These PCs are the linear relationships between your variables with the smallest variance possible. Thus they detect "exact" or close to exact multivariate relationships in your data.

A plot of the PC scores for the last PC will show outliers not easily detectable by looking individually at each variable. One example is for height and weight - some who has "above average" height and "below average" weight would be detected by the last PC of height and weight (assuming these are positively correlated), even if their height and weight were not "extreme" individually (e.g. someone who was 180cm and 60kg).

It may be an overshoot, but you may train an unsupervised Random Forest on the data and use the object proximity measure to detect outliers. More details here.

For moderate dimensions, like 3, then some sort of kernel cross-validation technique as suggested elsewhere seems reasonable and is the best I can come up with.

For higher dimensions, I'm not sure that the problem is solvable; it lands pretty squarely into 'curse-of-dimensionality' territory. The issue is that distance functions tend to converge to very large values very quickly as you increase dimensionality, including distances derived from distributions. If you're defining an outlier as "a point with a comparatively large distance function relative to the others", and all your distance functions are beginning to converge because you're in a high-dimensional space, well, you're in trouble.

Without some sort of distributional assumption that will let you turn it into a probabilistic classification problem, or at least some rotation that lets you separate your space into "noise dimensions" and "informative dimensions", I think that the geometry of high-dimensional spaces is going to prohibit any easy -- or at least robust -- identification of outliers.

• The procedure I already gave can be applied in large dimension, as I said, using a gaussian assumption. If the dimension is really large with respect to the sample size (i.e. p>>n) then you can make some sparcity assumption (assume that the parameters of your gaussian distribution lie in a low dimensional space for example) and use a thresholding estimation procedure for the estimation of the parameters... Jul 22, 2010 at 13:16
• Very insightful that curse of dimensionality precludes a solution without distribution assumptions (unless you simply have insane amounts of data or the dimension is small) Sep 11, 2012 at 19:15

I'm not sure what you mean when you say you aren't thinking of a regression problem but of "true multivariate data". My initial response would be to calculate the Mahalanobis distance since it doesn't require that you specify a particular IV or DV, but at its core (as far as I understand it) it is related to a leverage statistic.

• Yes, I could compute the Mahalanobis distance of each observation from the mean of the data. But the observations with the greatest distance from the mean are not necessarily multivariate outliers. Think of a bivariate scatterplot with high correlation. An outlier can be outside the main cloud of points but not that far from the mean. Jul 20, 2010 at 8:14

I'm not aware that anyone is doing this, but I generally like to try dimensionality reduction when I have a problem like this. You might look into a method from manifold learning or non-linear dimensionality reduction.

An example would be a Kohonen map. A good reference for R is "Self- and Super-organizing Maps in R: The kohonen Package".

My first response would be that if you can do multivariate regression on the data, then to use the residuals from that regression to spot outliers. (I know you said it's not a regression problem, so this might not help you, sorry !)

I'm copying some of this from a Stackoverflow question I've previously answered which has some example R code

First, we'll create some data, and then taint it with an outlier;

> testout<-data.frame(X1=rnorm(50,mean=50,sd=10),X2=rnorm(50,mean=5,sd=1.5),Y=rnorm(50,mean=200,sd=25))
> #Taint the Data
> testout$X1[10]<-5 > testout$X2[10]<-5
> testout$Y[10]<-530 > testout X1 X2 Y 1 44.20043 1.5259458 169.3296 2 40.46721 5.8437076 200.9038 3 48.20571 3.8243373 189.4652 4 60.09808 4.6609190 177.5159 5 50.23627 2.6193455 210.4360 6 43.50972 5.8212863 203.8361 7 44.95626 7.8368405 236.5821 8 66.14391 3.6828843 171.9624 9 45.53040 4.8311616 187.0553 10 5.00000 5.0000000 530.0000 11 64.71719 6.4007245 164.8052 12 54.43665 7.8695891 192.8824 13 45.78278 4.9921489 182.2957 14 49.59998 4.7716099 146.3090 <snip> 48 26.55487 5.8082497 189.7901 49 45.28317 5.0219647 208.1318 50 44.84145 3.6252663 251.5620  It's often most usefull to examine the data graphically (you're brain is much better at spotting outliers than maths is) > #Use Boxplot to Review the Data > boxplot(testout$X1, ylab="X1")
> boxplot(testout$X2, ylab="X2") > boxplot(testout$Y, ylab="Y")


You can then use stats to calculate critical cut off values, here using the Lund Test (See Lund, R. E. 1975, "Tables for An Approximate Test for Outliers in Linear Models", Technometrics, vol. 17, no. 4, pp. 473-476. and Prescott, P. 1975, "An Approximate Test for Outliers in Linear Models", Technometrics, vol. 17, no. 1, pp. 129-132.)

> #Alternative approach using Lund Test
> lundcrit<-function(a, n, q) {
+ # Calculates a Critical value for Outlier Test according to Lund
+ # See Lund, R. E. 1975, "Tables for An Approximate Test for Outliers in Linear Models", Technometrics, vol. 17, no. 4, pp. 473-476.
+ # and Prescott, P. 1975, "An Approximate Test for Outliers in Linear Models", Technometrics, vol. 17, no. 1, pp. 129-132.
+ # a = alpha
+ # n = Number of data elements
+ # q = Number of independent Variables (including intercept)
+ F<-qf(c(1-(a/n)),df1=1,df2=n-q-1,lower.tail=TRUE)
+ crit<-((n-q)*F/(n-q-1+F))^0.5
+ crit
+ }

> testoutlm<-lm(Y~X1+X2,data=testout)

> testout$fitted<-fitted(testoutlm) > testout$residual<-residuals(testoutlm)

> testout$standardresid<-rstandard(testoutlm) > n<-nrow(testout) > q<-length(testoutlm$coefficients)

> crit<-lundcrit(0.1,n,q)

> testout$Ynew<-ifelse(testout$standardresid>crit,NA,testout\$Y)

> testout
X1         X2        Y    newX1   fitted    residual standardresid
1  44.20043  1.5259458 169.3296 44.20043 209.8467 -40.5171222  -1.009507695
2  40.46721  5.8437076 200.9038 40.46721 231.9221 -31.0183107  -0.747624895
3  48.20571  3.8243373 189.4652 48.20571 203.4786 -14.0134646  -0.335955648
4  60.09808  4.6609190 177.5159 60.09808 169.6108   7.9050960   0.190908291
5  50.23627  2.6193455 210.4360 50.23627 194.3285  16.1075799   0.391537883
6  43.50972  5.8212863 203.8361 43.50972 222.6667 -18.8306252  -0.452070155
7  44.95626  7.8368405 236.5821 44.95626 223.3287  13.2534226   0.326339981
8  66.14391  3.6828843 171.9624 66.14391 148.8870  23.0754677   0.568829360
9  45.53040  4.8311616 187.0553 45.53040 214.0832 -27.0279262  -0.646090667
10  5.00000  5.0000000 530.0000       NA 337.0535 192.9465135   5.714275585
11 64.71719  6.4007245 164.8052 64.71719 159.9911   4.8141018   0.118618011
12 54.43665  7.8695891 192.8824 54.43665 194.7454  -1.8630426  -0.046004311
13 45.78278  4.9921489 182.2957 45.78278 213.7223 -31.4266180  -0.751115595
14 49.59998  4.7716099 146.3090 49.59998 201.6296 -55.3205552  -1.321042392
15 45.07720  4.2355525 192.9041 45.07720 213.9655 -21.0613819  -0.504406009
16 62.27717  7.1518606 186.6482 62.27717 169.2455  17.4027250   0.430262983
17 48.50446  3.0712422 228.3253 48.50446 200.6938  27.6314695   0.667366651
18 65.49983  5.4609713 184.8983 65.49983 155.2768  29.6214506   0.726319931
19 44.38387  4.9305222 213.9378 44.38387 217.7981  -3.8603382  -0.092354925
20 43.52883  8.3777627 203.5657 43.52883 228.9961 -25.4303732  -0.634725264
<snip>
49 45.28317  5.0219647 208.1318 45.28317 215.3075  -7.1756966  -0.171560291
50 44.84145  3.6252663 251.5620 44.84145 213.1535  38.4084869   0.923804784
Ynew
1  169.3296
2  200.9038
3  189.4652
4  177.5159
5  210.4360
6  203.8361
7  236.5821
8  171.9624
9  187.0553
10       NA
11 164.8052
12 192.8824
13 182.2957
14 146.3090
15 192.9041
16 186.6482
17 228.3253
18 184.8983
19 213.9378
20 203.5657
<snip>
49 208.1318
50 251.5620


Obviosuly there are other outlier tests than the Lund test (Grubbs springs to mind), but I'm not sure which are better suited to multivariate data.

• Univariate boxplots are useful for spotting univariate outliers. But they can completely miss multivariate outliers. The regression idea is ok if I had a Y and a bunch of X variables. But as I said in the question, there is no Y so regression is inappropriate. Jul 20, 2010 at 8:12
• @RobHyndman - one may fix Y (make it a constant) and try to model a multi-regression without an intercept. It works beautifully in practice. In essence this transforms a fundamentally unsupervised problem, which a typical anomaly-detection is, into a supervised one. May 28, 2016 at 2:48
• @arielf Interesting. Do you have a reference for that? May 28, 2016 at 3:30
• @RobHyndman not an academic paper, but perhaps better, you can try it yourself. See my github page on vw-top-errors @ goo.gl/l7SLlB (note that the examples there happen to have a Y, but I used the same technique, very successfully on unsupervised problems too by fixing Y.) May 28, 2016 at 3:41

One of the above answers touched in mahalanobis distances.... perhaps anpther step further and calculating simultaneous confidence intervals would help detect outliers!