# Hierarchical Variable Importance with Statistical Control for Medical Data-Based Prediction

Joseph Paillard, Antoine Collas, Denis A Engemann, Bertrand Thirion

► To cite this version:

Joseph Paillard, Antoine Collas, Denis A Engemann, Bertrand Thirion. Hierarchical Variable Im-
portance with Statistical Control for Medical Data-Based Prediction. IPMI 2025 - Information Pro-
cessing in Medical Imaging, May 2025, Kos Island, Greece. pp.79-93, 10.1007/978-3-031-96625-5_6.
hal-05207600

HAL Id: hal-05207600

https://hal.science/hal-05207600v1

Submitted on 12 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hierarchical Variable Importance with Statistical Control for Medical Data-Based Prediction

Joseph Paillard¹,², Antoine Collas², Denis A. Engemann¹, and Bertrand Thirion²

¹ Roche Pharma Research & Early Development, F. Hoffmann-La Roche Ltd, Basel, Switzerland

² Université Paris-Saclay, Inria, CEA, Paris, Palaiseau, France
Correspondance: joseph.paillard@roche.com

**Abstract.** Recent advances in machine learning have greatly expanded the repertoire of predictive methods for medical imaging. However, the interpretability of complex models remains a challenge, which limits their utility in medical applications. Recently, model-agnostic methods have been proposed to measure conditional variable importance and accommodate complex non-linear models. However, they often lack power when dealing with highly correlated data, a common problem in medical imaging. We introduce Hierarchical-CPI, a model-agnostic variable importance measure that frames the inference problem as the discovery of groups of variables that are jointly predictive of the outcome. By exploring subgroups along a hierarchical tree, it remains computationally tractable, yet also enjoys explicit family-wise error rate control. Moreover, we address the issue of vanishing conditional importance under high correlation with a tree-based importance allocation mechanism. We benchmarked Hierarchical-CPI against state-of-the-art variable importance methods. Its effectiveness is demonstrated in two neuroimaging datasets: classifying dementia diagnoses from MRI data (ADNI dataset) and analyzing the Berger effect on EEG data (TDBRAIN dataset), identifying biologically plausible variables.

**Keywords:** Statistics, neuroimaging, interpretable machine learning

## 1 Introduction

Within the field of medical imaging, machine learning holds great promise to facilitate prediction of clinical outcomes, see e.g. [1, 2, 3, 4, 5, 6]. However, these advances have also opened major interpretability challenges. A key issue is how to infer the importance of features from prediction models going beyond ordinary least squares to accommodate a large number of predictors and represent non-linear associations between features and outcomes. Therefore, developing methods to measure variable importance in a model-agnostic manner is critical in order to obtain clinical insights and develop biomarkers, for instance, using brain images for the diagnosis of Alzheimer Disease (AD) based on existing co-horts, or data from clinical trials [7, 8, 2, 6]. However, to develop trustworthy
methods, it is essential to understand their theoretical guarantees, particularly concerning the risk of making false discoveries, which can be captured by the Family-Wise Error Rate (FWER) (see e.g. [9, 10]). Only few variable importance methods give access to such guarantees. Moreover, we focus here on *conditional importance*, meaning the importance measure whether a variable is *directly predictive* of the outcome, without being explained away by other variables [11, 12]. Such conditional importance is needed to establish that a marker carries independent information about the outcome, rather than merely reflecting distributed factors that are also present in other variables. Conditional importance analysis is particularly difficult in datasets that exhibit strong correlation structures such as image- or genomics-based biomarkers, or health data that reflect common latent factors [13].

We assess the face validity of the approach with two tasks that have been extensively studied in the literature—the effect of AD on structural MRI and the Berger effect on electroencephalography (EEG)—to allow for a form of confirmation, addressing the challenge posed by the absence of ground truth in variable importance methods.

## 1.1 Related Work

This work focuses *global* variable importance, as opposed to local variable importance methods such as *LIME* [14] or *SHAP* [15]. Global variable importance is estimated using methods such as global sensitivity analysis [11] or the popular Leave One Covariate Out (LOCO) approach [16, 17]. These methods can accommodate different types of learners, taking advantage of advances in machine learning to measure importance in complex and nonlinear models [17, 18]. Similarly, conditional permutation approaches have been shown to estimate a quantity equivalent to LOCO at a cheaper computational cost and with a faster convergence rate [19]. These methods have in common to provide a good control of the type-1 error rate, that is considering a null variable (or group) as important. However, all approaches suffer from an inherent limitation: conditional importance decreases as correlation increases. For instance, considering two random normal variables $X_1, X_2$ with correlation $\rho$ in a simple linear model $y = \beta_1X_1 + \beta_2X_2$ the importance of $X_1$ decreases proportionally to $(1 - \rho^2)$.

