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

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.

• If the distribution has a long right tail, and you want a representative sample, then you can hope the sample has 'outliers' to the right. Why do you need to identify them? Surely, not to eliminate them, thus making a non-representative sample--I hope. – BruceET May 25 at 17:52

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