The correlation coefficient is:

$$ r = \frac{\sum_k \frac{(x_k - \bar{x}) (y_k - \bar{y_k})}{s_x s_y}}{n-1} $$

The sample mean and the sample standard deviation are sensitive to outliers.

As well, the mechanism where,

$$ r = \frac{\sum_k \text{stuff}_k}{n -1} $$

is sort of like a mean as well and maybe there might be a variation on that which is less sensitive to variation.

The sample mean is:

$$ \bar{x} = \frac{\sum_k x_k}{n} $$

The sample standard deviation is:

$$ s_x = \sqrt{\frac{\sum_k (x_k - \bar{x})^2}{n -1}} $$

I think I want

The median:

$$ \text{Median}[x]$$

The median absolute deviation:

$$ \text{Median}[\lvert x - \text{Median}[x]\rvert] $$

And for the correlation:

$$ \text{Median}\left[\frac{(x -\text{Median}[x])(y-\text{Median}[y]) }{\text{Median}[\lvert x - \text{Median}[x]\rvert]\text{Median}[\lvert y - \text{Median}[y]\rvert]}\right] $$

I tried this with some random numbers but got results greater than 1 which seems wrong. See the following R code.

 x<- c(237, 241, 251, 254, 263)
 y<- c(216, 218, 227, 234, 235)

 median.x <- median(x)
 median.y <- median(y)

 mad.x <- median(abs(x - median.x))
 mad.y <- median(abs(y - median.y))

 r <- median((((x - median.x) * (y - median.y)) / (mad.x * mad.y)))

 ## Prints 1.125

  • 1
    $\begingroup$ I'm not sure what your actual question is, unless you mean your title? If so, the Spearman correlation is a correlation that is less sensitive to outliers. It's basically a Pearson correlation of the ranks. $\endgroup$
    – Ashe
    Nov 14, 2016 at 21:25
  • 8
    $\begingroup$ Are you asking for a robust estimator of the usual correlation or for an alternative measure of co-variation that happens to be robust? $\endgroup$
    – whuber
    Nov 14, 2016 at 22:20
  • $\begingroup$ A related question with answers: stats.stackexchange.com/questions/381194/… $\endgroup$ Mar 6, 2019 at 13:41

4 Answers 4


I think you want a rank correlation. Those are generally more robust to outliers, although it's worth recognizing that they are measuring the monotonic association, not the straight line association. The most commonly known rank correlation is Spearman's correlation. It is just Pearson's product moment correlation of the ranks of the data.

I wouldn't go down the path you're taking with getting the differences of each datum from the median. The median of the distribution of X can be an entirely different point from the median of the distribution of Y, for example. That strikes me as likely to cause instability in the calculation.


Another answer for discrete as opposed to continuous variables, e.g., integers versus reals, is the Kendall rank correlation. In contrast to the Spearman rank correlation, the Kendall correlation is not affected by how far from each other ranks are but only by whether the ranks between observations are equal or not.

The Kendall τ coefficient is defined as:

$\tau = \frac{(\text{number of concordant pairs}) - (\text{number of discordant pairs})}{n (n-1) /2}$

The Kendall rank coefficient is often used as a test statistic in a statistical hypothesis test to establish whether two variables may be regarded as statistically dependent. This test is non-parametric, as it does not rely on any assumptions on the distributions of $X$ or $Y$ or the distribution of $(X,Y)$.

The treatment of ties for the Kendall correlation is, however, problematic as indicated by the existence of no less than 3 methods of dealing with ties. A tie for a pair {(xiyi), (xjyj)} is when xi = xj or yi = yj; a tied pair is neither concordant nor discordant.


This is a solution which works well for the data and problem proposed by IrishStat.


The idea is to replace the sample variance of $Y$ by the predicted variance $$\sigma_Y^2=a^2\sigma_x^2+\sigma_e^2$$. so that the formula for the correlation becomes $$ r=\sqrt{\frac{a^2\sigma^2_x}{a^2\sigma_x^2+\sigma_e^2}}$$ Now the reason that the correlation is underestimated is that the outlier causes the estimate for $\sigma_e^2$ to be inflated. To deal with this replace the assumption of normally distributed errors in the regression with a normal mixture $$\frac{0.95}{\sqrt{2\pi} \sigma} \exp(-\frac{e^2}{2\sigma^2}) +\frac{0.05}{\sqrt{2\pi} 3\sigma} \exp(-\frac{e^2}{18\sigma^2}) $$ I first saw this distribution used for robustness in Hubers book, Robust Statistics. This is "moderately" robust and works well for this example. It also has the property that if there are no outliers it produces parameter estimates almost identical to the usual least squares ones. So this procedure implicitly removes the influence of the outlier without having to modify the data. Fitting the data produces a correlation estimate of 0.944812.

  • $\begingroup$ Your .94 is uncannily close to the .94 I computed when I reversed y and x . Is this by chance ? $\endgroup$
    – IrishStat
    Nov 19, 2016 at 19:25
  • $\begingroup$ I think it is just by chance. $\endgroup$ Nov 19, 2016 at 20:33

My answer premises that the OP does not already know what observations are outliers because if the OP did then data adjustments would be obvious. Thus part of my answer deals with identification of the outlier(s)

