Counterfactual Fairness

Building counterfactually fair models

Data

To evaluate counterfactual fairness we will be using the “law school” dataset (McIntyre and Simkovic 2018).

The Law School Admission Council conducted a survey across 163 law schools in the United States. It contains information on 21,790 law students such as their entrance exam scores (LSAT), their grade-point average (GPA) collected prior to law school, and their first year average grade (FYA). Given this data, a school may wish to predict if an applicant will have a high FYA. The school would also like to make sure these predictions are not biased by an individual’s race and sex. However, the LSAT, GPA, and FYA scores, may be biased due to social factors.

We start by importing the data into a Pandas DataFrame.

import pandas as pd

df = pd.read_csv("_data/law_data.csv", index_col=0)
df.head()

Pre-processing

We now pre-process the data. We start by creating categorical “dummy” variables according to the race variable.

df = pd.get_dummies(df, columns=["race"], prefix="", prefix_sep="")
df.iloc[:, : 7].head()

We also want to expand the sex variable into male / female categorical variables and remove the original.

df["male"] = df["sex"].map(lambda x: 1 if x == 2 else 0)
df["female"] = df["sex"].map(lambda x: 1 if x == 1 else 0)
df = df.drop(axis=1, columns=["sex"])
df.iloc[:, 0:7].head()

We will also convert the entrance exam scores (LSAT) to a discrete variable.

df["LSAT"] = df["LSAT"].astype(int)
df.iloc[:, :6].head()

Protected attributes

Counterfactual fairness enforces that a distribution over possible predictions for an individual should remain unchanged in a world where an individual’s protected attributes $A$ had been different in a causal sense. Let’s start by defining the /protected attributes/. Obvious candidates are the different categorical variables for ethnicity (Asian, White, Black, etc) and gender (male, female).

A = [
    "Amerindian",
    "Asian",
    "Black",
    "Hispanic",
    "Mexican",
    "Other",
    "Puertorican",
    "White",
    "male",
    "female",
]

Training and testing subsets

We will now divide the dataset into training and testing subsets. We will use the same ratio as in (Kusner et al. 2017), that is 20%.

from sklearn.model_selection import train_test_split

df_train, df_test = train_test_split(df, random_state=23, test_size=0.2);

Models

Unfair model

As detailed in (Kusner et al. 2017), the concept of counterfactual fairness holds under three levels of assumptions of increasing strength.

The first of such levels is Level 1, where $\hat{Y}$ is built using only the observable non-descendants of $A$. This only requires partial causal ordering and no further causal assumptions, but in many problems there will be few, if any, observables which are not descendants of protected demographic factors.

For this dataset, since LSAT, GPA, and FYA are all biased by ethnicity and gender, we cannot use any observed features to construct a Level 1 counterfactually fair predictor as described in Level 1.

Instead (and in order to compare the performance with Level 2 and 3 models) we will build two unfair baselines.

• A Full model, which will be trained with the totality of the variables
• An Unaware model (FTU), which will be trained will all the variables, except the protected attributes $A$.

Let’s proceed with calculating the Full model.

Full model

As mentioned previously, the full model will be a simple linear regression in order to predict ZFYA using all of the variables.

from sklearn.linear_model import LinearRegression

linreg_unfair = LinearRegression()

The inputs will then be the totality of the variabes (protected variables $A$, as well as UGPA and LSAT).

import numpy as np

X = np.hstack(
    (
        df_train[A],
        np.array(df_train["UGPA"]).reshape(-1, 1),
        np.array(df_train["LSAT"]).reshape(-1, 1),
    )
)
X

As for our target, we are trying to predict _ZFYA (first year average grade).

y = df_train["ZFYA"]
y[:10]

We fit the model:

linreg_unfair = linreg_unfair.fit(X, y)

And perform some predictions on the test subset.

X_test = np.hstack(
    (
        df_test[A],
        np.array(df_test["UGPA"]).reshape(-1, 1),
        np.array(df_test["LSAT"]).reshape(-1, 1),
    )
)
X_test

predictions_unfair = linreg_unfair.predict(X_test)
predictions_unfair

We will also calculate the /unfair model/ score for future use.

score_unfair = linreg_unfair.score(X_test, df_test["ZFYA"])
score_unfair

from sklearn.metrics import mean_squared_error

RMSE_unfair = np.sqrt(mean_squared_error(df_test["ZFYA"], predictions_unfair))
float(RMSE_unfair)

Fairness through unawareness (FTU)

As also mentioned in (Kusner et al. 2017), the second baseline we will use is an Unaware model (FTU), which will be trained will all the variables, except the protected attributes $A$.

linreg_ftu = LinearRegression()

We will create the inputs as previously, but without using the protected attributes, $A$.

X_ftu = np.hstack(
    (
        np.array(df_train["UGPA"]).reshape(-1, 1),
        np.array(df_train["LSAT"]).reshape(-1, 1),
    )
)
X_ftu

