Outlier Detection on skewed Distributions

Question

Under a classical definition of an outlier as a data point outide the 1.5* IQR from the upper or lower quartile, there is an assumption of a non-skewed distribution. For skewed distributions (Exponential, Poisson, Geometric, etc) is the best way to detect an outlier by analyzing a transform of the original function?

For example, distributions loosely governed by an exponential distribution, could be transformed with a log function - at which point is it acceptable to look for outliers based on the same IQR definition?

Answer 1

Under a classical definition of an outlier as a data point outide the 1.5* IQR from the upper or lower quartile,

This is the rule for identifying points outside the ends of the whiskers in a boxplot. Tukey himself would no doubt object to calling them outliers on this basis (he didn't necessarily regard points outside those limits as outliers). These would rather be points which - if your data was expected to be from a distribution somewhat similar to a normal distribution - one might subject to further investigation (such as checking you didn't transpose two digits, for example) -- at most these could be potential outliers. As Nick Cox points out in comments under this answer, a tail of many such points would be taken more as a indicator that a re-expression might be suitable than an indication of the need to regard the points as outliers.

there is an assumption of a non-skewed distribution.

I assumed by 'non-skewed' you mean symmetric. The two are not synonyms, but if you really mean non-skewed, we'd have to discuss with care what measure of skewness you're using, and we'd get tied up with issues of asymmetric cases with zero skewness. At least the symmetric case is not reliant on a particular choice of measurement and simpler to discuss. Then the assumption is more than just that it's symmetric. A heavy-tailed but symmetric distribution might have many points outside the bounds on that rule.

For skewed distributions (Exponential, Poisson, Geometric, etc) is the best way to detect an outlier by analyzing a transform of the original function?

That depends on what constitutes an outlier for your purposes. There's no single definition that's suitable for each purpose - indeed, generally you're probably better off doing other things that (say) picking outliers and omitting them.

For the exponential or geometric, you might do a similar calculation to that for a boxplot, but which would identify a similar fraction in the right tail only (you won't have low-end points identified in an exponential or geometric)$^{\dagger}$ ... or you might do something else.

$\dagger$ In large samples from a normal distribution, the boxplot marks about 0.35% of points at each end, or about 0.7% in total. For an exponential you might mark some multiple of the median, for example. If you wanted to tag roughly 0.7% of points in total for an actual exponential, that would suggest marking points beyond about 7.1 times the median.

Marking points above 7.1 times the median for n=1000 will typically hit between 0.4% to 1.1% of values:

ae <- rexp(1000)
table( ae > 7.1*median(ae) )

FALSE  TRUE 
  993     7 

For example, distributions loosely governed by an exponential distribution, could be transformed with a log function - at which point is it acceptable to look for outliers based on the same IQR definition?

That totally depends on what you mean by "acceptable". Note, however that -

i) the resulting distribution (from taking logs of an exponential variable) isn't actually symmetric, but distinctly left-skew.

As a result, you'll usually only mark points in the left end (i.e. close to zero, where you expect exponential values to be anyway) rather than in the right (where the "outliers" might be), unless they're really extreme.

ii) suitability of such a rule is going to be heavily dependent on what you're doing.

If you're concerned about the odd strange value affecting your inference, in general, you're probably better off using robust procedures than formally identifying outliers.

If you really do want to use a normal-based rule for transformed exponential or Poisson data, I'd at least suggest applying it to the square root$^{\ddagger}$ for a Poisson (as long as the mean isn't too small, it should be roughly normalish) and to cube root or even fourth root for the exponential (and perhaps, by extension, the geometric).

$\ddagger$ or perhaps $\sqrt{X+\frac{3}{8}}$, as in the Anscombe transform

For an exponential, in large samples the cube-root approach will tend to mark points only in the upper tail (at roughly the same rate it marks them in the upper tail for a normal) and the fourth-root approach marks points in both tails (slightly more in the lower tail, in total at something near 40% of the rate it does so for a normal). Of the possibilities, the cube root makes more sense to me than the other two, but I wouldn't necessarily advise using this as some hard and fast rule.

