Identifying Outliers with Quantile Regression and Python

Question

I am trying to identify outliers in a dataset using the 5th and 95th percentiles of a regression line so I'm using quantile regression in Python with statsmodel, matplotlib and pandas. Based on this answer from blokeley, I can create a scatterplot of my data and show the best fit line and the lines for the 5th and 95th percentile based on quantile regression. But how do I identify those points that fall above and below those lines and then save them out to a pandas dataframe?

My data looks like this (there are 95 values in total):

    Month   Year    LST     NDVI
0   June    1984    310.550975  0.344335
1   June    1985    310.495331  0.320504
2   June    1986    306.820900  0.369494
3   June    1987    308.945602  0.369946
4   June    1988    308.694022  0.31863

2

and the script I have so far is this:

import pandas as pd
excel = my_excel
df = pd.read_excel(excel)
df.head()

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

model = smf.quantreg('NDVI ~ LST',df)
quantiles = [0.05,0.95]
fits = [model.fit(q=q) for q in quantiles]
figure,axes = plt.subplots()
x = df['LST']
y = df['NDVI']
axes.scatter(x,df['NDVI'],c='green',alpha=0.3,label='data point')
fit = np.polyfit(x, y, deg=1)
axes.plot(x, fit[0] * x + fit[1], color='grey',label='best fit')
_x = np.linspace(x.min(),x.max())
for index, quantile in enumerate(quantiles):
    _y = fits[index].params['LST'] * _x + fits[index].params['Intercept']
    axes.plot(_x, _y, label=quantile)

title = 'LST/NDVI Jun-Aug'
plt.title(title)
axes.legend()
axes.set_xticks(np.arange(298,320,4))
axes.set_yticks(np.arange(0.25,0.5,.05))
axes.set_xlabel('LST')
axes.set_ylabel('NDVI');

And the chart I get out of that is this:

So I can definitely see data points above the 95th line and below the 5th line that I would classify as outliers, but I want to identify those in my original dataframe and maybe plot them on the cart or highlight them in some way to show them as "outliers".

I am searching on a method but coming up empty and could use some help.

Answer 1

You need to figure out if certain point are above the 95% quantile line or below the 5% quantile line. This you can do using the cross product, see this answer for a straightforward implementation.

In your example, you would need to combine the points above and below the quantile lines, possibly in a mask.

Here's is an example:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

df = pd.DataFrame(np.random.normal(0, 1, (100, 2)))
df.columns = ['LST', 'NDVI']

model = smf.quantreg('NDVI ~ LST', df)
quantiles = [0.05, 0.95]
fits = [model.fit(q=q) for q in quantiles]
figure, axes = plt.subplots()
x = df['LST']
y = df['NDVI']

fit = np.polyfit(x, y, deg=1)
_x = np.linspace(x.min(), x.max(), num=len(y))

# fit lines
_y_005 = fits[0].params['LST'] * _x + fits[0].params['Intercept']
_y_095 = fits[1].params['LST'] * _x + fits[1].params['Intercept']

# start and end coordinates of fit lines
p = np.column_stack((x, y))
a = np.array([_x[0], _y_005[0]]) #first point of 0.05 quantile fit line
b = np.array([_x[-1], _y_005[-1]]) #last point of 0.05 quantile fit line

a_ = np.array([_x[0], _y_095[0]])
b_ = np.array([_x[-1], _y_095[-1]])

#mask based on if coordinates are above 0.95 or below 0.05 quantile fitlines using cross product
mask = lambda p, a, b, a_, b_: (np.cross(p-a, b-a) > 0) | (np.cross(p-a_, b_-a_) < 0)
mask = mask(p, a, b, a_, b_)

axes.scatter(x[mask], df['NDVI'][mask], facecolor='r', edgecolor='none', alpha=0.3, label='data point')
axes.scatter(x[~mask], df['NDVI'][~mask], facecolor='g', edgecolor='none', alpha=0.3, label='data point')

axes.plot(x, fit[0] * x + fit[1], label='best fit', c='lightgrey')
axes.plot(_x, _y_095, label=quantiles[1], c='orange')
axes.plot(_x, _y_005, label=quantiles[0], c='lightblue')

axes.legend()
axes.set_xlabel('LST')
axes.set_ylabel('NDVI')

plt.show()