To mitigate this issue, methods based on variables grouping have been proposed to identify groups of highly-correlated, hence indistinguishable variables that predict the outcome [12, 13]. Variable grouping can be performed based on prior knowledge about the data or by using clustering techniques. While this effectively increases the statistical power by averaging correlated variables, it also reduces the precision in the sense that error control only holds at the group level. When performing the grouping in a data-driven way, choosing the clustering scheme and parameters has a critical impact yet has no obvious solution. A line of work relying on linear models and agglomerative clustering has been proposed in that direction with applications to genomics data [9, 20, 10]. Agglomerative clustering offers a compelling solution as it naturally explores different groupings at various resolutions along the hierarchical tree learned from
the data. However, this approach relied on Lasso regression and would conse-
quently limit the user in the choice of the model used to predict the outcome of
interest from the variables.

Another popular model explainability approach is Shapley Additive Global importance (SAGE) [21]. Based on Shapley values, this approach estimates an importance score for a given variable by conditioning on all subgroups that include this variable. This procedure provides a more nuanced view than strict conditional importance, because it decomposes additively the variance explained by the model into variable importance. However, it suffers from two main limitations. First, as an aggregated statistic, it obscures the role of variables in the prediction function [18], and is unable to distinguish between a predictive variable and another, non-predictive variable yet correlated with a predictive one. Second, the exploration of all submodels comes with an exponential explosion of computation cost, making this approach intractable. While implementations rely on Monte-Carlo sampling instead of exhausting the full combinatorial sum, the number of sampling steps needed to obtain accurate estimates still leads to intractable computation costs [18]. This two limitations are clearly visible in our experiments in Figure 2.

Our contributions are *i*) to introduce Hierarchical-CPI, a model-agnostic vari-
able importance measure that improves FWER control; it explores subgroups
in a tree-guided manner, using agglomerative clustering to provide more infor-
mation than variable-level importance while remaining tractable; *ii*) to present
an approach that enforces importance conservation through downstream impor-
tance allocation strategy, addressing the issue of vanishing importance under
high correlation.

2 Methods

*Notations:* We denote $X$ as the variables, $y$ as the outcome, and $\mu$ as the predictive model. $X_G$ represents the set of variables belonging to a group $G$, and $X_{-G}$ denotes the set of variables in the complement of $G$. The importance of a group is denoted as $\psi_G$. We use $S^*$ to denote the support (or set of active variables) and $S_0$ for the set of null variables. In the hierarchical tree defined by the clustering, $P$ refers to a parent node, and $L$ and $R$ refer to its left and right child nodes, respectively.

## 2.1 Hierarchical-CPI

We present a method for measuring variable importance while conditioning on
others, with conditioning sets taken in a tree-organized hierarchical representa-
tion of the variables. It balances *precision* and *statistical power*; precision refers
to extracting the information located in groups of variables that are as small
as possible; statistical power is achieved by considering condition sets different
from the set of all variables. The proposed method builds on top of Conditional
**Algorithm 1 Hierarchical CPI**

**Input:** K: number of folds, μ: predictive model, ν: imputation model, (X, y): data

1: tree ← Fit hierarchical clustering on X

2: **for** k in [1, ..., K] **do**

3: &mu;&#x0302;<sub>k</sub> ← Fit using (X<sub>train</sub>, y<sub>train</sub>) // fit the full model

4: **for** node in tree **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;G ← traversal(node) // search variables belonging to the node

6: &nbsp;&nbsp;&nbsp;&nbsp;&#x0302;<sub>G</sub><sup>k</sup> ← E[X<sub>G</sub><sup>train</sup>|X<sub>-G</sub><sup>train</sup>] // estimate the conditional distribution

7: &nbsp;&nbsp;&nbsp;&nbsp;&#x0302;<sub>G</sub><sup>test</sup> &sim; &#x0302;<sub>G</sub><sup>k</sup>(X<sub>-G</sub><sup>test</sup>) // sample from the conditional distribution

8: &nbsp;&nbsp;&nbsp;&nbsp;&#x0302;<sub>G</sub><sup>k</sup> ← L(y, μ(<span style="text-decoration: overline;">&#x0302;<sub>G</sub><sup>test</sup></span>)) - L(y, μ(X<sup>test</sup>)) // compute variable importance

9: **end for**

10: **end for**

11: p<sub>G</sub> ← t-test(<span style="text-decoration: overline;">ψ<sub>G</sub><sup>1</sup></span>, ..., <span style="text-decoration: overline;">ψ<sub>G</sub><sup>K</sup></span>) // compute p-value over folds

12: p<sub>G</sub><sup>h</sup> ← max<sub>G&isin;D</sub> p<sub>D</sub><sup>h</sup> // hierarchical adjustment

13: **return** p<sub>G<sub>i</sub></sub><sup>h</sup> for i = 1, ..., 2p - 1

