Explaining a Decision, and Its Limits

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(3)
n = 8000

# A protected attribute, NOT given to the model.
group = rng.choice([0, 1], n, p=[0.7, 0.3])

# A feature strongly associated with it - a postcode, a school, a
# device. This is the proxy.
postcode_score = rng.normal(0, 1, n) + 1.5 * group

# The true outcome depends only on the honest feature.
ability = rng.normal(0, 1, n)
repaid = ((1.4 * ability + rng.normal(0, 0.8, n)) > 0).astype(int)

# But the HISTORICAL decisions were biased: group 1 was approved less
# often for the same ability, so the training labels encode that.
historical_approval = ((1.4 * ability - 0.9 * group + rng.normal(0, 0.7, n)) > 0).astype(int)

frame = pd.DataFrame({'ability': ability, 'postcode_score': postcode_score})
X_tr, X_te, y_tr, y_te, g_tr, g_te = train_test_split(
    frame, historical_approval, group, test_size=0.3, random_state=0)

model = HistGradientBoostingClassifier(random_state=0).fit(X_tr, y_tr)
predictions = model.predict(X_te)

print('the model never saw the protected attribute. Its predictions:')
rates = pd.DataFrame({'group': g_te, 'approved': predictions}).groupby('group')['approved'].agg(
    ['mean', 'size'])
rates.columns = ['approval_rate', 'n']
print(rates.round(3).to_string())
print()
print(f'approval rate gap: '
      f'{rates.loc[0, "approval_rate"] - rates.loc[1, "approval_rate"]:+.1%}')
print()
print(f'correlation between postcode_score and group: '
      f'{np.corrcoef(postcode_score, group)[0, 1]:.2f}')
print()
print('Dropping the protected attribute does not remove the bias. The')
print('proxy carries it, and the model finds the proxy.')

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(7)
n = 10_000

group = rng.choice([0, 1], n, p=[.65, .35])
ability = rng.normal(0, 1, n) - 0.35 * group        # a genuine difference in the data
proxy = ability + 0.8 * group + rng.normal(0, 0.5, n)
truth = ((1.5 * ability + rng.normal(0, 0.9, n)) > 0).astype(int)

X = pd.DataFrame({'proxy': proxy, 'noise': rng.normal(0, 1, n)})
X_tr, X_te, y_tr, y_te, g_tr, g_te = train_test_split(
    X, truth, group, test_size=0.3, random_state=0, stratify=truth)

model = HistGradientBoostingClassifier(random_state=0).fit(X_tr, y_tr)
scores = model.predict_proba(X_te)[:, 1]
predictions = (scores >= 0.5).astype(int)

rows = []
for g in [0, 1]:
    mask = g_te == g
    tp = int(((predictions == 1) & (y_te == 1) & mask).sum())
    fp = int(((predictions == 1) & (y_te == 0) & mask).sum())
    fn = int(((predictions == 0) & (y_te == 1) & mask).sum())
    tn = int(((predictions == 0) & (y_te == 0) & mask).sum())
    rows.append({
        'group': g, 'n': int(mask.sum()),
        'selection_rate': round((tp + fp) / mask.sum(), 3),
        'true_positive_rate': round(tp / max(tp + fn, 1), 3),
        'false_positive_rate': round(fp / max(fp + tn, 1), 3),
        'precision': round(tp / max(tp + fp, 1), 3),
        'base_rate': round(float(y_te[mask].mean()), 3),
    })
table = pd.DataFrame(rows)
print(table.to_string(index=False))
print()
print(f'demographic parity gap (selection rates): '
      f'{table["selection_rate"].iloc[0] - table["selection_rate"].iloc[1]:+.3f}')
print(f'equal opportunity gap (true positive rates): '
      f'{table["true_positive_rate"].iloc[0] - table["true_positive_rate"].iloc[1]:+.3f}')
print(f'predictive parity gap (precision): '
      f'{table["precision"].iloc[0] - table["precision"].iloc[1]:+.3f}')

import numpy as np
import pandas as pd

# Two groups with genuinely different base rates. Satisfying any two
# of the three definitions forces the third to be violated.
print(pd.DataFrame([
    {'base rates': 'equal', 'can satisfy all three': 'yes'},
    {'base rates': 'different', 'can satisfy all three': 'no - proved impossible'},
]).to_string(index=False))
print()
print('This is a theorem, not an engineering limitation. You must')
print('choose which definition matters for your decision, and say so.')

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.inspection import permutation_importance, partial_dependence
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(11)
n = 4000
frame = pd.DataFrame({
    'income': rng.lognormal(10.2, 0.5, n).round(-2),
    'loan': rng.lognormal(9.0, 0.6, n).round(-2),
    'years_employed': rng.gamma(2.2, 2.4, n).round(1),
    'prior_defaults': rng.poisson(0.25, n),
})
frame['ratio'] = frame['loan'] / frame['income']
risk = (-1.0 + 3.0 * frame['ratio'] + 0.8 * frame['prior_defaults']
        - 0.08 * frame['years_employed'] + rng.normal(0, 0.8, n))
y = (risk > 0).astype(int)

