Responsible AI: Fairness, Bias, and Interpretability Questions

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

Whether an observed subgroup performance gap reflects systemic bias or ordinary sampling variability is answered with two complementary tools: a permutation test, which builds an explicit null distribution by shuffling the group label and asks "how often would a gap this large appear if group membership had no real relationship to the outcome," and a bootstrap, which estimates how much the observed gap's SIZE would vary on repeated sampling and gives an interval rather than a single p-value. When many subgroups are tested at once, the resulting p-values must be corrected for multiple comparisons (Benjamini-Hochberg false-discovery-rate control is the usual choice) before any individual gap is called significant, and the entire analysis plan, including which comparisons will be tested and how they will be corrected, should be fixed before looking at the data that will be used for the final test, specifically to avoid p-hacking through cherry-picking which comparison to report.

Structured elaboration

The permutation test. Under the null hypothesis that group membership has no real relationship to the outcome, the observed labels and predictions would look the same no matter how group membership had been (counterfactually) assigned. Operationalize this by repeatedly shuffling the group label among the same individuals, recomputing the gap statistic (e.g. callback-rate difference) on each shuffle, and building an empirical null distribution of "gaps you'd see from noise alone." This test makes no assumption about the underlying distribution and is exact under the null by construction, which matters when a gap is measured from proportions or rates that are not well approximated by a normal distribution, especially at small counts.

Computing the p-value. Given $B$ permutations and an observed statistic $T_{obs}$, the two-sided permutation p-value is

where $T_b$ is the statistic recomputed on the $b$-th shuffled dataset. The "+1" in both numerator and denominator (counting the observed statistic itself as one of the $B+1$ possible outcomes under the null) guarantees the p-value is never reported as exactly zero, which would overstate certainty no matter how many permutations were run.