Permutation Importance (CPI) [12] which, given a group of variables G, a model μ and loss L estimates the conditional importance

$$ \psi_G = \mathcal{L}(y, \mu(\tilde{X}_G)) - \mathcal{L}(y, \mu(X)), \qquad (1) $$

where $\tilde{X}_G$ is obtained by substituting into the group $X_G$ variables sampled from the conditional distribution ($X_G|X_{-G}$) and leaving the $X_{-G}$ variables unchanged. In brief, CPI quantifies the loss increase when conditioning on all other variables than those in G. This approach estimates the well known total Sobol index [11]. In addition, hierarchical-CPI leverages Ward's minimum variance method for agglomerative clustering to learn the hierarchical group structure [22]. The proposed method is presented in Algorithm 1, for a problem with p variables, it consists in estimating the conditional permutation importance of each group of variables within the hierarchical structure. Empirical importance values are obtained in a K-fold cross-validation scheme, yielding K estimates per group, ($\psi_G^1, ..., \psi_G^K$). A p-value $p_G$ is then derived based on a one-sample t-test. Finally, the node-level p-values $p_G$ are hierarchically adjusted by,

$$ p_G^h = \max_{G \subseteq D} p_D, \qquad (2) $$

to enforce that the p-value of a node is larger than the p-value of its parent.

## 2.2 Hierarchical CPI achieves FWER control

In this section, we demonstrate that the hierarchical-CPI approach controls the FWER under assumptions of estimator optimality and regularity. The assumptions (A.1, A.2, B.1, B.2) from [17] stipulate that the estimator μ must be optimal
and exhibit sufficient regularity. These assumptions have been validated by independent work and are considered not too restrictive [18, 17, 19]. We refer to a tree cut as a set of non-overlapping nodes within a hierarchical tree. Let $S_0$ denote the set of groups that only contain null variables, and let $\hat{S}_\alpha = \{G | p_G \le \alpha\}$ be the estimated set of active variables for a given significance level $\alpha \in [0, 1]$. The following result holds.

**Theorem 1.** Under the assumption that the conditions (A.1, A.2, B.1, B.2) stated in [17] on $\mu$, for any significance level $\alpha \in [0, 1]$ the multiplicity corrected p-values $\tilde{p}_G^h = \min(1, C \cdot p_G^h)$, with $C = p$ control the family-wise error rate at level $\alpha$, i.e. $\mathbb{P}(S_0 \cap \hat{S}_\alpha \neq \emptyset) \le \alpha$. Where G is a node of a tree cut.

*Proof.*

$$ \mathbb{P}(S_0 \cap \hat{S}_\alpha \neq \emptyset) = \mathbb{P}\left(\min_{G \subseteq S_0} \tilde{p}_G^h \le \alpha\right) = \mathbb{P}\left(\bigcup_{G \subseteq S_0} p_G^h \le \frac{\alpha}{C}\right) $$

Then, given Boole's inequality,

$$ \begin{aligned} \mathbb{P}(S_0 \cap \hat{S}_\alpha \neq \emptyset) &\le \sum_{G \subseteq S_0} \mathbb{P}\left(p_G^h \le \frac{\alpha}{C}\right) \\ &\le C \cdot \mathbb{P}\left(p_G^h \le \frac{\alpha}{C}\right) \end{aligned} $$

Since a tree cut contains less than $p$ nodes, $\mathrm{card}(\{G_i | G_i \subseteq S_0\}_{i \in [1, C]}) \le C$. Furthermore, given that $p_G^h = \max_{G \subseteq D}\{p_D\}$, where the maximum is taken over ancestor nodes, it comes that $\mathbb{P}(p_G^h \le \alpha) \le \mathbb{P}(p_G \le \alpha)$. Finally, under assumptions (A.1, A.2, B.1, B.2), it has been shown in [19] that, $\forall G \subseteq S_0$, $\mathbb{P}(p_G \le \alpha) \le \alpha$. We then have, $\mathbb{P}(S_0 \cap \hat{S}_\alpha \neq \emptyset) \le \alpha$ which completes the proof.

While this result holds when considering inference at the variable level, it is more general and applies to any node of the tree. Hierarchical CPI allows to learn a tree structure from the data and make inference at different levels.

## 2.3 Importance conservation to prevent importance vanishing