As noted in some of my answers to other questions, a Weibull distribution is close to symmetric for shape parameter in a broad interval near 3.6 (exactly where it's closest depends on how you measure it), and fairly close to a normal. This suggests that rather than the cube root or fourth root with exponential variavles, something in the vicinity the 3.6-th root should very near symmetric, but cube roots and fourth roots are still quite reasonable.

Answer 2

I will answer your questions in the opposite order in which you asked them, so that the exposition proceeds from the specific to the general.

First, let us consider a situation where you can assume that except for a minority of outliers, the bulk of your data can be well described by a known distribution (in your case the exponential).

If $x$ has pdf:

$$p_X(x)=\sigma^{-1}\mbox{exp}\left(\frac{-(x-\theta)}{\sigma}\right),\;x>0;\sigma>0$$

then $x$ is said to follow an exponential distribution (the special case where we set $\theta=0$ is called the one-parameter or standard exponential distribution).

The usual MLE estimator of the parameters are [0,p 506]:

$$\hat{\theta}=\min_i x_i$$

and

$$\hat{\sigma}=\mbox{ave}_ix_i-\min_i x_i$$

Here is an example in R:

n<-100
theta<-1
sigma<-2
set.seed(123) #for reproducibility
x<-rexp(n,rate=1/sigma)+theta
mean(x)-min(x)

the MLE of $\sigma$ is $\approx2.08$.

Unfortunately, the MLE estimates are very sensitive to the presence of outliers. For example, if I corrupt the sample by replacing 20% of the $x_i$'s by $-x_i$:

m<-floor(0.2*n)
y<-x
y[1:m]<--y[1:m]
mean(y)-min(y)

the MLE of $\sigma$ based on the corrupted sample is now $\approx11.12$(!). As a second example, if I corrupt the sample by replacing 20% of the $x_i$'s by $100x_i$ (say if the decimal place was accidentally misplaced):

m<-floor(0.2*n)
z<-x
z[1:m]<-100*z[1:m]
mean(z)-min(z)

the MLE of $\sigma$ based on this second corrupted sample is now $\approx54$(!).

An alternative to the raw MLE is to (a) find the outliers using a robust outlier identification rule, (b) set them aside as spurious data and (c) compute the MLE on the non spurious part of the sample.

The most well known of these robust outlier identification rule is the med/mad rule proposed by Hampel[3] who attributed it to Gauss (I illustrated this rule here). In the med/mad rule, the rejection threshold are based on the assumption that the genuine observations in your sample are well approximated by a normal distribution.

Of course, if you have extra information (such as knowing that the distribution of the genuine observations is well approximated by a poisson distribution as in this example) there is nothing to prevent you from transforming your data and using the baseline outlier rejection rule (the med/mad) but this strikes me as a bit awkward to transform the data to preserve what is after all an ad-hoc rule.

It seems much more logical to me to preserve the data but adapt the rejection rules. Then, you would still use the 3 step procedure I described in the first link above, but with rejection threshold adapted to the distribution you suspect the good part of the data has. Below, I give the rejection rule in situations where the genuine observations are well fitted by an exponential distribution. In this case, you can construct good rejection thresholds using the following rule:

1) estimate $\theta$ using [1]:

$$\hat{\theta}'=\mbox{med}_ix_i-3.476\mbox{Qn}(x)\ln2$$

The Qn is a robust estimate of scatter that is not geared towards symmetric data. It is widely implemented, for example in the R package robustbase. For exponential distributed data, the Qn is multiplied by consistency factor of $\approx3.476$, see [1] for more details.

2) reject as spurious all observations outside of [2,p 188]

$$[\hat{\theta}',9(1+2/n)\mbox{med}_ix_i+\hat{\theta}']$$

