python
import numpy as np
from scipy.stats import mannwhitneyu
from statsmodels.stats.multitest import multipletests
def bootstrap_mean_diff(x, y, n_boot=1000, seed=None):
rng = np.random.RandomState(seed)
diffs = []
nx, ny = len(x), len(y)
for _ in range(n_boot):
bx = rng.choice(x, size=nx, replace=True)
by = rng.choice(y, size=ny, replace=True)
diffs.append(np.mean(bx) - np.mean(by))
diffs = np.array(diffs)
return np.mean(diffs), np.percentile(diffs, [2.5, 97.5])
def compare_shap(shap_A, shap_B, feature_meta, n_boot=1000, alpha=0.05):
# shap_A / shap_B: dict feature -> array of per-sample shap contributions (signed or abs as chosen)
pvals, effects, cis, names = [], [], [], []
for feat, x in shap_A.items():
y = shap_B[feat]
names.append(feat)
# If categorical, feature_meta[feat]['type']=='cat' and meta provides categories array per sample
if feature_meta.get(feat, {}).get('type') == 'cat':
# test per-category, combine using max-statistic p-value (conservative) OR weighted mean
cats = feature_meta[feat]['categories'] # list of categories
cat_p = []
cat_effects = []
for c in cats:
maskA = feature_meta[feat]['values_A'] == c
maskB = feature_meta[feat]['values_B'] == c
if maskA.sum() < 5 or maskB.sum() < 5:
continue
_, ci = bootstrap_mean_diff(x[maskA], y[maskB], n_boot=n_boot)
cat_effects.append((np.mean(x[maskA]) - np.mean(y[maskB]), ci))
# nonparametric p-value
p = mannwhitneyu(x[maskA], y[maskB], alternative='two-sided').pvalue
cat_p.append(p)
pval = max(cat_p) if cat_p else 1.0
effect = np.mean([e[0] for e in cat_effects]) if cat_effects else 0.0
ci = np.mean([e[1] for e in cat_effects], axis=0) if cat_effects else (0,0)
else:
# numeric feature: compare distributions robustly
# use Mann-Whitney U for p-value and bootstrap for effect & CI
pval = mannwhitneyu(x, y, alternative='two-sided').pvalue
effect, ci = bootstrap_mean_diff(np.abs(x), np.abs(y), n_boot=n_boot)
pvals.append(pval); effects.append(effect); cis.append(ci)
# Multiple-testing correction (FDR)
rej, pvals_corr, _, _ = multipletests(pvals, alpha=alpha, method='fdr_bh')
results = []
for i, name in enumerate(names):
results.append({
'feature': name,
'pvalue': pvals[i],
'pvalue_adj': pvals_corr[i],
'significant': bool(rej[i]),
'effect_mean_diff': effects[i],
'ci_95': tuple(cis[i])
})
return results