A common pitfall of conditional importance is that it vanishes under high correlation: For a parent node $P$ with two strongly correlated children nodes $L$ and $R$, then we have that $\psi_L^k + \psi_R^k \ll \psi_P^k$ for all $k$ in $[1..K]$. This effect is illustrated in Figure 1 a and c. Hierarchical CPI can infer that $P$ is important, but will give little to no importance to $L$ and $R$ (and all downstream groups). To obtain a Shapley-like additive decomposition of the model fit, inspired by variance partitioning ideas [23], we introduce a transfer mechanism that ensures additivity at the node level, meaning that the importance of a parent node is split into the sum of the importance of its children. This is meant to allow a more refined allocation of the importance budget compared to traditional clustering methods.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Importance conservation prevents the importance from vanishing as the group size decreases in a high-correlation setting. Example using simulated data with $n = 300$ samples and $d = 24$ variables. The data is generated by blocks, each corresponding to an AR(1) with autocorrelation parameter $\rho_{max} = 0.95$. a and b show the same dendrogram obtained through Ward's clustering. Each node's color encodes the conditional importance of the variables it contains. Without importance conservation, importance quickly vanishes down the tree. c and d present the p-value distributions, demonstrating that both methods accurately rank important variables.</figcaption></figure>

For a node $R$, let $\tilde{\psi}_R^k$ be the corrected value for $\psi_R^k$ that ensures importance conservation through the hierarchical structure. Then, $1_R(\epsilon)$ is the indicator function equal to one when $\psi_R/\hat{\sigma}_R \ge \epsilon > 0$, where $\hat{\sigma}_R$ is the standard deviation of the right node importance estimated over the k-folds and zero otherwise. Importance conservation aims at satisfying the equation $\tilde{\psi}_P^k = \psi_L^k + \psi_R^k$. This condition ensures that the importance of the top node, which measures the full model's importance, is allocated to its children nodes. The allocation mechanism for the child node $L$ with sibling $R$ and parent $P$ proceeds as follows:

$$ \tilde{\psi}_L^k = \begin{cases} \psi_L^k + 1_R(\epsilon) \frac{\tilde{\psi}_P^k - \psi_L^k - \psi_R^k}{2} & \text{if } 1_L(\epsilon) = 1 \\ \tilde{\psi}_P^k \frac{\psi_L^k}{\psi_R^k + \psi_L^k} (1 - 1_R(\epsilon)) + 1_R(\epsilon) (\tilde{\psi}_P^k - \psi_R^k) & \text{if } 1_L(\epsilon) = 0 \end{cases} \quad (3) $$

This equation distributes the parent's importance proportionally to the children's importance when both nodes' importance values are greater than a threshold $\epsilon$. If the importance $\psi_L$ of node $L$ is smaller than $\epsilon$, it remains unchanged, avoiding false positives (not all children of an important parent are important). When both children have sub-threshold importance, indicating that their correlation leads to mutual importance cancellation, the importance of the parent is allocated equally between the two nodes.

# 3 Results

## 3.1 Control of Family-Wise Error Rate on Simulated Data

This experiment benchmarks the hierarchical-CPI approach described in Algorithm 1 with other state of the art variable importance methods on simulated
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Hierarchical CPI empirically controls the FWER while high statistical power. Results from simulated data with $n = 400$ samples and $p = 124$ variables sampled from a normal distribution with block correlations described in subsection 3.1. The results present a summary of 100 repetitions of the simulation. **a** and **d** present the AUC for important/non-important variable classification as a function of correlation and SNR. Error bars represent 95% confidence interval. **b** and **c** show the evolution of the FWER at level $\alpha = 0.05$ with Bonferroni correction. The top row explores varying $\rho_{max}$ values at a fixed SNR=2, while the bottom row examines varying SNRs at a fixed $\rho_{max} = 0.9$. **c** and **f** shows the average computation time taken by each method over the simulations for the linear and non-linear scenario respectively.</figcaption></figure>

data. The data is generated by blocks, each corresponding to an AR(1) with autocorrelation parameter $\rho_{max}$. The outcome is modeled using two different scenarios. The first is a linear model with additive noise $y = X\beta + \sigma_N\varepsilon$, represented with dotted lines in Figure 2. The support $S^* = \{j; |\beta_j \neq 0\}$ is kept sparse, with $|S^*|$ set to either 5 or 10, and coefficients values sampled from the set $\beta_j \in \{-2, -1, 1, 2\}$ with uniform probability. The second is a non linear scenario presented in [19]: $y = X_{j1} + 2\log(1 + 2X_{j2}^2 + (X_{j3} + 1)^2) + X_{j4}X_{j5} + \sigma_N\varepsilon$. The support $S^* = \{j_1, \dots, j_5\} \in [0, p]^5$ is randomly sampled at each simulation run. In both cases, additive noise $\varepsilon \sim N(0, I_n)$, controls the signal-to-noise ratio (SNR), defined as $SNR = ||y^*||_2^2/\sigma_N^2||\varepsilon||_2^2$, where $y^*$ is the noiseless outcome. The SNR is a simulation parameter. In the experiments shown in Figure 2, we used $p = 124$ variables grouped into five blocks of correlated features with respective sizes 4, 8, 16, 32, and 64. The number of samples was fixed at $n = 400$ to match the dimensionality commonly found in medical imaging applications. To solve the non-linear regression task, a multilayer perceptron (MLP) with 100 hidden units was used. It was trained using Adam optimizer for 400 epochs with early stopping (patience 10 epochs). For the linear scenario, a Ridge-regularized model was used. Its regularization parameter was learned via nested cross-validation. For CPI and HCPI, the loss used in Equation 1 is the root mean squared error.
The top row (a, b, c) explores the influence of $\rho_{max}$ while the bottom row the influence of the SNR. We considered two metrics. First, the Area Under the Receiver Operating Characteristic Curve (AUC): This compares the predicted importance to the true importance. It aims to assess each method's ability to recover the true support. Second, the FWER: This measures the probability of making at least one false discovery, estimated over 100 simulation repetitions. We compare three methods: CPI, Hierarchical-CPI (HCPI) and SAGE. For CPI and HCPI, the predicted importance correspond to 1-p-value and the estimated support is $\hat{S}_{\alpha} = \{j | p_j \le \alpha\}$, we considered a level $\alpha = 0.05$. Regarding SAGE, we used a publicly available implementation<fnref n="3" /> which provides estimated standard deviation from which 95% confidence interval were derived. The estimated support for SAGE consisted of variables for which the confidence interval did not include 0.