When you construct an OLS model ($y$ versus $x$), you get a regression coefficient and subsequently the correlation coefficient I think it may be inherently dangerous not to challenge the "givens" . In this way you understand that the regression coefficient and its sibling are premised on no outliers/unusual values. Now if you identify an outlier and add an appropriate 0/1 predictor to your regression model the resultant regression coefficient for the $x$ is now robustified to the outlier/anomaly. This regression coefficient for the $x$ is then "truer" than the original regression coefficient as it is uncontaminated by the identified outlier. Note that no observations get permanently "thrown away"; it's just that an adjustment for the $y$ value is implicit for the point of the anomaly. This new coefficient for the $x$ can then be converted to a robust $r$.

An alternative view of this is just to take the adjusted $y$ value and replace the original $y$ value with this "smoothed value" and then run a simple correlation.

This process would have to be done repetitively until no outlier is found.

I hope this clarification helps the down-voters to understand the suggested procedure . Thanks to whuber for pushing me for clarification. If anyone still needs help with this one can always simulate a $y, x$ data set and inject an outlier at any particular x and follow the suggested steps to obtain a better estimate of $r$.

I welcome any comments on this as if it is "incorrect" I would sincerely like to know why hopefully supported by a numerical counter-example.


A small example will suffice to illustrate the proposed/transparent method of “obtaining of a version of r that is less sensitive to outliers” which is the direct question of the OP. This is an easy to follow script using standard ols and some simple arithmetic . Recall that B the ols regression coefficient is equal to r*[sigmay/sigmax).

Consider the following 10 pairs of observations.

enter image description here

And graphically

enter image description here

The simple correlation coefficient is .75 with sigmay = 18.41 and sigmax=.38

Now we compute a regression between y and x and obtain the following

enter image description here

Where 36.538 = .75*[18.41/.38] = r*[sigmay/sigmax]

The actual/fit table suggests an initial estimate of an outlier at observation 5 with value of 32.799 . enter image description here

If we exclude the 5th point we obtain the following regression result

enter image description here

Which yields a prediction of 173.31 using the x value 13.61 . This prediction then suggests a refined estimate of the outlier to be as follows ; 209-173.31 = 35.69 .

If we now restore the original 10 values but replace the value of y at period 5 (209) by the estimated/cleansed value 173.31 we obtain enter image description here

and enter image description here

Recomputed r we get the value .98 from the regression equation

r= B*[sigmax/sigmay] .98 = [37.4792]*[ .38/14.71]

Thus we now have a version or r (r =.98) that is less sensitive to an identified outlier at observation 5 . N.B. that the sigmay used above (14.71) is based on the adjusted y at period 5 and not the original contaminated sigmay (18.41). The effect of the outlier is large due to it's estimated size and the sample size. What we had was 9 pairs of readings (1-4;6-10) that were highly correlated but the standard r was obfuscated/distorted by the outlier at obervation 5.

There is a less transparent but nore powerfiul approach to resolving this and that is to use the TSAY procedure http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html to search for and resolve any and all outliers in one pass. For example enter image description here suggsts that the outlier value is 36.4481 thus the adjusted value (one-sided) is 172.5419 . Similar output would generate an actual/cleansed graph or table. enter image description here . Tsay's procedure actually iterativel checks each and every point for " statistical importance" and then selects the best point requiring adjustment. Time series solutions are immediately applicable if there is no time structure evidented or potentially assumed in the data. What I did was to supress the incorporation of any time series filter as I had domain knowledge/"knew" that it was captured in a cross-sectional i.e.non-longitudinal manner.

  • 2
    $\begingroup$ What does correlation have to do with time series, "pulses," "level shifts", and "seasonal pulses"? $\endgroup$
    – whuber
    Nov 14, 2016 at 22:58
  • 4
    $\begingroup$ Since time is not involved in regression in general, even something as simple as an autocorrelation coefficient isn't even defined. You cannot make every statistical problem look like a time series analysis! $\endgroup$
    – whuber
    Nov 14, 2016 at 23:07
  • 4
    $\begingroup$ @Engr I'm afraid this answer begs the question. It has several problems, of which the largest is that it provides no procedure to identify an "outlier." Another is that the proposal to iterate the procedure is invalid--for many outlier detection procedures, it will reduce the dataset to just a pair of points. $\endgroup$
    – whuber
    Nov 15, 2016 at 16:52
  • 4
    $\begingroup$ I fear that the present proposal is inherently dangerous, especially to naive or inexperienced users, for at least the following reasons (1) how to identify outliers objectively (2) the likely outcome is too complicated models based on ad hoc decisions (3) the procedure may not converge, or not converge nicely. Beginners usually over-identify outliers and make too little use of transformations and/or non-identity link functions as ways of taming them. $\endgroup$
    – Nick Cox
    Nov 15, 2016 at 17:35
  • 3
    $\begingroup$ No offence intended, @Carl, but you're in a mood to rant, and I am not and I am trying to disengage here. If it's the other way round, and it can be, I am not surprised if people ignore me. If I appear to be implying that transformation solves all problems, then be assured that I do not mean that. $\endgroup$
    – Nick Cox
    Nov 15, 2016 at 20:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.