And we fit the model:

linreg_ftu = linreg_ftu.fit(X_ftu, y)

Again, let’s perform some predictions on the test subset.

X_ftu_test = np.hstack(
    (np.array(df_test["UGPA"]).reshape(-1, 1), np.array(df_test["LSAT"]).reshape(-1, 1))
)
X_ftu_test

predictions_ftu = linreg_ftu.predict(X_ftu_test)
predictions_ftu

As previously, let’s calculate this model’s score.

ftu_score = linreg_ftu.score(X_ftu_test, df_test["ZFYA"])
print(ftu_score)

RMSE_ftu = np.sqrt(mean_squared_error(df_test["ZFYA"], predictions_ftu))
print(RMSE_ftu)

Latent variable model

Still according to (Kusner et al. 2017), a Level 2 approach will model latent ‘fair’ variables which are parents of observed variables.

If we consider a predictor parameterised by $\theta$, such as:

0

with $X_{\nsucc A} \subseteq X$ are non-descendants of $A$. Assuming a loss function $l(.,.)$ and training data $\mathcal{D}\equiv\{(A^{(i), X^{(i)}, Y^{(i)}})\}$, for $i=1,2\dots,n$, the empirical loss is defined as

1

which has to be minimised in order to $\theta$. Each $n$ expectation is with respect to random variable $U^{(i)}$ such that

2

where $P_{\mathcal{M}}(U|x,a)$ is the conditional distribution of the background variables as given by a causal model $M$ that is available by assumption.

If this expectation cannot be calculated analytically, Markov chain Monte Carlo (MCMC) can be used to approximate it as in the following algorithm.

We will follow the model specified in the original paper, where the latent variable considered is $K$, which represents a student’s knowledge. $K$ will affect GPA, LSAT and the outcome, FYA. The model can be defined by:

3

The priors used will be:

4

import pymc as pm

K = len(A)

def MCMC(data, samples=100):  # Reduced default samples for faster execution

    N = len(data)
    # Convert to numeric array explicitly to avoid object dtype
    a = np.array(data[A], dtype=np.float64)

    model = pm.Model()

    with model:
        # Priors
        k = pm.Normal("k", mu=0, sigma=1, shape=(N,))
        gpa0 = pm.Normal("gpa0", mu=0, sigma=1)
        lsat0 = pm.Normal("lsat0", mu=0, sigma=1)
        w_k_gpa = pm.Normal("w_k_gpa", mu=0, sigma=1)
        w_k_lsat = pm.Normal("w_k_lsat", mu=0, sigma=1)
        w_k_zfya = pm.Normal("w_k_zfya", mu=0, sigma=1)

        w_a_gpa = pm.Normal("w_a_gpa", mu=np.zeros(K), sigma=np.ones(K), shape=K)
        w_a_lsat = pm.Normal("w_a_lsat", mu=np.zeros(K), sigma=np.ones(K), shape=K)
        w_a_zfya = pm.Normal("w_a_zfya", mu=np.zeros(K), sigma=np.ones(K), shape=K)

        sigma_gpa_2 = pm.InverseGamma("sigma_gpa_2", alpha=1, beta=1)

        mu = gpa0 + (w_k_gpa * k) + pm.math.dot(a, w_a_gpa)

        # Observed data
        gpa = pm.Normal(
            "gpa",
            mu=mu,
            sigma=pm.math.sqrt(sigma_gpa_2),
            observed=data["UGPA"].values,
        )
        lsat = pm.Poisson(
            "lsat",
            mu=pm.math.exp(lsat0 + w_k_lsat * k + pm.math.dot(a, w_a_lsat)),
            observed=data["LSAT"].values,
        )
        zfya = pm.Normal(
            "zfya",
            mu=w_k_zfya * k + pm.math.dot(a, w_a_zfya),
            sigma=1,
            observed=data["ZFYA"].values,
        )

        # Optimize sampling: reduce tuning steps, use parallel chains, faster target_accept
        trace = pm.sample(
            draws=samples,
            tune=100,  # Reduce tuning steps (default is 1000)
            cores=8,  # Use parallel chains
            chains=2,  # Use 2 chains instead of default 4
            target_accept=0.8,  # Slightly lower than default 0.95 for faster sampling
            progressbar=True,
            return_inferencedata=False,
        )

    return trace

Let’s plot a single trace for $k^{(i)}$.

train_k = np.mean(train_estimates["k"], axis=0).reshape(-1, 1)
train_k

We can now estimate $k$ using the test data:

#| eval: false
test_map_estimates = MCMC(df_test, samples=50)

test_k = np.mean(test_map_estimates["k"], axis=0).reshape(-1, 1)
test_k

We now build the Level 2 predictor, using $k$ as the input.

linreg_latent = LinearRegression()

linreg_latent = linreg_latent.fit(train_k, df_train["ZFYA"])

predictions_latent = linreg_latent.predict(test_k)
predictions_latent