As shown in Figure 2, the hierarchical approach effectively controls the FWER even in challenging simulation settings, e.g. with very high correlation or low SNR. This can be attributed to the hierarchical adjustment, described in Equation 2, that bounds a node's p-value below using the p-value of its parent. As shown in subpanels a and d of the Figure 2, this additional control does not decrease the power of the method when compared to CPI. While the exploration of the nodes entails an additional computation cost, it remains of the same order as CPI, which is a fast method. Indeed, for a hierarchical clustering of p variables, the total number of nodes is $2^p - 1$ which makes the computation scale linearly with the dimension instead of the exponential explosion inherent to SAGE [21]. This fact is illustrated on the panels c and f of Figure 2, where the logarithmic axis illustrates the untractable computation time of SAGE. In the non-linear case where a neural network is used, this trend becomes even more pronounced.

## 3.2 Hierarchical CPI Identifies Characteristic Markers of AD

This study explores image-based diagnosis using the ADNI dataset [7]. Cohort selection was based on the availability of T1-weighted images, similarly to [8]. A total of 1616 patients were included: 760 controls (CN), 529 diagnosed with Mild Cognitive Impairment (MCI), and 327 with AD. Gray Matter (GM) density maps were computed using the *sMRIPrep* pipeline [24], which is part of the widely used *fMRIPrep* pipeline [25]. The mean GM densities were extracted from 116 Regions of Interest (ROIs) defined by the Automated Anatomical Labeling Atlas 3 [26] using *Nilearn* [27] and used as features for classification tasks. The high correlation between ROIs (Pearson correlation ranging from 0.32 to 0.94) and the interest in locating the pathology's impact precisely motivated the proposed approach. The methodology for data processing, model optimization, and hyper-parameter tuning followed [8], to ensure result reproducibility. We used a Support Vector Classifier (SVC) as implemented in *libsvm* on GM densities in the 116 regions of the AAL atlas. For the HCPI method, importance was measured by the

<footnote n="3">³ https://github.com/iancovert/sage</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Hierarchical CPI discovers groups of characteristic markers of AD progression. Importance obtained for classifying AD and MCI subjects from the ADNI dataset using a support vector classifier and grey matter densities in the 116 AAL regions. **a** Signed importance values obtained using the SAGE method. Only values for which the 95% CI did not overlap with 0 are reported. **b** Important regions identified by hierarchical CPI at the $\alpha = 0.05$ level. **c** Dendrogram derived from the hierarchical CPI approach. The important nodes ($\alpha = 0.05$) of the tree learned by hierarchical clustering are colored in red. Regions identified as important are labeled in red. **d** Important regions identified by HCPI for classifying patients with AD versus controls. **e** Important regions for classifying controls versus patients with MCI.</figcaption></figure>

hinge loss difference in Equation 1. All results are reported using 10-fold cross-validation with stratification. Hyper-parameter tuning was performed using a nested cross-validation loop to avoid information leakage from the test set [28]. Using a linear or *rbf* kernel led to similar predictive performance and importance scores. The results were reported for a linear kernel. Three classification tasks were considered: MCI vs. AD, AD vs. CN, and MCI vs. CN. The average AUC
on the test set over 10 folds were 0.78, 0.93, and 0.74, respectively, which is consistent with existing literature [8].

