# Does Pearson correlation require removal of bivariate or univariate outliers?

Does the Pearson's correlation estimator require no bivariate outliers, or no outliers in each of two individual vectors of data? The answer will impact on how I winsorize outliers before calculating my Pearson estimate.

-
@whuber: it's interesting how the standards for what constitutes an answer has involved over the last two years. Many of the answers to that earlier question (including my own) would probably not pass the bar now (not enough examples, illustrations, citations, pointer to good implementations, ect...). – user603 Dec 27 '12 at 17:23
@user603 It's interesting to look back, isn't it? That was one of the first questions answered on this site. I think much of that evolution occurred within the first few weeks. It wouldn't hurt to keep improving or augmenting those early answers. – whuber Dec 27 '12 at 17:32

You will have to remove both. But that will not be enough. You will also have to remove those observations that are outlying on any projections of your data along any direction in $\mathbb{R}^2$. This is because of so-called multivariate outliers.

It is possible for outliers to depart significantly from the pattern of the majority of the data without necessarily standing out on any of your variables taken individually (As an example, consider the cluster of red dots in the plot attached to this answer).

Because of this, multi (in this case bi)-variate outliers can in general not be reliably detected using a coordinate-wise approaches. The only reliable way of finding them is in fact to use a multivariate trimming/winsorizing approach (i.e. one that considers univariate projections of your data along all directions). Many such method exists and you will find good implementations in most modern statistical packages (R, MATLAB, STATA, SAS,....).

It is important to recognize that multivariate outliers are equally effective at wrecking the Pearson correlation as their better known coordinate-wise cousins. An example below illustrates this.

Geometrically, coordinate-wise trimming/winsorizing amounts to drawing a rectangle around the majority of your data and considering any observations outside that rectangle as outlying. In contrast, multivariate trimming amounts to drawing an ellipse around the majority of your data and considering any point outside of that ellipse as an outlier.

The former approach can only detect outliers if they are outlying on at least one of the coordinates (or in other words along a direction parallel to an axis of a scatter-plot of your data). The second approach, in contrast, does not suffer from this limitation. In other words, multivariate trimming approaches can detect outliers regardless of the multivariate direction in which they are outlying (and this include directions parallel to the an axis of a scatter-plot of your data).

Consider this example (the code to reproduce it is below this post):

The Pearson correlation computed on the full data (that is for the black+red points considered together) is 0.5. In this case, none of the observations will be down-weighted by a coordinate-wise winsorizing approach since for all observations $1\leq i\leq n$, the distance

$$\max\left(\frac{|x_{i1}-\text{median}(x_1)|}{\text{mad}(x_1)},\frac{|x_{i2}-\text{median}(x_2)|}{\text{mad}(x_2)}\right).$$

is smaller than 3. Therefore, in this case, the Pearson correlation computed on the winsorized observations would be identical to the Pearson correlation computed on the original data.

Now, 0.5 is very far from the correlation of the good part of the data (e.g. considering only the black dots) which in this case is 84%.

Contrast this with the results obtained by carrying a bivariate trimming and estimating the correlation of the remaining (untrimmed) observations. In this case, the multivariate trimming was done using the FastMCD(1) algorithm, perhaps the most popular multivariate trimming algorithm. FastMCD correctly identifies the red dots as being too far from the ellipse enclosing the majority of the data and flags them as outliers. Then, the correlation estimated on the remaining observations is now 85%, which is close enough to the correct result.

(1) Rousseeuw P. J. and Van Driessen K. (1999). A Fast Algorithm for the Minimum Covariance Determinant Estimator. Technometrics, 41, 212--223.

library(MASS)
library(rrcov)
n<-100
p<-2
set.seed(123)
A<-matrix(rnorm((p+1)*p),p+1,p)
A<-eigen(var(A))$vector B<-A%*%diag(c(16,1))%*%t(A) C<-t(A)%*%diag(c(4,1))%*%A x<-mvrnorm(n,rep(0,p),B) y<-mvrnorm(n,rep(0,p),C) d<-which.max(mahalanobis(y,rep(0,p),B)) y<-mvrnorm(floor(n/5),y[d,],diag(2)/100) z<-rbind(x,y) plot(z,asp=1,type="n") points(x,col="black",pch=16) points(y,col="red",pch=16) cor(z) [,1] [,2] [1,] 1.0000000 0.5018708 [2,] 0.5018708 1.0000000 cov2cor(CovMcd(z)@cov) [,1] [,2] [1,] 1.0000000 0.8592597 [2,] 0.8592597 1.0000000 cor(x) [,1] [,2] [1,] 1.0000000 0.8485822 [2,] 0.8485822 1.0000000 d1<-which((abs(z[,1]-median(z[,1]))/mad(z[,1])<3) & (abs(z[,2]-median(z[,2]))/mad(z[,2])<3)) length(d1) [1] 120  - (+1) What is the basis for your winsorization rule$(x_j−median(x_j))/mad(x_j)<3$,$j={1,2}\$? Also you told me we have to use both univariate and bivariate winsorization but that's not what you decided to do? – Jase Dec 27 '12 at 11:12
@Jase: the conclusions would be the same for any univariate winsorization rule. If you describe which one you had in mind, I’d be happy to adapt my answer. I edited my answer to address your second question. Let me know if this is not clear. – user603 Dec 27 '12 at 16:15