The bootstrap, and how it differs from the permutation test. The permutation test answers "is this gap distinguishable from noise." The bootstrap answers a different question: "how large is this gap, and how much would that size vary if we re-sampled from the same population." Resample within each group (with replacement, matching each group's own size) many times, recompute the gap on each resample, and take the empirical 2.5th/97.5th percentiles as a 95% interval on the gap's magnitude. The two tools are complementary and should usually be reported together: the permutation p-value for "is this real," the bootstrap interval for "how big, with how much uncertainty."

Multiple-hypothesis correction. Testing $m$ subgroups (or subgroup-by-category combinations) simultaneously inflates the chance that at least one shows a "significant" gap purely by chance, even if none of the underlying gaps are real. Apply Benjamini-Hochberg false-discovery-rate (FDR) control: sort the $m$ p-values ascending, find the largest rank $k$ such that $p_{(k)} \le \frac{k}{m}\alpha$, and flag every hypothesis at or below that rank. FDR control is generally preferred over a stricter Bonferroni correction ($\alpha/m$ for every test) in this setting because it controls the EXPECTED PROPORTION of false flags among flagged subgroups rather than the probability of any single false flag, preserving more power to detect real disparities as the number of tested subgroups grows.

Avoiding p-hacking via a pre-registered split. The single most damaging way this analysis goes wrong in practice is using the SAME data both to choose which comparison looks most interesting and to test that comparison's significance; a comparison selected for looking extreme will, by construction, look more extreme than it truly is (a form of selection bias sometimes called the "winner's curse"), so testing it on the same data does not produce a valid p-value even though the arithmetic looks identical to a legitimate test. The fix is to split the data BEFORE any analysis into a discovery portion (used freely to explore, generate hypotheses, or decide which subgroups are worth testing) and a confirmatory portion (used ONLY to run the final, pre-specified set of hypothesis tests with the pre-specified correction method), and to report every test that was run on the confirmatory set, not only the ones that came out significant.

Worked example

A resume-screening callback-rate audit across 5 job categories, with a paired-resume design (identical resumes except for a name signaling group A or B), where category 2 carries a genuine 0.17 gap and the other four categories have none:

import numpy as np
import pandas as pd

rng = np.random.default_rng(88)
n_categories, true_gap_category, base_callback, true_gap, n_per_cell = 5, 2, 0.24, 0.17, 400
records = []
for cat in range(n_categories):
    for group in (0, 1):
        p = base_callback + (true_gap / 2 if (cat == true_gap_category and group == 0) else 0) \
                           - (true_gap / 2 if (cat == true_gap_category and group == 1) else 0)
        for c in rng.random(n_per_cell) < p:
            records.append((cat, group, int(c)))
df = pd.DataFrame(records, columns=["category", "group", "callback"])

# pre-registered split, BEFORE any analysis
rng_split = np.random.default_rng(1)
df["is_discovery"] = rng_split.random(len(df)) < 0.40
discovery, confirmatory = df[df.is_discovery], df[~df.is_discovery]

def observed_gap(data, cat):
    sub = data[data.category == cat]
    return sub.loc[sub.group == 0, "callback"].mean() - sub.loc[sub.group == 1, "callback"].mean()

def permutation_pvalue(data, cat, n_perm, rng_p):
    sub = data[data.category == cat]
    obs = observed_gap(data, cat)
    group_vals, callback_vals = sub["group"].values, sub["callback"].values
    null_gaps = np.empty(n_perm)
    for i in range(n_perm):
        shuffled = rng_p.permutation(group_vals)  # ONE shuffle per iteration, reused for both means
        null_gaps[i] = callback_vals[shuffled == 0].mean() - callback_vals[shuffled == 1].mean()
    return obs, (np.sum(np.abs(null_gaps) >= abs(obs)) + 1) / (n_perm + 1)

# THE P-HACKING TRAP: cherry-pick the best-looking category from discovery, test on that SAME data
discovery_gaps = {cat: observed_gap(discovery, cat) for cat in range(n_categories)}
cherry = max(discovery_gaps, key=lambda c: abs(discovery_gaps[c]))
_, p_cherry = permutation_pvalue(discovery, cherry, 5000, np.random.default_rng(2))
print(f"cherry-picked category {cherry}, p on the SAME discovery data (circular): {p_cherry:.4f}")

# THE CORRECT APPROACH: pre-registered, all 5 categories, tested on the INDEPENDENT confirmatory set
results = [(cat, *permutation_pvalue(confirmatory, cat, 5000, np.random.default_rng(100+cat)))
           for cat in range(n_categories)]
res_df = pd.DataFrame(results, columns=["category", "observed_gap", "pvalue"]).sort_values("pvalue")
m = len(res_df)
res_df["bh_threshold"] = (np.arange(1, m + 1) / m) * 0.05
passes = res_df["pvalue"].values <= res_df["bh_threshold"].values
k = np.where(passes)[0].max() if passes.any() else -1
res_df["flagged"] = False
if k >= 0:
    res_df.iloc[:k+1, res_df.columns.get_loc("flagged")] = True
print(res_df.round(4).to_string(index=False))

# bootstrap 95% CI for the flagged category's gap (permutation test answers "is it real",
# bootstrap answers "how big"), computed on the SAME independent confirmatory set
flagged_cat = int(res_df.loc[res_df.flagged, "category"].iloc[0])
sub = confirmatory[confirmatory.category == flagged_cat]
idx_a, idx_b = sub.index[sub.group == 0], sub.index[sub.group == 1]
rng_boot = np.random.default_rng(500)
boot_gaps = [
    sub["callback"].loc[rng_boot.choice(idx_a, len(idx_a), replace=True)].mean() -
    sub["callback"].loc[rng_boot.choice(idx_b, len(idx_b), replace=True)].mean()
    for _ in range(5000)
]
boot_lo, boot_hi = np.percentile(boot_gaps, [2.5, 97.5])
print(f"bootstrap 95% CI for category {flagged_cat}'s gap: [{boot_lo:.4f}, {boot_hi:.4f}] "
      f"(point estimate {observed_gap(confirmatory, flagged_cat):.4f})")

Executed output:

cherry-picked category 2, p on the SAME discovery data (circular): 0.0004
 category  observed_gap  pvalue  bh_threshold  flagged
        2        0.1093  0.0052          0.01     True
        3        0.0107  0.8284          0.02    False
        0        0.0095  0.8294          0.03    False
        4        0.0093  0.8398          0.04    False
        1       -0.0066  0.9234          0.05    False
bootstrap 95% CI for category 2's gap: [0.0358, 0.1812] (point estimate 0.1093)

The permutation test correctly identifies category 2 as the only genuinely significant gap after Benjamini-Hochberg correction across all 5 pre-registered categories, and the bootstrap CI on that same category's gap, [0.0358, 0.1812], is consistent with the 0.17 true value used to generate the data, giving a sense of the gap's likely size on top of the permutation test's yes/no significance call.

Quantifying the p-hacking risk directly. To make the danger of the cherry-picking approach concrete rather than anecdotal, the same pipeline was run 300 times under a scenario with NO true effect anywhere (every category's true gap set to exactly zero), comparing the naive "pick the best-looking category in discovery, then test it on that same discovery data" approach against the proper "pre-register all categories, test on the independent confirmatory set, apply BH-FDR" approach:

import numpy as np
import pandas as pd

n_categories, n_per_cell, base_callback = 5, 400, 0.24  # NULL: true_gap = 0 in EVERY category

def build_null_dataset(seed):
    rng = np.random.default_rng(seed)
    records = []
    for cat in range(n_categories):
        for group in (0, 1):
            for c in rng.random(n_per_cell) < base_callback:
                records.append((cat, group, int(c)))
    return pd.DataFrame(records, columns=["category", "group", "callback"])

def observed_gap(data, cat):
    sub = data[data.category == cat]
    return sub.loc[sub.group == 0, "callback"].mean() - sub.loc[sub.group == 1, "callback"].mean()

def permutation_pvalue(data, cat, n_perm, rng_p):
    sub = data[data.category == cat]
    obs = observed_gap(data, cat)
    group_vals, callback_vals = sub["group"].values, sub["callback"].values
    null_gaps = np.empty(n_perm)
    for i in range(n_perm):
        shuffled = rng_p.permutation(group_vals)
        null_gaps[i] = callback_vals[shuffled == 0].mean() - callback_vals[shuffled == 1].mean()
    return obs, (np.sum(np.abs(null_gaps) >= abs(obs)) + 1) / (n_perm + 1)

def bh_flag(pvals, alpha=0.05):
    order = np.argsort(pvals)
    m = len(pvals)
    sorted_p = np.array(pvals)[order]
    thresh = (np.arange(1, m + 1) / m) * alpha
    passes = sorted_p <= thresh
    k = np.where(passes)[0].max() if passes.any() else -1
    flagged = np.zeros(m, dtype=bool)
    if k >= 0:
        flagged[order[:k + 1]] = True
    return flagged

N_REPEATS, N_PERM = 300, 1000
naive_false_positive, proper_false_positive = 0, 0
for rep in range(N_REPEATS):
    df = build_null_dataset(seed=10_000 + rep)
    rng_split = np.random.default_rng(20_000 + rep)
    df["is_discovery"] = rng_split.random(len(df)) < 0.40
    discovery, confirmatory = df[df.is_discovery], df[~df.is_discovery]

    gaps = {cat: observed_gap(discovery, cat) for cat in range(n_categories)}
    cherry = max(gaps, key=lambda c: abs(gaps[c]))
    _, p_naive = permutation_pvalue(discovery, cherry, N_PERM, np.random.default_rng(30_000 + rep))
    if p_naive < 0.05:
        naive_false_positive += 1

    rng_p2 = np.random.default_rng(40_000 + rep)
    pvals = [permutation_pvalue(confirmatory, cat, N_PERM, rng_p2)[1] for cat in range(n_categories)]
    if bh_flag(pvals, alpha=0.05).any():
        proper_false_positive += 1

print(f"naive cherry-pick-then-test-on-same-data false-positive rate: "
      f"{naive_false_positive}/{N_REPEATS} = {naive_false_positive/N_REPEATS:.4f}")
print(f"pre-registered split + confirmatory BH-FDR false-positive rate: "
      f"{proper_false_positive}/{N_REPEATS} = {proper_false_positive/N_REPEATS:.4f}")

Executed output:

naive cherry-pick-then-test-on-same-data false-positive rate: 66/300 = 0.2200
pre-registered split + confirmatory BH-FDR false-positive rate: 16/300 = 0.0533

Under a true null where NO category has a real effect, the naive cherry-pick-and-test-on-the-same-data approach falsely declares significance 22.00% of the time, more than four times the nominal 5% rate, purely because selecting "the most extreme of 5 categories" and then testing that same selection on the same data is a circular, invalid procedure. The pre-registered split-and-correct approach comes in at 5.33%, within ordinary Monte Carlo noise of the intended 5% false-positive rate, confirming it is properly calibrated.

Trade-offs and pitfalls

The single most damaging wrong turn, demonstrated numerically above, is testing a hypothesis on the same data used to select it; a p-value computed this way is not a valid p-value no matter how correctly the permutation or bootstrap machinery is implemented underneath it, because the selection step itself already used up some of the "surprise" that a valid p-value is supposed to measure. A second pitfall is treating Bonferroni correction as strictly safer than Benjamini-Hochberg; Bonferroni does control the probability of ANY false flag more strictly, but at a real cost in power, and for an exploratory fairness sweep across many subgroups, BH-FDR's expected-proportion guarantee is usually the more appropriate and commonly used standard. Third, a permutation test and a bootstrap CI can disagree in edge cases if the underlying statistic's sampling distribution is unusual (heavily skewed, or bounded near 0 or 1 at very small counts); when they disagree, trust the permutation test's exactness for the yes/no significance call and use the bootstrap for characterizing size and uncertainty, not for the significance decision itself. Finally, discipline about the discovery/confirmatory split has to be genuine: if the "confirmatory" set is quietly re-used after an unsatisfying result to try a different set of categories or a different statistic, the entire protection the split was meant to provide evaporates, so the pre-registration (which comparisons, which statistic, which correction method) needs to be written down and fixed BEFORE the confirmatory data is touched, not merely intended.

Direct answer

For a model predicting one of $K$ mutually exclusive classes, extend each binary fairness definition by treating every class in turn as the "positive" outcome in a one-versus-rest (OvR) framing: OvR selection-rate parity checks whether each specific class is assigned to different groups at different rates, and a per-class equalized-odds extension checks whether each class's recall (true-positive rate for that class) and false-positive rate are equal across groups. The one property that makes multi-class fairness genuinely different from a set of independent binary checks is that a single example gets exactly one predicted class, so the OvR selection rates across all $K$ classes must sum to 1 for every group, which means the per-class gaps are mathematically linked: a disparity favoring a group on one class must be offset by disparities against that group spread across the other classes. A concise dashboard should show the full class-by-group matrix as a heatmap, but lead with the worst single-class gap and a prevalence-weighted summary, not a single collapsed score.

Structured elaboration

One-versus-rest (OvR) parity checks. For each class $k \in \{1, \ldots, K\}$, define the binary indicator "was this example assigned to class $k$" and compute its selection rate per group, exactly as you would for demographic parity in a binary setting. Repeat for all $K$ classes to get a full class-by-group selection-rate table. The key structural fact to keep in mind: because classes are mutually exclusive, $\sum_{k=1}^{K} \text{SelectionRate}_g(k) = 1$ for every group $g$, so $\sum_{k=1}^{K} \big(\text{SelectionRate}_1(k) - \text{SelectionRate}_0(k)\big) = 0$: the per-class gaps are forced to sum to zero across classes. This is the sharpest distinction from a multi-task fairness setting with independent binary labels, where per-task gaps have no such constraint; in the multi-class case, "fair on every class" is a much stronger and more internally consistent requirement, since every gap you close on one class necessarily closes an offsetting gap somewhere else.

Per-class equalized-odds extension. For each class $k$, treat "is the true class $k$" as the binary positive label and compute, per group: recall for class $k$ (of examples whose TRUE class is $k$, what fraction the model correctly predicted as $k$) and the false-positive rate for class $k$ (of examples whose true class is NOT $k$, what fraction the model incorrectly predicted as $k$). Equalized odds for class $k$ is satisfied when both of these match across groups. Unlike the OvR selection rates, per-class recall and FPR values are NOT constrained to sum to anything in particular across classes, so this check does not inherit the same "offsetting" structure; it can reveal genuinely independent problems at each tier even after the selection-rate view has been accounted for.

Presenting multi-class fairness results in a concise dashboard. A full $K \times G \times \text{(several metrics)}$ table is too much to put in front of a stakeholder directly. Structure the dashboard in layers: a heatmap with classes as rows, groups (or the gap between them) as columns, one for selection-rate gaps and one for recall/FPR gaps, as the primary visual (a heatmap makes the worst cell immediately visible without reading numbers); a headline pair of numbers, the single WORST per-class gap (protects against averaging hiding a severe single-tier problem, exactly analogous to the worst-task summary in a multi-task setting) and a PREVALENCE-weighted average absolute gap (a business-relevant summary that weights each class by how often it actually occurs, rather than treating a rare tier and a common tier as equally important by default); and the full table available on drill-down for anyone auditing a specific class. Never report only the OvR selection-rate gaps without the per-class equalized-odds numbers alongside them, since the sum-to-zero constraint on selection rates means a model can show a deceptively small overall selection-rate profile while still failing badly on recall or FPR at a specific tier.

Worked example

A four-tier risk-scoring model (low, medium, high, very-high risk) where the TRUE risk distribution is identical across two groups, but the model systematically pushes group 1 up by one tier far more often than group 0 (a realistic miscalibration failure mode, not a difference in underlying risk):

import numpy as np
import pandas as pd

rng = np.random.default_rng(19)
n_per_group = 4000
classes = ["low_risk", "medium_risk", "high_risk", "very_high_risk"]
K = len(classes)
true_class_probs = [0.40, 0.30, 0.20, 0.10]

records = []
for g in (0, 1):
    true_class = rng.choice(K, size=n_per_group, p=true_class_probs)
    push_up_prob = 0.06 if g == 0 else 0.22
    push_up = rng.random(n_per_group) < push_up_prob
    pred_class = np.minimum(true_class + push_up.astype(int), K - 1)
    noise = rng.random(n_per_group) < 0.05
    pred_class = np.where(noise, rng.integers(0, K, n_per_group), pred_class)
    for t, p in zip(true_class, pred_class):
        records.append((g, t, p))
df = pd.DataFrame(records, columns=["group", "true_class", "pred_class"])

ovr_rows = [(classes[k], g, (df[df.group == g].pred_class == k).mean())
            for k in range(K) for g in (0, 1)]
ovr_pivot = pd.DataFrame(ovr_rows, columns=["class", "group", "selection_rate"]) \
    .pivot(index="class", columns="group", values="selection_rate").loc[classes]
ovr_pivot["gap"] = ovr_pivot[1] - ovr_pivot[0]
print(ovr_pivot.round(4).to_string())
print(f"sum of gaps across classes: {ovr_pivot['gap'].sum():.6f}")

eo_rows = []
for k in range(K):
    for g in (0, 1):
        sub = df[df.group == g]
        actual_k, pred_k = sub.true_class == k, sub.pred_class == k
        recall_k = (actual_k & pred_k).sum() / max(actual_k.sum(), 1)
        fpr_k = ((~actual_k) & pred_k).sum() / max((~actual_k).sum(), 1)
        eo_rows.append((classes[k], g, recall_k, fpr_k))
eo_pivot = pd.DataFrame(eo_rows, columns=["class", "group", "recall", "fpr"]) \
    .pivot(index="class", columns="group", values=["recall", "fpr"]).loc[classes]
eo_pivot[("recall", "gap")] = eo_pivot[("recall", 1)] - eo_pivot[("recall", 0)]
eo_pivot[("fpr", "gap")] = eo_pivot[("fpr", 1)] - eo_pivot[("fpr", 0)]
print(eo_pivot.round(4).to_string())

prevalence = dict(zip(classes, true_class_probs))
weighted_recall_gap = sum(prevalence[c] * abs(eo_pivot.loc[c, ("recall", "gap")]) for c in classes)
weighted_fpr_gap = sum(prevalence[c] * abs(eo_pivot.loc[c, ("fpr", "gap")]) for c in classes)
print(f"worst recall gap: {eo_pivot[('recall','gap')].abs().max():.4f}")
print(f"prevalence-weighted average |recall gap|: {weighted_recall_gap:.4f}")
print(f"prevalence-weighted average |FPR gap|: {weighted_fpr_gap:.4f}")

Executed output:

group                0       1     gap
class
low_risk        0.3638  0.3085 -0.0553
medium_risk     0.3088  0.3165  0.0077
high_risk       0.2068  0.2275  0.0208
very_high_risk  0.1208  0.1475  0.0267
sum of gaps across classes: -0.000000

                recall             fpr          recall     fpr
group                0       1       0       1     gap     gap
class
low_risk        0.9033  0.7644  0.0144  0.0112 -0.1389 -0.0033
medium_risk     0.9016  0.7550  0.0489  0.1268 -0.1466  0.0779
high_risk       0.8881  0.7647  0.0354  0.0898 -0.1234  0.0545
very_high_risk  0.9457  0.9647  0.0278  0.0575  0.0191  0.0296

worst recall gap: 0.1466
prevalence-weighted average |recall gap|: 0.1261
prevalence-weighted average |FPR gap|: 0.0385

The OvR selection-rate gaps sum to exactly zero across the four classes (-0.000000), confirming the structural constraint: group 1 is under-selected into low_risk (-0.0553) and over-selected into every higher tier, and those movements exactly balance. But the per-class equalized-odds numbers show the real harm the selection-rate view alone would understate: group 1's recall is 0.10 to 0.15 points WORSE than group 0's at every tier except the top one (0.7644 vs 0.9033 at low_risk, 0.7550 vs 0.9016 at medium_risk), meaning group 1 individuals who genuinely belong in a lower-risk tier are substantially less likely to be correctly recognized there, while group 1's false-positive rate into medium_risk and high_risk is markedly elevated (0.1268 vs 0.0489, and 0.0898 vs 0.0354), directly reflecting the push-up miscalibration. The worst single-class recall gap (0.1466, medium_risk) and the prevalence-weighted average (0.1261) both correctly signal a serious, broad problem that a dashboard showing only the (zero-summing) OvR selection-rate gaps would have completely missed.

Trade-offs and pitfalls

The most damaging wrong turn is reporting only the OvR selection-rate table and treating a small or near-zero sum as evidence of fairness; the worked example shows the model has a severe recall and false-positive-rate problem at every non-top tier even while its selection-rate gaps individually look modest and necessarily sum to zero. A second pitfall is applying a per-task-style naive averaging instinct (like a straight, unweighted mean of the per-class gaps) without recognizing that the OvR gaps are constrained to cancel by construction, so an unweighted average of them is close to structurally guaranteed to look small regardless of the model's real behavior; always pair it with the per-class equalized-odds numbers, which have no such constraint. Third, a prevalence-weighted summary can under-represent a severe problem concentrated in a rare but high-stakes class (the very-high-risk tier here, at only 10% prevalence, would be down-weighted in the summary even if it carried the most serious real-world consequences), so the weighting scheme itself, and whether it should reflect volume versus severity of consequence, needs to be a deliberate, disclosed choice rather than a default. Finally, resist collapsing the class-by-group matrix into one dashboard score; a heatmap with a worst-cell callout preserves the information a stakeholder needs to know WHERE the problem is, which a single number by construction throws away.

Direct answer

A from-scratch partial dependence function needs three pieces: a grid-builder that produces evenly spaced percentile points for continuous features but the feature's own observed unique values for categorical features, an outer loop over that grid that overwrites the column of interest across the whole reference dataset and averages model.predict(), and (for two features) a nested version of the same loop producing a 2D grid. Naively this costs one predict() call per grid point; the practical speedup, especially valuable when features are correlated and the reference set is large, is to stack every grid point's full copy of the dataset into ONE big array and issue a single vectorized predict() call instead of many small ones.

Structured elaboration

API design. The function needs to support 1 or 2 features (by column index), a grid_resolution for continuous features, and a categorical_features set telling the grid-builder which columns should use unique observed values instead of a percentile grid. This mirrors the exact distinction scikit-learn's own partial_dependence makes, without depending on it.

Grid construction. For a continuous feature, use the 5th to 95th percentile range (not min/max, to avoid extrapolating into sparse tails where the model was barely trained) split into grid_resolution evenly spaced points. For a categorical feature, the grid IS the feature's finite set of observed values; there is no "resolution" concept, since interpolating between category codes 1 and 2 is meaningless.

Core PDP computation. For each grid value $v$, replace every row's value for the feature(s) of interest with $v$, call model.predict() on the whole (modified) reference set, and average the predictions. This is exactly the Monte Carlo estimate of $\hat{f}_S(x_S) = \frac{1}{n}\sum_{i=1}^n \hat{f}(x_S, x_{C}^{(i)})$ introduced in the PDP/ICE definitional answer: fix the feature(s) of interest, keep everything else at its real observed value, average over the reference rows.

Runtime complexity. For a 1-feature PDP with grid resolution $G$ over $n$ reference rows: $O(G \cdot n)$ total predictions, issued as $G$ separate predict() calls of $n$ rows each in the naive loop version. For a 2-feature PDP with resolutions $G_0, G_1$: $O(G_0 \cdot G_1 \cdot n)$ predictions, issued as $G_0 \cdot G_1$ separate calls. This is the dominant cost; the grid-construction step is $O(n)$ or $O(n \log n)$ (percentile computation) and negligible by comparison.

Speeding it up, especially for correlated features.

1. Vectorize the grid loop. Instead of $G$ (or $G_0 \cdot G_1$) separate predict() calls, np.tile the reference set once per grid point into one big stacked array of shape $(G \cdot n, p)$, overwrite the relevant column(s) per block with np.repeat/np.tile, and issue exactly ONE predict() call. Most models (tree ensembles, neural nets, linear models) are far more efficient per-row when called on one large batch than on many small ones, because they amortize whatever fixed per-call overhead the framework has (Python dispatch, JIT warm-up, batch-normalization mode checks) across every row instead of paying it $G$ times. The trade-off is memory: the stacked array uses $G \cdot n \cdot p$ floats, so for a very large reference set or very fine grid you may need to sub-sample the reference set instead of using the full training data.
2. Sub-sample the reference set. Since the PDP is already a Monte Carlo AVERAGE over the reference rows, using a random sample of a few thousand rows instead of the full dataset gives a nearly identical curve with far lower cost per grid point, trading a small amount of Monte Carlo noise for a large constant-factor speedup; this matters more, not less, as the grid gets finer.
3. Correlated-feature-specific speedup: restrict the grid to on-manifold combinations. When the feature(s) of interest are correlated with features you are holding fixed, a naive percentile grid on the marginal distribution of $x_S$ generates rows that never occur in the joint distribution (for example sweeping "house size" to its 95th percentile while every held-fixed row still has "number of bedrooms" at its own, now-inconsistent, observed value). Rather than evaluating the model at every combination, restrict the grid to the OBSERVED joint values actually present in a nearest-neighbor ball around each reference row (an ICE-like conditional grid), which both avoids extrapolating into unrealistic regions and, because you no longer sweep the full percentile range for every row, cuts the number of predictions needed per row.

Worked example (implementation, without sklearn's built-in PDP)

import numpy as np

def partial_dependence(model, X, features, grid_resolution=15, categorical_features=()):
    """From-scratch partial dependence. features: tuple of 1 or 2 column indices."""
    features = tuple(features)
    if len(features) not in (1, 2):
        raise ValueError("supports exactly 1 or 2 features")

    def make_grid(col_idx):
        col = X[:, col_idx]
        if col_idx in categorical_features:
            return np.unique(col)
        lo, hi = np.percentile(col, [5, 95])
        return np.linspace(lo, hi, grid_resolution)

    grids = [make_grid(f) for f in features]
    if len(features) == 1:
        f0 = features[0]
        pdp = np.empty(len(grids[0]))
        Xg = X.copy()
        for i, v in enumerate(grids[0]):
            Xg[:, f0] = v
            pdp[i] = model.predict(Xg).mean()
        return grids, pdp

    f0, f1 = features
    pdp = np.empty((len(grids[0]), len(grids[1])))
    Xg = X.copy()
    for i, v0 in enumerate(grids[0]):
        Xg[:, f0] = v0
        for j, v1 in enumerate(grids[1]):
            Xg[:, f1] = v1
            pdp[i, j] = model.predict(Xg).mean()
    return grids, pdp

def partial_dependence_vectorized(model, X, features, grid_resolution=15, categorical_features=()):
    """Same semantics, ONE predict() call per feature count (1 or 2)."""
    features = tuple(features)
    def make_grid(col_idx):
        col = X[:, col_idx]
        if col_idx in categorical_features:
            return np.unique(col)
        lo, hi = np.percentile(col, [5, 95])
        return np.linspace(lo, hi, grid_resolution)
    n = X.shape[0]
    if len(features) == 1:
        f0 = features[0]
        grid0 = make_grid(f0)
        G = len(grid0)
        Xrep = np.tile(X, (G, 1))
        Xrep[:, f0] = np.repeat(grid0, n)
        preds = model.predict(Xrep).reshape(G, n)
        return [grid0], preds.mean(axis=1)
    f0, f1 = features
    grid0, grid1 = make_grid(f0), make_grid(f1)
    G0, G1 = len(grid0), len(grid1)
    Xrep = np.tile(X, (G0 * G1, 1))
    v0 = np.repeat(grid0, G1 * n)
    v1 = np.tile(np.repeat(grid1, n), G0)
    Xrep[:, f0] = v0
    Xrep[:, f1] = v1
    preds = model.predict(Xrep).reshape(G0, G1, n)
    return [grid0, grid1], preds.mean(axis=2)

# ---- correctness check against a hand-computable case ----
from sklearn.ensemble import RandomForestRegressor
rng = np.random.default_rng(11)
n = 800
x0 = rng.uniform(-2, 2, n)
x1 = rng.uniform(-2, 2, n)
cat = rng.integers(0, 3, n)
y = x0 + 0.5 * cat + rng.normal(0, 0.05, n)  # additive, no interaction: PDP should match ground truth
X = np.column_stack([x0, x1, cat.astype(float)])
model = RandomForestRegressor(n_estimators=200, max_depth=8, random_state=0).fit(X, y)
print("Train R2:", round(model.score(X, y), 4))

grids, pdp0 = partial_dependence(model, X, features=(0,), grid_resolution=9)
slope = np.polyfit(grids[0], pdp0, 1)[0]
print("Fitted slope from PDP curve:", round(slope, 3), "(true slope: 1.0)")

grids_cat, pdp_cat = partial_dependence(model, X, features=(2,), categorical_features={2})
print("Categorical PDP (true effect +0.5/level):", np.round(pdp_cat, 3))

grids2, pdp2 = partial_dependence(model, X, features=(0, 1), grid_resolution=6)
print("2-feature PDP grid shape:", pdp2.shape, " row-std across x1 (no true effect):", np.round(pdp2.std(axis=1), 3))

grids_v, pdp0_v = partial_dependence_vectorized(model, X, features=(0,), grid_resolution=9)
print("Loop vs vectorized match:", np.allclose(pdp0, pdp0_v, atol=1e-9))

Executed output:

Train R2: 0.9993
Fitted slope from PDP curve: 1.001 (true slope: 1.0)
Categorical PDP (true effect +0.5/level): [-0.068  0.426  0.928]
2-feature PDP grid shape: (6, 6)  row-std across x1 (no true effect): [0.007 0.003 0.003 0.004 0.009 0.004]
Loop vs vectorized match: True

The fitted slope of the recovered PDP curve (1.001) matches the true generating slope (1.0) to three decimal places, the categorical PDP recovers the +0.5-per-level ground truth (steps of about 0.494 and 0.502 between levels), the 2-feature grid has the expected (6, 6) shape with near-zero row-wise standard deviation across the no-effect axis (all under 0.01, confirming $x_1$ correctly shows no influence), and the vectorized implementation is numerically identical to the naive loop version (allclose returns True). For a 2-feature PDP at grid_resolution=12 (144 grid cells) on this 800-row reference set, the loop version issues 144 separate predict() calls of 800 rows each, while the vectorized version issues exactly 1 predict() call of 115,200 rows: the same total prediction work, restructured into one call instead of 144. Wall-clock numbers are hardware-dependent and are deliberately not reported; the call-count and total-row-count structure above is what actually determines relative cost and is verified directly from the code.

Edge cases and complexity summary

• Edge cases handled: categorical features use their own observed values as the grid (never interpolated); continuous features use percentile bounds (5th to 95th) rather than min/max to avoid extrapolating into sparse tails; the reference dataset is copied (.copy()) before each grid-point overwrite so the original X is never mutated between iterations.
• Complexity: $O(G \cdot n)$ for 1 feature, $O(G_0 \cdot G_1 \cdot n)$ for 2 features, in prediction-row count; the vectorized version has the same total FLOP count but restructures it into $O(1)$ predict calls at the cost of $O(G \cdot n \cdot p)$ peak memory for the stacked array.

Trade-offs and pitfalls

The vectorized version trades memory for call count: for a very fine grid ($G$ in the hundreds) on a very large reference set, materializing the full stacked array can exceed available memory before it ever helps, so a middle ground (batching a few dozen grid points per predict() call) is often the practical choice, not the two extremes shown here. A second pitfall is grid extrapolation: using the full min-max range instead of a percentile-trimmed range can push the sweep into regions with almost no training data, where the model's predictions are unreliable extrapolations rather than genuine learned behavior, and the PDP curve will silently reflect that unreliability as if it were signal. Finally, none of the speedups here address the CORRELATED-FEATURES problem itself (a flat, on-manifold-restricted grid is a mitigation, not a full fix); a PDP computed this way over strongly correlated features is still evaluating the model at some combinations it never saw during training unless the on-manifold restriction described above is actually implemented.

Direct answer

Responsible machine learning is the practice of engineering a model so that its behavior in the real world is safe, equitable, explainable, and accountable, not just statistically accurate on a held-out test set. It rests on four dimensions that a model can score well on accuracy while still failing badly on: privacy, fairness, interpretability, and governance. A model can be 95% accurate and still leak the training data it memorized, systematically disadvantage a protected group, be un-auditable when a regulator asks why it denied someone, or have no owner who catches any of that before it ships. Responsible ML is the set of practices that closes those gaps deliberately rather than by accident.

Structured elaboration

The four dimensions, defined:

1. Privacy. Protecting the personal or sensitive data a model is trained on and queried with from unauthorized access, leakage, or re-identification. This covers data minimization and anonymization at training time, and defending against inference-time attacks such as membership inference (an attacker determines whether a specific record was in the training set) and model inversion (an attacker reconstructs approximate training examples from model outputs).
2. Fairness. Ensuring the model's decisions or error rates do not unjustifiably disadvantage people because of a protected characteristic, measured against an explicitly chosen fairness criterion (demographic parity, equalized odds, calibration within groups, and so on, each of which encodes a different notion of "fair" and can conflict with the others).
3. Interpretability. The degree to which a human, whether a developer debugging the model, an auditor checking compliance, or an affected end user, can understand why the model produced a specific output or how it behaves in general. This is delivered either by choosing an inherently interpretable model class (a shallow decision tree, a linear model) or by applying post-hoc explanation techniques (SHAP, LIME, feature-attribution methods) to a more complex model.
4. Governance. The organizational infrastructure that makes the first three dimensions durable rather than a one-time engineering effort: documented model ownership, a pre-deployment review process (a model card, a sign-off gate), incident response procedures for when something goes wrong in production, and audit trails that let someone reconstruct what the model did and why months later.

Two concrete harms responsible ML aims to avoid:

1. An allocative fairness harm from training on biased historical data. A well-documented real case: an internal resume-screening tool was found to systematically penalize resumes containing the word "women's" (as in "women's chess club captain") and downgrade graduates of women's colleges, because it had been trained on a decade of past hiring decisions in a male-dominated field and had learned to treat markers of female gender as negative signal. The tool was scrapped before wide deployment, but it illustrates the core failure mode: the model was accurately learning the pattern in its training data, and the pattern itself was the harm.
2. A privacy harm from insufficient anonymization. Researchers re-identified individuals in a supposedly anonymized public movie-ratings dataset by cross-referencing rating patterns and dates against public reviews on another site, showing that removing direct identifiers (names, account IDs) is not sufficient to prevent re-identification when enough auxiliary structure survives in the data. The harm is not hypothetical: several individuals in that dataset were identifiable well enough to reveal information (such as apparent sexual orientation, inferred from rental patterns) they had not intended to disclose.

Three KPIs to track in production, one per operational concern:

1. Demographic parity gap or disparate-impact ratio, computed per protected group, tracked continuously on live traffic. This is a fairness KPI: it directly measures whether the model's positive-decision rate differs across groups by more than an agreed tolerance, and it is cheap enough to compute on every batch of scored traffic that it can alert automatically rather than waiting for a periodic audit.
2. Differential-privacy budget consumption or PII (personally identifiable information) leakage-scan hit rate on logged model inputs and outputs. This is a privacy KPI: for systems using formal differential privacy it is the cumulative privacy budget spent to date against an approved ceiling; for systems without formal DP guarantees it is the fraction of logged requests or model outputs that an automated PII scanner flags as containing unredacted sensitive fields, which should trend toward zero and alert on any increase.
3. Human-override rate on predictions above an escalation or high-uncertainty threshold. This is a joint interpretability and governance KPI: it measures whether the human-in-the-loop review step that governance policy requires is actually functioning as a check, rather than becoming a rubber stamp. A rate that drifts to near zero for a sustained period is itself a signal worth investigating, since it can mean either that the model has become excellent, or that reviewers have stopped meaningfully engaging with the explanations they are shown.

Worked example

The demographic parity gap and disparate-impact ratio for KPI 1, computed on a synthetic loan-approval scenario where two groups' underlying scores are deliberately offset (standing in for the kind of historically-encoded pattern the two harm examples above describe):

import numpy as np

rng = np.random.default_rng(42)
n = 5000

group = rng.integers(0, 2, n)  # 0 = group A, 1 = group B
base_score = rng.normal(0, 1, n)
score = np.where(group == 0, base_score + 0.35, base_score - 0.15)
approve = (score > 0.3).astype(int)

p_a = approve[group == 0].mean()
p_b = approve[group == 1].mean()
dp_gap = abs(p_a - p_b)
di_ratio = min(p_a, p_b) / max(p_a, p_b)

print(f"P(approve | group A) = {p_a:.4f}")
print(f"P(approve | group B) = {p_b:.4f}")
print(f"demographic parity gap = {dp_gap:.4f}")
print(f"disparate impact ratio (min/max) = {di_ratio:.4f}")
print(f"four-fifths rule (0.80 threshold) violated: {di_ratio < 0.80}")

Executed output:

P(approve | group A) = 0.5289
P(approve | group B) = 0.3259
demographic parity gap = 0.2029
disparate impact ratio (min/max) = 0.6163
four-fifths rule (0.80 threshold) violated: True

The two KPI forms are:

$DP_{gap} = \left| P(\hat{y}=1 \mid A=a) - P(\hat{y}=1 \mid A=b) \right|$

$DI_{ratio} = \frac{\min\big(P(\hat{y}=1\mid A=a),\, P(\hat{y}=1\mid A=b)\big)}{\max\big(P(\hat{y}=1\mid A=a),\, P(\hat{y}=1\mid A=b)\big)}$

In this synthetic run, group A is approved 52.89% of the time versus 32.59% for group B, a gap of 0.2029 and a disparate-impact ratio of 0.6163, below the commonly cited 0.80 "four-fifths" threshold used in U.S. employment-selection adverse-impact analysis. A production monitoring job would compute exactly this pair of numbers on a rolling window of live decisions per group and alert when the ratio crosses the chosen threshold, rather than waiting for a quarterly fairness audit to notice.

Trade-offs and pitfalls

Treating any single KPI, including the demographic-parity gap above, as a complete fairness guarantee is the most common wrong turn: it captures only one fairness criterion (statistical independence between the prediction and the protected attribute) and says nothing about calibration within groups or equalized error rates, which can move in the opposite direction as parity is optimized. The four-fifths rule itself is a legal heuristic from a specific regulatory context (adverse-impact analysis in U.S. employment selection), not a universal fairness threshold, and applying it uncritically to, say, a medical-risk model can be actively wrong. Governance is the dimension most likely to be under-resourced in practice, because unlike an accuracy improvement it has no obvious short-term metric of its own and requires organizational buy-in (a review board, a sign-off gate) that a single engineer cannot manufacture unilaterally. Finally, a privacy KPI that only measures data-at-rest encryption or access-control compliance misses the inference-time attack surface entirely; membership-inference and re-identification risk, as in the movie-ratings example above, exists even when every access-control checkbox is green.