Figure 3 presents the importance maps at the individual feature resolution, computed using the SAGE method (a) and the proposed hierarchical-CPI (b) for classifying patients with AD and MCI. Similar to the previous section, for both methods, results are reported at a significance threshold of $\alpha = 0.05$. While SAGE identifies more regions as important, it is likely that many of these are only marginally, not conditionally, associated with the outcome. By contrast, HCPI identifies fewer regions that summarize the specific markers of the disease. Importantly, all regions identified by HCPI are also identified by SAGE, revealing a form of consistency. These regions include the hippocampi, and the orbitofrontal cortex (rectus in the AAL) which have been extensively described in the literature as areas where atrophy is substantial [29, 30]. The putamen and thalamus were also identified as predictive, consistently with published work documenting the association between AD and decreased global GM in these regions [31]. Moreover, HCPI allows to learn clusters of predictive variables with varying resolutions. This is depicted in c, which presents the dendrogram learned by agglomerative clustering, with nodes having a p-value below the significance threshold highlighted in red. This information complements panel b by highlighting the importance of selected subgroups. For instance, the node including (Caudate_L, Caudate_R) is important whereas individual variables are not, due to the high correlation between these two regions, (Pearson correlation of 0.84) leading to a cancellation of their conditional importance. In addition to these results Figure 3 presents the importance map for two additional tasks: AD vs CN (d) and MCI vs CN (e). Similar to Figure 3, the HCPI approach identifies hallmarks of AD pathology, such as the gray matter density in the hippocampi, which has been extensively described in the literature [29].

## 3.3 Importance Conservation Enables Inference on Highly Correlated EEG Data

The importance conservation approach was then applied to the EEG data from the TDBRAIN dataset to characterize the known Berger effect [32]. EEG data is known to exhibit high correlation due to latent sources spreading across the scalp as a result of field spreads. Resting state EEG were acquired from 1234 healthy subjects who were asked to open and close their eyes during labeled periods. The dataset was preprocessed using the pipeline presented in [33] in order to remove artifacts generated by non-brain sources. Specifically, independent component analysis was applied in order to remove eye-movement artifacts which would make the task trivial. The power at each of the 26 electrodes was computed across 17 logarithmically spaced frequency bands, ranging from 1 to 64 Hz, using Morlet wavelets. The 442 resulting features present a very high correlation structure, with minimum Pearson correlation above 0.9. We considered the task of classifying the eyes status (closed vs open) using a pipeline consisting of a logarithm computation followed by logistic regression. Similarly to the previous section, 10-fold nested cross validation was used. The loss used to measure the
<figure><img src="image_5.png" /><figcaption>**Fig. 4. Importance conservation enables variable importance inference in high-dimensional settings with very high correlation.** Comparison of the variable importance obtained using hierarchical CPI with importance conservation (HCPI-IC), without (HCPI) and using SAGE. Absolute SAGE values are represented for readability. **a** The distribution of importance over frequencies for significant variables at the $α = 0.05$ threshold. Each point represents the sum of important channels at a given frequency. **b, c** The distribution of importance over scalp topography for significant variables at the $α = 0.05$ threshold. At a given channel, the sum is taken over all important frequencies. **d** Presents the performance of sub-models that use only a fraction of variables identified as significant at a level $α$, with $α = 1$ corresponding to the full model. The boxes represent the distribution over 10-fold cross-validation. The dotted line represents the chance level.</figcaption></figure>

conditional importance in Equation 1 is the cross-entropy loss. The transmission threshold $ϵ$ was set to the 95% quantile of the normal distribution.

As illustrated by Figure 4 **a**, the high correlation in the data causes the conditional importance to vanish, resulting in no significant discoveries at a threshold of $α = 0.05$ (HCPI in red). To mitigate this issue and increase statistical power, the importance conservation mechanism (HCPI-IC, in blue) introduced in subsection 2.3 is applied. The importance however remains more focal than with SAGE (orange), which spreads the importance over a wide range of frequencies. Panels **a**, **b** and **c** show the distribution of importance in the frequency and sensor spaces among significant variables at a threshold of $α = 0.05$. The pattern observed, with most of the importance located around 10Hz at occipital electrodes, corresponds to the well-studied Berger effect, characterized by increased occipital activity in the alpha-band [34]. The pattern is precisely identified by HCPI-IC, while SAGE distributes more broadly the importance over electrodes and frequencies. Finally, panel **d** demonstrates the performance of submodels using only significant variables at thresholds $α < 0.05$, $α < 0.1$, and 1 (all variables). It reveals that at the strictest threshold ($α = 0.05$), the procedure selects 55 variables out of 442, recovering 96% of the full model's performance.

# 4 Discussion and conclusion

The HCPI approach was motivated by the challenge of making inference on high-dimensional and highly correlated neuroimaging data. To achieve this, it frames
the inference problem as the discovery of groups of variables that are jointly predictive of the outcome. It can recover statistical control in high correlations regimes where standard methods lose consistency. Statistical guarantees were empirically validated on simulated data, and the method was applied to two neuroimaging modalities using publicly available datasets. Its effectiveness was demonstrated on both classification and regression tasks. By successfully testing different tasks, models and losses, we proved the practical utility of this model-agnostic approach. HCPI flexibility exceeds that of existing methods relying on linear models or Lasso-based knockoffs [9, 13, 20].

