MinShap: A Modified Shapley Value Approach for Feature Selection Chenghui Zheng ∗
Garvesh Raskutti ∗
arXiv:2604.15107v1 [stat.ML] 16 Apr 2026
Abstract Feature selection is a classical problem in statistics and machine learning, and it continues to remain an extremely challenging problem especially in the context of unknown non-linear relationships with dependent features. On the other hand, Shapley values are a classic solution concept from cooperative game theory that is widely used for feature attribution in general non-linear models with highly-dependent features. However, Shapley values are not naturally suited for feature selection since they tend to capture both direct effects from each feature to the response and indirect effects through other features. In this paper, we combine the advantages of Shapley values and adapt them to feature selection by proposing MinShap, a modification of the Shapley value framework along with a suite of other related algorithms. In particular for MinShap, instead of taking the average marginal contributions over permutations of features, considers the minimum marginal contribution across permutations. We provide a theoretical foundation motivated by the faithfulness assumption in DAG (directed acyclic graphical models), a guarantee for the Type I error of MinShap, and show through numerical simulations and real data experiments that MinShap tends to outperform state-of-the-art feature selection algorithms such as LOCO, GCM and Lasso in terms of both accuracy and stability. We also introduce a suite of algorithms related to MinShap by using the multiple testing/p-value perspective that improves performance in lower-sample settings and provide supporting theoretical guarantees. Keywords: conditional independence; feature importance; feature selection; Shapley value.
1
Introduction
Shapley values are a feature attribution method widely used to interpret black-box prediction methods [43, 16, 7, 33]. One of the drawbacks of Shapley values and other feature attribution methods is that unlike feature importance methods, Shapley values are not amenable to feature selection (see e.g. [22, 14]) as they incorporate both direct and indirect effects. From an interpretability perspective, this makes them slightly limited. On the other hand, existing feature importance and selection methods are known to often be inaccurate and unstable when features are highly dependent (see e.g. [50, 19, 49]), especially for non-linear models. One of the arguments for using Shapley values over standard feature importance methods is that they provide a more meaningful interpretation when features are correlated. In this paper, we adapt the existing Shapley value framework which may be interpreted as averaging feature importance over models under different permutations [6] by introducing the MinShap algorithm which simply replaces the average feature importance by the minimum feature importance. By exploiting connections with the faithfulness assumption for directed acyclic graphical (DAG) models, our MinShap algorithm naturally performs feature selection. Further, since MinShap incorporates the advantages of Shapley values which combine multiple feature importances, as we display through our simulation study, MinShap tends to significantly out-perform existing feature selection algorithms in terms of both accuracy and stability. Computationally, our MinShap algorithm is as costly as existing approximate Shapley value approaches (see e.g. [6]). ∗
Department of Statistics, University of Wisconsin - Madison, Madison, WI 53706
1
Furthermore, by exploiting asymptotic normality of feature importance scores [52], we connect the MinShap algorithm to its parallel Max-p value test, and further extend to partial conjunction hypothesis testing [2] to offer higher power in dense and finite-sample regime. We also empirically validate our algorithms through simulation studies and real data analysis.
1.1
Our Contributions
The main contributions of our work are summarized as follows: • Develop the MinShap algorithm based on the causal faithfulness condition for DAG models. • Provide theoretical guarantees for Type I error control for our MinShap algorithm. • Connect MinShap to multiple hypothesis testing and p-value approaches that lead to a suite of methods that are more suitable in a lower sample size, higher noise settings. • Demonstrate through simulation and real data the superior performance of MinShap and related p-value algorithms for permutation approaches with comparison to state-of-the-art methods such as GCM, LOCO and Lasso.
1.2
Related Work
Feature selection Feature selection is widely used in many machine learning settings to enhance prediction performance, improve interpretability and reduce computation. It refers to the procedure of selecting relevant input features to predict target values in a given dataset. There are three main categories: filter-based, wrapper-based, and embedded methods [18]. Filter methods select the variable of interest in a pre-processing step, which is independent of any machine learning algorithm. They are computationally cheap, but they can underperform model-based selection when features are correlated or have interactions. Wrapper methods are model-agnostic: they treat the learning algorithm as a black box and perform variable selection by searching over feature subsets, scoring/ranking each subset by predictive performance. There are several global model-agnostic interpretability methods assigning an importance score for each predictor. The most prominent one is leave-one-covariate-out (LOCO, [23]). It assesses a feature’s importance by comparing the prediction error of a full model to that of a retraining model with the targeted feature removed. The larger difference indicate the greater importance of the feature. Another recently popular method related to regression-based conditional independence is generalized covariance measure (GCM, [38]). It assesses a feature’s relevance by testing whether the residual of the target is conditionally uncorrelated with the residual of that feature after regressing both on the remaining covariates. However, both methods struggled with highly dependent and non-additive non-linear relationship [27, 49]. Embedded methods integrate feature selection into model training by directly computing feature importance score from the fitted model. Lasso [25] is a common example which is only suitable for linear model. Embedded methods are generally more efficient than wrappers but are not model-agnostic. Shapley values and feature selection The Shapley value, originally derived from cooperative game theory [6], offers a principled framework for attributing a model’s performance to its input features. Shapley value has been studied in terms of features importance in different contexts [51, 47, 10, 34]. A popular example where Shapley values are applied to feature selection is to rank features by their mean absolute SHAP values [24] via a user-specified threshold [46, 28, 17]. It is modelagnostic but lacks statistical inferential guarantee for pre-defined threshold and does not address the fundamental issue that Shapley values capture both direct and indirect effects. A following variation of this idea is BorutaSHAP [20]. It adapts the Boruta algorithm by generating permuted "shadow" features and identifying a feature significant only if its SHAP importance consistently exceeds the 2
best shadow. It resolves the instability induced by a user-specified SHAP threshold. However, it is computationally intensive due to iteratively retraining and SHAP recomputation on an expanded feature set, and it is not model-agnostic, as it is typically tied to tree-based learners. A relevant recent work, SHAP-XRT [45], combines Shapley values with randomization-based conditional independence tests to generate p-values for feature relevance. However, it relies on accurate modeling of conditional distribution of the data, and is more computationally demanding than standard SHAP due to repeated model evaluation on randomized samples to construct a null distribution for each feature. Our work is different from the existing methods. Although we do not directly aim to address the computational limits of Shapley values, our method does not incur extra cost beyond standard Shapley value estimation. Moreover, MinShap algorithm and associated p-value test avoid the model-averaging issue inherent in Shapley value and provides statistical inferential guarantees. Multiple testing for feature selection Our MinShap algorithm can be viewed as a multiple testing algorithm since a minimum feature importance is taken over feature importance across multiple feature lists. In section 4 we focus on the connection to multiple testing approaches. A related approach is stability selection [29, 39] which creates multiple feature lists under different data perturbation. Our MinShap algorithm, instead of using resampling strategy, uses multiple hypothesis testing across different orderings or conditioning subsets. Unlike stability selection, by exploiting the faithfulness condition we are able to construct a single test statistic and theoretically justified threshold for feature selection.
2
Background and Preliminary
2.1
Feature selection
In this section, we will introduce key notation, and then define feature selection from the conditional independence perspective. First, define covariates X = (X1 , . . . , Xp ) ∈ Rp and response Y ∈ R. Suppose each observation O consists of (X, Y ), and observations O1 , ..., On are independently drawn from true distribution P0 and O1 , ..., On combine to form the empirical distribution Pn . For notation convenience, for any population functional F we use the notation F0 = F (P0 ) and for it’s empirical counterpart Fn = F (Pn ). Denote S ⊂ {1, ..., p} as index subset. Let XS ∈ R|S| be the sub-vector of X indexed by S, and let X−S ∈ Rp−|S| be the complementary subvector with indices Sc. Definition 2.1. The optimal subset of selected features is defined as: Se0 ∈ arg min Y ⊥ ⊥ X−S |XS |S|
where Y ⊥ X−S |XS refers to the conditional independence statement Y independent of X−S conditioned on XS and conditional independence is defined in terms of the distribution P0 . Note that this is a standard definition (see e.g. [1]) but is extremely computationally expensive as it requires searching over 2p subsets. Thus from a computational perspective, it is often natural to consider the feature-wise definition: Definition 2.2. The selected subset of features is defined as: S0∗ = {j : Y ̸⊥ ⊥ Xj |X−j } While Se0 and S0∗ are not always equivalent, under conditions such as e.g. [42] Se0 = S0∗ and by using the feature-wise definition, we reduce feature selection to determine whether each feature satisfies the conditional independence statement Y ⊥ ⊥ Xj |X−j or not. For the remainder of the paper, we adopt this definition for feature selection. 3
2.2
Feature importance
Feature-wise selection or conditional independence testing naturally lends itself to hypothesis testing and the key is which statistics to use. A natural framework involves test statistics based on feature importance. We use the definitions and framework of feature importance adopted in Williams et al.[52]. First, define feature importance as a population parameter in terms of oracle predictive power. Let V (f, P0 ) be a measure of the predictive power of a prediction function f of the model, where f ∈ F, F is a rich function class. A natural candidate prediction function would be a population maximizer f0 over the class F: f0 ∈ argmaxf ∈F V (f, P0 ). Similarly, define prediction function f0,−S be the maximizer over function class F−S where F−S is a collection of function f ∈ F whose evaluation ignores entries of input X with index in S. We construct estimators fn , fn,−S for f0 , f0,S on the observed samples receptively. In the population setting, for the subset of features XS , importance is estimated by the difference in the predictive power by excluding features in S. Thus, the feature importance score of the subgroup XS is defined as: V I0,S := V (f0 , P0 ) − V (f0,−S , P0 ). We assume that V I0,S is non-negative, and larger value indicates higher predictive power of XS . Examples of measures V (f, P0 ) include R2 , deviance, and classification accuracy, etc. Throughout this paper we specifically use mean-squared error (MSE) as predictive measure: V (f0 , P0 ) = E0 [(Y − f0 (X1 , ..., Xp ))2 ]. The corresponding estimator of variable importance score V I0,S is V In,S = V (fn , Pn ) − V (fn,−S , Pn ).
2.3
Shapley value (feature attribution)
Feature importance is a global, population-level notion used to summarize how relevant a variable is to overall prediction or inference, and is used for tasks such as feature ranking or selection. In contrast, Shapley value is a local feature attribution method [43, 16, 7, 33], because it explains individual predictions by performing a local, instance-specific computation. In this paper, we adapt Shapley values for feature selection. Instead of defining the Shapley value in terms of the game theory perspective [40], we use the global population version of Shapley value for feature importance. Specifically, the Shapley value ϕ0,j (resp. ϕ0,j (V )) for feature Xj , j ∈ {1, ..., p} is defined as X
ϕ0,j =
S⊂{1,...,p}\{j}
1 p − 1 −1 {V (f0,S , P0 ) − V (f0,S∪j , P0 )} p |S|
(1)
In other words, the Shapley value for the j th feature is the weighted average gain in oracle predictive power by including Xj over all possible subsets. One of the benefits of Shapley values as a global feature importance is that it is more stable when features are highly dependent because it tends to split/attribute importance. The Shapley value can alternatively be expressed in terms of all possible order of permutations, which was employed for estimation from game theory perspective by Castro et al. [6]. We will re-express it in terms of all possible order of features from a feature importance perspective. ϕ0,j =
1 X {V (f0,[π]j−1 , P0 ) − V (f0,[π]j−1 ∪{j} , P0 )} p! π∈Π(p)
=
(2)
1 X π V I0,j p! π∈Π(p)
4
where Π(p) is all possible permutations of the feature set {1, ..., p}, and [π]j−1 represents the π := set of predecessors of the feature j in the ordering π. The Shapley value summand V I0,j V (f0,[π]j−1 , P0 ) − V (f0,[π]j−1 ∪{j} , P0 ) is the marginal contribution for j in a given order π ∈ Π(p). The Shapley value for feature j is the average of these marginal contributions over the possible different permutations of the feature set p. To see the equivalence between (1) and (2), we can first expand the Equation (1) as ϕ0,j =
X S⊂{1,...,p}\{j}
|S|!(p − |S| − 1)! {V (f0,S , P0 ) − V (f0,S∪j , P0 )}. p!
Then, consider that there are |S|! distinct orderings in which the features in subset S precede feature j in the permutation, whereas (p − |S| − 1)! denotes the number of unique orderings of the remaining features in {1, ..., p} \ S ∪ {j} that appear after feature j. This paper focuses on the definition from a permutation perspective in Equation (2), but the algorithm and related p-value tests in following sections of the main paper are also applicable to Shapley value framework based on combination formulation. Properties The Shapley value expressed in terms of variable importance also inherits the following properties: • Non-negativity: by construction, ϕ0,j ≥ 0. P • Efficiency: pj=1 ϕ0,j = V (f0,[π]p , P0 ) − V (f0,[π]∅ , P0 ). • Symmetry: if Xi = Xj , then ϕ0,i = ϕ0,j . • Null feature: for a feature Xj , if V (f0,[π]j−1 , P0 ) = V (f0,[π]j−1 ∪{j} , P0 ) for all π ∈ Π(p), then ϕ0,j = 0. • Linearity: given two variable importance measures V and V ′ , if V ′ = αV , then the Shapley value of feature j is ϕ0,j (V ′ ) = αϕ0,j (V ). The cost of computing the Shapley value grows exponentially with p. To estimate the Shapley value efficiently, we adopt the widely used Monte Carlo strategy first proposed by Castro et al. [6], which approximates the value by averaging marginal contributions over K randomly sampled permutations, achieving polynomial-time complexity on average. Let Π = {π1 , . . . , πK } ⊆ Π(p) denote the sampled permutations. For each πk ∈ Π, the Equation (2) and the predictive measure V (e.g., MSE) suggest using the plug-in estimator πk V In,j := V (fn,[πk ]j−1 , Pn ) − V (fn,[πk ]j−1 ∪{j} , Pn )
where fn,[πk ]j−1 and fn,[πk ]j−1 ∪{j} are estimators of the population optimizer f0,[πk ]j−1 and f0,[πk ]j−1 ∪{j} . The Shapley value plug-in estimator is ϕn,j :=
1 X πk V In,j . K πk ∈Π
The adaptive full estimation procedure of ϕn,j for the variable importance measure is given in Appendix B Algorithm 2.
5
2.4
Limitations of Shapley value
One of the limitations of Shapley values from a feature selection perspective is that it is possible that: Y ⊥ ⊥ Xj |X−j but ϕ0,j ̸= 0. For example, consider a DAG model Example 2.1 based on variable importance with MSE measure, clearly it follows that Y ⊥ ⊥ X1 |X2 and Y ⊥ ⊥ X1 |X3 , but Y ̸⊥ ⊥ X1 . Then it follows that a non-zero Shapley value does not indicate a variable has direct influence on the target. Thus, the main limitation of using Shapley value for feature selection is that ϕ0,j = 0 =⇒ Y ⊥ ⊥ Xj |X−j , but ϕ0,j ̸= 0 =⇒ ̸ Y ̸⊥ ⊥ Xj |X−j Example 2.1. DAG: X1 → X2 → X3 → Y , X2 = X1 + γ; X3 = X2 + δ; Y = X3 + ϵ, where Xj ∼ N (0, 1), j ∈ {1, 2, 3} and ϵ, δ, γ ∼ N (0, 1) are independent of X. Table 1: Marginal contributions for features in chain DAG example. Permutation
π V I0,1
π V I0,2
π V I0,3
[X1 , X2 , X3 ] [X1 , X3 , X2 ] [X2 , X1 , X3 ] [X2 , X3 , X1 ] [X3 , X1 , X2 ] [X3 , X2 , X1 ]
1 1 0 0 0 0
1 0 2 2 0 0
1 2 1 1 3 3
ϕ0,3 =
11 . 6
The mean Shapley values for (X1 , X2 , X3 ) are: 2 ϕ0,1 = , 6
5 ϕ0,2 = , 6
Besides, using Shapley value for feature selection can also result in suboptimal performance since the top-ranked variables with the largest Shapley value score does not imply they are non-redundant. See Markov boundary example discussed in [26, 12] where the Markov boundary members all have smaller Shapley values than the non-Markov boundary member. Given that rank selection based on Shapley values tend to overselect insignificant variables and underselect significant variables, it leads to inflated Type I error and reduced power.
3
Our Algorithm: MinShap
In this section we introduce our MinShap algorithm. First we motivate the algorithm by linking to directed acyclic graphical models (DAGs) and their connections to conditional independence.
3.1
Motivating MinShap via DAGs and faithfulness
A DAG represents directed relationships among variables, where nodes/vertices correspond to variables (e.g. features or outcome). If there is a directed path Xj → ... → Xk , then Xk is called a descendant of Xj , and Xj is an ancestor of Xk . Our motivation for the new algorithm is based on the faithfulness assumption in DAGs (see e.g. [41]). To see the significance of faithfulness, let’s first introduce its definition tailored for feature selection. Definition 3.1. (d-separation). Two nodes Xj (feature) and Y (outcome) are d-separated by a set of nodes Z if all of the paths between Xj and Y are blocked by Z.
6
The d-separation rules can be used to determine whether two variables are conditionally independent based only on the arrows in DAGs. In other words, if two nodes Xj and Y are d-separated by a set of nodes Z in the graph, then X and Y are conditionally independent given Z (Z = XS ) in the data. Assumption 1. (Faithfulness). The probability distribution P0 is faithful to DAGs, such that the conditional independence relationships that hold in P0 are exactly those implied by the DAGs via d-separation. As a consequence of the faithfulness assumption, we are able to claim that any independence relations in the data are caused by the underlying structure of the graph that generated it, so there is no accidental conditional independence arising from random coincidence or special cancellations. Without this assumption, a Shapley summand could be zero because of accidental cancellation, even though there is still an active path in the graph. Under this setting, for any feature Xj , it is uninformative to Y if and only if Y ⊥ ⊥ Xj |X−j . In other words, conditioned on X−j , the variable Xj carries no additional information about Y . Furthermore, to see how the MinShap algorithm is related to conditional independence testing, we need a no reverse causation assumption: Assumption 2. (No reverse causation). The outcome Y in DAGs does not causally influence any covariate, such that Y is not ancestor of any vertex Xj in DAGs. Lemma 1. Under assumption 2 that there is no reverse causation/data leakage: Y ⊥ ⊥ Xj |XS ⇒ Y ⊥ ⊥ Xj |X−j for any S ⊂ {1, ..., p}\{j}. Proof. The proof follows from the fact that since S ⊂ {1, ..., p}\{j}, using the language of DAG models [35, 36], conditioning on an additional node Xk can only create a path to Y if it is a collider Xj → Xk and Y → Xk . But under the no data leakage assumption, there can be no edge from Y to the covariates. Therefore conditioning on an additional node cannot lead to conditional dependence. The no reverse causation is standard in machine learning because we assume that the response Y is a function of the features (X1 , X2 , ..., Xp ) and not the other way around. π is the Shapley value summand for feature j captures the For an ordering π ∈ Π(p), V I0,j additional contribution of feature j to the target given ordering set [π]j−1 . That is to say, under a true ordering π ∗ , π∗ V I0,j =0⇔Y ⊥ ⊥ Xj |X[π∗ ]j−1 . Since π ∗ is typically unknown, if we assume (X1 , ..., Xp , Y ) follow a recursive structural DAG model where Y is a terminal node in the DAG [35, 36], then we have the following result: Theorem 1. If assumptions 1 and 2 hold in DAG modeling, then we have (min)
ϕ0,j
π := min V I0,j =0⇔Y ⊥ ⊥ Xj |X−j π∈Π(p)
or equivalently π Y ̸⊥ ⊥ Xj |X−j if and only if V I0,j > 0, for all π ∈ Π(p). π is non-negative and each permutation π can also be seen Proof. Since each Shapley summand V I0,j π = 0 is equivalent to the statement that there exists a valid subset S as a subset S, minπ∈Π(p) V I0,j such that Y ⊥ ⊥ Xj |XS . Thus, the ⇒ direction holds due to Assumption 1 and Lemma 1, and the ⇐ direction holds by choosing the subset S = {1, ..., p}\{j}.
7
Hence taking the minimum Shapley Value summand over permutations ensures the desired feature selection guarantee. As a simple illustration, consider Example 2.1, where X1 and X2 provide no extra information once X3 is observed. Yet the conventional Shapley values for X2 and X3 are both nonzero, violating the nature axiom for feature selection, which should identify a null feature when there is no direct edge to the target. However, MinShap focuses on the existence of null importance of feature. The MinShap for X1 and X2 are zeros, and thus only X3 will be selected. Note that while we use the DAG modeling framework to motivate the MinShap algorithm, the approach applies whether the data-generating process is a DAG model or not.
3.2
Theory for setting the threshold
Now, let’s present the main result about how to select the significant features by MinShap πK π1 algorithm. Denote VIn,j = (V In,j , ..., V In,j ) ∈ RK as the vector of K feature importance statistics for sampled permutations, and the order statistics (1)
(K)
V In,j ≤ ... ≤ V In,j . (min)
(1)
We define the MinShap test statistic as ϕn,j := V In,j . The MinShap algorithm resides in the idea π > 0 for all π, we reject the null that X and Y are conditionally independent. Define that if V I0,j j 2,(k)
2,π 2,π σn,j as the estimator of σ0,j := Var(ϵ2[π]j−1 ∪{j} − ϵ2[π]j−1 ), where ϵ[π] := Y − f0 (X[π] ). σn,j associated order statistics, k ∈ [K]. We now formally state the validity of this test.
is the
π = Theorem 2. If assumptions 1 and 2 hold, then under the null hypothesis H0 : minπ∈Π(p) V I0,j (1)
0⇔Y ⊥ ⊥ Xj |X−j , the MinShap V In,j estimated with sufficiently large K = O(|S0∗ |) (|S0∗ | is the size (1)
of optimal subset/ the number of significant features) permutations is valid, ie: P(V In,j > tj |H0 ) < α, q 2,(1) for all α ∈ [0, 1] with threshold tj > −2 ln(α)σn,j . Remark 1. The MinShap algorithm works well in filtering out the insignificant feature j when |S0∗ |, the number of significant features, is relatively small. Consider the worst case when feature j carries partial information of all significant features, the chances of selecting the correct true permutation order where j is placed after all significant features will be high if |S0∗ | is small. Remark 2. The Type II error is driven by the number of permutation K due to the independent multiple hypothesis testing procedure. Thus, fewer permutation K maximizes the power of MinShap. For further details, see Appendix A.2.
3.3
MinShap algorithm
Based on the threshold we setup in the prior section 3.2, we state Algorithm 1, which implements the MinShap for feature selection.
3.4
Computational complexity
The computation complexity of MinShap is the same as Shapley value except that we also compute the variance term for each feature j under permutation π for further determination of threshold. There are a number of ways to lower the computation complexity of computing MinShap. First, for different feature j and different permutation π ∈ Π, the computational jobs can be run in parallel and aggregate results at the end. In this way, the computation complexity can be reduced to O(p) at the price of higher storage space. Second, the number of permutation order set we used to approximate Shapley value can be proportional to the number of significant features |S0∗ | not the number of features p, as we discuss in Theorem 2. The classic Shapley value approximation 8
Algorithm 1 MinShap algorithm for feature selection Require: Dataset (X, Y ); predictiveness measure V ; number of features p; number of permutations K; significance level α 1: Fit null model fn,∅ ; set V∅ ← V (fn,∅ , Pn ) and residuals e2∅ 2,πk πk 2: Initialize storage arrays {ϕn,j } and {σn,j } for all j ∈ [p], k ∈ [K] 3: for k = 1, . . . , K do 4: Sample a random permutation πk 5: Set Vcur ← V∅ and e2cur ← e2∅ 6: for j = 1, . . . , p do πk 7: Set predecessor set Pj−1 ← { i ∈ [p] : πk (i) < πk (j) } 8: Fit model fn,P πk ∪{j} and compute Vnew ← V (fn,P πk ∪{j} , Pn ) and residuals e2new j−1
πk V In,j ← Vcur − Vnew 2,πk πk πk ; σn,j ← V ar(e2cur − e2new )/n 10: ϕn,j ← V In,j 2 11: Vcur ← Vnew ; ecur ← e2new 12: end for 13: end for 14: for j = 1, . . . , p do (min) k 15: ϕn,j ← mink∈[K] ϕπn,j q 2,(1) 16: tj ← −2 log(α) · σn,j
j−1
9:
(min)
Reject H0,j if ϕn,j ≥ tj 18: end for (min) 19: Return ϕn,j , tj , and decisions for all j ∈ [p]
17:
via Monte Carlo sampling [6, 37] requires O(p2 ) in time complexity where we have two nested loops: outer over permutations p, and inner over features in the permutation p. Now, the time complexity can be reduced to O(|S0∗ |p). Third, to save time from retraining models, we can either use Dropout method where the j th feature is replaced by its corresponding marginal mean and then reinserted into the pre-trained full model f0 , or the Lazy-training [13] by gradient computation in neural network set up to estimate reduced model parameter.
4
Additional Approaches Based on p-values and Multiple Testing
MinShap can naturally be connected to the maximum p-value in multiple conditional independence tests. The construction of p-value test for feature importance is derived from the framework of Williamson et al.[52] that the estimator V (fn , Pn ) of value V (f0 , P0 ) is asymptotically normal and efficient if several standard conditions related to smoothness, identifiability, and finite moments hold. π by continuous Thus, the result can immediately be applied to variable importance estimator V In,j mapping theorem. Theorem 3 (Re-formulated from [52]). If assumptions in [52] hold, then each estimator of Shapley π for j ∈ [p] under permutation order π, has the following result: value summand V In,j d √ π 2,π π n V In,j − V I0,j → N (0, σ0,j ). With Theorem 3, for each j, we can get a p-value for each permutation order π with the test π = V I π /σ π under the hypothesis that H π statistic T0,j 0π : V I0,j = 0 where π ∈ Π(p), and we have 0,j 0,j (max)
p0,j
:= max pπ0,j < α ⇔ Y ̸⊥ ⊥ Xj |X−j . π∈Π(p)
9
k Denote pπn,j as estimated p-value for pπ0,j of feature j under permutation πk ∈ Π with null hypothesis
(1)
(K)
k H0k , and pn,j ≤ ... ≤ pn,j as the order statistics of {pπn,j : k = 1, ..., K}. From a multiple testing perspective, each permutation/ordering generates a feature importance statistic and a corresponding p-value for feature j; we then aggregate these p-values to make the best decision. MinShap modifies the Shapley value by taking the minimum marginal contribution across orderings, which targets the direct effect. The link between them is that MinShap operates on the same collection of ordering statistics, which we map to p-values via a normal approximation and perform multiple-testing procedures. Both MinShap and Max-p test the following hypothesis:
K [
πk H0 : Xj ⊥ ⊥ Y |X−j ⇐⇒ ∃ k ∈ [K] s.t. V I0,j = 0 ⇐⇒
πk {V I0,j = 0};
k=1 πk Ha : Xj ̸⊥ ⊥ Y |X−j ⇐⇒ ∀ k ∈ [K], V I0,j ̸= 0 ⇐⇒
K \
πk {V I0,j ̸= 0}.
k=1
Corollary 4. If assumptions 1 and 2 hold in DAG modeling, then under the null hypothesis (K) H0 : Xj ⊥ ⊥ Y |X−j , the Max-p value pn,j estimated with sufficiently large K = O(|S0∗ |) permutations (K)
is valid, ie: P(pn,j < α|H0 ) < α, for all α ∈ [0, 1]. The corresponding proof is available in AppendixTA.3. Note that, in contrast to common multiple πk hypothesis testing - where the global null is H0 : K k=1 {V I0,j = 0} and is rejected if at least one S (k) πk πk V I0,j = ̸ 0 - our conjunction null is H0 : K k=1 {V I0,j = 0}, which is rejected if V I0,j ̸= 0 for all k ∈ [K].
4.1
Extension to partial conjunction hypothesis testing
Due to the possibility that a single p-value is highly significant while all other p-values show no evidence for the alternatives, Max-p might be conservative for testing. Further, the faithfulness assumption tends to often fail in the finite sample setting due to near cancellations [48] which means we may need to allow some error. Benjamini and Heller [2] proposed partial conjunction hypothesis testing (PCHT), which is less stringent that at least u out of K tested hypotheses are false. When u = 1, this is well-known global null hypothesis testing problem where we reject the null if there is at least one false null. Our Max-p value corresponds to the case when u = K, where we reject the null if all are false null. To define PCHT, let’s consider null hypothesis H01 , ..., H0K with corresponding p-value p1 , p2 ..., pK , then u/K
H0
: at least K − u + 1 nulls are true
u/K
H1
: at least u out of K are nonull
where u is a data–driven value selected on a tuning set on a range of u’s with the minimum loss/highest performance. To obtain a p-value for PCHT, we can choose any one of three adjusted u/K p-value as proposed by Benjamini and Heller [2] to provide a valid p-value, pu/K , for H0 : 1. General dependency p-value (Bonferroni): u/K (u) pn,j = (K − u + 1)pn,j ; 2. Independent p-value (Stouffer): P u/K
pn,j = 2 ∗ (1 − Φ(
K−u+1 (k) z k=1 √ )); K−u+1
3. Independent p-value (Fisher): P u/K (k) pn,j = P(χ22(K−u+1) ≥ −2 K k=u log pn,j ). 10
a Shapley values are approximated Figure 1: Overview of the MinShap algorithm and variations. ○: by averaging marginal contributions/feature importance over K randomly sampled permutations; b MinShap only tests the minimum contribution under faithfulness in DAGs; ○: c each feature ○: d connect MinShap to Max-p via multiple test importance test statistic yields a p-value test; ○: e PCHT, a resolution for conservativeness of Max-p. perspective; ○: ∗u/K
Then apply the Holm adjustment to get pn,j 2, ..., K. We have:
∗u/K
Theorem 5 (Re-formulated from [2]). pn,j
∗1/K
where pn,j
(1)
∗u/K
= pn,j , pn,j
∗ u−1
u/K
= max{pn,jK , pn,j }, u =
is a valid p-value for the partial conjunction null
u/K ∗u/K u/K hypothesis H0 , ie: P (pn,j < α|H0 ) < α, for all α ∈ [0, 1].
We remark that Bonferroni is more robust since it is valid under arbitrary dependence among p-values, whereas Stouffer and Fisher require independence. If there is any underlying structure of the selected permutation order set, then Bonferroni method is recommended. In our case, permutation sampled i.i.d. with replacement yields independent p-values across repetitions, so we can apply any of the three adjustment methods. The Bonferroni method is more conservative since it looks only at a single order statistic: the (K − u + 1)th largest p-value, whereas Fisher and Stouffer methods aggregate all evidence across the tail of largest p-values. The procedure of computing adjusted p-values does not require any extra model retraining, so the computation complexity is the same as MinShap and Max-p. When screening for different levels of u, one might worry about multiple testing, Holm adjustment to all above three methods are used to control overall false discovery rate (FDR) [2]. The procedure of the computation for Max-p and adjusted p-value is summarized in Algorithm 3. Fig.1 provides an overview of our procedure in this paper.
4.2
Practical guidelines for the choice of methods
We now outline some practical guidance on which methods we could choose over others in different settings. When the Type II error is close to 0, MinShap and Max-p perform similarly to the adjusted p-value methods, and it is therefore sufficient to use MinShap and Max-p alone. A potential concern is the conservativeness of MinShap and Max-p, that is, the Type II error is not close to 0. In this case, one may instead apply the adjusted p-value methods, which can control the inflated Type II error without incurring a substantial additional computational cost. 11
In a data application, where the significant feature set is unknown, a practical strategy is to apply all methods initially. If their results are largely consistent, then MinShap and Max-p are preferred, since the choice u = K in the adjusted p-value approach yields a procedure equivalent to them. If, however, the results differ substantially, this may indicate that the Type II error of MinShap and Max-p is non-negligible; in that case, the adjusted p-value methods can be applied.
4.3
Connection and comparison to conditional independence testing methods
As pointed out earlier, MinShap (and by extension Max p-value) may be viewed as resolving the conditional independence problem: (min)
ϕ0,j
π := min V I0,j =0⇔Y ⊥ ⊥ Xj |X−j . π∈Π(p)
Thus it is natural to ask the question, how does MinShap compare as a feature selection approach to other conditional independence testing approaches? (see e.g. [38, 23, 4, 44]) The main distinction between these standard approaches and MinShap as pointed out in the previous section, involves performing multiple conditional independence tests (corresponding to different permutation orderings) and using the connection to DAG models to perform a single test. Many conditional independence testing approaches (e.g. GCM [38], LOCO [23] e.t.c.) perform a single test. We demonstrate in the next section the advantage in terms of both accuracy and stability of using MinShap and the multiple testing approaches.
5
Simulations
In this section, we will present simulation results validating our MinShap algorithm together with Max-p for feature selection, and further provide guidance on adjusted p-value methods. We will also compare them with other state-of-the-art methods such as GCM, LOCO and Lasso. We assess their accuracy and stability for the linear model, non-linear non-additive model, conditional interaction model and logistic model. Specifically, we consider 20 covariates X ∼ N (0, Σ20×20 ) with diagonal equals 1, and the following different settings are used: (a) Linear model: Y ∼ 4X1 + 4X2 + 3X3 ∗ X4 + 3X5 + 2X6 + 2X5 ∗ X6 + X7 + X8 + ϵ, where Corr(X3 , X4 ) = 0.5, ϵ ∼ N (0, 1). (b) Non-linear non-additive model : Y ∼ 2 sin(X1 ) + 2 log(|X2 | + 1) + X1 X2 + 3 cos(X3 + X4 ) + max(0, X5 ) + X6 X7 X8 + ϵ, where Σ is block diagonal matrix, and within-block correlations are (0, 0.2, 0.5, 0.8), ϵ ∼ N (0, 1); (c) Conditional interaction model : Y = 1.5X1 X2 ∗ I(X3 > 0) + X4 X5 ∗ I(X3 < 0) + 3X6 X7 ∗ I(X8 > 0) + X9 X10 ∗ I(X8 < 0) + ϵ, where Corr(X1 , X2 ) = Corr(X6 , X7 ) = 0.9; Corr(X4 , X5 ) = Corr(X9 , X10 ) = 0.5, ϵ ∼ N (0, 1); (d) Logistic Model : Y ∼ f (2.5X1 + 2.5X2 + 2X3 X4 + 1.5X5 + 1.5X6 + X72 + X83 ) + ϵ, where f is sigmoid function, Corr(X1 , X2 ) = 0.5, ϵ ∼ N (0, 0.12 ).
5.1
Feature selection comparison across methods
For all models, we use 3000 samples over 100 simulations. The number of permutations K is 50 for MinShap and Max-p value approaches. Our algorithms are model-agnostic: we evaluate them both using XGBoost and feed-forward neural networks for comparison and obtain consistent 12
Conditional Interaction
Logistic
0.54
0.59
0.59
0.5
0.59
0.52
0.52
0.53
0.52
O
in M
C
ax
LO
C G
M
in M
La s
O C
ax
LO
M
M C G
in M
La s
O C
ax
LO
M
C G
La s
Sh ap
0.53
−p
0.54
M
0.54
so
0.46
Sh ap
0.51
−p
0.46
so
0.46
Sh ap
0.51
−p
0.5
M
0.51
so
0.51
M
in
Sh ap
O
−p
C
Non−Linear
ax
LO
M
G
C
M
1.00 0.75 0.50 0.25 0.00
so
0.6 0.4 0.2 0.0
La s
Accuracy F1 Recall
1.00 0.75 0.50 0.25 0.00
FDR
Performance
1.00 0.75 0.50 0.25 0.00
Type.I.Error
Linear 1.0 0.8 0.6 0.4
Method Method
GCM
Lasso
LOCO
Max−p
MinShap
Figure 2: Average performance measure with standard deviation for model (a)-(d) for feature selection comparison using XGBoost. We see that MinShap and Max-p outperform other feature selection methods. results. Specifically, we first mainly discuss results based on XGBoost, and for results in feed-forward neural networks, see Appendix C.1. The overall performance comparison is presented in Fig.2 and the dashed line is the empirical threshold from random fair coin-flipping. From Fig.2, we see that although Lasso requires the least training time, it only works relatively well in linear model and its accuracy are low in all other models. GCM tends to under-select significant features, resulting in low power, whereas LOCO tends to over-select insignificant features, leading to high Type I error. LOCO can be inefficient under strong feature dependence: correlated features can substitute for each other, inducing correlation bias [49]. GCM alleviates the effect of feature correlation by considering the covariance of residuals from two regressions, but it suffers from a natural limitation: a covariance of zero does not necessarily imply independence. For example, the scenario occurs when feature Xj follows a symmetric distribution. GCM and LOCO do not perform well when the model are not additive, as we can see that their corresponding F1 scores are low in Fig.2. Across different models (a)-(d), the performance of MinShap and Max-p out-perform those three baseline methods with higher accuracy and F1 score, and maintain Type I error control. We assess stability of the selected features via the Jaccard index [31], defined as the average pairwise similarity of selected feature sets: N −1
N
X X 2 JN = Jij . N (N − 1) i=1 j=i+1
where N is number of simulation, Jij ∈ [0, 1] is the proportion of commonly selected features in two simulations out of the union of all features selected in two simulations. A higher Jaccard index indicates a higher overall similarity and lower dispersion among the selected feature sets, and thus the algorithm is more stable [30]. From Table 2, we find that MinShap and Max-p are more stable than LOCO and GCM. While Lasso is stable, its accuracy is completely compromised in all nonlinear settings.
13
Table 2: Jaccard Index for the stability of feature selection methods on model (a)–(d) using XGBoost. Feature Selection Method
5.2
Model
MinShap
Max-p
LOCO
GCM
Lasso
(a) (b) (c) (d)
0.98 0.78 0.91 0.98
0.99 0.77 0.89 0.87
0.71 0.74 0.71 0.67
0.84 0.51 0.54 0.72
0.82 1.00 0.94 1.00
Comparison to stability selection
One of the motivations of our MinShap and Max-p algorithms is that they increase stability. In this section, we will compare the performance in terms of various accuracy matrics as well as run-time with the stability selection approach developed in [29, 39]. Stability selection is a resampling-based feature selection method that repeatedly runs a base feature selection method on different sub-samples of the data and retains variables with high selection frequency across splits. Specifically, we will apply the stability selection to LOCO, GCM and Lasso to model (a)-(d) over 100 simulations with a 50% sub-sampling rate and a selection threshold of 0.8. Then, we compare their performance with MinShap and Max-p, which, as discussed in the previous section 4, can be regarded as the multiple testing based feature selection method. The setup for MinShap and Max-p remain the same as previous section 5.1 with XGBoost method. Our results both shows how our method has similar performance to LOCO-stability but also significantly shorter run-time.
0.51
0.46
0.54
0.54
0.53
0.54
0.59
0.59
0.5
0.59
0.52
0.52
0.53
0.52
M C
ss
G
La
−S
o−
M C
ss
G
La
−S M C G
M C G
La
−S ta b o− ility LO St ab C il O −S ity ta bi lit y M ax −p M in Sh ap
0.46
ta bi lit y LO St ab C ilit O y −S ta bi lit y M ax −p M in Sh ap
0.46
ta bi lit y LO St ab C il O −S ity ta bi lit y M ax −p M in Sh ap
0.51
−S
0.5
o−
400 300 200 100 0
Logistic
0.51
ss
1.00 0.75 0.50 0.25 0.00
Conditional Interaction
0.51
La
FDR
0.6 0.4 0.2 0.0
Non−Linear
t ss abi lit o− y LO St ab C ilit O y −S ta bi lit y M ax −p M in Sh ap
Accuracy F1 Recall
1.00 0.75 0.50 0.25 0.00
Type.I.Error
Performance & Time
1.00 0.75 0.50 0.25 0.00
Time (s)
Linear 1.0 0.8 0.6 0.4
Method Method
GCM−Stability
Lasso−Stability
LOCO−Stability
Max−p
MinShap
Figure 3: Average performance measure and time with standard deviation for model (a)-(d) for multiple testing based and stability selection based feature selection comparison using XGBoost. We see that MinShap and Max-p are more efficient than stability-based feature selection methods. From the Jaccard index in Table 3, we observe that the stability of our method and stability selection are similar especially for LOCO-stability. Fig. 3 shows that after applying stability selection, the performance of Lasso and GCM remains roughly unchanged or is slightly improved in terms of accuracy and F1 score. The performance of LOCO improves significantly. Moreover, the Type I error and FDR are both well controlled for all three baseline feature selection methods after applying 14
stability selection. However, despite the improvements brought by stability selection, our algorithms achieve similar performance in models (a), (c), and (d), and outperforms LOCO with stability selection in model (b) while always having significantly faster run-time. Therefore, MinShap and Max-p are more efficient. Table 3: Jaccard Index for the stability comparison of multiple testing based and stability selection based feature selection methods on model (a)–(d) using XGBoost. Feature Selection Method
5.3
Model
MinShap
Max-p
LOCO-Stability
GCM-Stability
Lasso-Stability
(a) (b) (c) (d)
0.98 0.78 0.91 0.98
0.99 0.77 0.89 0.87
1.00 0.72 0.99 0.93
1.00 0.92 1.00 1.00
0.85 0.98 1.00 1.00
Strength and limitation of MinShap/Max-p and the application of other adjusted p-value methods
In this section, we will further discuss the strength and limitation of MinShap and Max-p, and provide practical guidance for when we use adjusted p-value methods. We will mainly focus on the first two models provided in section 5: linear and non-linear models. We will discuss the relationship among the number of permutation orders K, the proportion of significant features |S0∗ |, the number of features p and further extend to high-dimensional case. Type I Error 0.8
0.8
0.6
0.6
0.4
Type II Error
0.8
Type I Error
1.0
0.4
0.6 0.4 0.2
0.0
0.0
0.0
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.4 0.2
Type II Error
0.2
Type I Error
0.2
F1
Non-linear
Type II Error
1.0
F1
Linear
F1 1.0
0.4 0.2
0.0 20
30
40
50
Number of permutations
60
0.4 0.2
0.0 10
0.6
0.0 10
20
30
40
50
Number of permutations
Bonferroni Fisher GCM
LOCO Lasso Stouffer
60 max_p min
10
20
30
40
50
Number of permutations
60
Figure 4: Average performance measure with standard deviation of different number of permutation for linear and non-linear model with sample size 3000. We see that when Type II error is nearly 0, more permutations drive higher performance of MinShap and Max-p. We first assess the effect of the number of permutation K to MinShap and Max-p value methods. For Fig. 4, the setup of the model is the same as section 5, except that we use different number permutations that is proportional to |S0∗ | = 8, with K range from 8 to 64. The screening range u ∈ [K − 7, K] to avoid overlap with the previous selecting interval. When the sample size is large enough, and the Type II error is extremely close to 0, the more permutations we perform, the higher F1 score MinShap and Max-p value will achieve. The higher number of permutation K drives the Type I error down. Next, we discuss the relation between the number of permutation K and when 15
Type I Error 0.8
0.8
0.6
0.6
0.6
0.4
Type II Error
0.8 Type I Error
1.0
0.4
0.4 0.2
0.0
0.0
0.0
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4 0.2
Type II Error
0.2
Type I Error
0.2
F1
Non-linear
Type II Error
1.0
F1
Linear
F1 1.0
0.4 0.2
0.0
0.2
0.0 10
20
30 40 50 Number of Permutations
60
0.4
0.0 10
20
30 40 50 60 Number of Permutations Bonferroni LOCO max_p Fisher Lasso min GCM Stouffer
10
20
30 40 50 Number of Permutations
60
Figure 5: Average performance measure with standard deviation of different number of permutation for linear and non-linear model with sample size 1000. We see that when Type II error is between 0 and 1, more permutations drive higher Type II error of MinShap and Max-p. Type II error is between 0 and 1. We still compare the same linear and non-linear model but with sample size 1000. The number of permutation K still ranges from 8 to 64 with the screening range u ∈ [K − 7, K] . From Fig. 5, we still see that more permutations always drive the Type I error down for all methods. When the Type II error is between 0 and 1, more permutations make the Type II error of MinShap and Max-p increase. Thus, when Type II error is not close to 1, the adjusted p-value methods can be applied. Type I Error 0.8
0.8
0.6
0.6
0.4
Type II Error
0.8
Type I Error
1.0
0.4
0.6 0.4 0.2
0.0
0.0
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0 0.2
0.3
0.4
0.5
0.6
0.7
Proportion of Significant Features
0.8
Type II Error
0.2
0.0
Type I Error
0.2
F1
Non-linear
Type II Error
1.0
F1
Linear
F1 1.0
0.6 0.4 0.2 0.0
0.2
0.3
0.4
0.5
0.6
0.7
Proportion of Significant Features
Bonferroni Fisher GCM
LOCO Lasso Stouffer
0.8 max_p min
0.2
0.3
0.4
0.5
0.6
0.7
Proportion of Significant Features
0.8
Figure 6: Average performance with standard deviation of varied proportion of significant features for linear and non-linear model. We see the advantage of the adjusted p-value tests in dense regime. Next, we explore how the proportion of significant feature, |S0∗ |, affect the performance of MinShap and Max-p value methods. We fix p = 56, and vary the signal proportion by repeating the same pattern of significant features; each step adds 8 additional significant features by repeating the linear and nonlinear pattern. We set the number of permutations to K = 2|S0∗ |. The best u for 16
Type I Error 0.8
0.8
0.6
0.6
0.4
Type II Error
0.8
Type I Error
1.0
0.4
0.6 0.4 0.2
0.0
0.0
0.0
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.4 0.2
Type II Error
0.2
Type I Error
0.2
F1
Non-linear
Type II Error
1.0
F1
Linear
F1 1.0
0.4 0.2
0.0 20
30
40
Number of Total Features
50
60
0.4 0.2
0.0 10
0.6
0.0 10
20
30
40
Number of Total Features
Bonferroni Fisher GCM
LOCO Lasso Stouffer
50
60
10
max_p min
20
30
40
Number of Total Features
50
60
Figure 7: Average performance measure with standard deviation of fixed number of significant features with increasing number of total features for linear and non-linear model. We see that MinShap and Max-p work well under sparsity regime.
F1
0.6
0.6
0.6
Type II Error
0.8
Type I Error
0.8
0.4
0.4
0.4 0.2
0.0
0.0
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.2
Type II Error
0.2
0.0
Type I Error
0.2
0.4
0.4 0.2
0.0
0 0 0 100 150 200 Number of samples
0
250
0
300
0.4 0.2
0.0
100 300 500
Type II Error
1.0
0.8
F1
Non-linear
Type I Error
1.0
F1
Linear
1.0
0.0
0 0 0 100 150 200 Number of samples Bonferroni LOCO Fisher Lasso GCM Stouffer
100 300 500
0
250
max_p min
0
300
100 300 500
0 0 0 100 150 200 Number of samples
0
250
0
300
Figure 8: Average performance with standard deviation of different sample size for linear and non-linear model. We see the advantage of the suite of adjusted p-value (Bonferroni, Stouffer, Fisher) tests in lower-sample setting.
17
adjusted p-value is chosen from range [0.7K, K] with the highest F1 score. In Fig.6, we find that the MinShap and Max-p value approach work well when the model is sparse. As |S0∗ | increases, the F1 score for both methods gradually decrease whereas F1 score of LOCO and GCM increase. When the sparsity structure of the model does not hold, the adjusted p-value methods can be used to alleviate the conservative effect of MinShap and Max-p caused by dense signal in the model. Slightly different from previous setup, we now fix the number of significant features at |S0∗ | = 8, but increase the total number of features p from 10 to 60, the number of permutation is K = 2|S0∗ |, and the screening range u ∈ [0.7K, K]. From the Fig.7, in terms of both F1 score and Type I error, we see that MinShap and Max-p work well compared to other baseline feature selection methods when the sparsity holds, or equivalently, when the noise the in data is getting higher, and the adjusted p-value methods do not help with large improvement in this case. Then, we compare how sample size n affect their performance. We keep the models’ pattern fixed with varied sample size from 100 to 3000, K = 2|S0∗ |, and the best u is chosen from range [0.7K, K]. In Fig.8, we find that when sample size is small, GCM and LOCO have higher F1 score than MinShap and Max-p, the Type II error of MinShap and Max-p is large and our adjustment methods such as Bonferroni, Stouffer and Fisher help reduce Type II error. To summarize, the MinShap and Max-p coincide with adjusted p-value methods when the Type II error is close to 1. When the Type I error lies between 0 and 1, increasing the number of permutations decreases the power of our MinShap and Max-p. Also, MinShap and Max-p work well when the model is sparse and the sample size is relatively large. As the number features increase or the sample size decreases, the Bonferroni, Stouffer and Fisher adjusted p-value methods improve performance. In particular, Bonferroni is relatively more conservative, and Fisher and Stouffer yield comparable power.
5.4
High-dimensional setting
Next, we focus on high-dimensional setting. We set the total number of features p = 200, and the pattern of significant features repeats once (ie: |S0∗ | = 16 in two models), the number of permutation K = 1.5|S0∗ |, the screening level of u in [0.6K, K] and the sample size increases from 200 to 500. From Fig.9, we see that Lasso works well in sparse linear model, but fails to keep its performance in non-linear model. In high-dimensional case, when the sample size is small, the performance of MinShap and Max-p is not ideal even if the Type I error is still under controlled. The Type II error of MinShap and Max-p is affected by the sample size. As sample size increases, the Type II error decreases. In linear model, the adjusted Stouffer p-value method performs slightly better than LOCO and GCM. In non-linear model, the Stouffer and Fisher adjusted p-value methods are consistently better than other baseline methods. Thus, in the high-dimensional setting where p is relatively close to n, the adjusted p-value methods are favored.
6
Real Data Analysis
In this section, we apply the variable selection methods MinShap, Max-p, GCM and LOCO to two real and well-studied datasets, the well-known wine quality and California housing datasets. Each dataset is analyzed using 5-fold cross-validation. Variable selection is carried out within each training fold, and the mean squared prediction error on the corresponding test fold is used to evaluate predictive performance. In addition, the stability of the feature selections is measured by Jaccard index. Additional real-data experiment is provided in the Appendix D.
6.1
Wine quality data
With increasing interest in wine, quality certification has become essential, relying on sensory evaluations by human experts. Our goal is to predict wine quality using objective analytical 18
Type I Error 0.8
0.8
0.6
0.6
0.4
Type II Error
0.8
Type I Error
1.0
0.4
0.6 0.4 0.2
0.0
0.0
0.0
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0 200
250
300
350
Sample Size
400
450
Type II Error
0.2
Type I Error
0.2
F1
Non-linear
Type II Error
1.0
F1
Linear
F1 1.0
500
0.6 0.4 0.2 0.0
200
250
300
Bonferroni Fisher GCM
350
Sample Size
400
LOCO Lasso Stouffer
450 max_p min
500
200
250
300
350
Sample Size
400
450
500
Figure 9: Average performance with standard deviation of high-dimensional setting for linear and non-linear model. We see the advantage of the suite of adjusted p-value tests in in high-dimensional setting. measurements collected during the certification process. The Wine Quality dataset, introduced by Cortez et al. [9], consists of n = 1599 red wine samples with p = 11 physicochemical features (e.g., acidity, alcohol content). The response is a wine quality score on a 0–10 scale, reflecting the perceived quality of each wine. From Table 4, we find that features selected by LOCO are much less stable than other methods. The most significant features frequently selected by both MinShap and Max-p using XGBoost and random forest(ie: "alcohol", "sulphates" and "volatile acidity") suggest that physio-chemical properties play an important role in predicting wine quality. These findings are consistent with prior literature on wine quality analysis [32]. GCM always select both significant and moderate significant features, such as "pH", "density", etc. [15, 53]. Table 4: Comparison of methods in terms of prediction error and selection stability on wine quality data under two models. Model
XGBoost
Random forest
Method MinShap
MSE 0.4045
se(MSE) 0.0152
Jaccard index 1.00
Max-p
0.4048
0.01495
0.92
LOCO
0.4027
0.01580
0.69
GCM
0.4039
0.01548
0.83
MinShap Max-p LOCO GCM
0.4288 0.4288 0.4285 0.4244
0.0156 0.0156 0.0156 0.0163
1.00 1.00 0.69 0.83
19
Selected Features [‘alcohol’, ‘sulphates’, ‘total sulfur dioxide’, ‘volatile acidity’] [‘alcohol’, ‘sulphates’, ‘total sulfur dioxide’, ‘volatile acidity’] [‘alcohol’, ‘sulphates’, ‘total sulfur dioxide’, ‘volatile acidity’] [‘alcohol’, ‘density’, ‘pH’, ‘sulphates’, ‘total sulfur dioxide’, ‘volatile acidity’] [‘alcohol’, ‘sulphates’, ‘volatile acidity’] [‘alcohol’, ‘sulphates’, ‘volatile acidity’] [‘alcohol’, ‘sulphates’, ‘volatile acidity’] [‘alcohol’, ‘chlorides’, ‘density’, ‘fixed acidity’, ‘pH’, ‘sulphates’, ‘total sulfur dioxide’, ‘volatile acidity’]
6.2
California housing data
Accurate prediction of house prices requires accounting for many factors. We are interested in reducing pricing discrepancies by feature selection for predicting real estate prices. The California housing dataset, originally constructed by Pace and Barry (1997) [21] using data derived from the 1990 U.S. Census, contains n = 20,640 observations of California districts with p = 8 predictive features. These features describe demographic and housing characteristics, including median income (MedInc), median house age (HouseAge), average number of rooms and bedrooms (AveRooms, AveBedrms), population, average number of people per household in that district(AveOccup), and geographic location (latitude and longitude). The response variable is the median house value in each district in units of $100,000. From Table 5, we see that both GCM and LOCO tend to select more moderate or less significant features. The feature selection performed by MinShap and Max-p are stable and always choosing the most significant features such as household income and the location of the property, which is aligned with prior studies [5, 8]. Table 5: Comparison of methods in terms of prediction error and selection stability on California housing data under two models. Model
XGBoost
Random forest
7
Method MinShap Max-p LOCO
MSE 0.3281 0.3283 0.3262
se(MSE) 0.0079 0.0080 0.0083
Jaccard index 0.92 0.82 0.77
GCM
0.3257
0.0081
0.84
MinShap Max-p LOCO
0.3456 0.3450 0.3204
0.0075 0.0072 0.0076
1.00 0.90 0.88
GCM
0.3162
0.0068
0.95
Selected Features [‘AveOccup’, ‘Latitude’, ‘Longitude’, ‘MedInc’] [‘AveOccup’, ‘Latitude’, ‘Longitude’, ‘MedInc’] [‘AveOccup’, ‘AveRooms’, ‘HouseAge’, ‘Latitude’, ‘Longitude’, ‘MedInc’, ‘Population’] [‘AveBedrms’, ‘AveOccup’, ‘HouseAge’, ‘Latitude’, ‘Longitude’, ‘MedInc’, ‘Population’] [‘Latitude’, ‘Longitude’, ‘MedInc’] [‘Latitude’, ‘Longitude’, ‘MedInc’] [‘AveOccup’, ‘AveBedrms’, ‘HouseAge’, ‘Latitude’, ‘Longitude’, ‘MedInc’, ‘Population’] [‘AveBedrms’, ‘AveOccup’, ‘HouseAge’, ‘Latitude’, ‘Longitude’, ‘MedInc’, ‘Population’]
Conclusion
In this paper, we make a connection between a modified Shapley value and feature selection via conditional independence testing. We propose the MinShap algorithm inspired from DAG model’s perspective, establish its connection to the Max-p algorithm, and further extend it with adjusted p-values to mitigate the conservativeness of these earlier procedures in certain circumstances. We present numerical experiments on both simulated and real data sets of varying complexity to validate the theoretical results in practice, and also demonstrate its superior performance than other feature selection methods, such as GCM, LOCO and Lasso. Future work for this study could focus on reducing the computational cost of the proposed procedures. In particular, it would be of interest to develop faster approximations, more efficient sampling or permutation schemes for MinShap and Max-p.
References [1] Constantin F. Aliferis, Alexander Statnikov, Ioannis Tsamardinos, Subramani Mani, and Xenofon D. Koutsoukos. Local causal and Markov blanket induction for causal discovery and
20
feature selection for classification part I: algorithms and empirical evaluation. Journal of Machine Learning Research, 11(1), 2010. [2] Yoav Benjamini and Ruth Heller. Screening for Partial Conjunction Hypotheses. Biometrics, 64(4):1215–1222, 2008. [3] Danilo Bzdok, Denis Engemann, and Bertrand Thirion. Inference and Prediction Diverge in Biomedicine. Patterns (New York, N.Y.), 1(8):100119, 2020. [4] Emmanuel Candes, Yingying Fan, Lucas Janson, and Jinchi Lv. Panning for gold:‘modelX’knockoffs for high dimensional controlled variable selection. Journal of the Royal Statistical Society Series B: Statistical Methodology, 80(3):551–577, 2018. [5] Tomás Capdevielle and Santiago Cifuentes. Feature Relevancy, Necessity and Usefulness: Complexity and Algorithms. 2025. [6] Javier Castro, Daniel Gómez, and Juan Tejada. Polynomial calculation of the Shapley value based on sampling. Computers & Operations Research, 36(5):1726–1730, 2009. [7] Hugh Chen, Ian C. Covert, Scott M. Lundberg, and Su-In Lee. Algorithms to estimate Shapley value feature attributions. Nature Machine Intelligence, 5(6):590–601, 2023. [8] Yucong Chen. Analysis and Forecasting of California Housing. Highlights in Business, Economics and Management, 3:128–135, 2023. [9] Paulo Cortez, António Cerdeira, Fernando Almeida, Telmo Matos, and José Reis. Modeling wine preferences by data mining from physicochemical properties. Decision Support Systems, 47(4):547–553, 2009. [10] Ian Covert, Scott M. Lundberg, and Su-In Lee. Understanding global feature contributions with additive importance measures. Advances in neural information processing systems, 33:17212– 17223, 2020. [11] Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani. Least angle regression. The Annals of Statistics, 32(2):407–499, 2004. [12] Daniel Fryer, Inga Strumke, and Hien Nguyen. Shapley Values for Feature Selection: The Good, the Bad, and the Axioms. IEEE Access, 9:144352–144360, 2021. [13] Yue Gao, Abby Stevens, Garvesh Raskutti, and Rebecca Willett. Lazy estimation of variable importance for large neural networks. In International Conference on Machine Learning, pages 7122–7143. PMLR, 2022. [14] Florian Huber and Volker Steinhage. Conditional Feature Selection: Evaluating Model Averaging When Selecting Features with Shapley Values. Geomatics, 4(3):286–310, September 2024. [15] Khushboo Jain, Keshav Kaushik, Sachin Kumar Gupta, Shubham Mahajan, and Seifedine Kadry. Machine learning-based predictive modelling for the enhancement of wine quality. Scientific Reports, 13(1):17042, 2023. [16] Dominik Janzing, Lenon Minorics, and Patrick Blöbaum. Feature relevance quantification in explainable AI: A causal problem. In International Conference on artificial intelligence and statistics, pages 2907–2916. PMLR, 2020. [17] Neesha Jothi, Wahidah Husain, and Nur’Aini Abdul Rashid. Predicting generalized anxiety disorder among women using Shapley value. Journal of infection and public health, 14(1):103–108, 2021. 21
[18] Alan Jović, Karla Brkić, and Nikola Bogunović. A review of feature selection methods with applications. In 38th international convention on information and communication technology, electronics and microelectronics (MIPRO), pages 1200–1205, 2015. [19] Hiromasa Kaneko. Cross-validated permutation feature importance considering correlation between features. Analytical Science Advances, 3(9-10):278–287, September 2022. [20] Eoghan Keany. BorutaShap: A wrapper feature selection method which combines the Boruta feature selection algorithm with Shapley values. Zenodo, 2020. [21] R. Kelley Pace and Ronald Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(3):291–297, 1997. [22] I. Elizabeth Kumar, Suresh Venkatasubramanian, Carlos Scheidegger, and Sorelle Friedler. Problems with Shapley-value-based explanations as feature importance measures. In International conference on machine learning, pages 5491–5500. PMLR, 2020. [23] Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J. Tibshirani, and Larry Wasserman. Distribution-Free Predictive Inference for Regression. Journal of the American Statistical Association, 113(523):1094–1111, 2018. [24] Scott M. Lundberg and Su-In Lee. A unified approach to interpreting model predictions. In Proceedings of the 31st International Conference on Neural Information Processing Systems, 2017. [25] Shuangge Ma and Jian Huang. Penalized feature selection and classification in bioinformatics. Briefings in bioinformatics, 9(5):392–403, 2008. [26] Sisi Ma and Roshan Tourani. Predictive and Causal Implications of using Shapley Value for Model Interpretation, 2020. [27] Felipe Maia Polo, Yuekai Sun, and Moulinath Banerjee. Conditional independence testing under misspecified inductive biases. Advances in Neural Information Processing Systems, 36:58577–58612, 2023. [28] 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, 2020. [29] Nicolai Meinshausen and Peter Bühlmann. Stability selection. Journal of the Royal Statistical Society Series B: Statistical Methodology, 72(4):417–473, 2010. [30] Majid Mohammadi, Hossein Sharifi Noghabi, Ghosheh Abed Hodtani, and Habib Rajabi Mashhadi. Robust and stable gene selection via maximum–minimum correntropy criterion. Genomics, 107(2-3):83–87, 2016. [31] C. P. Mohana and Karuppasamy Perumal. A survey on feature selection stability measures. International Journal of Computer and Information Technology, 5(1):98–103, 2016. [32] Christoph Molnar, Timo Freiesleben, Gunnar König, Julia Herbinger, Tim Reisinger, Giuseppe Casalicchio, Marvin N. Wright, and Bernd Bischl. Relating the Partial Dependence Plot and Permutation Feature Importance to the Data Generating Process. In Explainable Artificial Intelligence, pages 456–479, 2023. [33] Mohd Amril Nazir, Edmund Evangelista, Syed M Salman Bukhari, and Ravishankar Sharma. A survey of feature attribution techniques in explainable AI: taxonomy, analysis and comparison. Annals of Mathematics and Computer Science, 28:115–126, 2025. 22
[34] Roma Patel, Marta Garnelo, Ian Gemp, Chris Dyer, and Yoram Bachrach. Game-theoretic vocabulary selection via the shapley value and banzhaf index. In Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, pages 2789–2798, 2021. [35] Judea Pearl. Causal inference in statistics: An overview. Statistics Survey, 2009. [36] Jonas Peters, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. 2014. [37] Benedek Rozemberczki, Lauren Watson, Péter Bayer, Hao-Tsung Yang, Olivér Kiss, Sebastian Nilsson, and Rik Sarkar. The shapley value in machine learning. pages 5572–5579. International Joint Conferences on Artificial Intelligence Organization, 2022. [38] Rajen D. Shah and Jonas Peters. The Hardness of Conditional Independence Testing and the Generalised Covariance Measure. The Annals of Statistics, 48(3), 2020. [39] Rajen D. Shah and Richard J. Samworth. Variable selection with error control: another look at stability selection. Journal of the Royal Statistical Society Series B: Statistical Methodology, 75(1):55–80, 2013. [40] Lloyd S. Shapley. A value for n-person games. Contributions to the Theory of Games, 2(28):307– 317, 1953. [41] Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search. 2001. [42] Alexander Statnikov, Jan Lemeir, and Constantin F. Aliferis. Algorithms for discovery of multiple Markov boundaries. The Journal of Machine Learning Research, 14(1):499–566, 2013. [43] Mukund Sundararajan and Amir Najmi. The many Shapley values for model explanation. In International conference on machine learning, pages 9269–9278. PMLR, 2020. [44] Wesley Tansey, Victor Veitch, Haoran Zhang, Raul Rabadan, and David M. Blei. The Holdout Randomization Test for Feature Selection in Black Box Models. Journal of Computational and Graphical Statistics, 31(1):151–162, 2022. [45] Jacopo Teneggi, Beepul Bharti, Yaniv Romano, and Jeremias Sulam. SHAP-XRT: The Shapley Value Meets Conditional Independence Testing. Transactions on Machine Learning Research, 2023. [46] Sandhya Tripathi, N. Hemachandra, and Prashant Trivedi. Interpretable feature subset selection: A Shapley value based approach. In 2020 IEEE International Conference on Big Data (Big Data), pages 5463–5472, 2020. [47] Sandhya Tripathi, N. Hemachandra, and Prashant Trivedi. On feature interactions identified by Shapley values of binary classification games. stat, 1050:12, 2020. [48] Caroline Uhler, Garvesh Raskutti, Peter Bühlmann, and Bin Yu. Geometry of the faithfulness assumption in causal inference. The Annals of Statistics, pages 436–463, 2013. [49] Isabella Verdinelli and Larry Wasserman. Decorrelated variable importance. Journal of Machine Learning Research, 25(7):1–27, 2024. [50] Haohan Wang, Benjamin J Lengerich, Bryon Aragam, and Eric P Xing. Precision Lasso: accounting for correlations and linear dependencies in high-dimensional genomic data. Bioinformatics, 35(7):1181–1187, 2019. 23
[51] Brian Williamson and Jean Feng. Efficient nonparametric statistical inference on population feature importance using Shapley values. In International conference on machine learning, pages 10282–10291. PMLR, 2020. [52] Brian D. Williamson, Peter B. Gilbert, Noah R. Simon, and Marco Carone. A general framework for inference on algorithm-agnostic variable importance. Journal of the American Statistical Association, 118(543):1645–1658, 2023. [53] Siphendulwe Zaza, Marcellin Atemkeng, and Sisipho Hamlomo. Wine Feature Importance and Quality Prediction: A Comparative Study of Machine Learning Algorithms with Unbalanced Data. volume 566, pages 308–327, 2024.
24
A
Appendix - Theorem Related Proofs
A.1
Proof of Theorem 2
Proof. By Theorem 1, we know that under the faithfulness assumption and there is no data leakage assumption, we are performing the following hypothesis testing: π π H0 : min V I0,j =0⇔Y ⊥ ⊥ Xj |X−j ; Ha : min V I0,j > 0 ⇔ Y ̸⊥ ⊥ Xj |X−j π∈Π(p)
π∈Π(p)
Let’s randomly select K permutations, and recall that we define Π = (π1 , π2 , ..., πK ), and πk V In,j , k ∈ [K] as the estimated variable importance of feature j under permutation πk . Define tj as the thresholds we use to decide if we reject the null hypothesis. Define an event A = πk e ={Among the K randomly selected {For a given order of permutation πk : V In,j ≤ tj }, and event A πk order of permutations, there exists at least one order πk , k ∈ [K] s.t. V In,j ≤ tj }. e = {B1 , B2 , ..., BI } be the Suppose we have p features, and there are |S0∗ | significant features. Let B collection of minimal feature subset of feature j, where each Bi ⊂ [p]\{j} s.t. V (Bi ∪{j})−V (Bi ) = 0. Then, define the event E(Bi ) = {feature j is placed after all features in Bi }, we have p 1 |Bi |+1 · |Bi |! = P(E(Bi )) = p |Bi | + 1 |B |+1 · |Bi + 1|! i
e = P( P(A) = P(E(B))
I [
E(Bi )) ≥ max P(E(Bi )) = Bi ∈B̃
i=1
1 1 ≥ ∗ minBi ∈B̃ |Bi | + 1 |S0 | + 1
πk Then, the probability that at least one out of K independent permutations yielding V In,j = 0, k ∈ [K] is: 1 e = 1 − (1 − P(A))K ≥ 1 − (1 − P(A) )K ∗ |S0 | + 1
Thus, for any fixed ϵ ∈ (0, 1), if K is sufficiently large with order K = O(|S0∗ |) then we have e > 1 − ϵ. Finally, let’s justify that the test statistic V I (1) is under Type I error control: P(A) n,j (1) (1) e · P(A) e P(V In,j > tj |H0 ) = P(V In,j > t|H0 , A) 1 2(1) ≤ exp( λ2 σn,j − λt), by Chernoff Inq 2
To give a tight bound, we optimize λ by minimizing the quadratic equation, and we get λ̂ = Then, we have (1)
P(V In,j > tj |H0 ) ≤ exp
1 t2 t2 − 2 σ 2,(1) σ 2,(1) n,j
= exp −
n,j
t2
1 <α 2 σ 2,(1) n,j
q 2,(1) Thus, we have threshold tj > −2 ln(α)σn,j .
25
t 2,(1) . σn,j
A.2
Type II error control
To provide some further intuition of maximizing power, let’s define an event C = {For a given order πk of permutation πk : V In,j > tj }. (1)
Type II error = 1 − P(V In,j > tj |Ha ) = 1 − P(
K \
πk V In,j > t|Ha )
k=1
=1−
K Y
πk P(V In,j > t|Ha )
k=1
= 1 − P(C|Ha )K This result suggests that more permutations drive the Type II error down when P (C|Ha ) is close to 1. If P (C|Ha ) ∈ (0, 1), fewer permutations K yields better Type II error control.
A.3
Proof of Corollary 4
The MinShap uses concentration inequality to generate threshold whereas Max-p uses normal approximation and multiple hypothesis testing. Proof. Under the union of K null hypothesis, there is at least one true null. The proof for the equivalence that Max-p is also under Type I error control is easy to see, since (K)
(1)
P(pn,j < α|H0 ) = P(|Zn,j | > t) (k)
≤ P(|Zn,j | > t) = α (k)
where Zn,j =
B
(k)
V In,j (k) σn,j
, k ∈ [K] and t = Φ−1 (1 − α/2).
Appendix - Algorithms
This section includes two algorithms: the first one is for Monte carol permutation sampling approximation of Shapley value based on feature importance context (see Algorithm 2), and the second one is the Max-p related p-value algorithm (see Algorithm 3). Algorithm 2 Monte Carlo Permutation Sampling Approximation of Variable-Importance Shapley values Require: Dataset (X, Y ); measure V ; number of features p; number of permutations K 1: Fit null model fn,∅ ; set V∅ ← V (fn,∅ , Pn ) 2: Initialize ϕn,j ← 0 for all j ∈ [p] 3: for k = 1, . . . , K do 4: Sample a random permutation πk 5: Set Vcur ← V∅ 6: for j = 1, . . . , p do πk 7: Set predecessor set: Pj−1 ← { i ∈ [p] : πk (i) < πk (j) } 8: Fit model fn,P πk ∪{j} and compute Vnew ← V fn,P πk ∪{j} , Pn j−1
j−1
πk V In,j ← Vcur − Vnew πk 10: ϕn,j ← ϕn,j + V In,j 11: Vcur ← Vnew 12: end for 13: end for 14: Return ϕn,j ← ϕn,j /K for all j ∈ [p]
9:
26
Algorithm 3 Max-p and adjusted p-value testing for feature selection Require: Dataset (X, Y ); measure V ; features p; number of permutations K; significance level α Fisher , p̃Stou and pMax for each feature j Ensure: Adjusted p-values p̃Bon u , p̃u u 1: Fit null model fn,∅ ; set V∅ ← V (fn,∅ , Pn ) and residuals e2∅ 2,πk πk 2: Allocate arrays ϕn,j and σn,j for j ∈ [p], k ∈ [K] 3: for k = 1, . . . , K do 4: Sample permutation πk 5: Set Vcur ← V∅ and e2cur ← e2∅ 6: for j = 1, . . . , p do πk 7: Pj−1 ← { i ∈ [p] : πk (i) < πk (j) } 8: Fit fn,P πk ∪{j} ; compute Vnew and residuals e2new j−1
k ϕπn,j ← Vcur − Vnew 2,πk 10: σn,j ← V ar(e2cur − e2new )/n 11: Vcur ← Vnew ; e2cur ← e2new 12: end for 13: end for 14: for j = 1, . . . , p do 2,πk K k k K 15: Convert {ϕπn,j , σn,j }k=1 to z-scores and p-values {pπn,j }k=1
9:
16:
(1)
(K)
(1)
(K)
Sort pn,j ≤ · · · ≤ pn,j and |z|n,j ≤ · · · ≤ |z|n,j
(u)
Bonferroni-Holm: for u = 1, . . . , K, set pBon ← (K − u + 1)pn,j u Bon Bon Bon Bon ; for u = 2, . . . , K set p̃u ← min max(p̃Bon 18: p̃1 ← p1 u−1 , pu ), 1 P (k) Fisher ← 1 − 19: Fisher-Holm: for u = 1, . . . , K, compute TuFisher ← −2 K k=u log pn,j and pu (TuFisher ) Fχ2 17:
2(K−u+1)
20:
Fisher ) p̃Fisher ← pFisher ; for u = 2, . . . , K set p̃Fisher ← max(p̃Fisher u 1 1 u−1 , pu
21:
Stouffer-Holm: for u = 1, . . . , K, set zuStou ←
22: 23: 24:
PK−u+1 (k) |z|n,j k=1 √ ← 2 1 − Φ(|zuStou |) and pStou u K−u+1 Stou ) p̃Stou ← pStou ; for u = 2, . . . , K set p̃Stou ← max(p̃Stou u 1 1 u−1 , pu (K) Max-p: pMax ← pn,j Decision: reject H0,j if the selected p-value < α
25: end for
C
Appendix - Simulation
C.1
Experiment on neural network
Our modified Shapley value based variable selection procedure is model-agnostic and thus not restricted to tree-based algorithm. In addition to XGBoost, we also evaluate the method using a two-layer ReLU network as the predictive model for comparison. Since neural network training is more computationally intensive, we run these experiments at a reduced scale, decreasing the feature dimension from p = 20 to p = 10 with only the first 5 features significant in all 4 models, and the sample size from n = 3000 to n = 1000. The number of simulation is 50 and the number of permutation K is 10 for MinShap and Max-p value approaches. The overall performance comparison is presented in Fig.10. MinShap and Max-p value approaches are much better than GCM, achieving higher accuracy and F1 score in all four models. Although the F1 score and accuracy for LOCO are similar to that of MinShap and Max-p in model (a) and (d), LOCO fails to control false positives; empirically, its Type I error is worse than that of a random selector with probability 0.5, and false discovery rate is higher than all other methods across all models. From Jaccard index Table 6, we see that the stability of MinShap and Max-p are higher 27
Conditional Interaction
Logistic
0.54
0.39
0.39
0.39
0.39
0.49
0.49
0.49
0.49
O
M in
O
C
ax
LO
M
C G
La s
M in
C
ax
LO
M
M C G
M in
La s
O C
ax
LO
M
C G
Sh ap
0.54
−p
0.54
M
0.54
so
0.56
Sh ap
0.56
−p
0.56
so
0.56
Sh ap
0.53
−p
0.53
M
0.53
so
0.53
La s
Sh ap
O
−p
M in
C
Non−Linear
ax
LO
M
G
C
M
1.00 0.75 0.50 0.25 0.00
so
0.6 0.4 0.2 0.0
La s
Accuracy F1 Recall
1.00 0.75 0.50 0.25 0.00
FDR
Performance
1.00 0.75 0.50 0.25 0.00
Type.I.Error
Linear 1.0 0.8 0.6 0.4
Method Method
GCM
Lasso
LOCO
Max−p
MinShap
Figure 10: Average performance measure with standard deviation for model (a)-(d) for feature selection comparison using Neural Network. We see that MinShap and Max-p outperform other feature selection methods and have well-controlled Type I error in the meanwhile. Table 6: Jaccard Index for the stability of feature selection methods on model (a)–(d) using Neural Network. Feature Selection Method Model
MinShap
Max-p
LOCO
GCM
Lasso
(a) (b) (c) (d)
0.99 0.91 0.57 0.89
0.98 0.89 0.49 0.87
0.73 0.72 0.94 0.72
0.79 0.72 0.38 0.80
1.0 0.96 1.0 1.0
than LOCO and GCM in models (a), (b) and (d). The selection stability of LOCO in model (c) is high, but its inflated false positives undermine reliability. Similarly, Lasso remains stable, but it is not reliable in all nonlinear settings with low F1 score.
D
Additional Real Data Analysis
In this section, we still mainly apply the variable selection methods of MinShap, Max-p, GCM and LOCO to additional real dataset. The dataset is analyzed using 5-fold cross-validation. Variable selection is carried out within each training fold, and the mean squared prediction error on the corresponding test fold is used to evaluate predictive performance. In addition, the stability of the feature selections is measured by Jaccard index.
D.1
Diabete data
Understanding the link between clinical and biochemical measurements to disease progression is essential to medical decision-making. We aim to predict how patients’ characteristics can be used to predict health outcomes with diabete dataset which is available in the sklearn Python library. The data has ten base features, age, sex, body mass index (BMI), average blood pressure (BP), and six blood serum measurements. Data is obtained for each of 442 diabetes patients, as well as 28
Table 7: Comparison of methods in terms of prediction error and selection stability on diabetes data under two models. Model
XGBoost
Random forest
Method MinShap Max-p LOCO GCM MinShap Max-p LOCO GCM
MSE 3455.04 3431.36 3310.99 3228.06 3708.00 3429.36 3437.57 3294.31
se(MSE) 112.99 124.33 163.90 130.00 228.36 137.94 178.67 133.62
Jaccard index 0.80 0.73 0.63 0.81 0.70 0.80 0.47 0.85
Selected Features [‘bmi’, ‘s5’] [‘bmi’, ‘s5’] [‘bmi’, ‘bp’, ‘s5’, ‘sex’] [‘bmi’, ‘bp’, ‘s5’, ‘sex’] [‘bmi’, ‘s5’] [‘bmi’, ‘s5’] [‘bmi’, ‘bp’] [‘bmi’, ‘bp’, ‘s2’, ‘s3’, ‘s5’, ‘sex’]
the response of interest, a quantitative measure of disease progression one year after the time of measurement of the base features. From Table 7, we see that GCM always selects more variables than other methods, and the selection stability of LOCO is inferior to others. Features selected by MinShap and Max-p, "bmi" (body mass index) and "s5" (log of serum triglycerides level), using XGBoost and random forest are the most commonly selected in previous studies [3, 11].
29