Responsible AI: Fairness, Bias, and Interpretability Questions

Direct answer

Vanilla gradient saliency measures the LOCAL slope of the target class score at the exact input point (one backward pass); Integrated Gradients (IG) instead accumulates the gradient along the whole straight-line path from a baseline to the input, then multiplies by the total input-minus-baseline difference. Saliency is cheap but fragile: it can saturate to near-zero on confident predictions and jump sharply from an imperceptible input change if the decision surface has a nearby kink. IG costs steps times as many backward passes but is far smoother under small perturbations and satisfies a completeness guarantee saliency does not. In production, choose IG whenever confident, saturated predictions are common or whenever you need attributions stable enough to compare across similar images, and accept plain gradient saliency only when its steps-times cost is the deciding constraint (very high request volume, real-time latency budgets) and some fragility is tolerable.

Structured elaboration

Algorithmic idea, side by side.

• Vanilla gradient saliency: $S_i = \partial f / \partial x_i$ (or $x_i \cdot \partial f/\partial x_i$ for a gradient-times-input variant), a single backward pass evaluated at exactly the observed input.
• Integrated Gradients: $IG_i(x) = (x_i - x_i') \times \int_{\alpha=0}^{1} \frac{\partial F(x' + \alpha(x-x'))}{\partial x_i} \, d\alpha$, where $x'$ is a baseline (commonly all-zero or all-black). In practice the integral is approximated by a Riemann or trapezoidal sum over a fixed number of interpolation steps between baseline and input, each requiring its own forward-and-backward pass.

Failure modes.

• Saliency: gradient saturation (a confidently correct prediction can sit in a locally flat region of the score function, so the gradient collapses toward zero even though the input pixels are entirely responsible for that confidence); and generically, being a purely LOCAL measurement, it can miss the pixel's actual overall contribution when the function is highly nonlinear between the baseline "absence of signal" state and the actual input.
• IG: choice-of-baseline sensitivity (a black-image baseline versus a blurred-image baseline versus a random-noise baseline can all give materially different attributions for the same input, since IG's answer is explicitly relative to "the path from baseline to input"); and if steps is too low, the discretized integral can meaningfully deviate from the true continuous integral (an approximation-error failure mode saliency does not have, since saliency is exact by definition, just exact for the wrong, purely local, question).

Sensitivity to noise. Saliency's local-slope nature makes it exactly as smooth or as jumpy as the decision surface is AT THAT ONE POINT: if a tiny input perturbation crosses a nearby kink in a piecewise-linear network (a ReLU threshold flipping), the gradient can change sharply even though the prediction itself barely moves. IG's path-integral average over the WHOLE path from baseline to input is far less sensitive to any single kink near the evaluation point, because one additional or missing kink along a long path is a small fraction of the path's average gradient.

Practical validation approaches for both. Insertion/deletion curves (rank pixels by the map's magnitude, then progressively insert/delete them and track the class probability) apply unchanged to either method. Model-randomization sanity checks (recompute the map after resetting the model to random weights and check the map changes substantially) also apply to both, and are specifically important for IG because IG's completeness axiom guarantees the attributions SUM to the right total regardless of whether the model has learned anything meaningful, so completeness alone is not evidence of a faithful explanation and the randomization check is the check that actually probes that.

When to select IG over plain gradient saliency in production. Prefer IG when: predictions are frequently confident/saturated (common in well-trained classifiers, exactly where plain saliency is least reliable); attributions need to be compared ACROSS similar inputs (IG's smoother, completeness-anchored attributions are more comparable image-to-image than saliency's locally jumpy ones); or the explanation feeds a downstream decision (a compliance report, a user-facing "why was this flagged" feature) where a single misleading near-zero saliency map on a confidently correct prediction would be actively wrong, not just noisy. Prefer plain gradient saliency when the volume of explanation requests is very high and latency-sensitive (real-time, per-request explanation on every inference), since IG's steps extra backward passes multiply cost linearly, and a rough, occasionally-jumpy map for a debugging or exploratory use case is an acceptable trade for a fraction of the cost.

Worked example

A single scalar-input ReLU network with 200 densely packed kinks in $[0,1]$ (a deliberately non-smooth decision surface, standing in for a realistic CNN's genuinely piecewise-linear behavior), comparing how much each method's attribution changes for a tiny input perturbation ($\epsilon = 0.005$) at 40 random evaluation points:

import numpy as np
import torch
import torch.nn as nn

torch.manual_seed(4)
rng = np.random.default_rng(4)

def vanilla_gradient(model, x, target=0):
    x = x.clone().detach().requires_grad_(True)
    out = model(x); out[0, target].backward()
    return x.grad[0].item()

def integrated_gradients_1d(model, x, baseline, target=0, steps=200):
    alphas = torch.linspace(0, 1, steps + 1)
    diff = x - baseline
    total = 0.0
    for a in alphas:
        pt = (baseline + a * diff).clone().detach().requires_grad_(True)
        model(pt)[0, target].backward()
        total += pt.grad[0].item()
    return diff[0].item() * (total / (steps + 1))

n_units = 200
thresholds = torch.linspace(0.02, 0.98, n_units)
net = nn.Sequential(nn.Linear(1, n_units), nn.ReLU(), nn.Linear(n_units, 1))
with torch.no_grad():
    net[0].weight.copy_(torch.full((n_units,), 80.0).view(-1, 1))
    net[0].bias.copy_(-torch.full((n_units,), 80.0) * thresholds)
    net[2].weight.copy_(torch.sin(torch.linspace(0, 40, n_units)).view(1, -1) * 0.3)

baseline, epsilon = torch.tensor([[0.0]]), 0.005
grad_diffs, ig_diffs = [], []
for _ in range(40):
    xc = float(rng.uniform(0.1, 0.9))
    x_center, x_pert = torch.tensor([[xc]]), torch.tensor([[xc + epsilon]])
    grad_diffs.append(abs(vanilla_gradient(net, x_pert) - vanilla_gradient(net, x_center)))
    ig_diffs.append(abs(integrated_gradients_1d(net, x_pert, baseline) - integrated_gradients_1d(net, x_center, baseline)))

grad_diffs, ig_diffs = np.array(grad_diffs), np.array(ig_diffs)
print("vanilla gradient: mean|change|=", round(grad_diffs.mean(), 4), " max=", round(grad_diffs.max(), 4))
print("Integrated Gradients: mean|change|=", round(ig_diffs.mean(), 4), " max=", round(ig_diffs.max(), 4))
print("ratio (vanilla/IG):", round(grad_diffs.mean() / ig_diffs.mean(), 2), "x")
print("frac > 1.0 (vanilla, IG):", (grad_diffs > 1.0).mean(), (ig_diffs > 1.0).mean())

Executed output (literal stdout of the code block above, run as one fresh module in a clean interpreter):

vanilla gradient: mean|change|= 16.7895  max= 41.676
Integrated Gradients: mean|change|= 0.6583  max= 1.2563
ratio (vanilla/IG): 25.5 x
frac > 1.0 (vanilla, IG): 0.95 0.275

For the exact same 40 evaluation points and the exact same tiny perturbation (0.005), vanilla gradient's attribution changes by a mean of 16.7895 (and as much as 41.676 at its worst), while Integrated Gradients changes by a mean of only 0.6583 (max 1.2563): a 25.5x difference in noise sensitivity for the identical input nudge. The printed "frac > 1.0" line shows 95% of vanilla-gradient trials had a change greater than 1.0 in this run, versus 27.5% for IG, confirming the qualitative claim (saliency is exactly as jumpy as the nearest kink; IG's path average smooths that out) with an actual measured ratio rather than an assertion.

Trade-offs and pitfalls

The most common wrong turn is treating IG as a strict, always-better upgrade and never checking its baseline sensitivity: an all-zero baseline is only a neutral "absence of signal" choice for domains where zero genuinely means absence (raw pixel intensity), and using it uncritically on standardized tabular features (where zero is the training-set mean after standardization, not an "absent" value) produces attributions relative to an arbitrary reference point dressed up as a principled one. A second pitfall is under-provisioning steps: too few interpolation points leaves a visible gap between the attributions' sum and the actual completeness target ($F(x) - F(\text{baseline})$), and that gap should be checked and reported, not assumed away, whenever IG is deployed. Finally, neither method is a causal explanation of the model's decision process in a deeper mechanistic sense; both describe how the SCORE responds to small (saliency) or path-integrated (IG) changes to the input, which is a genuinely useful and different question from "what feature did the model's internal computation actually use," and conflating the two in a user-facing explanation risks overclaiming what either method actually establishes.

Direct answer

Define faithfulness as: rank the input's features by the local explanation's attribution magnitude, progressively replace the top-$k$ features with baseline values, and measure how much the model's output moves away from its original prediction toward the fully-baseline prediction as $k$ grows. A faithful explanation should produce a steep drop early (the features it called important really do move the output), while an unfaithful one produces a flat curve no better than a random feature order. Averaging that drop across all prefix lengths $k=0,\dots,p$ gives a single number, the normalized Area Over the Perturbation Curve (AOPC), and dividing by the total gap between the real prediction and the fully-baseline prediction is what makes the score comparable across inputs that start at very different confidence levels.

Structured elaboration

The metric, precisely. Let $f$ be the model's scalar output of interest (a predicted probability, a logit, a regression value), $x$ the input, and $\phi_1, \dots, \phi_p$ the local explanation's per-feature attribution scores for $x$ (from SHAP, LIME, or any other local method). Sort features into a ranking $\pi$ by $|\phi_i|$ descending, and let $x^{(k)}$ be $x$ with the top-$k$ ranked features replaced by baseline values (and the rest left at their original values in $x$). Define the per-$k$ drop:

$D_k(x) = f(x) - f\big(x^{(k)}\big), \qquad D_0(x) = 0$

Aggregate over the whole ranking into a single scalar, the raw AOPC:

$\text{AOPC}(x) = \frac{1}{p+1}\sum_{k=0}^{p} D_k(x)$

and normalize by the maximum achievable gap, the difference between the real prediction and the prediction when every feature is replaced by the baseline:

$\text{Faithfulness}(x) = \frac{\text{AOPC}(x)}{f(x) - f(x^{(p)})}$

A value near 1 means removing the explanation's top features accounts for essentially the whole gap between "the model's real answer" and "the model with no case-specific information," which is exactly what a faithful explanation should do. A value near 0 means the top-ranked features barely moved the output at all, meaning the explanation pointed at the wrong things. Values can go negative if removing the "important" features somehow pushes the output further from the baseline (a real, informative failure mode: it flags a non-monotonic interaction the explanation missed).

Algorithm.

1. Compute the local attribution $\phi$ for $x$ using whatever explanation method is being evaluated.
2. Rank features by $|\phi_i|$.
3. For $k = 0, 1, \dots, p$: build $x^{(k)}$ by substituting the top-$k$ ranked features with baseline values, call the model, record $f(x^{(k)})$.
4. Average the drops to get raw AOPC, divide by the full-perturbation gap to normalize.
5. Repeat for every audited input and report the mean (or median, for robustness to outlier cases where the denominator is near zero) normalized faithfulness across the dataset.

Complexity analysis. Let $p$ be the number of features, $N$ the number of audited inputs, and $M$ the number of background rows used if the baseline is a sampled expectation rather than a single fixed vector.

• Fixed baseline vector (a single mean-imputed or zero vector): computing the full curve for one input costs $p$ model forward passes (one per $k$ from 1 to $p$; $k=0$ is free, it is just $f(x)$ already on hand). Across a dataset, that is $O(N \cdot p)$ model calls.
• Sampled-background (expectation) baseline: for each $k$, the baseline value used for the perturbed features is averaged over $M$ background rows rather than a single vector, so each $k$ costs $M$ model calls (or one batched call of size $M$). That is $O(N \cdot p \cdot M)$ model calls total.
• This cost is on top of, and usually much cheaper than, generating the attributions themselves: TreeSHAP for a tree ensemble costs $O(T L D^2)$ per instance ($T$ trees, $L$ max leaves, $D$ max depth), so for tree models the faithfulness evaluation is typically the cheap half of the pipeline; for a model-agnostic KernelSHAP or LIME explanation, generating $\phi$ itself can dominate, since those already cost many perturbation-and-refit calls per instance before faithfulness evaluation even starts.
• If the model itself is expensive to call (a large neural network), the dominant cost is model inference calls, not the bookkeeping around them, so batching all $k$ (and all $M$ background rows, if used) into one forward pass over a matrix of perturbed rows, rather than $p$ or $p \cdot M$ separate calls, is the practical optimization that matters most.

Normalization and baseline choices, and why both are load-bearing for comparability across inputs.

• Why normalize at all. Two inputs can have very different "room to move": one input the model is already 99% confident about has almost no gap left to explain, so even a perfect explanation produces a small raw drop; an input near the decision boundary has a much larger gap. Comparing RAW drops across such inputs, or across a whole dataset, conflates "how confident was the model to begin with" with "how faithful was the explanation." Dividing by $f(x) - f(x^{(p)})$, the total achievable gap for that specific input, puts every input on the same roughly $[0,1]$ scale regardless of its starting confidence.
• Baseline choice changes the number, not just its scale. A fixed mean-imputation vector is cheap and deterministic but can be off the true data manifold if features are correlated (the "average customer" may not resemble any real customer). A sampled-background (interventional) expectation baseline, averaging the model's output over real background rows substituted in for the removed features, better respects the marginal distribution of each feature and is the choice most consistent with how interventional SHAP itself defines its value function, at the cost of $M\times$ more model calls. A conditional baseline, sampling replacement values from the distribution of the removed features CONDITIONED on the retained features, is the most faithful to the data manifold but requires a model of that conditional distribution and is rarely used outside of research settings because of that extra machinery. An off-manifold constant baseline (for example, a fixed extreme value chosen without reference to the training distribution) is the cheapest option and the most dangerous: since the model was never trained near it, its behavior there is unconstrained extrapolation, and the resulting "drop" reflects model behavior in a region the model was never trained to be sensible in, not a genuine test of the explanation. Any reported faithfulness number is only comparable to another faithfulness number computed with the SAME baseline policy; comparing a paper's off-manifold-baseline faithfulness score to another's sampled-background score is comparing two different metrics that happen to share a name.

Worked example

A synthetic classification problem with two truly predictive features ($x_1, x_2$) and four pure-noise features, so the ground truth of "what the model should attribute weight to" is known, lets the faithfulness metric be sanity-checked against a real explanation (TreeSHAP) and a deliberately bad one (a random feature ranking used as a negative control), plus three baseline choices for the same explanation:

import numpy as np
from sklearn.ensemble import RandomForestClassifier
import shap

rng = np.random.default_rng(11)
n = 3000
x1 = rng.normal(0, 1, n)
x2 = rng.normal(0, 1, n)
x3 = rng.normal(0, 1, n)  # noise
x4 = rng.normal(0, 1, n)  # noise
x5 = rng.normal(0, 1, n)  # noise
x6 = rng.normal(0, 1, n)  # noise
logit = 3.0 * x1 - 2.0 * x2
p_true = 1 / (1 + np.exp(-logit))
y = (rng.uniform(0, 1, n) < p_true).astype(int)
X = np.column_stack([x1, x2, x3, x4, x5, x6])
feature_names = ["x1", "x2", "x3", "x4", "x5", "x6"]

X_train, y_train = X[:2500], y[:2500]
X_bg, y_bg = X[2500:2600], y[2500:2600]      # background/reference set for baselines (100 rows, avoids SHAP auto-subsampling)
X_test, y_test = X[2800:], y[2800:]

model = RandomForestClassifier(n_estimators=200, max_depth=4, random_state=0)
model.fit(X_train, y_train)

def predict_proba1(Xrows):
    # Model output used by the faithfulness metric: P(y=1).
    return model.predict_proba(Xrows)[:, 1]

# ---- 1. Local explanation via TreeSHAP ----
explainer = shap.TreeExplainer(model, data=X_bg, feature_perturbation="interventional")
x0 = X_test[0]
shap_values = explainer.shap_values(x0.reshape(1, -1))[0][:, 1]  # class-1 (P(y=1)) column
print("Instance x0:", dict(zip(feature_names, np.round(x0, 3))))
print("SHAP values:", dict(zip(feature_names, np.round(shap_values, 4))))

ranking_shap = np.argsort(-np.abs(shap_values))
print("SHAP-based importance ranking (most attributed first):", [feature_names[i] for i in ranking_shap])

# ---- 2. Faithfulness metric: normalized AOPC under top-k deletion ----
def faithfulness_curve(x, ranking, baseline_vec, predict_fn, p_features):
    # Return the drop-from-original curve f(x) - f(x with top-k perturbed), k=0..p_features,
    # using a single fixed baseline VECTOR (mean-imputation baseline).
    f_x = predict_fn(x.reshape(1, -1))[0]
    drops = [0.0]  # k=0: nothing perturbed, drop is 0 by definition
    x_pert = x.copy()
    for k in range(1, p_features + 1):
        idx = ranking[:k]
        x_pert = x.copy()
        x_pert[idx] = baseline_vec[idx]
        f_k = predict_fn(x_pert.reshape(1, -1))[0]
        drops.append(f_x - f_k)
    return np.array(drops), f_x

def faithfulness_curve_sampled_baseline(x, ranking, background, predict_fn, p_features, seed=0):
    # Same curve, but the baseline for each perturbed feature set is the AVERAGE model
    # output over M background rows substituted in for the top-k features (interventional
    # expectation baseline), instead of a single fixed mean vector.
    f_x = predict_fn(x.reshape(1, -1))[0]
    drops = [0.0]
    M = background.shape[0]
    for k in range(1, p_features + 1):
        idx = ranking[:k]
        Xp = np.tile(x, (M, 1))
        Xp[:, idx] = background[:, idx]
        f_k = predict_fn(Xp).mean()
        drops.append(f_x - f_k)
    return np.array(drops), f_x

p_features = X.shape[1]
mean_baseline = X_bg.mean(axis=0)
f_full_baseline_meanvec = predict_proba1(mean_baseline.reshape(1, -1))[0]

drops_shap_mean, f_x0 = faithfulness_curve(x0, ranking_shap, mean_baseline, predict_proba1, p_features)
aopc_shap_mean_raw = drops_shap_mean[1:].mean()
norm_denominator = f_x0 - f_full_baseline_meanvec
aopc_shap_mean_norm = aopc_shap_mean_raw / norm_denominator

print("\n--- Mean-vector baseline ---")
print("f(x0) =", round(f_x0, 4), " f(full mean-baseline) =", round(f_full_baseline_meanvec, 4))
print("Drop curve k=0..6 (SHAP ranking):", np.round(drops_shap_mean, 4))
print("Raw AOPC (SHAP ranking):", round(aopc_shap_mean_raw, 4))
print("Normalized AOPC (SHAP ranking):", round(aopc_shap_mean_norm, 4))

# Negative control: a RANDOM feature ranking should be less faithful (smaller AOPC)
rng2 = np.random.default_rng(99)
ranking_random = rng2.permutation(p_features)
drops_random_mean, _ = faithfulness_curve(x0, ranking_random, mean_baseline, predict_proba1, p_features)
aopc_random_raw = drops_random_mean[1:].mean()
aopc_random_norm = aopc_random_raw / norm_denominator
print("\nRandom ranking used as negative control:", [feature_names[i] for i in ranking_random])
print("Drop curve k=0..6 (random ranking):", np.round(drops_random_mean, 4))
print("Normalized AOPC (random ranking):", round(aopc_random_norm, 4))

# ---- 3. Baseline-choice sensitivity: sampled-background baseline instead of mean vector ----
drops_shap_sampled, _ = faithfulness_curve_sampled_baseline(x0, ranking_shap, X_bg, predict_proba1, p_features)
aopc_shap_sampled_raw = drops_shap_sampled[1:].mean()
f_full_baseline_sampled = predict_proba1(X_bg).mean()
norm_denom_sampled = f_x0 - f_full_baseline_sampled
aopc_shap_sampled_norm = aopc_shap_sampled_raw / norm_denom_sampled

print("\n--- Sampled-background baseline (M =", X_bg.shape[0], "rows) ---")
print("Drop curve k=0..6 (SHAP ranking, sampled baseline):", np.round(drops_shap_sampled, 4))
print("Normalized AOPC (SHAP ranking, sampled baseline):", round(aopc_shap_sampled_norm, 4))
print("\nBaseline-choice gap: normalized AOPC differs by",
      round(abs(aopc_shap_mean_norm - aopc_shap_sampled_norm), 4),
      "between the mean-vector baseline and the sampled-background baseline for the SAME explanation.")

# ---- 4. Off-manifold baseline: an extreme out-of-distribution constant, to show
# why the baseline must be documented for cross-method / cross-paper comparability ----
offmanifold_baseline = np.full(p_features, 4.0)   # features are ~N(0,1); +4 is far off-manifold
drops_shap_offmanifold, _ = faithfulness_curve(x0, ranking_shap, offmanifold_baseline, predict_proba1, p_features)
f_full_offmanifold = predict_proba1(offmanifold_baseline.reshape(1, -1))[0]
norm_denom_offmanifold = f_x0 - f_full_offmanifold
aopc_shap_offmanifold_raw = drops_shap_offmanifold[1:].mean()
aopc_shap_offmanifold_norm = aopc_shap_offmanifold_raw / norm_denom_offmanifold

print("\n--- Off-manifold constant baseline (all features = +4.0, vs training range ~[-3,3]) ---")
print("f(full off-manifold baseline) =", round(f_full_offmanifold, 4))
print("Drop curve k=0..6 (SHAP ranking, off-manifold baseline):", np.round(drops_shap_offmanifold, 4))
print("Normalized AOPC (SHAP ranking, off-manifold baseline):", round(aopc_shap_offmanifold_norm, 4))
print("\nSame explanation, same ranking, three baselines -> normalized AOPC of",
      round(aopc_shap_mean_norm, 4), "(mean vector),", round(aopc_shap_sampled_norm, 4),
      "(sampled background), and", round(aopc_shap_offmanifold_norm, 4), "(off-manifold constant):",
      "the off-manifold choice moves the score the most, because the model was never trained on inputs near it.")

Executed output:

Instance x0: {'x1': np.float64(-0.187), 'x2': np.float64(1.31), 'x3': np.float64(0.923), 'x4': np.float64(0.194), 'x5': np.float64(-0.556), 'x6': np.float64(0.462)}
SHAP values: {'x1': np.float64(-0.1192), 'x2': np.float64(-0.1524), 'x3': np.float64(-0.0003), 'x4': np.float64(0.0061), 'x5': np.float64(0.0019), 'x6': np.float64(0.0002)}
SHAP-based importance ranking (most attributed first): ['x2', 'x1', 'x4', 'x5', 'x3', 'x6']

--- Mean-vector baseline ---
f(x0) = 0.2056  f(full mean-baseline) = 0.3686
Drop curve k=0..6 (SHAP ranking): [ 0.     -0.0494 -0.1729 -0.1716 -0.1713 -0.1694 -0.163 ]
Raw AOPC (SHAP ranking): -0.1496
Normalized AOPC (SHAP ranking): 0.9179

Random ranking used as negative control: ['x1', 'x4', 'x5', 'x6', 'x3', 'x2']
Drop curve k=0..6 (random ranking): [ 0.     -0.0713 -0.07   -0.069  -0.0637 -0.0602 -0.163 ]
Normalized AOPC (random ranking): 0.5085

--- Sampled-background baseline (M = 100 rows) ---
Drop curve k=0..6 (SHAP ranking, sampled baseline): [ 0.     -0.1579 -0.273  -0.2654 -0.2627 -0.2634 -0.2637]
Normalized AOPC (SHAP ranking, sampled baseline): 0.9393

Baseline-choice gap: normalized AOPC differs by 0.0215 between the mean-vector baseline and the sampled-background baseline for the SAME explanation.

--- Off-manifold constant baseline (all features = +4.0, vs training range ~[-3,3]) ---
f(full off-manifold baseline) = 0.5031
Drop curve k=0..6 (SHAP ranking, off-manifold baseline): [ 0.      0.0099 -0.3331 -0.3601 -0.3182 -0.2944 -0.2975]
Normalized AOPC (SHAP ranking, off-manifold baseline): 0.8928

Same explanation, same ranking, three baselines -> normalized AOPC of 0.9179 (mean vector), 0.9393 (sampled background), and 0.8928 (off-manifold constant): the off-manifold choice moves the score the most, because the model was never trained on inputs near it.

Two things are worth reading closely in this output. First, the sanity check works: SHAP's ranking (which correctly puts $x_1$ and $x_2$, the two truly predictive features, first) reaches a normalized faithfulness of 0.9179, while a random ranking on the exact same instance and baseline reaches only 0.5085, confirming the metric actually distinguishes a good explanation from a bad one rather than always reading near 1. Second, the baseline sensitivity is real even here: the identical SHAP ranking scores 0.9179 with a mean-vector baseline, 0.9393 with a sampled-background baseline, and 0.8928 with an off-manifold constant, and the off-manifold curve even shows a small POSITIVE drop at $k=1$ (+0.0099), meaning removing just the top feature into an unseen region briefly made the prediction move the wrong way, a symptom of the model extrapolating unpredictably rather than the explanation being wrong.

Trade-offs and pitfalls

The single most common mistake is comparing faithfulness scores computed under different baseline policies as if they were the same metric; a paper or internal report claiming "our explanation method scores 0.95 faithfulness" is uninterpretable without stating the baseline, since the example above shows almost a 5-point swing on the identical ranking just from swapping baselines. A second pitfall is skipping the negative-control check: without confirming that a deliberately bad (random) ranking scores meaningfully lower than the real explanation on the same data, there is no evidence the metric is measuring anything beyond how sensitive the model is to perturbation in general, rather than how good the specific ranking is. A third is using an off-manifold baseline for convenience (a global constant is trivial to implement) without acknowledging that the resulting score partly measures the model's un-validated extrapolation behavior rather than the explanation's quality, which is exactly the failure mode this worked example surfaced at $k=1$. Finally, this metric evaluates whether the TOP-$k$ features are influential, which is necessary but not sufficient for faithfulness: an explanation could correctly rank the top few features while badly misattributing the relative weight among them, or while being unstable across reruns (a LIME-style concern); a full faithfulness audit typically pairs this deletion-based AOPC with a stability check (does the same instance get the same ranking across repeated runs) rather than relying on AOPC alone.

Direct answer

For labeled image data, the recurring sources of annotation bias are annotator demographic mismatch with the depicted subjects, ambiguous or culturally loaded labeling instructions, fatigue/anchoring effects that make later images in a batch drift toward whatever label came before, and image-quality confounds (lighting, camera angle, resolution) that correlate with a demographic attribute and get misread as the target label. The practical detection process is to measure inter-annotator agreement overall and then AGAIN split by subgroup, and to run a stratified error analysis comparing each annotator's label rate against a trusted reference (or against the group average) split by subgroup, because an annotator who is well-calibrated in aggregate can still be quietly biased on one subgroup, which an overall agreement number will hide.

Structured elaboration

Typical sources of annotation bias in labeled image datasets.

1. Annotator-subject demographic mismatch. Annotators unfamiliar with a subgroup's visual presentation (skin tone, hairstyle, cultural dress, age-related appearance) make systematically more labeling errors or apply harsher/more permissive judgment for that subgroup, especially on subjective labels (attractiveness, emotion, professionalism) rather than purely objective ones (object category).
2. Ambiguous or under-specified labeling instructions. When the labeling guideline leaves room for judgment (what counts as "provocative attire," what emotion a facial expression conveys), annotators fill the gap with their own cultural priors, and those priors are not uniform across subgroups of the labeled subjects.
3. Fatigue and anchoring/order effects. Long annotation sessions produce label drift: an annotator's threshold shifts over a shift, or a label is unconsciously anchored to the previous few images. If certain subgroups are clustered together in the annotation queue (a common artifact of how images are batched by source or by scraping order), this becomes a subgroup-correlated bias rather than pure noise.
4. Image-quality confounds correlated with the sensitive attribute. Some subgroups are systematically photographed under different lighting, camera equipment, or resolution in the source dataset (older public datasets undercount certain skin tones under standard exposure settings), and annotators (or, more insidiously, the downstream model) can pick up on those quality artifacts as if they were signal for the label.
5. Historical or embedded stereotype bias. For subjective categorical labels (occupation guesses from an image, perceived competence), annotators bring pre-existing stereotypes into ambiguous cases, producing a systematic skew that is not random noise and does not average out with more annotators unless the annotator pool itself is diverse.

Practical detection process.

1. Assign redundant annotation. Every image gets labeled by multiple (at least 2, ideally 3+) independent annotators, and every annotator's assignment is stratified so no annotator only ever sees one subgroup (which would make it impossible to separate an annotator's per-subgroup behavior from a confound in which images they happened to get).
2. Compute agreement metrics overall. Use Cohen's kappa for two annotators or Fleiss' kappa for three or more, which correct for chance agreement (unlike raw percent agreement, which looks artificially high whenever one label is common).
3. Compute the SAME agreement metric split by subgroup. This is the step that actually surfaces annotation bias: an annotator whose overall kappa looks fine can have a much lower kappa specifically within one subgroup, revealing that their disagreements are not spread evenly but concentrated exactly where a demographic-mismatch or stereotype effect would predict.
4. Stratified error analysis against a trusted reference. Where a small, carefully-adjudicated gold-standard subset exists (ideally labeled by a demographically diverse expert panel), compare each annotator's positive-label rate against that reference, split by subgroup, using a two-proportion test so a rate difference is judged against its own sampling noise rather than eyeballed.
5. Look for order/fatigue confounds. Check whether an annotator's label rate for a subgroup correlates with that subgroup's position in the annotation queue or session length; if a subgroup is disproportionately annotated late in long sessions, a fatigue effect will look like a demographic effect unless queue position is controlled for.

Corrective actions.

• Rebalance annotator-to-image assignment so every annotator sees a representative mix of subgroups in every session, breaking the confound between "which annotator" and "which subgroup."
• Tighten the labeling guideline for the specific ambiguous cases the stratified analysis flagged, with example images from each affected subgroup, and re-annotate a sample to confirm the tightened guideline actually closes the gap.
• Recruit a more demographically diverse annotator pool, particularly when the flagged bias tracks a demographic mismatch between annotators and the labeled population.
• Route flagged, high-disagreement, or subgroup-skewed images to an adjudication step (a senior reviewer or a majority vote among a larger annotator set) rather than trusting a single annotator's label for exactly the cases the audit shows are least reliable.
• Down-weight or exclude a systematically biased annotator's affected labels during model training once the pattern is confirmed statistically, rather than silently keeping labels known to carry a directional error.

Worked example

I built a synthetic labeling scenario to make the detection step concrete rather than asserted: 800 images with a true label independent of a binary subgroup attribute, labeled by four simulated annotators. Three annotators (ann_1, ann_2, ann_4) label with ordinary random noise (5-7% flip rate) uncorrelated with subgroup; one annotator (ann_3_biased) has an extra 22% chance of flipping a true-negative image to positive specifically when the image is in subgroup 1, simulating a real stereotype-driven annotation bias.

import numpy as np
from sklearn.metrics import cohen_kappa_score
from scipy.stats import norm

rng = np.random.default_rng(17)

n_images = 800
subgroup = rng.binomial(1, 0.5, size=n_images)
true_label = rng.binomial(1, 0.4, size=n_images)

def annotate(true_label, subgroup, noise_rate, biased_extra_flip_for_subgroup1):
    labels = true_label.copy()
    noise_flip = rng.binomial(1, noise_rate, size=len(true_label)).astype(bool)
    labels[noise_flip] = 1 - labels[noise_flip]
    if biased_extra_flip_for_subgroup1 > 0:
        target = (subgroup == 1) & (true_label == 0)
        extra_flip = rng.binomial(1, biased_extra_flip_for_subgroup1, size=len(true_label)).astype(bool)
        labels[target & extra_flip] = 1
    return labels

annotators = {
    "ann_1": annotate(true_label, subgroup, 0.05, 0.0),
    "ann_2": annotate(true_label, subgroup, 0.06, 0.0),
    "ann_3_biased": annotate(true_label, subgroup, 0.05, 0.22),
    "ann_4": annotate(true_label, subgroup, 0.07, 0.0),
}

def two_proportion_z(p1, n1, p2, n2):
    p_pool = (p1 * n1 + p2 * n2) / (n1 + n2)
    se = np.sqrt(p_pool * (1 - p_pool) * (1 / n1 + 1 / n2))
    z = (p1 - p2) / se
    return z, 2 * (1 - norm.cdf(abs(z)))

for name, labels in annotators.items():
    p0, p1 = labels[subgroup == 0].mean(), labels[subgroup == 1].mean()
    n0, n1 = int((subgroup == 0).sum()), int((subgroup == 1).sum())
    z, p_value = two_proportion_z(p1, n1, p0, n0)
    print(f"{name}: rate|sub=0={p0:.3f} rate|sub=1={p1:.3f} z={z:.2f} p={p_value:.5f}")

names = list(annotators.keys())
for i in range(len(names)):
    for j in range(i + 1, len(names)):
        k_all = cohen_kappa_score(annotators[names[i]], annotators[names[j]])
        m1 = subgroup == 1
        k_sub1 = cohen_kappa_score(annotators[names[i]][m1], annotators[names[j]][m1])
        print(f"kappa({names[i]},{names[j]}) overall={k_all:.3f} subgroup1_only={k_sub1:.3f}")

Running this prints (seed 17, numpy 2.5.1):

ann_1: rate|sub=0=0.409 rate|sub=1=0.380 z=-0.84 p=0.40027
ann_2: rate|sub=0=0.412 rate|sub=1=0.378 z=-0.99 p=0.32423
ann_3_biased: rate|sub=0=0.419 rate|sub=1=0.511 z=2.61 p=0.00911
ann_4: rate|sub=0=0.422 rate|sub=1=0.380 z=-1.20 p=0.23136

kappa(ann_1,ann_2) overall=0.759 subgroup1_only=0.770
kappa(ann_1,ann_3_biased) overall=0.701 subgroup1_only=0.579
kappa(ann_1,ann_4) overall=0.752 subgroup1_only=0.743
kappa(ann_2,ann_3_biased) overall=0.701 subgroup1_only=0.604
kappa(ann_2,ann_4) overall=0.726 subgroup1_only=0.727
kappa(ann_3_biased,ann_4) overall=0.683 subgroup1_only=0.569

The two-proportion test on ann_3_biased's own label rate correctly flags it (p = 0.009, below a 0.01 threshold) since its positive-label rate jumps from 41.9% to 51.1% specifically in subgroup 1, while the other three annotators show no significant subgroup gap. The kappa breakdown adds independent confirmation: every pair INVOLVING ann_3_biased loses noticeably more agreement within subgroup 1 specifically (for example ann_1 vs ann_3_biased drops from kappa 0.701 overall to 0.579 within subgroup 1) than pairs excluding ann_3_biased (ann_1 vs ann_4 barely moves, 0.752 to 0.743). Both signals independently point at the same annotator and the same subgroup, which is exactly the corroborating evidence a real audit needs before accusing a specific annotator of biased labeling.

Trade-offs and pitfalls

• Checking overall agreement alone is the single most common miss: ann_3_biased's overall kappa with the other annotators (0.68-0.70) looks unremarkable next to the other pairs (0.73-0.76), and only the subgroup-stratified breakdown exposes the concentrated disagreement.
• A stratified analysis needs enough per-subgroup, per-annotator sample size to be trustworthy; with a small subgroup or a small per-annotator image count, both the two-proportion test and the subgroup-restricted kappa become noisy, and a real audit should report a confidence interval, not just a point estimate. Even in the worked example above, where subgroup 1 has a comfortable 397 images, the 95% confidence interval around ann_3_biased's 9.2-percentage-point gap (computed the same Wald normal-approximation way as the two-proportion z-test above: gap +/- 1.96 * standard error) works out to roughly [2.3, 16.1] points, not a tight band around 9.2; with a genuinely small subgroup, for example 30-40 images instead of 397, that same style of calculation produces an interval wide enough to include zero, meaning the audit could no longer distinguish the biased annotator from ordinary noise at all. A reviewer acting on the bare point estimate, without that interval, is treating a number that is inherently noisy as more precise than it actually is.
• Do not treat a flagged annotator as proof of malicious bias; the same statistical signature can come from a genuinely ambiguous subset of images that happens to correlate with the subgroup (an image-quality confound), so the corrective step should include reviewing a sample of the flagged annotator's specific disagreements before acting, not just the aggregate numbers.
• Rebalancing annotator-to-image assignment after the fact does not retroactively fix already-collected biased labels; a real remediation plan needs both a forward-looking process fix and a backward-looking re-annotation or down-weighting pass on the affected historical data.

Direct answer

With a 7-day label delay, any fairness metric computed "right now" from only the labels that happen to have already resolved is exposed to censoring bias: if HOW FAST a label resolves is correlated with the outcome itself (a common and realistic pattern, an expedited-review process resolves the model's caught frauds fast while a missed fraud only surfaces later through a customer chargeback), the currently-observed labels are a systematically skewed sample of the truth. The fix is to treat time-to-label-resolution as a survival problem: estimate the resolution-time distribution with a Kaplan-Meier curve, stratified by whatever is known at PREDICTION time (never by the future label itself), and reweight every currently-observed label by the inverse of its estimated probability of having resolved by now (inverse-probability-of-censoring weighting, IPCW). This produces a streaming estimate that is both lower-latency than waiting the full 7 days and, done correctly, less biased than a naive same-day computation.

Structured elaboration

Why naive streaming metrics are biased. At any moment, transactions younger than 7 days have a mix of resolved and not-yet-resolved labels. If resolution speed is INDEPENDENT of the outcome, restricting a metric to currently-resolved labels only adds noise (a smaller effective sample), which is manageable. But if resolution speed correlates with the outcome, for example because a transaction the model correctly flagged gets fast-tracked through an expedited investigation while a transaction the model missed only gets confirmed later via a customer dispute, then at any given moment the resolved sample OVER-represents "easy," fast-resolving cases and UNDER-represents the harder, slow-resolving ones. A naive TPR computed from today's resolved labels alone is then biased in a specific, predictable direction, not just noisier.

Kaplan-Meier for the resolution-time distribution. Fit a survival curve on HISTORICAL, fully-resolved data (transactions old enough, at least 7 days, that every one of them is guaranteed resolved, so there is no censoring in this particular fitting set) to estimate $F(t) = P(\text{label resolved by day } t)$, stratified by every variable known at PREDICTION time that plausibly affects resolution speed. Critically, the model's own real-time decision (flagged or not flagged) is exactly this kind of variable: it is known immediately, for every transaction, whether or not its label has resolved yet, so it is a valid stratification variable, unlike the (not-yet-known) true label itself.

Inverse-probability-of-censoring weighting. For a currently-observed (resolved) label at age $a$ days, weight it by $w = 1 / \hat{F}_{\text{stratum}}(a)$, the inverse of the estimated probability that a label in this stratum would have resolved by age $a$. A label that resolved unusually fast for its stratum (small $a$, so $\hat{F}(a)$ is small) gets a large weight, because it is standing in, on average, for the many same-stratum cases that have NOT yet resolved at that age; a label observed only once most of its stratum has already resolved ($\hat{F}(a)$ close to 1) gets a weight close to 1, since barely any of its peers were still unresolved at that point.

Pseudocode:

# offline, refreshed periodically on fully-resolved historical data:
for each stratum s in (group, predicted_decision):
    fit KaplanMeier(delays of historical transactions in stratum s) -> F_hat[s]

# streaming, run continuously as "today" advances:
for each transaction i with age a_i <= 7 and label observed (delay_i <= a_i):
    s_i = stratum(group_i, predicted_decision_i)
    weight_i = 1 / max(F_hat[s_i](min(a_i, 7)), epsilon)

metric_g = weighted_average( predicted_decision_i == true_label_i
                              for i in observed transactions of group g,
                              weights = weight_i )

Variance-versus-bias trade-off. Waiting the full 7 days before computing anything (a batch process) is unbiased by construction, since there is no censoring left to correct once every included transaction is guaranteed resolved, but it has maximum latency and, if you only look at a narrow recent window, a smaller effective sample (only transactions old enough to have matured). The IPCW streaming estimator uses every currently-resolved label, including very recent ones, which lowers latency and can even increase the effective sample size relative to a narrow mature-only window, but the weights themselves are estimates from a fitted survival curve, so any error in that curve (wrong stratification, an unmodeled dependency between resolution speed and the outcome) leaks directly into the metric. In practice this trade-off is not always "IPCW costs precision for speed": if the naive alternative to full-wait is a much smaller mature-only sample, IPCW's much larger effective sample can win on both bias AND variance simultaneously, provided the stratification used for the survival curve is genuinely correctly specified.

Worked example

A fraud-detection model with a real TPR gap (majority 0.80, minority 0.55), where resolution delay depends on the model's own real-time decision (flagged transactions resolve fast via expedited review; unflagged transactions resolve slowly via customer disputes) plus a group-level operational slowdown:

import numpy as np
import pandas as pd

rng = np.random.default_rng(303)
MAX_DELAY = 7
group_share = {"majority": 0.75, "minority": 0.25}
true_fraud_rate = {"majority": 0.05, "minority": 0.05}
model_tpr = {"majority": 0.80, "minority": 0.55}
delay_probs = {
    ("majority", 1): np.array([0.55, 0.25, 0.10, 0.05, 0.03, 0.01, 0.01]),
    ("majority", 0): np.array([0.08, 0.12, 0.15, 0.20, 0.20, 0.15, 0.10]),
    ("minority", 1): np.array([0.35, 0.28, 0.15, 0.10, 0.06, 0.04, 0.02]),
    ("minority", 0): np.array([0.03, 0.05, 0.08, 0.14, 0.20, 0.25, 0.25]),
}

def simulate_transactions(n, seed):
    r = np.random.default_rng(seed)
    groups = r.choice(list(group_share), size=n, p=list(group_share.values()))
    is_fraud = np.array([r.random() < true_fraud_rate[g] for g in groups])
    predicted_fraud = np.array([(r.random() < model_tpr[g]) if f else (r.random() < 0.03)
                                 for g, f in zip(groups, is_fraud)]).astype(int)
    delay = np.array([r.choice(np.arange(1, MAX_DELAY + 1), p=delay_probs[(g, p)])
                       for g, p in zip(groups, predicted_fraud)])
    return pd.DataFrame({"group": groups, "is_fraud": is_fraud.astype(int),
                          "predicted_fraud": predicted_fraud, "delay": delay})

def kaplan_meier_resolution_cdf(delays, max_t=MAX_DELAY):
    n_at_risk, survival, cdf = len(delays), 1.0, {}
    for t in range(1, max_t + 1):
        events_at_t = int((delays == t).sum())
        if n_at_risk > 0:
            survival *= (1 - events_at_t / n_at_risk)
        cdf[t] = 1 - survival
        n_at_risk -= events_at_t
    return cdf

historical = simulate_transactions(80_000, seed=1)
strata = [(g, p) for g in group_share for p in (0, 1)]
resolution_cdf = {st: kaplan_meier_resolution_cdf(
    historical.loc[(historical.group == st[0]) & (historical.predicted_fraud == st[1]), "delay"].values)
    for st in strata}

cohort = simulate_transactions(300_000, seed=2)
cohort["age"] = rng.integers(0, MAX_DELAY + 1, len(cohort))
cohort["observed"] = cohort["delay"] <= cohort["age"]

def tpr(df):
    fraud = df[df.is_fraud == 1]
    return (fraud.predicted_fraud == 1).mean() if len(fraud) else np.nan

naive_tpr = {g: tpr(cohort[(cohort.group == g) & cohort.observed]) for g in group_share}
mature = cohort[cohort.age >= MAX_DELAY]
fullwait_tpr = {g: tpr(mature[mature.group == g]) for g in group_share}

def ipcw_tpr(df):
    fraud = df[(df.is_fraud == 1) & (df.age >= 1) & df.observed].copy()
    fraud["weight"] = [1 / max(resolution_cdf[(g, p)][min(a, MAX_DELAY)], 1e-6)
                        for g, p, a in zip(fraud.group, fraud.predicted_fraud, fraud.age)]
    return np.average(fraud.predicted_fraud, weights=fraud.weight) if len(fraud) else np.nan

ipcw_tpr_result = {g: ipcw_tpr(cohort[cohort.group == g]) for g in group_share}

naive_gap = naive_tpr["majority"] - naive_tpr["minority"]
fullwait_gap = fullwait_tpr["majority"] - fullwait_tpr["minority"]
ipcw_gap = ipcw_tpr_result["majority"] - ipcw_tpr_result["minority"]
true_gap = model_tpr["majority"] - model_tpr["minority"]
print(f"TPR gap: true={true_gap:.4f}, naive={naive_gap:.4f}, full-wait={fullwait_gap:.4f}, IPCW={ipcw_gap:.4f}")
print(f"absolute error vs true gap: naive={abs(naive_gap - true_gap):.4f}, full-wait={abs(fullwait_gap - true_gap):.4f}, IPCW={abs(ipcw_gap - true_gap):.4f}")

rng_boot = np.random.default_rng(9)
def bootstrap_gap(estimator_fn, pool, n_boot=400):
    n = len(pool)
    return np.std([estimator_fn(pool.iloc[rng_boot.integers(0, n, n)]) for _ in range(n_boot)])

ipcw_std = bootstrap_gap(lambda d: ipcw_tpr(d[d.group == "majority"]) - ipcw_tpr(d[d.group == "minority"]), cohort)
fullwait_std = bootstrap_gap(
    lambda d: tpr(d[(d.group == "majority") & (d.age >= MAX_DELAY)]) -
              tpr(d[(d.group == "minority") & (d.age >= MAX_DELAY)]), cohort)
print(f"bootstrap std: IPCW={ipcw_std:.4f}, full-wait={fullwait_std:.4f}")

Executed output:

TPR gap: true=0.2500, naive=0.1608, full-wait=0.2382, IPCW=0.2559
absolute error vs true gap: naive=0.0892, full-wait=0.0118, IPCW=0.0059
bootstrap std: IPCW=0.0216, full-wait=0.0238

The naive same-day estimate is off by 0.0892 from the true 0.2500 gap, badly understating the disparity because fast-resolving (expedited, correctly-flagged) transactions dominate the currently-observed sample. Waiting the full 7 days brings the error down to 0.0118, but the IPCW streaming estimator, using ALL of today's resolved labels reweighted by the Kaplan-Meier-estimated resolution probability, gets even closer (0.0059) while never requiring the 7-day wait, and its bootstrap standard deviation (0.0216) is actually slightly SMALLER than the full-wait estimator's (0.0238), because IPCW draws on a much larger effective sample (all currently-resolved transactions, not just the narrow age-7 slice) even after accounting for the extra variance the weights introduce.

Trade-offs and pitfalls

The most consequential mistake is stratifying (or not stratifying) the Kaplan-Meier fit incorrectly: it must condition only on variables known at PREDICTION time (the model's own decision, the group, other features present at scoring time), never on the eventual true label, since the true label is exactly the thing being censored and using it to build the weights would be circular. A second pitfall is applying a single, pooled resolution-time curve across all strata when resolution speed genuinely differs by stratum, as in the worked example; a pooled curve would systematically under-weight the slow-resolving, unflagged transactions relative to what a correctly stratified curve does, reintroducing much of the bias IPCW is meant to remove. Third, weights near the youngest ages (very small $\hat{F}(a)$) can become extremely large, so a practical implementation should floor the minimum age included (exclude transactions younger than 1 day, as done above) or cap the maximum weight, trading a small amount of remaining bias for materially lower variance. Finally, none of this substitutes for periodically validating the IPCW estimate against the eventual full-wait ground truth once it becomes available; if the two diverge by more than expected, that is the signal the survival model's assumed stratification has missed something, not a reason to distrust survival-analysis-based streaming metrics in general.

Direct answer

Detecting proxy features is a three-layer procedure: a fast statistical screen (correlation and mutual information against the sensitive attribute), a model-level explainability pass (SHAP or LIME feature importance, checked for overlap with the flagged statistical candidates), and a confirmatory check (can a simple classifier predict the sensitive attribute from the candidate feature alone, and does a counterfactual perturbation of the candidate flip the model's decision). Remediation then ranges from cheapest-and-weakest (drop the feature) to more robust but more expensive (adversarial debiasing, fair representation learning), with post-processing as a fallback when the model cannot be retrained.

Structured elaboration

Step 1: statistical screen. For every feature, compute mutual information (or Pearson/Spearman correlation for continuous pairs) against the sensitive attribute. This is cheap and catches most linear or simple nonlinear associations; it does NOT catch a proxy that only becomes predictive in combination with other features (an interaction-only proxy), which is why step 2 matters.

Step 2: explainability pass. Fit the actual production model (with the sensitive attribute excluded, as it should be for a compliant serving model) and compute global SHAP importance. A feature that ranks high in mean $|\text{SHAP}|$ AND was flagged in step 1 is a strong proxy candidate; a feature that is important in SHAP but was NOT flagged in step 1 is worth a second look for an interaction-driven proxy effect (the correlation screen alone would have missed it).

Step 3: confirmatory check. Train a small classifier (logistic regression is enough) that tries to predict the sensitive attribute FROM the candidate feature alone, cross-validated. An AUC meaningfully above 0.5, in practice greater than roughly 0.75, confirms the feature functions as a de facto label for the sensitive attribute even if no one intended it to. This step is what turns "looks suspicious" into "confirmed."

Remediation, in increasing order of engineering cost and typically effectiveness:

• Feature removal or orthogonalization: drop the feature, or regress it on the sensitive attribute and keep only the residual (the part of its signal NOT explained by group membership). Cheapest, but orthogonalization against a nonlinear downstream model can leave more of the original signal intact than intended: a linear residual still preserves the feature's rank order almost perfectly, and a tree ensemble can split on that rank order exactly as before, so this transform must be verified against the actual disparity metric post-transform, never assumed to work from the linear-algebra description alone.
• Adversarial debiasing: train the predictor jointly with an adversary that tries to recover the sensitive attribute from the model's internal representation, penalizing the main model when the adversary succeeds. Preserves more predictive signal than blunt removal, at the cost of a harder, less stable training loop.
• Fair representation learning: learn an intermediate representation explicitly constrained (via a distributional-distance penalty or an adversary) to be independent of the sensitive attribute, then train the final model on that representation. More modular and reusable across models, more engineering investment to build and validate.
• Post-processing: if retraining is not an option, calibrate per-group decision thresholds after the fact. Requires having (or being able to estimate) the sensitive attribute at evaluation time, but works on ANY already-deployed model with no retrain.

Worked example

Four candidate features, one true strong proxy (zip code, correlated with group at 90%), one legitimate predictor (income), one legitimate but group-independent predictor (years of experience), and one irrelevant control (random noise). The screen should flag zip code and only zip code:

import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.feature_selection import mutual_info_classif
import shap

np.random.seed(7)
n = 4000

group = np.random.binomial(1, 0.5, n)                       # protected attribute (held out of X)
zip_code = np.where(np.random.rand(n) < 0.9, group, 1 - group)  # strong proxy
income = np.random.normal(50000, 15000, n)                  # legitimate, ~independent of group
years_experience = np.clip(np.random.normal(8, 4, n), 0, None)  # legitimate, ~independent of group
random_noise = np.random.normal(0, 1, n)                    # irrelevant control feature

logit = -1.0 + 0.00005 * income + 0.05 * years_experience - 1.3 * group + np.random.normal(0, 1, n)
label = np.random.binomial(1, 1 / (1 + np.exp(-logit)))

feature_names = ["income", "years_experience", "zip_code", "random_noise"]
X = np.column_stack([income, years_experience, zip_code, random_noise])

# Step 1: correlation / mutual information scan against the sensitive attribute
mi_scores = mutual_info_classif(X, group, discrete_features=[False, False, True, False], random_state=0)
print("Mutual information with protected attribute (group):")
for name, mi in zip(feature_names, mi_scores):
    flag = "FLAGGED (proxy candidate)" if mi > 0.05 else ""
    print(f"  {name:18s} MI={mi:.4f} {flag}")

# Step 2: fit the actual production-style model on the candidate features (group excluded)
model = GradientBoostingClassifier(n_estimators=100, max_depth=3, random_state=0)
model.fit(X, label)

# Step 3: global explainability - rank features by mean |SHAP|
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X)
mean_abs_shap = np.abs(shap_values).mean(axis=0)
print("\nMean |SHAP| feature importance:")
order = np.argsort(-mean_abs_shap)
for i in order:
    print(f"  {feature_names[i]:18s} mean|SHAP|={mean_abs_shap[i]:.4f}")

# Step 4: confirm the top-ranked non-legitimate feature is actually predictive of the
# protected attribute (AUC of "predict group from this one feature alone")
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
zip_auc = cross_val_score(LogisticRegression(), zip_code.reshape(-1, 1), group, cv=5, scoring="roc_auc").mean()
noise_auc = cross_val_score(LogisticRegression(), random_noise.reshape(-1, 1), group, cv=5, scoring="roc_auc").mean()
print(f"\nAUC(zip_code -> group)     = {zip_auc:.3f}  (>0.75 threshold => confirmed proxy)")
print(f"AUC(random_noise -> group) = {noise_auc:.3f}  (near 0.5 => not a proxy)")

Executed output:

Mutual information with protected attribute (group):
  income             MI=0.0157 
  years_experience   MI=0.0000 
  zip_code           MI=0.3574 FLAGGED (proxy candidate)
  random_noise       MI=0.0000 

Mean |SHAP| feature importance:
  income             mean|SHAP|=0.5659
  zip_code           mean|SHAP|=0.3967
  years_experience   mean|SHAP|=0.1399
  random_noise       mean|SHAP|=0.1070

AUC(zip_code -> group)     = 0.895  (>0.75 threshold => confirmed proxy)
AUC(random_noise -> group) = 0.491  (near 0.5 => not a proxy)

Read top to bottom: the mutual-information screen immediately isolates zip_code (MI=0.357, everything else near zero). The SHAP pass on the trained model shows zip_code as the second-most-important feature overall (mean |SHAP|=0.397, just behind the legitimately predictive income at 0.566), which is exactly the "important AND flagged" combination that should trigger remediation. The confirmatory check nails it down: an AUC of 0.895 predicting group from zip_code alone versus 0.491 (indistinguishable from chance) for the random-noise control. years_experience never appears in either the correlation or the AUC screen, correctly clearing it as a legitimate, group-independent predictor despite also having real SHAP importance.

Trade-offs and pitfalls

The single most common mistake is running only the correlation screen and stopping there: it misses interaction-driven proxies where no single feature correlates with the sensitive attribute but a COMBINATION does (e.g. income conditional on zip code). Running the SHAP pass on the trained model, not just the raw feature table, is what catches those. A second pitfall is choosing remediation before confirming the proxy is real: dropping a feature that only weakly correlates with the sensitive attribute but carries substantial independent business signal can silently tank model performance for no measurable fairness gain, which is why the AUC confirmatory step exists, to make the trade-off decision on evidence rather than suspicion. Third, orthogonalization is often reached for as the "safe middle option" between removal and adversarial debiasing, but a LINEAR residualization only removes the LINEAR component of the correlation; if the downstream model is a nonlinear tree ensemble, it can often still recover most of the original signal from the residual's rank order alone, so orthogonalization needs to be validated post-hoc (recompute the disparity metric after transforming), never assumed to work.