By exploring subgroups within a learned hierarchical tree, HCPI balances precision and statistical power, allowing the identification of groups that are important, even if none of the individual variables is significant. It thus identifies the highest resolution at which importance can be narrowed down, without needing to optimize clustering parameters [13, 12, 20]. This information can easily be visualized using a dendrogram. Unlike additive methods like SAGE, which exhaustively explore all subgroups including a variable, eliminating many costly and useless evaluations. It provides a FWER control, thus contrasting with SAGE's known lack of type-1 error control [18]. Moreover, it remains tractable, requiring only $2p - 1$ importance evaluations. Finally, the importance conservation mechanism introduced in Equation 3 mitigates power loss due to vanishing importance, as shown with highly correlated EEG data features.

When tested on MRI and EEG data, our method identified biologically well-studied features consistent with existing literature such as hallmarks of AD in MRI and the Berger effect in EEG. Lastly, while we used Conditional Permutation Importance, because known to be more stable and efficient than LOCO, the latter could also be used as a drop-in replacement for estimating importance.

*Limitations:* We have explored scenarios where a single agglomerative clustering is performed, demonstrating that it can yield insightful learnings about data structure, with clusters of variables being predictive even if individual variables are not. However, this step can introduce randomness. For applications not requiring hierarchical tree learning, like voxel-level or applications to raw images—it may be beneficial to repeat the procedure and leverage p-value aggregation strategies or e-values [35, 36]. Future work could involve repeated agglomerative clustering on random data subsets, followed by aggregation to improve robustness. Another limitation concerns the theoretical guarantees of the importance conservation approach. While we empirically observed a type-1 error rate much lower than SAGE, a formal result remains to be established. The transmission threshold $\epsilon$ is critical in this context: it defines a threshold below which the importance of a parent node becomes indivisible because the contributions of its children nodes cancel each other out. We conjecture that it is possible to obtain guarantees for type-1 error control outside a neighborhood (which size depends on $\epsilon$) around the support.

The algorithm builds on open-source software available on Github<fnref n="4" />.

<footnote n="4">⁴ https://github.com/mind-inria/hidimstat</footnote>
**Acknowledgments.** This research has received funding from the H2020 Research In-
frastructures Grant EBRAIN-Health 101058516 and the VITE ANR-23-CE23-0016 and
PEPR Santé numérique, Brain health Trajectories ANR-22-PESN-0012 projects.
References

[1] Junhao Wen et al. "Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation". In: *Medical image analysis* 63 (2020), p. 101694.

[2] Duygu Tosun et al. "Identifying individuals with non-Alzheimer's disease co-pathologies: A precision medicine approach to clinical trials in sporadic Alzheimer's disease". In: *Alzheimer's & Dementia* 20.1 (2024), pp. 421-436.

[3] Lukas AW Gemein et al. "Machine-learning-based diagnostics of EEG pathology". In: *NeuroImage* 220 (2020), p. 117021.

[4] Rémi Cuingnet et al. "Automatic classification of patients with Alzheimer's disease from structural MRI: a comparison of ten methods using the ADNI database". In: *neuroimage* 56.2 (2011), pp. 766-781.

[5] Alexandre Abraham et al. "Machine learning for neuroimaging with scikit-learn". In: *Frontiers in neuroinformatics* 8 (2014), p. 71792.

[6] Adrian Tousignant et al. "Prediction of disease progression in multiple sclerosis patients using deep learning analysis of MRI data". In: *International conference on medical imaging with deep learning*. PMLR. 2019, pp. 483-492.

[7] Ronald Carl Petersen et al. "Alzheimer's disease Neuroimaging Initiative (ADNI) clinical characterization". In: *Neurology* 74.3 (2010), pp. 201-209.

[8] Jorge Samper-González et al. "Reproducible evaluation of classification methods in Alzheimer's disease: Framework and application to MRI and PET data". In: *NeuroImage* 183 (2018), pp. 504-521.

[9] Jacopo Mandozzi and Peter Bühlmann. "Hierarchical Testing in the High-Dimensional Setting With Correlated Variables". In: *Journal of the American Statistical Association* 111.513 (2016). ISSN: 0162-1459. DOI: 10.1080/01621459.2015.1007209.

[10] Bertrand Thirion Jérôme-Alexis Chevalier Joseph Salmon. "Statistical inference with ensemble of clustered desparsified lasso". In: *International Conference on Medical Image Computing and Computer-Assisted Intervention* (2018).

[11] I.M Sobol. "Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates". In: *Mathematics and Computers in Simulation* 55.1-3 (2001), pp. 271-280. DOI: 10.1016/s0378-4754(00)00270-6.