(the factor 9 in the rule above is obtained as the 7.1 in Glen_b's answer above, but using a higher cut-off. The factor (1+2/n) is small sample correction factor that was derived by simulations in [2]. For large enough sample sizes, it is essentially equal to 1).

3) use the MLE on the non spurious data to estimate $\sigma$:

$$\hat{\sigma}'=\mbox{ave}_{i\in H}x_i-\mbox{min}_{i\in H}x_i$$

where $H=\{i:\hat{\theta}'\leq x_i \leq 9(1+2/n)\mbox{med}_ix_i+\hat{\theta}'\}$.

using this rule on the previous examples, you would get:

library(robustbase)
theta<-median(x)-Qn(x,constant=3.476)*log(2)
clean<-which(x>=theta & x<=9*(1+2/n)*median(x)+theta)
mean(x[clean])-min(x[clean])

the robust estimate of $\sigma$ is now $\approx2.05$ (very close to the MLE value when the data is clean). On the second example:

theta<-median(y)-Qn(y,constant=3.476)*log(2)
clean<-which(y>=theta & y<=9*(1+2/n)*median(y)+theta)
mean(y[clean])-min(y[clean])

The robust estimate of $\sigma$ is now $\approx2.2$ (very close to the value we would have gotten without the outliers).

On the third example:

theta<-median(z)-Qn(z,constant=3.476)*log(2)
clean<-which(z>=theta & z<=9*(1+2/n)*median(z)+theta)
mean(z[clean])-min(z[clean])

The robust estimate of $\sigma$ is now $\approx2.2$ (very close to the value we would have gotten without the outliers).

A side benefit of this approach is that it yields a subset of indexes of suspect observations which should be set aside from the rest of the data, perhaps to be studied as object of interest in their own right (the members of $\{i:i\notin H\}$).

Now, for the general case where you do not have a good candidate distribution to fit the bulk of your observations beyond knowing that a symmetric distribution won't do, you can use the adjusted boxplot[4]. This is a generalization of the boxplot that takes into account a (non parametric and outlier robust) measure of skewness of your data (so that when the bulk of the data is symmetric is collapses down to the usual boxplot). You can also check this answer for an illustration.

• [0] Johnson N. L., Kotz S., Balakrishnan N. (1994). Continuous Univariate Distributions, Volume 1, 2nd Edition.
• [1] Rousseeuw P. J. and Croux C. (1993). Alternatives to the Median Absolute Deviation. Journal of the American Statistical Association, Vol. 88, No. 424, pp. 1273--1283.
• [2] J. K. Patel, C. H. Kapadia, and D. B. Owen, Dekker (1976). Handbook of statistical distributions.
• [3] Hampel (1974). The Influence Curve and Its Role in Robust Estimation. Journal of the American Statistical Association Vol. 69, No. 346 (Jun., 1974), pp. 383-393.
• [4] Vandervieren, E., Hubert, M. (2004) "An adjusted boxplot for skewed distributions". Computational Statistics & Data Analysis Volume 52, Issue 12, 15 August 2008, Pages 5186–5201.

Answer 3

First, I'd question the definition, classical or otherwise. An "outlier" is a surprising point. Using any particular rule (even for symmetric distributions) is a flawed idea, especially nowadays when there are so many huge data sets. In a data set of (say) one million observations (not all that big, in some fields), there will be many many cases beyond the 1.5 IQR limit you cite, even if the distribution is perfectly normal.

Second, I'd suggest looking for outliers on the original data. It will nearly always be more intuitive. For instance, with income data, it is quite common to take logs. But even here I'd look for outliers on the original scale (dollars or euros or whatever) because we have a better feel for such numbers. (If you do take logs, I'd suggest log base 10, at least for outlier detection, because it is at least a little intuitive).

Third, when looking for outliers, beware of masking.

Finally, I am currently researching the "forward search" algorithm proposed by Atkinson and Riani for various sorts of data and problems. This looks very promising.