RoSHAP: A Distributional Framework and Robust Metric for Stable Feature Attribution∗
arXiv:2605.15154v1 [stat.ML] 14 May 2026
Lanxin Xiang Department of Statistics Virginia Tech [email protected]
Liang Shi Department of Statistics Virginia Tech Transportation Institute Virginia Tech [email protected]
Boyu Jiang Department of Statistics Virginia Tech [email protected]
Youhui Ye Department of Statistics Virginia Tech [email protected]
Dawei Zhou Department of Computer Science Virginia Tech [email protected]
Feng Guo† Department of Statistics Virginia Tech Transportation Institute Virginia Tech [email protected]
Abstract Feature attribution analysis is critical for interpreting machine learning models and supporting reliable data-driven decisions. However, feature attribution measures often exhibit stochastic variation: different train–test splits, random seeds, or model-fitting procedures can produce substantially different attribution values and feature rankings. This paper proposes a framework for incorporating stochastic nature of feature attribution and a robust attribution metric, RoSHAP, for stable feature ranking based on the SHAP metric. The proposed framework models the distribution of feature attribution scores and estimates it through bootstrap resampling and kernel density estimation. We show that, under mild regularity conditions, the aggregated feature attribution score is asymptotically Gaussian, which greatly reduces the computational cost of distribution estimation. The RoSHAP summarizes the distribution of SHAP into a robust feature-ranking criterion that simultaneously rewards features that are active, strong, and stable. Through simulations and real-data experiments, the proposed framework and RoSHAP outperform standard single-run attribution measures in identifying signal features. In addition, models built using RoSHAP-selected features achieve predictive performance comparable to full-feature models while using substantially fewer predictors. The proposed RoSHAP approach improves the stability and interpretability of machine learning models, enabling reliable and consistent insights for analysis.
1
Introduction
Feature attribution provides a quantitative basis for interpreting machine learning (ML) models by assigning contribution scores to individual predictors. For many modern ML models, however, ∗ Code is available at https://github.com/Lanxin-Xiang/RobustSHAP. † Corresponding author.
Preprint.
these scores cannot be directly derived from the model structure and must instead be estimated using post hoc attribution methods. This task becomes especially challenging when the data are high-dimensional, noisy, and correlated, or when the predictor–response relationship is complex. In such settings, attribution results may be difficult to interpret and may not provide a stable basis for scientific or practical conclusions. A large body of work has developed feature attribution methods for explaining model predictions, including SHapley Additive exPlanations (SHAP) [1], Local Interpretable Model-agnostic Explanations (LIME) [2], DeepLIFT [3], and information-gain-based measures. These metrics capture different aspects of feature contribution, each with distinct theoretical properties and practical tradeoffs. SHAP, for instance, produces additive instance-level attribution scores grounded in cooperative game theory and is the unique additive attribution method satisfying local accuracy, missingness, and consistency [1]. These guarantees have made SHAP a popular choice across diverse domains, including medicine [4–7], natural language processing [8], and spatial data analysis [9]. While existing feature attribution methods are valuable, challenges remain in how attribution scores are estimated and interpreted. A standard workflow splits the data into training and test sets, fits a model, and computes attribution values from the trained model. Each step can introduce stochastic variation, including random data partitioning, random initialization, and stochastic optimization. Consequently, attribution scores are not fixed quantities; different train–test splits or random seeds can yield substantially different feature rankings [10]. This uncertainty is especially pronounced in high-dimensional settings such as genomics, where the number of features far exceeds the number of samples. Figure 1 illustrates this issue: the top-three SHAP-ranked genes for cancer classification vary across train–test splits, and their attribution magnitudes also differ substantially. Such randomness provides an unreliable basis for costly downstream gene testing.
(a) Round 1
(b) Round 2
(c) Round 3
Figure 1: Top SHAP-ranked features in the Golub data across different train–test sets. Although all three runs achieve perfect macro F1, the selected top features differ substantially. We propose a bootstrap-based framework for estimating the distribution of feature attribution scores. Based on this framework, we introduce RoSHAP (Robust SHAP), a robust attribution metric for stable feature ranking. The framework repeatedly resamples the data and refits the model to capture variability in feature attribution scores. The resulting attribution scores are summarized through kernel density estimation, yielding an estimated distribution that reflects both attribution magnitude and variability. We further derive a Gaussian approximation for the aggregated feature-level attribution score under a mild regularity condition, substantially reducing the computational cost of distribution estimation. RoSHAP summarizes the estimated attribution distribution through three components: activity, strength, and stability. Activity measures the proportion of nonzero attribution values, strength measures the median magnitude of the nonzero attribution values, and stability measures the signalto-noise ratio of the attribution distribution. By rewarding features that contribute frequently, have large attribution magnitudes, and remain stable across resampling runs, RoSHAP provides a robust basis for feature ranking and interpretation. Results on both simulated and real-world datasets show that RoSHAP yields stable interpretations across regression and classification tasks and across diverse data modalities. RoSHAP effectively identifies the most influential features, and models trained using only these selected features achieve performance comparable to that of models trained on the full feature set. The remainder of this paper is organized as follows. Section 2 reviews related work on SHAP and other feature contribution methods. Section 4 reports simulation studies under different settings. Section 5 presents three motivating examples and performance comparisons. Section 6 concludes with a discussion of the main findings and future directions. 2
2
Related Works
Feature attribution and additive attribution. Feature attribution is widely used to interpret predictive models. For classical linear models, p-values, t-tests, and F -tests provide feature-level inference under parametric assumptions. For ML models, common importance measures include gain-, impurity-, and entropy-based scores, as well as post hoc methods such as SHAP [1], LIME [2], and DeepLIFT [3]. SHAP is grounded in cooperative game theory, treating predictors as players and the model prediction as the payoff; each feature is assigned its average marginal contribution over all coalitions. This additive formulation makes SHAP broadly applicable across model classes. TreeSHAP [11] enables efficient computation for tree ensembles such as random forests, XGBoost [12], and LightGBM [13], while KernelSHAP and DeepSHAP extend SHAP to general black-box and deep learning models. SHAP stability and prediction uncertainty. Growing literature examines how feature attributions vary under data and model randomness. Shaikhina et al. [14] showed that predictive uncertainty can propagate to explanation uncertainty, producing unstable attribution scores, especially for out-ofdistribution samples. Other studies use repeated resampling or nested cross-validation to summarize SHAP values across folds through averages, selection frequencies, or rank distributions [15, 16, 10]. Related work explains predictive uncertainty: Watson et al. [17] extended SHAP with informationtheoretic value functions to attribute epistemic and aleatoric uncertainty. However, most approaches focus on averages or rankings rather than the full distributional structure of attribution scores. SHAP for feature selection. SHAP is also used for feature selection, typically by ranking variables by mean absolute SHAP values and retaining the top-ranked subset. This strategy has been adopted across application domains [18–20]. Comparative studies have evaluated SHAP-based feature selection against built-in model importance measures and other screening rules [21]. While often effective, these approaches remain largely magnitude-based, with limited attention to resampling instability, exact zeros, and strong predictor correlation. Overall, prior work has established SHAP as a useful tool for interpretation, uncertainty analysis, and feature selection. However, existing approaches mainly focus on average attribution magnitude, ranking stability, or predictive uncertainty, rather than the full empirical distribution of attribution scores across samples and resamples. This gap motivates the framework proposed in this paper.
3
Methodology
To account for the stochastic nature of feature attribution, we adopt a distributional view of feature attributions and introduce a new metric called RoSHAP. Since most additive attribution methods result from complex calculation procedure, the underlying distribution is not available in closed-form. We propose to use bootstrapping to approximate these distributions, and then estimate the resulting empirical distribution via kernel density estimation (KDE). Under mild regularity conditions, the attribution distribution can be approximated by a Gaussian distribution, reducing the computational burden associated with bootstrapping. RoSHAP then quantifies feature attribution by jointly accounting for inactivity, signal, and noise level. The overall pipeline is illustrated in Figure 2.
Figure 2: Overall framework for distributional feature attribution estimation.
3
Bootstrap-based Estimation of Feature At- Algorithm 1. Bootstrap-based SHAP distribution estimation tributions. Let Dn = {(x(i) , yi ) : i = 1, . . . , n} denote a labeled dataset with n Require: Data D; bootstrap runs B Ensure: Sample- and feature-level SHAP distributions observations and p features, where x(i) = 1: Preprocess D; optionally screen predictors ⊤ p (xi1 , . . . , xip ) ∈ R and yi is the response. 2: for b = 1, . . . , B do (b) Draw Dtrain from D Let fˆ(·) denote the fitted model, with prediction 3: (b) (b) Define the out-of-bag set Doob = D \ Dtrain fˆ(x(i) ) for observation i. For each observation 4: (b) Fit model fˆ(b) on Dtrain i and feature j, denote by Tij the contribution 5: (b) (i) 6: Compute OOB SHAP Tij of feature j at x . 7: end for
We estimate the distribution of attribution values 8: Summarize Uj(b) by Pi |Tij | by bootstrapping. By repeatedly resampling the 9: Estimate distribution for Tij and Uj training data, refitting the model, and recomputing attributions, we obtain multiple realizations of Tij , which serve as samples from the underlying distribution. Algorithm 1 summarizes the bootstrap-based distribution estimation procedure. Zero-Inflated Modeling and Gaussian Approximation of Feature-level Attributions. Our framework applies to any additive attribution method, and we use SHAP as a representative example throughout. In practice, SHAP values often contain exact zeros, especially in tree-based models, where features absent from a prediction path receive no attribution [1, 11]. To account for this behavior, we model the sample-feature SHAP value Tij as a zero-inflated mixture, Tij ∼ wij δ0 + (1 − wij ) Gij , where δ0 is a point mass at zero, wij = Pr(Tij = 0), and Gij denotes the nonzero component of Tij . We define the feature-level attribution as the aggregated magnitude Uj =
n X
|Tij | =
i=1
n X (1 − Bij )Hij . i=1
For asymptotic analysis, it is convenient to decompose |Tij | = (1 − Bij )Hij ,
Bij ∼ Bernoulli(wij ),
where Bij = 1 indicates Tij = 0, and Hij = |Gij | is the nonzero magnitude. This representation separates sparsity from magnitude while preserving the distribution of |Tij |. We estimate the distribution of Uj nonparametrically using KDE. When mild regularity conditions hold, the distribution of Uj can be well approximated by a Gaussian distribution and can therefore be summarized by its mean and standard deviation. This approximation is justified through a Lyapunov central limit theorem argument. The mean and variance of Uj are given by µj := E(Uj ) = s2j := Var(Uj ) =
n X
n X (1 − wij )E(Hij ), i=1
(1 − wij ) Var(Hij ) +
i=1
n X
wij (1 − wij ){E(Hij )}2 .
i=1
Assumption. For a fixed feature j, the pairs (Bij , Hij ), i = 1, . . . , n, are mutually independent, with Bij ⊥ Hij for each i. In addition, for some δ > 0, the Lyapunov condition holds: n 1 X
s2+δ j i=1
2+δ
E ||Tij | − E[|Tij |]|
→ 0.
This assumption ensures that the aggregate contribution Uj is governed by many small terms rather than a few dominant ones, which is crucial for the Gaussian approximation. Theorem. Under the above assumption, Uj − µj d − → N (0, 1). sj 4
This result justifies approximating the distribution of Uj by a Gaussian when n → ∞. In practice, it implies that for features with sufficiently stable aggregated contributions, uncertainty can be effectively summarized using only the mean and variance, substantially reducing the number of bootstrap resamples. A proof is provided in Appendix A.1. When the regularity conditions are approximately satisfied, we summarize Uj using its empirical median and standard deviation. Robust metric. Based on the estimated feature-level attribution distribution, we define several complementary metrics for feature evaluation. For feature j, define mj := Median(Uj | Uj > 0),
P0j := Pr(Uj = 0),
sj := SD(Uj ).
The robust feature attribution score, RoSHAP, is defined as RoSHAPj := (1 − P0j )
m2j . sj
Here, P0j measures inactivity, while RoSHAPj rewards features that are active, strong, and stable.
4
Simulation
This simulation study demonstrates that the proposed method provides more robust and interpretable feature contribution estimates. From a distributional perspective, the method helps identify stable signals and distinguish them from noise features selected by chance. The results verify that the feature contribution distributions are approximately Gaussian, allowing fewer bootstrap repetitions and substantially reducing computational cost. Simulation setup. A binary classification setting was considered with n observations and d predictors, where the first s predictors were true signals and the remaining d − s predictors were noise, with s < d. To mimic settings in which a feature is informative at the population level but may contribute little for some individual observations, predictors were generated from zero-inflated Gaussian distributions. For each signal feature j = 1, . . . , s, 2 Xj | Y = y ∼ πsignal δ0 + (1 − πsignal )N (µjy , σsignal ),
where Y ∈ {0, 1}, P (Y = 1) = 0.5, and µj1 = −µj0 . Thus, signal features differed between the two classes through their nonzero component means. For noise features j = s + 1, . . . , d, the class-conditional distributions had the same zero-inflated Gaussian form but with zero mean in both classes: 2 Xj | Y = y ∼ πnoise δ0 + (1 − πnoise )N (0, σnoise ).
Therefore, noise features had no class-dependent mean shift. Under this design, signal features may be locally inactive when they take the value zero, while remaining globally informative through their class-dependent nonzero means. In this experiment, n, d, and s were set to 600, 1000, and 10, respectively. For signal features, σsignal = 3 and πsignal = 0.3. For noise features, σnoise = 1 and πnoise = 0.2. The signal strengths were defined by µj1 = −µj0 , with µj1 evenly spaced from 1.5 to 0.4 across the ten signal features. An XGBoost classifier was then fitted using the same model parameters across all comparison settings. Results. Figure 3 shows the empirical distributions of the importance estimates over 1000 bootstrap runs. Among the selected predictors, the distributions for x1 and x3 are approximately Gaussian, while the distribution for x7 is slightly right-skewed. Overall, the nonzero portions of the distributions are close to Gaussian, supporting the use of Gaussian-based uncertainty summaries with fewer bootstrap repetitions. Three settings were compared: 1000-run bootstrap, 10-run bootstrap, and standard single train–test split SHAP. As shown in Table 1, the 1000-run bootstrap provides the most reliable ranking: eight signal features are recovered within the top 12, and the lowest-ranked signal feature, x9 , is placed 5
(a) x1
(b) x3
(c) x7
Figure 3: Empirical distributions of robust importance estimates over 1000 bootstrap runs with fitted Gaussian distributions for predictors x1 , x3 , and x7 . 85th among 1000 features based on the robust importance metric. The 10-run bootstrap also performs well, recovering eight signal features among the top 12, although the lowest-ranked signal feature drops to 317th among 1000 features. The zero percentage values P0 provide additional insight into stability. Strong signal features have near-zero P0 , indicating consistent contribution across bootstrap runs. In contrast, weaker or less stable features, such as x8 and x9 , exhibit larger P0 , suggesting more frequent zero contributions. Table 1: Comparison of feature rankings from bootstrap-based robust importance, single-round SHAP, information gain, XGBoost gain, and LIME screening metrics. 1000 Runs 10 Runs XGBoost SHAP Information Gain XGBoost Gain LIME Rank ID RoSHAP P0 Rank ID RoSHAP P0 Rank ID Score Rank ID Score Rank ID Score Rank ID Score 1 x1 512.66 0.00 1 x4 779.54 0.00 1 x2 35.24 1 x4 0.090 5 x1 1.826 1 x3 0.077 2 x3 347.21 0.00 2 x1 767.57 0.00 2 x6 25.86 48 x2 0.043 6 x2 1.782 2 x1 0.075 3 x4 339.45 0.00 3 x2 433.42 0.00 3 x4 15.38 187 x5 0.024 10 x3 1.696 3 x2 0.055 4 x2 270.75 0.00 4 x3 430.41 0.00 4 x3 15.14 235 x6 0.021 13 x4 1.450 4 x4 0.036 x6 271.18 0.00 12 x5 4.12 259 x1 0.019 26 x6 1.069 5 x6 0.035 5 x6 120.99 0.20 5 6 x5 80.65 0.40 6 x5 104.68 0.00 31 x10 0.14 281 x3 0.017 41 x7 0.887 7 x7 0.008 9 x7 14.73 7.20 7 x7 60.94 0.00 333 x7 0.00 292 x7 0.015 74 x10 0.686 15 x10 0.004 12 x10 10.65 11.00 10 x10 25.85 0.00 387 x8 0.00 329 x9 0.012 131 x5 0.555 58 x9 0.000 53 x8 1.92 27.60 95 x9 2.64 30.00 669 x1 0.00 658 x10 0.000 159 x8 0.494 960 x8 0.000 85 x9 1.08 37.00 317 x8 0.32 40.00 670 x9 0.00 744 x8 0.000 707 x9 0.063 999 x5 0.000
The other metric are substantially less stable. Only four signal features are recovered within the top 12, while four signal features are inactive, including x1 , which has the strongest signal by construction. This highlights that a single SHAP run may fail to identify important predictors, whereas the bootstrap-based robust importance more effectively captures stable signal features. The alternative metrics are substantially less effective. Single-round SHAP identifies only four signal features within the top 12 and assigns zero or near-zero importance to several true signals, including x1 , the strongest signal by construction. Information gain and gain-based importance also show unreliable rankings, placing several true signals far below noise variables. Despite recovering the first seven signal features reasonably well, LIME treats the remaining three signals as inactive. In contrast, RoSHAP identifies the signal features across repeated model fits more consistently. Computational comparison. TreeSHAP has computational complexity approximately O(T LD2 ), where T is the number of trees, L is the maximum number of leaves per tree, and D is the maximum tree depth [11]. With B bootstrap repetitions, the total computation is roughly B O(T LD2 ), so reducing B from 1000 to 10 reduces the SHAP computation by about 100 times. Each bootstrap sample covers about 1 − (1 − 1/n)n ≈ 63.2% of the original observations. Across B runs, the probability that an observation is selected at least once is approximately 1 − 0.368B , which is nearly 1 when B = 10. Thus, fewer bootstrap runs can still provide broad data coverage and useful Gaussian-based uncertainty summaries with much lower runtime.
5
Experiments
Four experiments demonstrate the proposed framework across data types. The Golub dataset [22] represents high-dimensional small-sample classification, Musk (Version 2) [23] represents moderatesized high-dimensional molecular classification, UJIIndoorLoc [24] represents high-dimensional regression, and CIFAR-10 represents 10-class image classification. 6
For tabular data, feature-level interpretability is evaluated across neural networks, tree-based models including LightGBM, XGBoost, random forest (RF), and CatBoost, and logistic regression (LR). Because the true important features are unknown, the proposed feature-selection results are compared with information gain (IG), model-based gain or LR coefficients, LIME, and SHAP. Classification performance is evaluated by accuracy, F1 score, average precision (AP), and AUC-ROC; regression performance is evaluated by root mean squared error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). For each method, the top k selected features are used to refit the same predictive model with fixed parameters. All comparisons use the same training–testing split, cross-validation procedure, and held-out test data, so performance differences reflect featuresubset quality rather than model or data-split differences. The main text reports results from one representative model, with other models provided in the appendix. For image data, the framework is applied at the sample level. In the CIFAR-10 experiment, a ViT-base classifier with patch-level attribution is used to identify image regions contributing most to the predicted class. 5.1
Classification: Golub
The Golub dataset [22] contains 72 samples and 7,129 gene expression probes. The response of interest is leukemia subtype: Acute Lymphoblastic Leukemia (ALL), including both allB and allT, versus Acute Myeloblastic Leukemia (AML). The dataset contains 47 ALL samples and 25 AML samples. This experiment studies gene-level feature attribution using an XGBoost classifier. The goal is to identify not only which genes contribute most strongly to prediction, but also which genes contribute most consistently across repeated model fits. The Golub dataset provides a representative high-dimensional, small-sample setting where interpretation is difficult and attribution instability is common. Figure 4 shows the distributions of the top 10 genes ranked by RoSHAP metric over 500 bootstrap runs. After the eighth-ranked gene, most genes show little or no contribution to cancer prediction. Compared with the seed-specific SHAP results in Figure 1, the proposed method provides a clearer summary of which genes are consistently important under the fitted XGBoost model. In this setting, a small number of observations can dominate the variance of, so the Lyapunov condition fails to hold. The distributions are therefore estimated nonparametrically using KDE. Figure 5 compares feature-selection performance across methods. Since experimental validation of candidate genes is often costly, Figure 4: Distribution of top 10 the evaluation focuses on the top 1–15 selected features from each XGBoost feature contribution esmethod. For each metric, bars show mean performance across timates for the Golub dataset. candidate feature-set sizes, and error bars show the corresponding variability. The proposed method, shown in purple, achieves strong accuracy and F1 score. RoSHAP also maintains comparable AP and AUC-ROC, suggesting stable and reliable feature selection.
Figure 5: Feature-selection performance on the Golub data using XGBoost. Methods are SHAP (blue), gain (orange), LIME (green), IG (red), and RoSHAP (purple). Bars show mean performance across candidate feature-set sizes; error bars show variation. Higher is better. 5.2 Classification: Musk (Version 2) The Musk (Version 2) dataset [23] contains 6,598 instances and 166 continuous features. The binary classification task is to distinguish musk from non-musk molecules, where each instance represents one molecular conformation described by physicochemical descriptors. The dataset is imbalanced, with 5,581 non-musk and 1,017 musk instances, making it useful for evaluating feature contribution robustness in high-dimensional molecular classification. 7
The Musk experiment studies feature-level attribution using a CatBoost classifier. Figure 6 shows the distributions of the top 15 CatBoost feature contribution estimates for the Musk dataset. Several features have similar contribution magnitudes and variability, indicating that the predictive signal is distributed across multiple features rather than dominated by a single feature. For variable selection, only the top 30 selected features are compared for Musk dataset, since prediction performance remains high beyond this point. Figure 7 shows the feature selection performance for the Musk dataset. Performance across selection metrics is largely comparable. The feature contribution distribution suggests that many features have similar levels of importance. As a result, differences in the selected feature subsets have limited Figure 6: Distribution of top 15 CatBoost feature attribution estiimpact on prediction performance. mates for the Musk dataset.
Figure 7: Feature-selection performance on the Musk (Version 2) data using CatBoost. Methods are SHAP (blue), gain (orange), LIME (green), IG (red), and RoSHAP (purple). Bars show mean performance across candidate feature-set sizes; error bars show variation. Higher is better. 5.3
Regression: UJIIndoorLoc
The UJIIndoorLoc dataset is used to evaluate the proposed framework in a regression task. It is a public indoor-localization benchmark collected at Universitat Jaume I for evaluating WLAN/WiFi fingerprint-based positioning systems [24]. The dataset contains 21,048 observations, including 19,937 training/reference records and 1,111 validation/test records. Each observation includes received signal strength intensity (RSSI) measurements from 520 wireless access points (WAPs), along with location and metadata variables. The RSSI values range from approximately −104 to 0 dBm, with 100 indicating an undetected access point. In this study, the 520 WAP signal-strength variables are used as predictors and longitude is used as the continuous response. The value 100 is recoded as −110 to represent a very weak or undetected signal. This experiment studies feature-level attribution using a LightGBM regression model. Figure 8 shows the distributions of the top 15 LightGBM feature attribution estimates for the UJIIndoorLoc dataset. The top two features have similar attribution patterns, while the remaining features exhibit different levels of variability: some distributions are more concentrated, whereas others are more dispersed. Figure 9 compares the feature-selection performance of different importance measures on the UJIIndoorLoc regression task using LightGBM. Overall, RoSHAP achieves competitive and stable performance across the evaluation metrics. 5.4
Classification: CIFAR-10
This experiment extends the proposed framework beyond tabular Figure 8: Distribution of top data to sample-level image interpretation. The CIFAR-10 dataset 15 LightGBM feature attribution contains 60,000 32 × 32 color images from 10 classes, with estimates for the UJIIndoorLoc 50,000 training and 10,000 test images [25]. We conduct 50 dataset. bootstrap-based SHAP analyses using a ViT-base model. For each test image, the SHAP array has dimension 50 × p × 10, where p = 24 × 24 is the number 8
Figure 9: Feature-selection performance on the UJIIndoorLoc data using LightGBM. Methods are SHAP (blue), gain (orange), LIME (green), IG (red), and the RoSHAP (purple). Bars show mean performance across candidate feature-set sizes; error bars show variation. Lower is better. of image patches and 10 is the number of classes. Results focus on the SHAP distribution for the predicted class. As an illustrative example, one ship image that was correctly predicted by most bootstrap-fitted models is selected. Figure 10 shows patch-level SHAP maps from three randomly selected bootstrap runs, where red indicates positive SHAP values, blue negative values, and white values near zero. Although the prediction is accurate, the contributing patches vary noticeably across runs, suggesting that single-run image explanations are sensitive to training randomness. For example, in Round 1 at For example, the largest SHAP value appears on the water region in Round 1 but shifts toward the ship body in Round 3. After aggregation, the robust attribution map in Figure 11a shows that the prediction is mainly attributed to patches on the ship body. The proposed framework provides additional insight into attribution patterns behind image predictions. Figure 11 illustrates two representative cases: a correct prediction with attribution focused on the target object, and an incorrect prediction in which attribution is concentrated on irrelevant regions. From Figure 10: CIFAR-10 ship example left to right, each row shows the original image, the median using 3 bootstrap rounds. SHAP map for the predicted class, the standard deviation of SHAP values across bootstrap runs, and the robust SHAP map. In the median SHAP map, red indicates positive contribution, blue negative contribution, and white near-zero contribution; in the standard deviation map, red indicates larger variation and white little variation; in the robust SHAP map, green indicates larger contribution and white no contribution. In Figure 11a, the model focuses mainly on the correct object region. In Figure 11b, however, the model focuses largely on the sky rather than the car, leading to an incorrect bird prediction. These examples show that the framework can help diagnose prediction errors by revealing whether the model attends to the appropriate image regions.
6
Discussion
This study addresses a key limitation of existing feature attribution methods: attribution scores are often treated as fixed quantities despite their stochastic nature. We propose a distributional framework for characterizing attribution uncertainty and introduce RoSHAP, a robust SHAP-based metric for stable feature ranking and contribution evaluation. RoSHAP prioritizes features that are consistently active, strong, and stable across repeated analyses. We extensively test our model on realistic datasets and demonstrate that the proposed robust metric improves interpretation stability, helps distinguish consistent signals from chance-selected features, and can be applied across different data types. Although the bootstrap procedure remains more costly than single-run attribution, the Gaussian approximation helps extensively reduce this burden when applicable. Future work can extend the framework from data-sampling randomness to model-level randomness. Overall, RoSHAP serves as a useful complement to existing feature attribution methods by supporting more stable and reliable interpretation.
9
(a) Correct prediction; attribution focuses on the ship.
(b) Wrong prediction; attribution does not focus on the car.
Figure 11: CIFAR-10 explanations after 50 ViT-base bootstrap runs. Columns show the original image, median SHAP map, SHAP standard deviation, and robust SHAP map. Red/blue denote positive/negative SHAP values, red in the standard deviation map denotes larger variability, and green denotes larger robust attribution.
10
References [1] Scott M Lundberg and Su-In Lee. A unified approach to interpreting model predictions. Advances in neural information processing systems, 30, 2017. [2] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. " why should i trust you?" explaining the predictions of any classifier. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 1135–1144, 2016. [3] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje. Learning important features through propagating activation differences. In International conference on machine learning, pages 3145–3153. PMlR, 2017. [4] Puskar Bhattarai, Deepa Singh Thakuri, Yuzheng Nie, and Ganesh B Chand. Explainable ai-based deepshap for mapping the multivariate relationships between regional neuroimaging biomarkers and cognition. European journal of radiology, 174:111403, 2024. [5] Ana Victoria Ponce-Bobadilla, Vanessa Schmitt, Corinna S Maier, Sven Mensing, and Sven Stodtmann. Practical guide to shap analysis: Explaining supervised machine learning model predictions in drug development. Clinical and translational science, 17(11):e70056, 2024. [6] Viswan Vimbi, Noushath Shaffi, and Mufti Mahmud. Interpreting artificial intelligence models: a systematic review on the application of lime and shap in alzheimer’s disease detection. Brain informatics, 11(1):10, 2024. [7] Yuchun Liu, Zhihui Liu, Xue Luo, and Hongjingtian Zhao. Diagnosis of parkinson’s disease based on shap value feature selection. Biocybernetics and Biomedical Engineering, 42(3):856–869, 2022. [8] Edoardo Mosca, Ferenc Szigeti, Stella Tragianni, Daniel Gallagher, and Georg Groh. Shap-based explanation methods: a review for nlp interpretability. In Proceedings of the 29th international conference on computational linguistics, pages 4593–4603, 2022. [9] Ziqi Li. Extracting spatial effects from machine learning model using local interpretation method: An example of shap and xgboost. Computers, Environment and Urban Systems, 96:101845, 2022. [10] Jeremy Goldwasser and Giles Hooker. Statistical significance of feature importance rankings. arXiv preprint arXiv:2401.15800, 2024. [11] Scott M Lundberg, Gabriel Erion, Hugh Chen, Alex DeGrave, Jordan M Prutkin, Bala Nair, Ronit Katz, Jonathan Himmelfarb, Nisha Bansal, and Su-In Lee. From local explanations to global understanding with explainable ai for trees. Nature machine intelligence, 2(1):56–67, 2020. [12] Tianqi Chen and Carlos Guestrin. Xgboost: A scalable tree boosting system. In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pages 785–794, 2016. [13] Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. Lightgbm: A highly efficient gradient boosting decision tree. Advances in neural information processing systems, 30, 2017. [14] Torgyn Shaikhina, Umang Bhatt, Roxanne Zhang, Konstantinos Georgatzis, Alice Xiang, and Adrian Weller. Effects of uncertainty on the quality of feature importance explanations. In AAAI workshop on explainable agency in artificial intelligence. AAAI Press Washington, DC, USA, 2021. [15] Riccardo Scheda and Stefano Diciotti. Explanations of machine learning models in repeated nested cross-validation: an application in age prediction using brain complexity features. Applied Sciences, 12 (13):6681, 2022. [16] Ramzi Halabi, Benoit H Mulsant, Mirkamal Tolend, Daniel M Blumberger, Alexandra DeShaw, Arend Hintze, Christina Gonzalez-Torres, Muhammad I Husain, Helena K Kim, Claire O’Donovan, et al. A systematic exploration of digital biomarkers for the detection of depressive episodes in bipolar disorder. npj Mental Health Research, 5(1):13, 2026. [17] David Watson, Joshua O’Hara, Niek Tax, Richard Mudd, and Ido Guy. Explaining predictive uncertainty with information theoretic shapley values. Advances in Neural Information Processing Systems, 36: 7330–7350, 2023. [18] Wilson E Marcílio and Danilo M Eler. From explanations to feature selection: assessing shap values as feature selection mechanism. In 2020 33rd SIBGRAPI conference on Graphics, Patterns and Images (SIBGRAPI), pages 340–347. Ieee, 2020.
11
[19] Yibrah Gebreyesus, Damian Dalton, Sebastian Nixon, Davide De Chiara, and Marta Chinnici. Machine learning for data center optimizations: feature selection using shapley additive explanation (shap). Future Internet, 15(3):88, 2023. [20] Anruo Shen, Jingnan Sun, Xiaogang Chen, and Xiaorong Gao. A data-centric and interpretable eeg framework for depression severity grading using shap-based insights. Journal of NeuroEngineering and Rehabilitation, 22(1):116, 2025. [21] Huanjing Wang, Qianxin Liang, John T Hancock, and Taghi M Khoshgoftaar. Feature selection strategies: a comparative analysis of shap-value and importance-based methods. Journal of Big Data, 11(1):44, 2024. [22] Todd R Golub, Donna K Slonim, Pablo Tamayo, Christine Huard, Michelle Gaasenbeek, Jill P Mesirov, Hilary Coller, Mignon L Loh, James R Downing, Mark A Caligiuri, et al. Molecular classification of cancer: class discovery and class prediction by gene expression monitoring. science, 286(5439):531–537, 1999. [23] David Chapman and Ajay Jain. Musk (Version 2). UCI Machine Learning Repository, 1994. [24] Joaquín Torres-Sospedra, Raúl Montoliu, Adolfo Martínez-Usó, Tomàs Arnau, and Joan P. Avariento. UJIIndoorLoc. UCI Machine Learning Repository, 2014. Accessed from the UCI Machine Learning Repository. [25] Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009. [26] Patrick Billingsley. Probability and Measure. Wiley, New York, 3 edition, 1995.
A
Appendix / supplemental material
A.1
Proof of Theorem
Lemma 1 (Lyapunov central limit theorem [26, Theorem 27.3]). Pn Let 2Zn1 , . . . , Znn be independent 2 random variables with E[Zni ] = 0, variances σni , and s2n = i=1 σni . If, for some δ > 0, n 1 X
s2+δ n i=1 then
E|Zni |2+δ → 0,
Pn
i=1 Zni d
sn
− → N (0, 1).
Proof. Since |Tij | = (1 − Bij )Hij , the independence of Bij and Hij gives E(|Tij |) = E(1 − Bij )E(Hij ) = (1 − wij )E(Hij ). Therefore, µj = E(Uj ) =
n X
n X E(|Tij |) = (1 − wij )E(Hij ).
i=1
i=1
Next, by the law of total variance, Var(|Tij |) = E(Var(|Tij | | Bij )) + Var(E(|Tij | | Bij )). Conditionally on Bij , E(|Tij | | Bij ) = (1 − Bij )E(Hij ),
Var(|Tij | | Bij ) = (1 − Bij ) Var(Hij ).
Hence
Var(|Tij |) = (1 − wij ) Var(Hij ) + wij (1 − wij ){E(Hij )}2 . Since the pairs (Bij , Hij ) are mutually independent, the variables |Tij | are mutually independent, and therefore n n n X X X 2 sj = Var(Uj ) = Var(|Tij |) = (1 − wij ) Var(Hij ) + wij (1 − wij ){E(Hij )}2 . i=1
i=1
i=1
12
Now define the centered summands Zni := |Tij | − E(|Tij |). By the stated Lyapunov condition with δ > 0, n 1 X 2+δ E ||Tij | − E(|Tij |)| → 0. s2+δ j i=1 Thus Lyapunov’s central limit theorem applies, giving Pn i=1 {|Tij | − E(|Tij |)} d − → N (0, 1). sj Pn Pn Since Uj = i=1 |Tij | and µj = i=1 E(|Tij |), this is exactly Uj − µj d − → N (0, 1). sj Remark 1 (On the Lyapunov condition). A natural choice is δ = 1, which reduces the condition to n 1 X 3 E |Tij | − E(|Tij |) → 0. 3 sj i=1 A simple sufficient condition is that the third moments are uniformly bounded, supi E|Tij |3 < ∞, and no single term dominates the variance, i.e., max Var(|Tij |) = o(s2j ). 1≤i≤n
Under these conditions, the numerator grows at most linearly in n, while s2j grows proportionally to the number of contributing terms, yielding n → 0, s3j and hence the Lyapunov condition holds. A.2
Additional Golub Data Results
We employed the same six classifiers as in the previous experiments, with hyperparameters re-tuned to suit the characteristics of the Golub dataset. The tree-based ensemble models were configured as follows: LightGBM (binary objective, binary logloss metric, learning rate 0.1, max depth 8, 31 leaves, minimum 15 samples per leaf, 80% feature and bagging fractions with bagging frequency of 1, 150 boosting rounds); XGBoost (binary logistic objective, logloss evaluation metric, learning rate 0.1, max depth 6, 80% subsample and 60% column sampling, 100 boosting rounds); and CatBoost (logloss metric, learning rate 0.1, depth 6, default iterations); Random Forest used 100 trees with max depth 6, while Gradient Boosting used 200 trees with max depth 3 at the default learning rate. As a linear baseline, Logistic Regression was trained with a maximum of 200 iterations using default regularization. All models were initialized with a fixed random seed to ensure reproducibility, and feature contribution was extracted using tree-based explainers for the ensemble models and a linear (coefficient-based) explainer for Logistic Regression. The XGBoost results are presented in the main text, while the feature selection performance for the remaining classifiers on the Golub dataset is shown in Figures 12, 13, 14, 15, and 16.
Figure 12: Golub gene selection performance comparison using CatBoost. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Higher better. 13
Figure 13: Golub gene selection performance comparison using Gradient Boosting. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Higher better.
Figure 14: Golub gene selection performance comparison using LightGBM. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Higher better.
Figure 15: Golub gene selection performance comparison using Logistic Regression. Feature selection methods include SHAP (blue), |coefficient| (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Higher better.
Figure 16: Golub gene selection performance comparison using Random Forest. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Higher better. A.3
Additional Musk (Version 2) Data Results
The same six classifiers in the previous experiments were used, with hyperparameters re-tuned to suit the characteristics of the Musk (Version 2) dataset. The tree-based ensemble models were configured as follows: LightGBM (binary objective, binary logloss metric, learning rate 0.1, max depth 8, 31 leaves, minimum 15 samples per leaf, 80% feature and bagging fractions with bagging frequency of 1, 150 boosting rounds); XGBoost (binary logistic objective, logloss evaluation metric, learning rate 0.1, max depth 6, 80% subsample and 60% column sampling, 100 boosting rounds); and CatBoost (logloss metric, learning rate 0.1, depth 6, default iterations); Random Forest used 100 trees with max depth 6, while Gradient Boosting used 200 trees with max depth 3 at the default learning rate. Logistic Regression was trained with a maximum of 200 iterations using default regularization. All models were initialized with a fixed random seed to ensure reproducibility, and feature contribution was extracted using tree-based explainers for the ensemble models and a linear (coefficient-based) 14
explainer for Logistic Regression. The CatBoost results are presented in the main text, while the feature selection performance for the remaining classifiers on the Musk Version 2 dataset is shown in Figures 17, 18, 19, 20, and 21.
Figure 17: Musk (Version 2) variable selection performance comparison using XGBoost. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple).
Figure 18: Musk (Version 2) variable selection performance comparison using Gradient Boosting. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple).
Figure 19: Musk (Version 2) variable selection performance comparison using LightGBM. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple).
Figure 20: Musk (Version 2) variable selection performance comparison using Logistic Regression. Feature selection methods include SHAP (blue), |coefficient| (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple).
Figure 21: Musk (Version 2) variable selection performance comparison using Random Forest. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). 15
A.4
Additional UJIIndoorLoc Results
For the UJIIndoorLoc regression experiment, except for the random forest model, same predictive models used in the previous experiments were considered, with hyperparameters re-tuned to suit the characteristics of the UJIIndoorLoc dataset. The response variable was longitude. The tree-based ensemble models were configured as follows: LightGBM used a regression objective with RMSE metric, learning rate 0.03, max depth 7, 63 leaves, minimum 30 samples per leaf, 70% feature fraction, 80% bagging fraction with bagging frequency of 1, ℓ1 regularization 0.01, ℓ2 regularization 1.0, and 300 boosting rounds; XGBoost used a squared-error regression objective with RMSE evaluation metric, learning rate 0.03, max depth 6, 80% subsample ratio, 70% column sampling, minimum child weight 10, ℓ1 regularization 0.01, ℓ2 regularization 2.0, histogram-based tree construction, and 300 boosting rounds; and CatBoost used RMSE loss and evaluation metric, learning rate 0.03, depth 6, 300 iterations, ℓ2 leaf regularization 5.0, and 70% random subspace sampling. Ridge Regression was included as a linear baseline with regularization parameter α = 10. All models were initialized with a fixed random seed to ensure reproducibility, and feature contribution was extracted using tree-based explainers for the ensemble models and a linear explainer for Ridge Regression. The LightGBM results are presented in the main text, while the feature selection performance for the remaining models on the UJIIndoorLoc dataset is shown in Figures 22, 23, 24, and 25.
Figure 22: UJIIndoorLoc variable selection performance comparison using CatBoost. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Bar heights represent the mean performance metric across different candidate feature-set sizes, and error bars indicate the corresponding variation. Lower better.
Figure 23: UJIIndoorLoc variable selection performance comparison using Gradient Boosting. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Bar heights represent the mean performance metric across different candidate feature-set sizes, and error bars indicate the corresponding variation. Lower better.
Figure 24: UJIIndoorLoc variable selection performance comparison using XGBoost. Feature selection methods include SHAP (blue), gain (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Bar heights represent the mean performance metric across different candidate feature-set sizes, and error bars indicate the corresponding variation. Lower better. 16
Figure 25: UJIIndoorLoc variable selection performance comparison using Linear Regressio (Ridge). Feature selection methods include SHAP (blue), |coef| (orange), LIME (green), information gain (IG; red), and the proposed robust method (purple). Bar heights represent the mean performance metric across different candidate feature-set sizes, and error bars indicate the corresponding variation. Lower better.
A.5
Additional CIFAR 10 Data Results
This experiment employ ViT-Base/16, pretrained on ImageNet-21k. The model takes 384×384 inputs divided into 16×16 patches (24×24 grid, 576 patches total), with 12 transformer layers, 12 attention heads, and a hidden dimension of 768 ( 86M parameters). The classification head is replaced with a linear layer mapping to 10 CIFAR-10 classes and fine-tuned for 10 epochs using Adam (lr=1e-4, cosine decay) with batch size 128.
Figure 26: Prediction and bootstrap results for image 8.
Figure 27: Prediction and bootstrap results for image 11.
Figure 28: Prediction and bootstrap results for image 20. 17
Figure 29: Prediction and bootstrap results for image 36.
Figure 30: Prediction and bootstrap results for image 41.
Figure 31: Prediction and bootstrap results for image 66.
Figure 32: Prediction and bootstrap results for image 85.
18