FEATURES = ['income', 'loan', 'years_employed', 'prior_defaults', 'ratio']
X_tr, X_te, y_tr, y_te = train_test_split(frame[FEATURES], y, test_size=0.3,
                                          random_state=0, stratify=y)
model = HistGradientBoostingClassifier(random_state=0).fit(X_tr, y_tr)

# GLOBAL: which features the model relies on overall.
importance = permutation_importance(model, X_te, y_te, n_repeats=10, random_state=0)
print('global importance:')
print(pd.DataFrame({'feature': FEATURES,
                    'importance': importance.importances_mean.round(4)}
                   ).sort_values('importance', ascending=False).to_string(index=False))
print()

# GLOBAL: the shape of the relationship, averaged over everyone.
result = partial_dependence(model, X_te, features=['ratio'], grid_resolution=6)
print('partial dependence on the loan-to-income ratio:')
for value, average in zip(result['grid_values'][0], result['average'][0]):
    print(f'  ratio {value:.2f} -> mean predicted risk {average:.3f}')
print()
print('Both of these describe the MODEL. Neither describes an')
print('individual decision, which is what a person who was refused')
print('actually wants to know.')

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(13)
n = 4000
frame = pd.DataFrame({
    'income': rng.lognormal(10.2, 0.5, n).round(-2),
    'loan': rng.lognormal(9.0, 0.6, n).round(-2),
    'years_employed': rng.gamma(2.2, 2.4, n).round(1),
    'prior_defaults': rng.poisson(0.25, n),
})
frame['ratio'] = (frame['loan'] / frame['income']).round(3)
risk = (-1.0 + 3.0 * frame['ratio'] + 0.8 * frame['prior_defaults']
        - 0.08 * frame['years_employed'] + rng.normal(0, 0.8, n))
y = (risk > 0).astype(int)

FEATURES = ['income', 'loan', 'years_employed', 'prior_defaults', 'ratio']
X_tr, X_te, y_tr, y_te = train_test_split(frame[FEATURES], y, test_size=0.3, random_state=0)
model = HistGradientBoostingClassifier(random_state=0).fit(X_tr, y_tr)

# A BORDERLINE refused applicant - the case where a counterfactual
# is actually useful. For someone at 0.99 the honest answer is "no
# small change would help", which is also worth being able to say.
scores_all = model.predict_proba(X_te)[:, 1]
refused = X_te.iloc[[int(np.argmin(np.abs(scores_all - 0.62)))]]
baseline = float(model.predict_proba(X_tr)[:, 1].mean())
score = float(model.predict_proba(refused)[0, 1])

print(f'this applicant scored {score:.3f}; the average is {baseline:.3f}')
print('their values:')
for feature in FEATURES:
    print(f'  {feature:<16} {float(refused[feature].iloc[0]):>12,.2f}')
print()

# A counterfactual: what single change would flip the decision? This
# is the explanation a person can act on.
print('counterfactual - the smallest change that moves the score below 0.5:')
found = False
for feature, direction in [('loan', 'down'), ('income', 'up'),
                           ('prior_defaults', 'down'), ('years_employed', 'up')]:
    original = float(refused[feature].iloc[0])
    factors = [0.9, 0.8, 0.7, 0.5, 0.3] if direction == 'down' else [1.2, 1.5, 2.0, 3.0]
    for factor in factors:
        probe = refused.copy()
        probe[feature] = original * factor
        probe['ratio'] = probe['loan'] / probe['income']
        new_score = float(model.predict_proba(probe)[0, 1])
        if new_score < 0.5:
            change = 'reduce' if factor < 1 else 'increase'
            print(f'  {change} {feature} from {original:,.2f} to '
                  f'{original * factor:,.2f}: score {score:.3f} -> {new_score:.3f}')
            found = True
            break
if not found:
    print('  no single-feature change of this size flips it - which is')
    print('  itself the honest answer to give')
print()

# And the applicant at the very top of the risk scale, for contrast.
worst = X_te.iloc[[int(np.argmax(scores_all))]]
print(f'for comparison, the highest-risk applicant scored '
      f'{float(model.predict_proba(worst)[0, 1]):.3f} with a ratio of '
      f'{float(worst["ratio"].iloc[0]):.2f} - no counterfactual of this')
print('kind exists for them, and saying so is better than inventing one.')
print()
print('"Your loan-to-income ratio was too high; borrowing 20% less')
print('would change the answer" is an explanation. "Feature importance')
print('ranks ratio first" is not.')

import pandas as pd

print(pd.DataFrame([
    {'obligation': 'measure disparities before deploying',
     'concretely': 'selection rate, TPR and precision by group, with intervals'},
    {'obligation': 'name the fairness definition you chose',
     'concretely': 'write it down; they conflict, so the choice is a decision'},
    {'obligation': 'record what the model was trained on',
     'concretely': 'period, population, exclusions, label source, version'},
    {'obligation': 'provide a per-decision explanation',
     'concretely': 'the factors, and what would change the answer'},
    {'obligation': 'provide a route to a human',
     'concretely': 'an appeal that a person actually reads'},
    {'obligation': 'monitor the disparities after deployment',
     'concretely': 'the same table, every month, with the same thresholds'},
    {'obligation': 'be able to turn it off',
     'concretely': 'a fallback rule, tested, that production can switch to'},
]).to_string(index=False))

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(seed=101)
n = 12_000