[12] Ahmad Chamma, Bertrand Thirion, and Denis Engemann. "Variable Importance in High-Dimensional Settings Requires Grouping". In: *Proceedings of the AAAI Conference on Artificial Intelligence* 38.10 (2024). DOI: 10.1609/aaai.v38i10.28997.

[13] Jérôme-Alexis Chevalier et al. "Decoding with confidence: Statistical control on decoder maps". In: *NeuroImage* 234 (2021), p. 117921. DOI: 10.1016/j.neuroimage.2021.117921.

[14] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. ""Why Should I Trust You?": Explaining the Predictions of Any Classifier". In: *Proceed-
ings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (2016). DOI: 10.48550/arxiv.1602.04938.

[15] Erik Štrumbelj and Igor Kononenko. "Explaining prediction models and individual predictions with feature contributions". In: *Knowledge and Information Systems* 41 (2014). DOI: 10.1007/s10115-013-0679-x.

[16] Toshimitsu Homma and Andrea Saltelli. "Importance measures in global sensitivity analysis of nonlinear models". In: *Reliability Engineering & System Safety* 52 (1996). DOI: 10.1016/0951-8320(96)00002-6.

[17] Brian D. Williamson et al. "A General Framework for Inference on Algorithm-Agnostic Variable Importance". In: *Journal of the American Statistical Association* 118 (2023). DOI: 10.1080/01621459.2021.2003200.

[18] Isabella Verdinelli and Larry Wasserman. "Feature importance: A closer look at shapley values and loco". In: *Statistical Science* 39.4 (2024), pp. 623-636.

[19] Ahmad Chamma, Denis A Engemann, and Bertrand Thirion. "Statistically Valid Variable Importance Assessment through Conditional Permutations". In: *Advances in Neural Information Processing Systems* 36 (2024).

[20] Tuan-Binh Nguyen, Jérôme-Alexis Chevalier, and Bertrand Thirion. "Ecko: Ensemble of clustered knockoffs for robust multivariate inference on fMRI data". In: *Information Processing in Medical Imaging: 26th International Conference, IPMI 2019* (2019).

[21] Ian Covert, Scott Lundberg, and Su-In Lee. "Understanding Global Feature Contributions With Additive Importance Measures". In: *Advances in Neural Information Processing Systems* 32 (2020).

[22] Joe H Ward Jr. "Hierarchical grouping to optimize an objective function". In: *Journal of the American statistical association* 58.301 (1963), pp. 236-244.

[23] Mark D Lescroart, Dustin E Stansbury, and Jack L Gallant. "Fourier power, subjective distance, and object categories all provide plausible models of BOLD responses in scene-selective visual areas". In: *Frontiers in Computational Neuroscience* 9.November (2015), p. 135.

[24] O Esteban et al. *sMRIPrep: structural MRI PREProcessing workflows*. 2021.

[25] Oscar Esteban et al. "fMRIPrep: a robust preprocessing pipeline for functional MRI". In: *Nature methods* 16.1 (2019), pp. 111-116.

[26] Edmund T Rolls et al. "Automated anatomical labelling atlas 3". In: *Neuroimage* 206 (2020), p. 116189.

[27] Nilearn contributors. *nilearn*. DOI: https://doi.org/10.5281/zenodo.8397156. URL: https://github.com/nilearn/nilearn.

[28] Fabian Pedregosa et al. "Scikit-learn: Machine learning in Python". In: *the Journal of machine Learning research* 12 (2011), pp. 2825-2830.

[29] David S Knopman et al. "Alzheimer disease". In: *Nature reviews Disease primers* 7.1 (2021), p. 33.
[30] Gary W Van Hoesen, Josef Parvizi, and Ching-Chiang Chu. "Orbitofrontal cortex pathology in Alzheimer's disease". In: *Cerebral Cortex* 10.3 (2000), pp. 243-251.

[31] Laura W de Jong et al. "Strongly reduced volumes of putamen and thalamus in Alzheimer's disease: an MRI study". In: *Brain* 131.12 (2008), pp. 3277-3285.

[32] Hanneke Van Dijk et al. "The two decades brainclinics research archive for insights in neurophysiology (TDBRAIN) database". In: *Scientific data* 9.1 (2022), p. 333.

[33] Philipp Bomatter et al. *Machine learning of brain-specific biomarkers from eeg. bioRxiv.* 2023.

[34] Wiremu Hohaia et al. "Occipital alpha-band brain waves when the eyes are closed are shaped by ongoing visual processes". In: *Scientific reports* 12.1 (2022), p. 1194.

[35] Nicolai Meinshausen, Lukas Meier, and Peter Bühlmann. "P-values for high-dimensional regression". In: *Journal of the American Statistical Association* 104.488 (2009), pp. 1671-1681.

[36] Vladimir Vovk and Ruodu Wang. "E-values: Calibration, combination and applications". In: *The Annals of Statistics* 49.3 (2021), pp. 1736-1754.