latent_score = linreg_latent.score(test_k, df_test["ZFYA"])
print(latent_score)

RMSE_latent = np.sqrt(mean_squared_error(df_test["ZFYA"], predictions_latent))
print(RMSE_latent)

Additive error model

Finally, in Level 3, we model GPA, LSAT, and FYA as continuous variables with additive error terms independent of race and sex (though these error terms may in turn be correlated with one-another).

This corresponds to

5

We estimate the error terms $\epsilon_{GPA}, \epsilon_{LSAT}$ by first fitting two models that each use race and sex to individually predict GPA and LSAT. We then compute the residuals of each model (e.g., $\epsilon_{GPA} =GPA−\hat{Y}_{GPA}(R, S)$). We use these residual estimates of $\epsilon_{GPA}, \epsilon_{LSAT}$ to predict $FYA$. In (Kusner et al. 2017) this is called Fair Add.

Since the process is similar for the individual predictions for GPA and LSAT, we will write a method to avoid repetion.

def calculate_epsilon(data, var_name, protected_attr):
    X = data[protected_attr]
    y = data[var_name]

    linreg = LinearRegression()
    linreg = linreg.fit(X, y)

    predictions = linreg.predict(X)

    return data[var_name] - predictions

Let’s apply it to each variable, individually. First we calculate $\epsilon_{GPA}$:

epsilons_gpa = calculate_epsilon(df, "UGPA", A)
epsilons_gpa

Next, we calculate $\epsilon_{LSAT}$:

epsilons_LSAT = calculate_epsilon(df, "LSAT", A)
epsilons_LSAT

Let’s visualise the $\epsilon$ distribution quickly:

We finally use the calculated $\epsilon$ to train a model in order to predict FYA. We start by getting the subset of the $\epsilon$ which match the training indices.

X = np.hstack(
    (
        np.array(epsilons_gpa[df_train.index]).reshape(-1, 1),
        np.array(epsilons_LSAT[df_train.index]).reshape(-1, 1),
    )
)
X

linreg_fair_add = LinearRegression()

linreg_fair_add = linreg_fair_add.fit(
    X,
    df_train["ZFYA"],
)

We now use this model to calculate the predictions

X_test = np.hstack(
    (
        np.array(epsilons_gpa[df_test.index]).reshape(-1, 1),
        np.array(epsilons_LSAT[df_test.index]).reshape(-1, 1),
    )
)

predictions_fair_add = linreg_fair_add.predict(X_test)
predictions_fair_add

And as previously, we calculate the model’s score:

fair_add_score = linreg_fair_add.score(X_test, df_test["ZFYA"])
print(fair_add_score)

RMSE_fair_add = np.sqrt(mean_squared_error(df_test["ZFYA"], predictions_fair_add))
print(RMSE_fair_add)

Comparison

The scores, so far, are:

print(f"Unfair score:\t{score_unfair}")
print(f"FTU score:\t{ftu_score}")
print(f"L2 score:\t{latent_score}")
print(f"Fair add score:\t{fair_add_score}")

print(f"Unfair RMSE:\t{RMSE_unfair}")
print(f"FTU RMSE:\t{RMSE_ftu}")
print(f"L2 RMSE:\t{RMSE_latent}")
print(f"Fair add RMSE:\t{RMSE_fair_add}")

Measuring counterfactual fairness

First, we will measure two quantities, the Statistical Parity Difference (SPD)11 See Statistical parity difference. and Disparate impact (DI)22 See Disparate Impact Ratio. .

Statistical Parity Difference / Disparate Impact

from fairlearn.metrics import demographic_parity_difference, demographic_parity_ratio

parities = []
impacts = []

for a in A:
    parity = demographic_parity_difference(df_train["ZFYA"], df_train["ZFYA"],
                                                sensitive_features = df_train[a])
    di = demographic_parity_ratio(df_train["ZFYA"], df_train["ZFYA"],
                                                sensitive_features = df_train[a])
    parities.append(parity)
    impacts.append(di)

df_parities = pd.DataFrame({'protected':A,'parity':parities,'impact':impacts})

Finding sensitive features

Typically a $SPD > 0.1$ and a $DI < 0.9$ might indicate discrimination on those features. In our dataset, we have two features (_Hispanic and _Mexican) which clearly fail the DI test (with $DI < 0.9$), indicating potential discrimination. The SPD values for these features are relatively small, suggesting that while there may be disparate impact, the absolute difference in proportions is not as pronounced.

for a in ["Mexican", "Hispanic"]:
    spd = demographic_parity_difference(y_true=df_train["ZFYA"],
                                        y_pred=df_train["ZFYA"],
                                        sensitive_features = df_train[a])
    print(f"SPD({a}) = {spd}")
    di = demographic_parity_ratio(y_true=df_train["ZFYA"],
                                  y_pred=df_train["ZFYA"],
                                  sensitive_features = df_train[a])
    print(f"DI({a}) = {di}")