# A hiring screen trained on historical decisions.
group = rng.choice(['A', 'B'], n, p=[.72, .28])
is_b = (group == 'B').astype(int)

# Genuine ability, equally distributed between the groups.
ability = rng.normal(0, 1, n)

# Two features. One is honest; one is a proxy that correlates with
# the group for reasons unrelated to ability.
test_score = ability + rng.normal(0, 0.6, n)
years_at_large_firms = ability * 0.4 - 0.9 * is_b + rng.normal(0, 0.7, n)

# The HISTORICAL label: past hiring decisions, which were biased.
hired_historically = ((1.1 * ability - 0.75 * is_b + rng.normal(0, 0.7, n)) > 0).astype(int)
# The outcome we actually care about, which is NOT biased.
would_succeed = ((1.3 * ability + rng.normal(0, 0.8, n)) > 0).astype(int)

frame = pd.DataFrame({'test_score': test_score,
                      'years_at_large_firms': years_at_large_firms})

print('the data')
print(f'  group B is {is_b.mean():.0%} of applicants')
print(f'  mean ability: A {ability[is_b == 0].mean():+.3f}, '
      f'B {ability[is_b == 1].mean():+.3f}  (equal by construction)')
print(f'  would succeed: A {would_succeed[is_b == 0].mean():.1%}, '
      f'B {would_succeed[is_b == 1].mean():.1%}')
print(f'  hired historically: A {hired_historically[is_b == 0].mean():.1%}, '
      f'B {hired_historically[is_b == 1].mean():.1%}  <- the bias')
print(f'  correlation of years_at_large_firms with group: '
      f'{np.corrcoef(years_at_large_firms, is_b)[0, 1]:.2f}')
print()

def audit(name, features, label):
    X_tr, X_te, y_tr, y_te, g_tr, g_te, s_tr, s_te = train_test_split(
        frame[features], label, is_b, would_succeed, test_size=0.3, random_state=0)
    model = HistGradientBoostingClassifier(random_state=0).fit(X_tr, y_tr)
    selected = model.predict(X_te)

    rows = []
    for g, gname in [(0, 'A'), (1, 'B')]:
        mask = g_te == g
        chosen = selected[mask].astype(bool)
        succeed = s_te[mask].astype(bool)
        rows.append({
            'group': gname, 'n': int(mask.sum()),
            'selection_rate': round(float(chosen.mean()), 3),
            'tpr_vs_success': round(float(chosen[succeed].mean()), 3),
            'precision_vs_success': round(float(succeed[chosen].mean()), 3)
                if chosen.sum() else float('nan'),
        })
    table = pd.DataFrame(rows)
    print(name)
    print(table.to_string(index=False))
    print(f'  selection gap {table["selection_rate"].iloc[0] - table["selection_rate"].iloc[1]:+.3f}   '
          f'opportunity gap {table["tpr_vs_success"].iloc[0] - table["tpr_vs_success"].iloc[1]:+.3f}')
    print()
    return table

# 1. Train on the biased historical label, with the proxy. The worst case.
audit('trained on historical hires, with the proxy feature:',
      ['test_score', 'years_at_large_firms'], hired_historically)

# 2. Drop the proxy. Some of the bias goes.
audit('trained on historical hires, proxy removed:',
      ['test_score'], hired_historically)

# 3. Train on the OUTCOME rather than the past decision. Most of it goes.
audit('trained on who actually succeeded, with the proxy:',
      ['test_score', 'years_at_large_firms'], would_succeed)

# 4. Both fixes together.
audit('trained on who actually succeeded, proxy removed:',
      ['test_score'], would_succeed)

print('Both fixes work, and in this data the proxy is the larger')
print('lever - because years_at_large_firms is correlated -0.44 with')
print('the group and carries almost the whole disparity on its own.')
print()
print('Note what the biased LABEL does even with the proxy removed:')
print('the model trained on historical hires selects fewer people')
print('overall than the one trained on outcomes, because it learned')
print('a decision rule that was stricter than the evidence warranted.')
print()
print('Which lever is larger depends on the data. What does not vary')
print('is that you have to measure both: auditing the feature list')
print('while training on a biased label leaves the bias in place.')

import pandas as pd

print(pd.DataFrame([
    {'factor': 'cost of a wrong decision to the person',
     'automate if': 'low and reversible', 'do not if': 'high or irreversible'},
    {'factor': 'volume',
     'automate if': 'thousands a day', 'do not if': 'a handful a week'},
    {'factor': 'quality of the labels',
     'automate if': 'the outcome is measured', 'do not if': 'the label is a past decision'},
    {'factor': 'explainability required',
     'automate if': 'none, or simple', 'do not if': 'a legal right to reasons'},
    {'factor': 'ability to appeal',
     'automate if': 'a human reviews on request', 'do not if': 'there is no route back'},
]).to_string(index=False))