v0.1 (21.05.2021)
(c) Ali Zaidi
A heuristic introduction to removing "unsavoury" features from datasets to avoid using them for making predictions.
A rigorous meta-analysis of the biases that arise in datasets, the sociology of data collection, a statement on fairness or how to measure it, ensure it, etc.
DISCLAIMER: I did not discover this bias in the dataset, nor did I invent this technique of cleaning data. A brief search on the internet will reveal the vast amounts of work related to this. I'm using this dataset merely for illustratives purposes. This work is original as far as the code and text are concerned. That is to say, I wrote this piece myself and haven't copied code, text or metrics from anywhere else. Howver, this post is inspired by Vincent Warmerdam's work.
One of the founding fathers of modern NeuralNets and AI, Prof. Yan LeCunn, raised a twitter storm a little over a year ago over a tweet on bias in algorithms. In essense, the statement was: algorithms themselves aren't biased, data is biased.
A major debate raged on regarding who is to blame and/or responsible for ensuring "fairness" in AI-based systems, however fairness is defined or understood.
It is a fact that biases are commonplace in our society, whether deliberate or accidental. In case you're interested, there's an extensive list on Wikipedia.
Since data collection might be subject to some of these biases, they naturally become part of the data. As a responsible data scientist, it is import to understand the variables we are working with, where they come from, what they mean, and how they affect our predictions.
As an example, we'll be checking out a very common dataset that most people come across very early on in their data science journeys: the Boston House price dataset.
Almost everyone who's played with regression would have come across this dataset, where the objective is to predict the price of a house based on various attributes regarding the house and neighbourhood.
We'll start by having a quick look at the input features...
# Some common preliminaries
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston # Get the dataset
data = load_boston() # Load the data
First things first, inspect the keys in the data object:
data.keys()
dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])
"feauture_names" seems promising. Let's have a quick look at the features.
data.feature_names
array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
'TAX', 'PTRATIO', 'B', 'LSTAT'], dtype='<U7')
Ok. Well that's an abbreviated list. Not very informative. Let's try something else...
"DESCR" seems like a description. Cool.
print(data.DESCR)
.. _boston_dataset:
Boston house prices dataset
---------------------------
**Data Set Characteristics:**
:Number of Instances: 506
:Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.
:Attribute Information (in order):
- CRIM per capita crime rate by town
- ZN proportion of residential land zoned for lots over 25,000 sq.ft.
- INDUS proportion of non-retail business acres per town
- CHAS Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
- NOX nitric oxides concentration (parts per 10 million)
- RM average number of rooms per dwelling
- AGE proportion of owner-occupied units built prior to 1940
- DIS weighted distances to five Boston employment centres
- RAD index of accessibility to radial highways
- TAX full-value property-tax rate per $10,000
- PTRATIO pupil-teacher ratio by town
- B 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
- LSTAT % lower status of the population
- MEDV Median value of owner-occupied homes in $1000's
:Missing Attribute Values: None
:Creator: Harrison, D. and Rubinfeld, D.L.
This is a copy of UCI ML housing dataset.
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/
This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.
The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic
prices and the demand for clean air', J. Environ. Economics & Management,
vol.5, 81-102, 1978. Used in Belsley, Kuh & Welsch, 'Regression diagnostics
...', Wiley, 1980. N.B. Various transformations are used in the table on
pages 244-261 of the latter.
The Boston house-price data has been used in many machine learning papers that address regression
problems.
.. topic:: References
- Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.
- Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.
Now that is informative!
I'd like to point your attention to the last two features, though: >
B: The proportion of blacks by town.
LSTAT: % Lower status of the population
Erm. That's a little awkward and uncomfortable. Should we really be using these as features for making predictions? Ideally we don't want to.
Maybe we can start by seeing how they affect our target variable.
But before that, a tiny bit preprocessing and cleaning of the data, to make life easier. :)
plt.hist((data.target - data.target.mean())/data.target.std(), 100, density=True);
# Linear regression is sensitive to outliers.
# Let's get rid of anything that may be z>2.5 from the target variable
# Find where the target z-score is above 2.5
max_val = data.target.mean() + (2.5 * data.target.std())
valid_targets = np.where(data.target<=max_val)
# filter out those values
data.data = data.data[valid_targets[0],:]
data.target = data.target[valid_targets[0]]
# Another good step is to normalize the columns.
# This will prevents headaches later on!
from sklearn.preprocessing import normalize
data_norm = normalize(data.data, norm='l2', axis=0, copy=True)
# Quick check of normalized cols
np.allclose(np.diagonal(data_norm.T @ data_norm),1)
True
Great!
Moving on to the issue in question: the questionable columns...
# Relationship of LSTAT to Target
from scipy.stats import pearsonr
r,p = pearsonr(data_norm[:,-1], data.target)
plt.scatter(data_norm[:,-1], data.target)
plt.xlabel(data.feature_names[-1])
plt.ylabel('Price')
plt.title(f'r={r:0.3f}, p={p:0.2E}');
The overall trend suggests that the price drops with increasing proportion of 'low-income' families.
Let's have a look at the relationship of 'B' with price
# Relationship of B to Target
r, p = pearsonr(data_norm[:,-2], data.target)
plt.scatter(data_norm[:,-2], data.target)
plt.xlabel(data.feature_names[-2])
plt.ylabel('Price')
plt.title(f'r={r:0.3f}, p={p:0.2E}');
A positive trend? Not, bad! The price seems to increase with 'B'. So does that mean that the higher the proportion of 'black' residents, the higher the price?
Nope! This variable is a little more cryptic. Let's revisit the description:
...
- B: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
...
So, they're taking the difference between a value and 0.63 and then squaring it.
Which means, the further the value deviates from 63%, the higher the score will be, irrespective of whether the actual value increases or decreases.
Meaning higher values could very well indicate a much lower proportion of black residents. This is what we sometimes call a proxy-variable.
These kinds of scenarios pop up often in real world data, and we need to be very careful while interpreting input variables.
We're now going to fit a simple linear model to our dataset, and use the mean squared error as a measure of 'loss'.
Good statical practices would demand a separation of trainig and test sets, to ensure proper validation. However, the purpose of this exercise is to study and evaluate bias, so we will be a little lapse regarding validation here and use the entire dataset. Bad practice in general, but since we're interested in another metric, it will be good enough for now.
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error as mse
from sklearn.utils import resample
# Fit a linear model and get the MSE
reg = LinearRegression().fit(data_norm, data.target)
mse_raw = mse(reg.predict(data_norm), data.target)
print(f'MSE={mse_raw:0.3f}')
MSE=12.588
So far so good. Any tutorial on linear regression would be wrapping up right about now, but this is where our journey actually begins.
While thinking about bias (or fairness), a difficult task is to identify a proper measure to quantify it. The right metric is a contentious subject in itself, and definitely beyond the scope of this tutorial.
Since we're short on time, and mildly irresponsible, we need to care about a measure of bias, instead of the measure of bias.
For example, one way to measure bias would be to look at the differences between the model predictions for houses with low prices versus high prices. If the model performance differs, we will consider our model biased.
As a proxy for the amount of bias, let's look at the difference between the predictions for low vs high priced houses. However, as we are going to be playing around with the features (by modifying or removing them), we need a way to measure the bias relative to the modified features. A simple way to accomplish that is by normalizing the difference of low vs high by the overall prediction.
In short, we're looking at whether higher priced houses have a predictive advantate over lower priced houses, and we're normalizing that to the overall prediction to enable comparisions across various methods of "tweaking" our dataset.
The ideal case is to ensure "fair" predictions by having the model perform equally well for both types of houses.
Why this metric? Ideally, you'd want to measure the effect of the "biased" features on your dataset. But then how do we do that if we intend to remove them? We therefore need to find an aspect within our target variable to be able to calculate the bias. The one feature we have is the price of the house, and if we think deeply, the price is a proxy for economic disparity, which we are using a proxy for bias here. Again, this is only for illustrative purposes.
# Measure fairness
# Let's combine the X and y arrays to enable easier indexing
d = np.hstack((data_norm, data.target[:,np.newaxis])) # keep shape tuple same
# Step 1:
# Fit the full model (we've done this before but just bear with me)
reg_full = LinearRegression().fit(d[:,:-1], d[:,-1])
s_full = mse(reg_full.predict(d[:,:-1]), d[:,-1])
# Step 2:
# Divide the target into equal lower and higher ranges based
d_l_idx = np.where(d[:, -1] < np.median(d[:,-1]))
d_h_idx = np.where(d[:, -1] >= np.median(d[:,-1]))
d_l = d[d_l_idx]
d_h = d[d_h_idx]
# Perform separate regressions for low vs high value houses
reg_l = LinearRegression().fit(d_l[:,:-1], d_l[:,-1])
s_l = mse(reg_l.predict(d_l[:,:-1]), d_l[:,-1])
reg_h = LinearRegression().fit(d_h[:,:-1], d_h[:,-1])
s_h = mse(reg_h.predict(d_h[:,:-1]), d_h[:,-1])
# Step 4:
# Calculate the bias (high - low / full)
score_diff = 100*(s_l - s_h)/s_full
print(f'Bias: {score_diff:0.3f}')
Bias: -49.313
Now that we have a measure for bias, we'd like to get some statistics in order to compare different methods of De-biasing.
To do that, we can do a simple bootstrap (sampling with replacement), and get a distribution of bias scores.
# Combine the data and targets for easier sampling
data_full = np.hstack((data_norm, data.target[:, np.newaxis]))
# array to accumulate bias score
score_diff = np.array([])
# Run bootstrap
for i in range(1000):
d = resample(data_full)
# Obtain MSE for the current bootstrapped dataset
reg_full = LinearRegression().fit(d[:,:-1], d[:,-1])
s_full = mse(reg_full.predict(d[:,:-1]), d[:,-1])
# Divide the target into equal lower and higher ranges based
d_l_idx = np.where(d[:, -1] < np.median(d[:,-1]))
d_h_idx = np.where(d[:, -1] >= np.median(d[:,-1]))
d_l = d[d_l_idx]
d_h = d[d_h_idx]
# Perform separate regressions for low vs high value houses
reg_l = LinearRegression().fit(d_l[:,:-1], d_l[:,-1])
s_l = mse(reg_l.predict(d_l[:,:-1]), d_l[:,-1])
reg_h = LinearRegression().fit(d_h[:,:-1], d_h[:,-1])
s_h = mse(reg_h.predict(d_h[:,:-1]), d_h[:,-1])
score_diff = np.append(score_diff, 100*(s_l - s_h)/s_full)
# Plot the results
plt.subplots(1,2, figsize=(10, 5))
plt.subplot(1,2,1)
plt.scatter(reg_full.predict(data_norm), data.target)
plt.title(f'MSE full fit = {mse(reg_full.predict(data_norm), data.target):0.3f}');
plt.subplot(1,2,2)
plt.hist(score_diff, density=True, alpha=0.3)
plt.xlabel("bias")
plt.ylabel("freq");
plt.title(f'bias mean={score_diff.mean():0.3f}, std={score_diff.std():0.3f}');
We'll be doing this exercise A LOT. Modifying the dataset, fitting a regression line, bootstrapping the bias and ploting the results.
Writing these methods in a class is going to make life super easy.
The class below has all the things we've done above: namely, fitting a line, calculating MSE, bootstrapping and plotting. It inherits from the sklearn LinearRegression class, and adds bootstrapping and plotting to the functionality.
# One Class to Rule Them All!
class linreg(LinearRegression):
def __init__(self, X, y, *args, **kwargs):
"""
Wrapper over the scikit-learn LinearRegression.
Init calls fit() and also calculates MSE (mean squared error).
Enables creation of multiple dataset/model objects.
"""
# Initialize the parent class
super().__init__()
# Set params
self.debug=False if not "debug" in kwargs else True # useful for debugging
self.nreps=1000 if not "nreps" in kwargs else kwargs['nreps']
# Fit the data
self.fit(X, y)
self.preds = self.predict(X)
self.X = X
self.y = y
# Calculate MSE
self.mse = self.mean_squared_error(X, y)
#Always good to have some housekeeping
if self.debug: print(f'mse: {self.mse:0.3f}')
def mean_squared_error(self, X, y):
"Calculate the mean squared error between prediction and observation"
return np.sum((self.predict(X) - y)**2)/len(y)
def bootstrap_bias(self, *args, **kwargs):
"Run a bootstrap over the data to get a distribution of bias scores"
score_diff = np.array([])
# Combine the data and targets for easier sampling
data = np.hstack((self.X, self.y[:, np.newaxis]))
# Run bootstrap
for i in range(self.nreps):
d = resample(data)
# Obtain MSE for the current bootstrapped dataset
reg_full = LinearRegression().fit(d[:,:-1], d[:,-1])
s_full = mse(reg_full.predict(d[:,:-1]), d[:,-1])
# Divide the data into lower and higher prices based on target value
d_l_idx = np.where(d[:, -1] < np.median(d[:,-1]))
d_h_idx = np.where(d[:, -1] >= np.median(d[:,-1]))
d_l = d[d_l_idx]
d_h = d[d_h_idx]
# Perform separate regressions for low vs high value houses
reg_l = LinearRegression().fit(d_l[:,:-1], d_l[:,-1])
s_l = mse(reg_l.predict(d_l[:,:-1]), d_l[:,-1])
reg_h = LinearRegression().fit(d_h[:,:-1], d_h[:,-1])
s_h = mse(reg_h.predict(d_h[:,:-1]), d_h[:,-1])
score_diff = np.append(score_diff, 100*(s_l - s_h)/s_full)
self.score_diff = score_diff
if self.debug: print(f'bias: {score_diff.mean():0.3f}')
return score_diff
def plot(self):
"Plot the fit and bias scores"
plt.subplots(1,2, figsize=(10, 5))
plt.subplot(1,2,1)
plt.scatter(self.preds, self.y)
plt.title(f'MSE full fit = {self.mse:0.3f}');
plt.subplot(1,2,2)
plt.hist(self.score_diff, density=True, alpha=0.3)
plt.xlabel("bias")
plt.ylabel("freq");
plt.title(f'bias mean={self.score_diff.mean():0.3f}, std={self.score_diff.std():0.3f}')
# Let's see this in action using the raw features and time it.
%%time
raw_reg = linreg(data_norm, data.target, nreps=1000, debug=True)
raw_bias = raw_reg.bootstrap_bias()
raw_reg.plot()
mse: 12.588 bias: -45.037 CPU times: user 3.07 s, sys: 13.5 ms, total: 3.08 s Wall time: 3.08 s
Awesome! Now we're equipped to evaluate the effects of our interventions on the training data.
Our first idea is not to use the questionable features for constructing our model. So we'll just drop them. No bias in, no bias out, right?
The class we wrote above will help us test and quantify our hypothesis.
# Create politically correct features
# Get 'clean' features
pol_corr_cols = [i for i in range(len(data.feature_names)) if data.feature_names[i]!='B' and data.feature_names[i]!='LSTAT']
# Filter the data
pol_corr_data = data_norm[:, pol_corr_cols]
# Model it
pol_corr_reg = linreg(pol_corr_data, data.target, debug=True)
pol_corr_bias = pol_corr_reg.bootstrap_bias()
pol_corr_reg.plot()
mse: 15.419 bias: -28.762
plt.figure()
plt.hist(raw_bias, alpha=0.3, density=True, label='Raw Bias');
plt.hist(pol_corr_bias, alpha=0.3, density=True, label='Dropped Columns');
plt.axvline(raw_bias.mean(), c='b')
plt.axvline(pol_corr_bias.mean(), c='r')
plt.legend();
We can see we've gained some error but lost some bias. That's good. But our bias is still kind of high.
The question is, when we removed the columns, did also remove the information from the columns?
In case the other columns are correlated to the ones we dropped, we're in trouble, and simply dropping them won't do us much good.
How do we check for that? A simple correlation analysis would do. A complicated way is to use Linear regression again! If we can express our dropped columns as a linear combination of our politically correct columns, that means the information is still there!
Have a look at the $R^2$ values. Theres definitely an information leak:
# Fit clean data to 'LSTAT'
reg_pc = LinearRegression().fit(pol_corr_data, data_norm[:,-1])
print(reg_pc.score(pol_corr_data, data_norm[:,-1]))
# Fit clean data to 'B'
reg_pc = LinearRegression().fit(pol_corr_data, data_norm[:,-2])
print(reg_pc.score(pol_corr_data, data_norm[:,-2]))
0.6709661975363947 0.24306629879425456
Someone familiar with linear algebra might say, "Use an SVD. That orthogonalizes the columns"
It decomposes the data matrix D into three matrices:
$$ D = U \Sigma V^{T}$$Where $U$ is an orthonormal set of vectors.
This should orthogonalize out the unwanted data, shouldn't it?
Here's a test
# Get the orthonormal basis from the SVD
U, _, _ = np.linalg.svd(pol_corr_data, full_matrices=False)
# Test if the columns are orthogonal to the unwanted columns
U.T @ data_norm[:,(-1,-2)]
array([[-0.92092404, -0.91490534],
[ 0.15816429, -0.25946357],
[ 0.00450424, 0.02901175],
[ 0.07869565, 0.14152012],
[ 0.02250988, 0.13558084],
[ 0.08881227, -0.00821001],
[ 0.00677911, 0.01736526],
[-0.10568813, 0.05822104],
[-0.04101328, 0.00384456],
[-0.0192633 , 0.02302542],
[ 0.12959063, 0.01288281]])
# Run a linreg just to see if the bias changes
svd_reg = linreg(U, data.target, debug=True)
svd_reg.bootstrap_bias()
svd_reg.plot()
mse: 15.419 bias: -28.455
Absolutely no change in our results!
Bummer!
Removing columns doesn't remove information when the columns are codependent. Even if we do an SVD, the information still remains with high fidelity.
Along with removing the features, we need to remove the information of the features from the dataset. Then can we truly call it 'clean'!
How do you remove the information of codependent features? Linear Algebra of course. A simple process called orthogonalization.
In essense, what you're doing is making the clean data orthogonal to the features you'd like to remove.
For a detailed explanation, check out the Wikipedia article here.
But for our purposes this explanation should suffice.
The process is pretty simple: Subtract the projections of the unwanted columns from the rest to make clean ones.
The formula for orthogonalizing two vector is as follows:
$$ v = a - (b \cdot \frac{<a,b>}{||a||\cdot||b||} )$$where $<a,b>$ is the inner product of a and b.
Let's break it down.
Project b onto a :
$$ Proj_{a}(b)= a \cdot \frac{<a,b>}{||a||\cdot||b||} $$We can then subtract this from b. This is effectively equivalent to walking along b, and then walking along $-Proj_{a}(b)$.
Here's an illustration to show the process:
Let's write a small script to orthogonalize our features in two simple steps.
For each column, remove the projection of 'B' and then remove the projection of 'LSTAT'.
# Create an matrix to store the values of the clean columns
orth_data = np.ones((data_norm.shape[0], data_norm.shape[1]-2))
# Loop over each column and remove the projections of the unwanted cols
for i in range(orth_data.shape[1]):
# Remove projection of 'B' from each column and re-normalize
orth_data[:,i] = data_norm[:,i] - data_norm[:,-2]*np.dot(data_norm[:,i], data_norm[:,-2])
orth_data[:,i] = orth_data[:,i]/np.sqrt(np.dot(orth_data[:,i], orth_data[:,i]))
# Remove projection of 'LSTAT' from each column
orth_data[:,i] = orth_data[:,i] - (np.dot(orth_data[:,i], data_norm[:,-1])*data_norm[:,1])
orth_data[:,i] = orth_data[:,i]/np.sqrt(np.dot(orth_data[:,i], orth_data[:,i]))
orth_data = normalize(orth_data, norm='l2', axis=0)
# Let's be sure the columns are orthonormal
np.allclose(np.diagonal(orth_data.T @ orth_data),1)
True
# And then let's ensure they are also orthogonal to the unwanted columns
orth_data.T @ data_norm[:,(-1,-2)]
array([[ 0.23840696, -0.13797649],
[-0.13066592, 0.07562206],
[ 0.28720424, -0.1662176 ],
[ 0.01376148, -0.00796437],
[ 0.30507012, -0.17655736],
[ 0.20102537, -0.11634213],
[ 0.2934696 , -0.16984364],
[-0.10581986, 0.06124256],
[ 0.27025785, -0.15640999],
[ 0.29945689, -0.17330874],
[ 0.28894943, -0.16722762]])
Wait, what?! Our code snippet above clearly removes the projections of two unwanted columns one after the other from the remaining columns.
# Could it be that the unwanted columns aren't normalized?
np.allclose(np.diagonal(data_norm[:,(-1,-2)].T @ data_norm[:,(-1,-2)]),1)
True
That's not it. There's something wrong with the orthogonalization.
Again, the devil is in the details.
We did the process sequentially. However the orthonormal vectors are supposed to be generated recursively. This is an extremely and often painfully common error during Gram-Schmidt. And we have to been very careful about this!
An illustration of the right and wrong methods is below:
In brief: Sequentially: wrong; Recursively: right
Correcing our code to enable orthogonalization the proper way, let's take another hack at it.
orth_data = np.ones((data_norm.shape[0], data_norm.shape[1]-2))
for i in range(orth_data.shape[1]):
# Remove the projection of 'B' from the cols
orth_data[:,i] = data_norm[:,i] - data_norm[:,-2]*np.dot(data_norm[:,i], data_norm[:,-2])
orth_data[:,i] = orth_data[:,i]/np.sqrt(np.dot(orth_data[:,i], orth_data[:,i]))
# Remove the projection of 'B' on 'LSTAT', let's call this C. (Step 2 in illustration)
proj_b_lstat = data_norm[:,-1] - (np.dot(data_norm[:,-2], data_norm[:,-1])*data_norm[:,-2])
proj_b_lstat = proj_b_lstat/np.sqrt(np.dot(proj_b_lstat, proj_b_lstat))
# Remove this projection of C from all cols (step 3 in illustration)
orth_data[:,i] = orth_data[:,i] - proj_b_lstat*np.dot(orth_data[:,i], proj_b_lstat)
orth_data[:,i] = orth_data[:,i]/np.sqrt(np.dot(orth_data[:,i], orth_data[:,i]))
orth_data = normalize(orth_data, norm='l2', axis=0)
# test for orthonormal cols
np.allclose(np.diagonal(orth_data.T @ orth_data),1)
True
# test for orthogonality to unwanted cols
np.allclose(orth_data.T @ data_norm[:,(-1,-2)],0)
True
Yippie! We've removed the unwanted information from our features!
Time to see how the regression model performs
ortho_reg = linreg(orth_data, data.target, debug=True)
ortho_bias = ortho_reg.bootstrap_bias()
ortho_reg.plot()
mse: 44.409 bias: -2.697
# A comparision of the three different datasets
plt.figure()
plt.hist(raw_bias, 10, alpha=0.3, density=True, label="raw");
plt.axvline(raw_bias.mean(), c='b')
plt.hist(pol_corr_bias, 10, alpha=0.3, density=True, label='drop cols');
plt.axvline(pol_corr_bias.mean(), c='r')
plt.hist(ortho_bias, 10, alpha=0.3, density=True, label="Gram-Schmidt");
plt.axvline(ortho_bias.mean(), c='g')
plt.xlabel('Bias');
plt.ylabel('Freq');
plt.legend()
<matplotlib.legend.Legend at 0x7f1542ab3d10>
Awesome! We've effectively removed the difference in the accuracy of predictions between low and high priced houses! So this goes to show that our input features did affect the bias of our model, and that it is possible to identify and remove those features with a little bit of effort. Ok maybe a slightly involved amount of effort. But it is possible.
However! We have sacrificed accuracy (MSE) to gain fairness. This is an important lesson: fairness has a price.
The obvious question is, can we find a balance between model peformance and fairness?
Of course! And there are many, many ways to do that. A super simple idea is to predict using both the biased and fair models, and take a weighted average of their predictions.
Let's have a look at how that might work:
bias_weights = np.arange(0,1.01,0.01)
mse_s = np.array([])
bias_s = np.array([])
for i in bias_weights:
# Take weighted average of models
weighted_preds = ortho_reg.predict(orth_data)*(1-i) + raw_reg.predict(data_norm)*i
# Compute MSE for weighted average
mse_s = np.append(mse_s, np.sum((weighted_preds - data.target)**2)/len(data.target))
# Calculate the bias (a quick and dirty hack)
idx_l = data.target<np.median(data.target)
idx_h = data.target>=np.median(data.target)
s_l = mse(weighted_preds[idx_l], data.target[idx_l])
s_h = mse(weighted_preds[idx_h], data.target[idx_h])
bias_s = np.append(bias_s, 100*(s_l - s_h)/mse_s[-1])
# NOTE: This is not the way we calculated bias. Earlier we fitted separate models
# So the values aren't exactly identical. This is a rough estimate for bias.
plt.plot(bias_weights, mse_s, label="MSE")
plt.plot(bias_weights, bias_s, label="Bias proxy")
plt.grid()
plt.xlabel('Weight of biased predictions')
plt.ylabel('MSE / QDBias') # QD: quick and dirty
plt.legend()
<matplotlib.legend.Legend at 0x7f1542b50450>
We can now find an acceptable balance between model peformance and fairness.