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. 
 A: 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

