How to find the upper outlier threshold in a right skewed distribution?

Question

Say we have a right skewed distribution (tail on the right side) like income. Above what value can I consider data points to be outliers?

One way could be Q3 (75th percentile) + 1.5 * Interquartile range (IQR)

The other way could be using box-cox transform to convert into a normal distribution and finding the value corresponding to z-score of 3.

However, I would like to know if there are more elegant ways of finding outlier thresholds.

Answer 1

The answer depends on what you mean by 'outlier' and why you want to identify outliers.

In R, I generate a random sample x of size 50 from an exponential distribution with mean 1. The mean 1.11 of the sample is not a bad estimate of the population mean.

set.seed(2020)
x = rexp(50)
mean(x)
[1] 1.117136

In R, boxplot.stats with $out lists outliers. So here are the boxplot outliers in x:

boxplot.stats(x)$out
[1] 5.867519 4.572054 5.645287 3.238821

I decide to eliminate the outliers in x, making the new dataset x1:

x1 = x[x < min(boxplot.stats(x)$out)]

The mean of the truncated sample is not such a good estimate of the population mean.

mean(x1)
[1] 0.794198

But wait! There's more! The truncated dataset has a boxplot outlier of its own. Where do you want to terminate this process? In some sense that matters to you, is 2.6436 also an outlier?

boxplot.stats(x1)$out
[1] 2.643577

Boxplots of the original (left) and truncated samples.

Note: Using R, I simulated $100\,000$ samples of size $n=50$ from an exponential distribution with $\mu = \sigma = 1.$ All but 977 had boxplot outliers.

For the original samples, sample means averaged 1.001 and sample SDs averaged 0.991. When boxplot outliers (almost 5 per sample on average) were removed, the sample means averaged 0.849 and the sample SDs averaged 0.713. The R program for the simulation is shown below.

set.seed(1234)
m = 10^5;  n = 100
a = s = a1 = s1 = n1 = numeric(m)
for(i in 1:m){
 x=rexp(n); a[i]=mean(x); s[i]=sd(x)
 x1 = x[x <= boxplot.stats(x)$stats[5]]
 n1[i] = length(x1)
 if(length(x1)==n){a1[i]=a[i];s1[i]=s[i]
 } else  {a1[i] = mean(x1);  s1[i] = sd(x1)}
}
mean(a); mean(s)
[1] 1.000513
[1] 0.990878
mean(n1); sum(n1==n)
[1] 95.1565
[1] 977
mean(a1); mean(s1)
[1] 0.848945
[1] 0.7130943