Measuring Variable Importance in Heterogeneous Treatment Effects with Confidence
Authors: Joseph Paillard
Angel David REYERO LOBO
Vitaliy Kolodyazhniy
Bertrand Thirion
Denis-Alexander Engemann
Venue: ICML 2025 poster (2025)
Keywords: Causal Learning
Interpretable ML
Primary area: general_machine_learning->causality
TL;DR: We show how to measure the importance of variables driving heterogeneity in treatment effects.

Abstract: Causal machine learning (ML) promises to provide powerful tools for estimating individual treatment effects. 
While causal methods have placed some emphasis on heterogeneity in treatment response, it is of paramount importance to clarify the nature of this heterogeneity, by highlighting which variables drive it.
We propose PermuCATE, an algorithm based on the Conditional Permutation Importance (CPI) method, for statistically rigorous global variable importance assessment in the estimation of the Conditional Average Treatment Effect (CATE).
Theoretical analysis of the finite sample regime and empirical studies show that PermuCATE has lower variance than the Leave-One-Covariate-Out (LOCO) method and provides a reliable measure of variable importance.
This property increases statistical power, which is crucial for causal inference applications with finite sample sizes.
We empirically demonstrate the benefits of PermuCATE in simulated and real datasets, including complex settings with high-dimensional, correlated variables.

Measuring Variable Importance in Heterogeneous Treatment Effects with Confidence
Joseph Paillard 12
Angel Reyero Lobo
Vitaliy Kolodyazhniy'
Bertrand Thirion Denís A. Engemann
Abstract
and development (Frank & Hargreaves, 2003; Hartl et al., 2021). Machine Learning (ML) promises improve the prediction of biomedical and clinical outcomes from high- dimensional, high-throughput data in various biomedical applications due to its ability to learn complex functions and representations from hetcrogencous and correlated data, such as those encountered in e.g. imaging, genomics or transcriptomics data (Le Goallec et al., 2022; Parrish et al., 2024; DeGroat al., 2023). Regardless, their notorious lack of interpretability hinders the applicability of ML models for generating scientific and clinical insights, motivating efforts for new methods for estimating variable importance within the field interpretable ML (Molnar, 2020; Biecek & Burzykowski, 2021). Most studies so far have focused classical MI models like random forests (Strobl et al., 2008), which well-suited for tabular data with few vari- ables (Grinsztajn et al., 2022). A more recent line of work, exposed to complex life-science data, has therefore em- phasised the importance of high-capacity MI. algorithms and proposed methods combining high-capacity ML algo- rithms with inference algorithms backed by statistical guar- antees (Mi et al., 2021; Williamson al., 2023; Lei et al., 2018; Chamma al., 2024).
Causal machine learning holds promise for esti- mating individual treatment effects from complex data. For successful real-world applications of machine learning methods, it is of paramount im- portance to obtain reliable insights into which variables drive heterogeneity in the response to treatment. We propose PermuCATE, an algo- rithm based on the Conditional Permutation Im- portance (CPI) method, for statistically rigorous global variable importance assessment in the es- timation of the Conditional Average Treatment Effect (CATE). Theoretical analysis of the finite sample regime and empirical studies show that PermuCATE has lower variance than the Leave- Tie-Covariate--O0 Out (LOCO) reference method and provides a reliable measure of variable im- portance. This property increases statistical power, which is crucial for causal inference in the limited-data regime common to biomedical applications. We empirically demonstrate the ben- efits of Pe PErTUAATE in simulated and real-world health datasets, including settings with up to hun- dreds of correlated variables.
Despite significant efforts and advances methods for in- terpretability, the utility of ML models remains limited as biomedical data are often of interventional nature. There- fore, such datasets are typically analyzed using biostatistical methods appropriate for causal inference and accepted by regulatory agencies. Over the past decade, findings from the causal inference and potential outcomes literature (Imbens & Rubin, 2015) have been adapted to propose the estimand framework (Akacha et al., 2017; Ratitch et al., 2020), flexible statistical-modellingg framework useful for causal analysis of clinical data. provides formal concepts and draws a conceptual distinction between the quantity of in- terest and the estimation method. This enables a formal language that can bridge the gap between classical biostatis- ties and ML, hence, pave the way for causal ML in biomed- ical applications (Feuerriegel et al., 2024; Sanchez et al., 2022; Doutreligne & Varoquaux, 2025). When studying heterogeneity treatment effects using ML, explainability with statistical guarantees crucial for uncovering useful predictive biomarkers ad-hoc methodology may increase
1. Introduction
The growing volume of high-dimensional biomedi- cal data holds promise to facilitate the discovery of novel biomarkers that can predict individual treatment success or risk of adverse events (Pang et al., 2023; Wu et al., 2020; Kong et al., 2022). Transforming complex biological sig- nals into actionable variables, therefore, plays a crucial role for increasing the efficiency of biomedical research
Roche Pharma Research & Early Development, Roche Innova- tion Center Basel, F. Hoffmann-La Roche Ltd, Basel, Switzer- land; UUniversité Paris-Saclay. Inria, CEA, Palaiseau, France: IIstiiut de Mathematiques de Toulouse, UMR5219 Univer- sité de Toulouse, France. Correspondence to: Joseph Paillard <joseph, paillard@roche.coo
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
knockoffs, both approaches measure the influence of each variable on the output of (possibly complex) model by plexity of multimodal sigh-dimensional biological data and support the research questions that arise from interventional studies with estimands such as average treatment effects (ATE) and conditional average treatment effects (CATE). In the causal setting, it remains unclear whether variable im- portance can be measured with intact statistical guarantees.
the false positives rate as well as the false negative rate.
Taken together, in the context of present and future biomed- ical applications, ML methods must address both the com-
estimating the total Sobol index. This well-known quantity can be interpreted as the loss increase when a variable is removed from the input of a predictive model. It has also been interpreted as a generalized form of ANOVA (Sobol, 2001; Williamson et al., 2021).
Prior research also studied the problem of variable impor- tance for CATE estimation in the case of causal forest esti- mators (Bénard & Josse, 2023) or Bayesian Additive Trees (BART) (Chipman et al., 2010). For biomarker detection in the CATE context, Sechidis et al. have investigated knockoff procedures with LASSO models and causal forests. How- ever, these approaches limit the choice of prediction models. Based on the behavior of random forests and knockoff's in standard prediction problems (Chamma et al., 2023), we would expect that both approaches encounter limits when it comes to nonlinear effects in high-dimensional data with high correlations. Another line of work focused on inter- pretability using attribution methods such as Shapley values or integrated gradients (Crabbé et al., 2022). While such methods are often good at ranking important variables, they do not provide type-1 error control. This is particularly true for Shapley values, a measure of importance based on an additivity axiom, meaning that estimated importance tends to be smeared over correlated variables. This property nec- essarily leads to false positives in the presence of correlation. Hines et al. (2022) have developed a LOCO procedure in concert with the potential-outeomes feasible risk, hence, ex- tending a statistically grounded variable-importancee method to the CATE context.
Contributions. In this work we study algorithms for esti- mating the importance of variables with statistical guaran- tees in the context of interventional data and CATE.
We propose a new variable importance inference pro- cedure for heterogeneous treatment effects termed PermuCATE, which generalizes Conditional Permu- tation Importance (CPI) (Chamma et al., 2023). We put it in perspective with the established leave-one- covariate-out (LOCO) method (Williamson al., 2023; Lei et al., 2018).
We develop theoretical analysis of the finite sample variance for both PermuCATE and LoCO, focusing on the impact of estimation error in finite samples.
Through experiments on simulated and real-world data, using various causal ML models, we show how Pe PermuATE achieves lower variance and higher sta- tistical power due the influence of finite sample estimation error, an effect that is more pronounced in small sample sizes or when using complex deep learning CATE estimators.
As pointed for standard prediction problems (Chamma et al., 2023; 2024), in contrast with permutation methods, LOCO is a refitting procedure requiring a full model fit per variable of interest. In the CATE context, the required amount of computation is therefore expected to be higher, depending the meta-learner and the choice of feasible risk. It is thus a question of interest if CPI can be further developed to support variable importance in CATE estimation. Finally, a neglected aspect of previous studies of LOCO and CP I is the impact of estimation error that stems from the refitting and the covariate estimation, respectively.
All proofs and additional experiments are given in appendix.
2. Related Work and Setting
2.1. Related Work
We first review prior art on variable importance in stan- dard prediction problems and, subsequently, their applica- tions to heterogeneous treatment effects. Previous works have analysed (Verdinelli & Wasserman, 2023) and bench- marked (Chamma et al., 2023) a wealth of procedures and statistical methods for variable importance estima- tion (Strobl et al., 2008; Nguyen et al., 2022; Watson & Wright, 2021; Covert et al., 2020; Candes et al., 2018; Louppe et al., 2013). Out of those approaches for stan- dard prediction problems, LOCO (Williamson et al., 2023) and CPI (Chamma al., 2023) satisfied the requirements of being model-agnostic (therefore compatible with high- capacity prediction models), offering statistical guarantees for controlling type-1 error, having a realistic computation cost and being interpretable (Chamma et al., 2023). In addi- tion to enabling variable selection, similar to methods like
The following sections present formal treatment of the related work and prepare the ground for our contributions.
2.2. Problem setting
Capital letters denote random variables, and lowercase let- ters denote their realizations. We consider the potential- outcomes framework with a set of observations Z; (X,, Y(A))) Here, X; E IR denotes the set of covari- ates for subject i, A E {0, 1} indicates a binary treatment assignment, and Y((A)) represents the observed outcome.
2

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
We focus on the estimation of the Conditional Average Treatment Effect r(x) = (1) Y captures the expected difference between the outcome under treatment (A, = 1) and control (A, = 0) for a subject with covariates X = x. Unlike classical MI. problems, this func- tion depends on potential outcomes, which are unobserved data since each subject is assigned to either the treatment or control group. This leads to fundamental challenges in the estimation procedure of this function as well as in scoring the candidate estimate since the ground truth is not readily available. For background on CATE estimation and feasi- ble risks, we recommend the work by Kennedy (2023) and by Doutreligne& & Varoquaux (2025), respectively.
propensity score. Alternatively, a recent review of causal risks (Doutreligne & Varoquaux, 2025) for model selection in causal inference suggested using the R-risk given by
zr., which
((Y - (0mmm(()))) - (A - x (X))
(3)
where m(X) m(X) (X) + x()) (X). Both risks are consistent with the oracle PEHE (proof in Ap- pendix A.2) when the ground-truth nuisance functions are used. Additionally, for variable importance, they empiri- cally produce the same results (Figure S1). We used the PO-risk for the rest of this work as it is easily comparable to the oracle PEHE, has been used in previous publication and is easy to compute.
Leave-One-Cove Oww-nnonrriite-Oni- Out (LOCO) for CATE: Unlike most efforts centered on the estimation of the CATE, in this work, we focus on measuring the importance of the differ- ent covariates X = 1, d for CATE prediction. In their work on variable importance for general ML problems, (Williamson et al., 2023) analysed in detail the Leave-One- Covariate-Our Out (LOCO) approach for estimating the impor- tance of a variable (or group of variables) as the decrease in performance between an estimator fitted on the full set of covariates X and the subset of covariates excluding the jth, that we denote -J This approach was then adapted to the causal framework (Hines et al., 2022), providing a variable importance measure (VIM) given a risk R that measures the predictiveness of a CATE estimator,
Conditional Permutation Importance (CPI): Another line of work on conditional VIM has introduced the Con- ditional Permutation Importance (CPI) measure (Chamma et al., 2023). Instead of re-fitting the estimator on a subset of covariates, this method involves sampling the given subset from the conditional distribution and evaluating the impact on the loss. This provides a second estimator of variable importance, denoted Verr, described in the case of CATE in equation 4. Similar to LOCO, this approach provides confidence intervals, is model-agnostic, and enjoys type-I error control. However, estimating the conditional density as required by CPI is often less error-prone than estimating sub-model in LOCO since in many practical situations, the relationships among the covariates are simpler than the mechanisms linking the covariates to the outcome as dis- cussed in Remark 3.3. Finally, both approaches present similarities, among which, converging asymptotically to the total Sobol index (Sobol, 2001).
V LLoco = R 2
- R (T,X,A, (1)
where 7ー being the CATE estimator fitted on the subset of covariates. In practice, CATE estimators are fitted on a training set, and ( Loco is evaluated on a disjoint test set.
The initial formula- tion of the variable importance measure by Williamson et al. 2023 allows for the choice of different metrics to evaluate the estimator. However, the extension of this framework encounters the fundamental problem in causal inference: the oracle CATE function is not accessible in practice (Hol- land, 1986). Thus, applying ML metrics commonly used for empirical risk minimization in prediction tasks, such as the precision of estimating heterogeneous effects (PEHE), which is the Mean-Squared Error (MSE) in estimating the CATE, cannot be computed. To overcome this hurdle, feasible risk that can be computed from real data should be used to provide practical methods. Hines et al. 2022 for instance used the PO-risk: the mean squared error between the CATE prediction and the pscudo-outcome given by
Feasible Risks for Causal Inference:
Unlike LOCO, CPI reuses the same estimator when predicting from the perturbed set of covariates. This reduces both the computational cost and the error that comes from stochastic models optimization. It is worth noting that in the work from Hines et al. 2022 adapting LOCO for CATE estimation, only the second stage of the De-learner is re-fitted on the subset of covariates: T-j = E[p(2) = r] By contrast, the pseudo- outcomes pox) (cf Equation 2). are obtained from the full set of covariates. Bypassing the intermediate estimation of the nuisance functions, this trick mitigates the extra com- putational cost incurred by LOCO when applied to causal estimators as opposed to regular ML problems and avoids vi- olating the unconfoundedness assumption when estimating the CATE with missing covariates. However, in scenarios where this final estimation step requires complex models, such as super-learners (Laan et al., 2007) or deep neural networks (Curth & Schaar, 2021) the computational burden can become intractable.
Computational efficiency:
(0())= = (y - Pe (a -aaaaaa--a T(x)) [())) + un (x) - -o(()), (2)
T (1x π (z))
where Ho and μ1 denote the response functions: (a =
EY|X- = x, A = a) and (x)) = E|A- = 1|X = x] is the
3

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
3. PermuCATE algorithm
estimating U, with consistent estimator 0, (see Proposition 3.3 from Reyero Lobo et al. (2025)). Then, the residual information rj X by can used to construct a perturbed version of the jdh covariate: 2 + +r where r, shuffle(ry) is permutation of the residual.
The method presented in this section aims at identifying important variable for CATE prediction. While similarly to LOCO, it is agnostic the meta-learning approach used to estimated the CATE, 2023). This choice is motivated by its favorable conver- gence properties and the fact that the scoring step requires deriving pseudo-outcomes explicitly, which limits the com- putational benefits of simpler and more direct meta-learning approaches such as the S- or T-learner.
used the DR-learner (Kennedy,
Proposition 3.2. Assume the consistency of the estimators i and Assumption 3.1. PermuCATE consistently estimates E V [xX) the quantity known as the total Sobol index.
The proof is provided in subsection A.3.
Remark 3.3. Assumption 3.1 amounts to assuming that the relationships between the covariates within X are simpler than the link between the CATE, 7 and the input X. Conse- quently, estimating the conditional density is simpler than estimating the outcome, akin to the 'model-X' knockoff framework presented in Candes et al. 2018. Numerous ex- amples illustrate the utility of this setting, such as in the study of biologically complex diseases using multiple vari- ables like single nucleotide polymorphisms (SNPs) where the relationship between covariates is much simpler than their association with clinical outcomes. Similarly, in clini- cal trials, the covariates are collected concurrently (e.g. at haseline) but the outcome of interest is measured at a later time. Furthermore, the amount of unlabelled data often ex- ceeds the number of available labels, especially in clinical settings where labelling can be costly.
Algorithm 1 Connditional Permutation Importance for CATE
Input: Dirains Preet two independent sets of observations Z, = (X,, A, X, E Rd
Parameter: Response functions: Ao; ii; Propensity: f;;
CATE: f; Covariate predictor: D(-), feasible risk: R(-),
Number of permutations:
Output: Importance estimate V E I&d
1: Using Dtrain
2: Fit Po. fr,, n. F // CV can be used for DR-learner
3: for j == 1, 1.... do
4: Fit (.)) + [X
5: end for
6: Using Dreat
7: for -,,... , d do
8:
ij + (3)
To check if Assumption 3.1 is realistic, we conducted com- putational experiments via simulations and analysis of real- world medical data. We focused on the IHDP benchmark which uses real-world data on
9:
rj + X
g+x' - f1j
10:
for k = 1..... do
11:
F, x shuffle(r,)
development. The dataset has Gon-Gausiian continuous features and imbal- anced categorical (up to four categories) features. Because the ground truth for the conditional distribution is not avail- able, the validity of Assumption 3.1 cannot be assessed. In addition to this benchmark, we propose additional experi- ment varying the complexity of the conditional distributions. For this purpose, we use simulations inspired by Hollmann et al. 2025 and sample the covariates using a latent variable model. We first sampled latent variables from mixtures of Gaussians and then generated observed covariates through non-linear transformation with interaction terms of the la- tent variables. This resulted in non non-Gaussian covariates with complex conditional distributions as shown in Figure S5.
12:
Xr. + [X1. , -x-111 + Fy.... ,Xa]
13:
J + (R((F(Xp)). Y, A) - R(F((), A))/2
14: end for
15: TT < P Ekmn 11
16: end for
17:
return ( U
LOCO estimates the importance of a given variable by refit- ting model subset of covariates, whereas permutation- based approaches reuse the same model to predict from perturbed design matrix. In the traditional permutation im- portance approach (Breiman, 2001), the perturbed matrix differs from the original by a permutation of the studied covariate. That approach is naive in the sense that it does not account for correlation, leading uncontrolled type-1 errors; CP I, refined version mitigates this issue (Chamma et al., 2023). The key idea is to sample from the conditional distribution p (i)) which requires the following as- sumption.
Adapting this idea to the context of CATE estimation, given causal risk Rt the following measure of variable impor- tance can be obtained:
Vepr = (R (T, A,Y) - R (T,,
(4)
2
Assumption 3.1. There exists a function 1/ such that X by (x(())) + with E IL x(1) and B[E,] = O.
where XPj = [X [XI.-...XX-11.14 + Fy. f.. Xa]. PermuCATE (Algorithm 1) describes the procedure to ob- tain V'pr-- Previous works focusing on predictive models have shown how a fross-fitting approach allows to obtain a
Sampling from the conditional can thus be achieved by
4

---

Mcasuring Variable Importance in Heterogencous Treatment Effeett with Confidence
valid estimate of the desired quantity with a control of the type-I error (Williamson al., 2023; Verdinelli & Wasser- man, 2023; Chamma et al., 2023; Reyero Lobo et al., 2025).
Dirain and Σ is the covariance matrix of X- (proof in Appendix A.8). Consequently, Equation (6) demonstrates that the variance of PermuCATE is mainly influenced by the error in estimating the conditional distribution via i but not on the error in +. Concerning the LOCO approach, the iterative refitting of the model on the subsets of covariates X leads to fluctuations due to the estimation of 7 and introduces bias term associated with the CATE prediction.
Properties The following paragraphs compare the proper- ties of the LOCC and PermuCATE TE estimators in the context of finite training sets Dtrain. to highlight variance-inducing error terms from finite sample estimations. A key advantage of PermuCATE is that it allows the use of an expressive model for the main function estimation, such as deep neu- ral network, without incurring the high cost for refitting the full model, provided that more parsimonious model is used for covariate prediction. However, in finite sample regimes, this estimation of ūg on a fixed train set Dtrain introduces variance. We denote || 1|a is the 12-norm given the fixed training set, Dirain-
Proposition 3.5. For a consistent CATE estimator 9, the finite sample estimation error of the importance for the jth covariate is
E T WioooolPene Perain
V7 Op (I[--) - f-1112 -1| -1I- =]1 (7)
Proof in Appendix A.5.
Proposition 3.4. Under Assumption 1, for a consistent CATE estimator i, the finite sample bias the estimation of the importance for the jth covariate is given by.
Denoting 28 and AB_; the difference between the esti- mated and true coefficients for respectively =(X) = 3X
and the sub-model T-j. In the linear scenario as presented above the finite sample variance of the LOCO estimator is then
E. Crri (CCriPeeraii Derain] Vj-
B[(+(Xp)) - {(xp)))² Dirain] + OP(TT| - f F(Derai
var E
V Loco (usiiu)
= var (114
where ; E [V [x(X)xx)-) is the targeted total Sobol index.
+ var (114811))
(8)
Remark: Although the first term involves the CATE esti- mate 4, it does not include the estimation error, thus mak- ing more robust to model misspecificatiios. The estima- tion error appears only in the second, linear term, instead of the mean- seaan- error (MSE) present for LOCO, The proof for that relies the fact that MSE terms cancel out because [((((Xpp)) - 7( p.j))² (Drriin]
In the above equations, which present the different bias and variances in linear setting assuming independent train- ing samples, a key distinction emerges between the two approaches. The variance of LOCO depends on the sum of the error in estimating the CATE f and its sub-model fy. In contrast, for PermuCATE it depends on the error in esti- mating the conditional distribution (through the error in it)). As discussed in Remark 3.3, the latter prediction problem is arguably easier to solve. Supplementary Figure $2 provides an illustration of this comparison.
[X[(T(X)-
=(X)))(Dtrain) since by construction Xpj iid. X. Corollary: Under the additional assumption that f is Lips- chitz continuous and è is consistent, the bias term becomes
Op(llo3 - -351 + Op - +))
(5)
Therefore, we investigated the variance behavior of both approaches using mamte- simulations and datasets. The results in the following section will show empirically that the error term in LOCO can lead to increased variance and reduced statistical power.
Proof in Appendix A.4.
In line with Proposition 3.2, it can be observed that, given a consistent estimator of vj. Cri is a consistent estimator of the total Sobol index. It also shows that the error comes from estimating 8, To better understand the impact of estimation error on variance, consider a linear setting with Gaussian covariates where the CATE is defined as r(X)= 3X,, with a finite training sample and an infinitely large testing sample, the variance of the PermuCATE estimator is given by
4. Experiments
Datasets: This work studies three different simulation sce- narios and real-world dataset. In the paper, we refer to high-dimensional data as it commonly understood in the field of medical research, assuming that prior dimensionality reduction has been performed, thereby limiting the number of variables to at most a few hundred. Scenarios involving much larger dimensions are not addressed here due to the high correlation between covariates, which leads to vanish- ing importance and a consequent lack of statistical power.
var B T CPPI Derain)) v (yll ,) (6)
where A = y - ? is the error when estimating the coeffi- cients of the function ()) = 7 -1 on the training set,
5

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
However, our proposed approach is compatible with group- ing strategies (Chamma et al., 2024), as illustrated in the Figure S3. This allows for inference the importance of groups of covariates rather than individual covariates, help- ing to mitigate the challenges posed by highly correlated variables and facilitating the handling of higher-dimensional data when they can be grouped. The following list presents an overview of the datasets, further details are provided in the Appendix A.10.
6
LOCO
4-
2
0
True positive
X1
6
CermuCATE
False negative
X2
4
X3
The Low Dimensional (LD) dataset is the data gen- erating process from the work of Hines et al. 2022. It is used as a benchmark to compare the proposed approach to their baseline method. It consists of six covariates, each pair ((X1,X2) (Xa, X4), (Xs, X6)) being correlated (p- = 0.5). Three of these covari- ates X1,X2 X3) are linked to the nuisance fune- tions and CATE through linear functions. The ana- lytical importance for these variables is respectively V 0.75, 2 3, 3 0.75.
3
2
0
10a
10a
sample size
Figure I. PermuCATE detects important variables with greater statistical power. Using the gatl-gennrating process proposed by (Hines et al., 2022), compared the estimates of variable im- portance using 1,000 (top subplot), the baseline proposed hy the authors, and PermuCATE. By computing p-values, observed different sample sizes whether each of the three important variables was correctly identified (true positive) or missed (false negative) significance level a = 0.05. For all three important variables, the minimum sample size (vertical dashed line) required to detect their importance was smaller with the PermuCATE approach. Both methods converged to the thcoretical importance score (horizontal dotted line). The shaded areas around the curves represent Monte Carlo estimates of the standard deviation over 100 random sam- ples.
The High Dimensional Linear (HL) dataset is an exten- sion of the work by Doutreligne& & Varoquaux 2025 to high dimensional settings. It consists of = 50, 100 covariates, 10 of which are linked the nuisance func- tions and CATE by a linear function. The loading of each important covariate is drawn from Rademacher distribution.
The High Dimensional Polynomial (HPP dataset is sim- ilar to HL but uses a linear combination of polynomial features of degree three, including interaction terms.
The Infant Health and Development Program (IHDP). The dataset consists 747 subjects with 25 real covari- ates, including 6 continuous and 19 binary variables, along with a simulated outcome that is both non-linear and noisy (Shalit et al., 2017).
To compute p-values over the different folds we used the Wadeau-Bengio corrected t-test (Nadeau & Bengio, 1999). True positives and false negative were defined with respect to a significance level of a = 0.05. Both methods provided a false positive rate below 0.05. As presented in Figure 1, both methods asymptotically converged to the analytical im- portance score represented by dotted lines. However, it can be seen that for FermuCATE, the standard deviation was smaller (represented by shaded areas). This larger variance of the LOCO estimator translated into decreased statistical power, as revealed by the minimum sample size needed to identify each of the truly important variables represented by vertical dashed lines. The following experiment aimed at studying its impact in more complex and high-dimensional settings.
Variable importance for low dimensional datasets We first compared both approaches using the previously pro- posed L.D benchmark (Hines al., 2022). Out of six covari- ates, the first three were important for the CATE prediction. To estimate the CATE, we used a DR- leanner (Kennedy, 2023) with regularized linear models for nuisances functions and the final regression step. For PermUCATE, we used the same regularized linear model for covariate prediction and used 50 permutations. The importance of variables was esti- mated using a nested cross-fitting scheme. In each split split, 20% of the data was left out for the importance estimation. The remaining 80% was used to fit the DR-learner using the cross-validation scheme presented in the work from Kennedy (2023). We used a five-fold cross-fitting strategy for both loops. Additionally, for each sample size, we re- peated these computations on 100 different random samples.
Variable importance in high dimensional settings To study the scaling behavior of both methods, we evaluated their capacity to identify important variables in datasets with higher dimension. We used two simulation scenarios, HL and HP, which respectively used linear and polynomial link functions between the covariates and the different nui-
6

---

Mcasuring Variable Importance in Heterogencous Treatment Effeett with Confidence
linear
non-linear
1.0
power 0.05)
10-1
Ridge
10
variance
PermuCATE LOCO
10-4
10a
0.5
C
d
0.0
CateNet
variance
10a
10-3
10a
power (a = 0.05)
1.0
b
109
PermuCATE
0.5
LOCO
variance
TabPFN 10a
d= 50
d 100
10-4
0.0
10
400
800
1600
3200
6400
12800
103
104
10a
10
sample size
sample size
sample size
Figure 3. The PermuCATE method identified more important variables in high-dimensional, linear, and complex scenarios. (a) Statistical power for detecting important variables as a func- tion of sample size the HP dataset (non-linear scenario) with Pe PermuCATE and LOCO methods. The CATE was estimated with DR-learner using and regularized linear models to estimate the nuisance functions. Results are shown for two scenarios with = and d = 100 covariates, respectively. In both cases, only 10 covariates were important. Error bars represent Monte Carlo estimates of the stan- dard deviation over 10 simulations. (b) Same (a) but using the HL Gata-generating scenario and only regularized linear models to estimate the nuisance functions.
Figure PermuCATE demonstrates lower variance mis- specified models and small sample sizes. The figure compares the variance of the CATE estimations when true CATE is lincar function of the covariates (left column) and polynomial function of degree three with interaction terms (right column). used: (a, DR-leamerr with Ridge to estimate the nuisances and pseudo-outcomes: (c, d) CateNet from Curth & Schaar 2021; and (e, f) DR-learner using TabPFN, Both simulation scenarios used (d=50) correlated Gaus- sian covariates with support size of 10. Variance was assessed for one predictive variable across 100 random simulations. Error bars represent bootstrap estimates of the variance.
stacking gradient boosting trees
Three different CATE estimators
higher robustness to model specification and finite sam- ple estimation error for PermuCATE are consistent with Equation (5) and Equation (7), which show that LOCO bias depends on the error in estimating the CATE, whereas Pe PermuCATE's does not.
sance functions and CATE. Both scenarios used d Gaussian covariates with correlation of 0.5. The support size was sparse using only 10 important variables. For this more complex setting, CATE was estimated using more pressive models: CateNet, a DL architecture proposed by Curth & Schaar 2021, TabPFN, a DR-learner using the pre- trained transformer architecture from Hollmann et al. 2025 to estimate the nuisances and pseudo-outcomes and finally baseline, Ridge as in the previous experiment. Figure 2 presents the variance when estimating the importance of predictive covariate using PermuCATE and LOCO for dif- ferent CATE estimators in linear and non non-linear scenarios. In the linear setting, where the eтror in estimating the CATE was similar to the error in estimating the conditional distri- bution, the variances were comparable, with smaller values in small sample sizes for Pe reruCATE and faster asymp- totic decrease for LOCO. In addition, PermuCATE shows systematic benefit for CateNer (panel e) where a stochastic optimization was used. However, in the non-linear sce- nario, PermuCATE demonstrated significantly lower vari- ance, sometimes by orders of magnitude. This effect was particularly pronounced in small sample sizes and tends to diminish asymptotically. These findings, which suggest
We then compared their statistical power at the significance level a = 0.05, defined by tp/(tp - fn), were tp is the number of true positive (variables that belong to the true support and have p-value S 0.05), and fn the number of false negatives. Figure presents the statistical power the HP dataset for different sample sizes. The CATE was computed using DR-learner for which the different nuisances and presudo-outcomes are estimated using stack- ing of gradient boosted trees and regularized linear model drawing inspiration from the popular Super Learner method (Laan et al., 2007) (further details in subsection A.11). The statistical power of PermuCATE was greater, especially in the non-linear setting (Figure 3, a). Besides this advantage, it is remarkable that PermuCATE suffered less from the increase of dimension as opposed to LOCO which presents a greater drop in power in the case where d = 100. This observation likely stems from the fact that the error dif- ference presented in Equation (7) fluctuates more as the

---

Mcasuring Variable Importance in Heterogencous Treatment Effects with Confidence
dimension increases. Although not represented on that plot, both methods demonstrated the desired type-1 error control. Figure 3b shows the same as (a) but in the linear scenario (HL). Here again, PermuCATE appeared to benefit from increased statistical power, leading to a better identification of important variables in small sample sizes.
a cr
CATENet
HH
TabPFN
HHH
-0.6
-0.4
0.2
0.4
0.6
negative PEHE
Power = 0.05)
b CF
d
PermuCATE
LOCO
Variable importance one the IHDP benchmark PermuCATE and LOCO were then compared on the estab- lished semi-symthetic IHDP benchmark (see datasets). For binary variables, sampling from the conditional distribution was achieved using the predicted probabilities of a logistic regression model. Three different learners were used to esti- mate the CATE: Causal Forest (CF) with 100 trees (Athey & Wager, 2019), deep neural networks (CATENett) from (Curth & Schaar, 2021) using the best model parameters reported in the publication, and a DR-learner incorporat- ing tabular foundation models (Hollmann et al.. 2025) to estimate nuisance functions and predict pseudo-outcomer Consistent with Curth & Schaar 2021, results were com- puted over 100 simulations using a single train-test split. The semi-synthetic nature of the dataset allowed measuring the Precision in Estimating Causal Effects (PEHE), which presented Figure 4 a. It further allowed quantifying the ability of each variable importance method to recover the true support, for each learner, as shown in Figure 4 b, d. p-values were computed using variable importance estimates derived from the 100 simulations. The type-1 er- ror corresponds to fp/(fp) + fn), where fp stands for false positive and tn for true negative. The statistical power is computed by tp/(tp + fn) where tp stands for true positive and tn for true negative.
CATENet-
HH
TabPFN
0.6
0.8 AUC
1.0
0,00
0.05
0.10
type-l error (a= 0.05)
Figure 4. Comparison of variable importance methods with three different learners on the IHDP benchmark. The CATE was estimated using a Causal Forest (CF, Athey & Wager 2019), deep neural network (CATENet, Curth & Schaar 2021), and pre-trained tabular foundation model (TabPEN, Hollmann et al. 2025), For cach learner, variable importance LOCO and PermuCATE. (a) Displays the negative Precision in Estimating Heterogeneous Effects (PEHE) across the models. (b) presents the AUC for estimating variable importance, which can be framed as a classification task. The predicted importance is rep- resented by -p-value, which is subsequently utilized to compute statistical power (c) and the type-1 error (d) at the significance level alpha=0.05.
estimated with
sample estimation error but also with optimization noise which is inherent to this learning paradigm. Finally, sub- panel d shows that both methods effectively controlled the type-1 error, except for LOCO when CATENet is used, likely due the aforementioned optimization noise. The experiment suggests that on real-world covariates, PermuCATE can achieve greater statistical power and better importance rank- ing than its LOCO counterpart, while allowing the use of ofate-of-the-art deep learning and foundation models.
These results indicate that Pe PErTUCATE, provided better es- timates of variable importance under model misspecification conditions, as evidenced by its performance on the IHDP benchmark. As shown in Figure 4 a, the learner based on TabFPN outperformed the popular CATENet, itself better than the Causal Forest in terms of PEHE for estimating the CATE. For all three models, FermuCATE achieved sig- nificantly better AUC and statistical power for identifying important variables. This suggests that PermuCATE per- forms faovrably when applied to complex models to ap- proximate an unknown generative mechanism. One would typically select the best-performing model to measure vari- able importance, as it is more likely to capture the data- generative model. However, subpanel revealed that the AUC when using LOCO combined with a CATENet was lower than when using the CF despite the former achieving a better PEHE. This effect is due to the highly stochastic nature of the optimization process for CATENets, as it re- lied on the Adam algorithm, which introduces additional variance. When using deep learning models, the successive per-covariate refitting operations not only come with a finite
Complexity of the conditional distribution estimation The experiment using the IHDP benchmark involves real-world covariates, for which the ground truth of the conditional distribution is unknown. We found that PermuCATE outperformed LoCO (Figure 4) in terms of statistical power and AUC for detecting important features while controlling type-1 eгтor on this dataset. These re- sults suggest that Assumption 3.1 is realistic. Moreover, in Appendix A.13 we systematically varied the complex- ity of conditional distributions by generating non- -au-Sauss,nn, multimodal covariates from latent-variable model (Fig- ure S5). Our analysis revealed that PermuCATE again exhibited lower variance compared to LOCO (Figure S6), Differences in variance trends for both methods resembled results from the Gaussian scenario shown in Figure 2. A possible explanation may be related to the fact that increased complexity of covariate distributions not only rendered con-
8

---

Mcasuring Variable Importance in Heterogeneous Treatment Effeets with Confidence
ditional estimation more challenging (preferentially impact- ing PermuCATE) but also CATE estimation (preferentially impacting LOCO). Both methods may thus have suffered in comparable way but for different reasons, leading to stable ranking favoring Pe PermuCATE for smaller sample sizes.
The present work responded to two current trends in research & development for life science and healthcare: (1) the unmet need for interpretable ML methods for multimodal corre- lated biomedical data (beyond tabular data), and (2) CATE estimation with MI. in interventional clinical studies, such as the discovery of predictive biomarkers. While this has shaped the perspective and scope of the present work, important to emphasize that the theoretical foundation of our study makes our results applicable beyond life sciences and medicine, including policy-impact, educational, and economic applications.
5. Discussion
In this work, studied two estimators of variable importance in the context of CATE estimation: the LOCO approach (Hines et al., 2022) and the proposed PermuCATE algorithm. Our study focused on the pre- asymptotic amite-sample regime, aiming to provide insights on how to select variable importance methods for real-world use cases. Through theoretical analysis and simulation studies covering a wide range of scenarios, we unveiled how noise originating from the finite sample model fitting impacts both algorithms in specific ways through their re- spective variance terms. We wish to emphasize that this consequential aspect was neglected by most previous work analysing variable importance estimation, which focused more on type-1 error whereas our findings concern type-2 error. Our results also reveals how the finite-sample vari- ance critically impacts real-world applications with small- to medium-sized datasets, characteristic of the biomedical field. For practical applications, we would, thus, recom- mend Pe PermuÇATE TE over LCCO as both methods present the same statistical guarantees for type-1 error control (as estab- lished for standard prediction problems by Chamma et al. 2023), while Pe PermuCATE combines shorter computation times with higher power for detecting important variables for predicting treatment success. In clinical trials, this could determine whether informative signals are detected.
Limitations and future work Validating variable impor- tance methods on real-world data is challenging due to the lack of known ground truth. We addressed this by us- ing semi-synthetic data that combines real-world covariates and simulated outcomes. This leaves room for further ex- ploratory analysis and validation in applications that could benefit from modern ML techniques to extract new insights from rich biomedical datasets. For example, it can be ap- plied to prevention by examining how lifestyle and habits influence genetic factors (considered as treatments) on var- ious outcomes like quality of life and survival. clinical trial data, this method can support the identification of safety biomarkers to reduce the risk of adverse events and biomark- ers to predict treatment response.
Impact Statement
This paper introduces a model inspection method aimed at enhancing our understanding of heterogeneity in treatment responses in interventional studies including randomized-controlled trials but also in observational contexts. By providing robust explanations for potentially complex models, this approach can facilitate the adop- tion of cutting-edge ML advancements in fields where interpretability is crucial. Notably, this includes medical applications where the increasing predictive power of deep learning models, combined with the growing availability of data, holds significant promise for accelerating treatment development and supporting physician decision-making
The proposed approach shares similarities with 'model-X" knockoffs (KO) (Candes et al., 2018) as both involve mod- elling the conditional distribution of covariates. Both ap- proaches are motivated by the assumption that estimating covariates is easier than estimating the quantity of interest (e.g., CATE) as discussed in Remark 3.3. A key differ- ence is that 'model-X" KO is designed for variable selection, whereas PermuCATE and LOCO focus on variable impor- tance estimation. The latter provides richer information by measuring importance (total Sobol index) rather than just binary selection. Additionally, KO handles multiple testing, while our method ensures type-I error control. Furthermore, constructing KO variables requires complete exchangeabil- ity, a stricter condition than the conditional independence needed for PermuCATE. While not estimating variable im- portance, the conditional randomization test (CRT) from Candes et al. 2018 is more comparable in terms of its statis- tical properties to our approaches for individual conditional independence testing. However, CRT is much more compu- tationally expensive.
Acknowledgments: J.P., V.K, and D.E. are full-time em- ployees of F. Hoffmann-La Roche Ltd. The employer of these authors had no influence on the design of this study nor on the interpretation of its findings. This research has received funding from the H2020 Research Infrastructures Grant EBRAINN-Halte 101058516 and the VITE ANR- 23-CE23-0016 and PEPR Santé numérique, Brain health Trajectories ANR-22- 2NR-22-PESN-AA projects
9

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
References
Doutreligne, M. and Varoquaux, G. How to select predictive models for decision-making causal inference. Giga- Science, 14, 2025.
Akacha, M., Bretz, F., and Ruberg, S. Estimands in clinical trials broadennog the perspective. Statistics in medicine,
36(1):5-19, 2017.
Feuerriegel. S., Frauen, D., Melnychuk, V., Schweisthal, J., Hess, K., Curth, A., Bauer, S., Kilbertus, N., Kohane, I. S., and van der Schaar, M. Causal machine learning for predicting treatment outcomes. Nature Medicine, 30: 958-968, 2024.
Athey, S. and Wager, S. Estimating treatment effects with causal forests: An application. Observational studies, 5 (2):37-51, 2019.
Biecek, P. and Burzykowski, T. Explanatory model anal- ysis: explore, explain, and examine predictive models. Chapman and Hall/CRC, 2021.
Frank, R. and Hargreaves, Clinical biomarkers in drug discovery and development, Nature reviews Drug discov- ery, 2:566-580, 2003.
Breiman, L.. Random forests. Machine learning, 45:5-32, 2001.
Grinsztajn, L., Oyallon, E., and Varoquaux, G. Why do tree- based models still outperform deep learning on typical tabular data? Advances Neural Information Processing Systems, 35, 2022.
Bénard, C. and Josse, Variable importance for causal forests: breaking down the heterogeneity of treatment effects. arXiv, 2023. doi:
Hartl, D., de Luca, V., Kostikova, A., Laramie, J., Kennedy, S., Ferrero, E., Siegel, R., Fink, M., Ahmed, S., Mill- holland, J., et al. Translational precision medicine: an industry perspective. Journal of translational medicine, 19:245, 2021.
Candes, E., Fan, Y., Janson, L.., and Lv, J. Panning for gold:'model-x' knoooo- for high dimensional controlled variable selection. Journal of the Royal Stafistical Society Series B: Statistical Methodology, 80, 2018.
Chamma, A., Engemann, D. A., and Thirion, B. Statistically Valid Variable Importance Assessment through Condi- tional Permutations. Advances in Neural Information Pro- cexsing Systems, 2023. doi: 10885000rrxi.2200.077500
Hines, O., Diaz-Ordaz, K., and Vansteelandt, S. Variable importance measures for heterogeneous causal effects. arXiv, 2022. doi: (0.485550/axxxx22244.00
Holland, W. Statistics and Causal Inference. Journal the American Statistical Association, (396):945, 1986. doi: 10. 10.2307/2280664.
Chamma, A., Thirion, B., and Engemann, D. Variable im- portance in high- -igh--iimnsionall settings requires grouping. In Proceedings of the AAAI Conference Artificial In- telligence, volume 38, pp. 11 11195-112000 2024.
Hollmann, N., Müller, S., Purucker, L., Krishnakumar, A., Körfer, M., Hoo, S. B., Schirrmeister, R. T., and Hutter, F. Accurate predictions small data with tabular foundation model. Nature, 63778045)3319-326 2025.
Chipman, H. A., George, F. and McCulloch, R. E. Bart: Bayesian additive regression trees. The Annals of Applied Statistics, 4, 2010.
doi: 10.
Covert, I., Lundberg, S. M., and Lee, S.-I. Understand- ing global feature contributions with additive importance measures. Advances in Neural Information Processing Systems, 33:17212-17223, 2020.
Imbens, W. and Rubin, D. B. Causal inference in statis- tics, social, and biomedical sciences. Cambridge univer- sity press, 2015.
Crabbé, J., Curth, A., Bica, I., and van der Schaar, M. Bench- marking heterogeneous treatment effect models through the lens of interpretability. Advances Neural Informa- tion Processing Systems, 35, 2022.
Kennedy, E. H. Towards optimal doubly robust estimation heterogeneous causal effects. Electronic Journal of Statistics, 2023, doi: 10.1214/23-ej22157.
Kong, J., Ha, D., Lee, J., Kim, I., Park, M., Im, S.-.., Shin, K., and Kim, S. Network-based machine learning approach to predict immunotherapy response in cancer patients. Nature communications, 13:3703, 2022.
Curth, A. and Schaar, M. d. Nonparametric Estimation of Heterogeneous Treatment Effects: From Theory to Learning Algorithms. Advances in Neural Information Processing Systems, 2021.
Laan, M. J. d., Polley, C., and Hubbard, A. E. Su- per Learner. Statistical Applications in Genetics and Molecular Biology, 2007. doi: 10.2202/1544--1155.1300
DeGroat, W., Mendhe, D., Bhusari, A., Abdelhalim, H., Zeeshan, S., and Ahmed, Z. Intelligenes: novel machine learning pipeline for biomarker discovery and predictive analysis using melti-genomic profiles. Bioinformatics, 39, 2023.
Le Goallec, A., Diai, S., Collin, S., Prost, J.-B., Vincent, T., and Patel, C. J. Using deep learning to predict abdominal
10

---

Mcasuring Variable Importance in Heterogencous Treatment Efffetts with Confidence
age from liver and pancreas magnetic resonance images. Nature Communications, 13(1):1979, 2022.
Sechidis, K., Kormaksson, M., and Ohlssen, D. Using
knockoffs for controlled predictive biomarker identifica- tion. Statistics in Medicine, 40.
Lei, J., G'Sell, M., Rinaldo, A., Tibshirani, R. J., and Wasserman, L. Distribution-fieee predictive inference for regression. Journal of the American Statistical Associa- tion, 113(523):109-1111, 2018.
Shalit, U., Johansson, D., and Sontag, Estimating individual treatment effect: generalization bounds and algorithms. pp. 3076-3085, 2017.
Louppe, G., Wehenkel, L.., Sutera, A., and Geurts, P. Under- standing variable importances in forests of randomized trees. Advances in Neural Information Processing Syx-
Sobol, I. Global sensitivity indices for nonlinear mathemati- cal models and their Monte Carlo estimates. Mathematics and Computers in Simulation, 55:271-280, 2001.
tems, 26, 2013.
Strobl, C., Boulesteix, A.-L., Kneib, T., Augustin, T., and Zeileis, A. Conditional variable importance for random forests. BMC bioinformatics, 9:1-11, 2008
Mi, X., Zou, B., Zou, F., and Hu, Permutation- iden- tification of important biomarkers for complex diseases via machine learning models. Nature communications, 12(1):3008, 2021.
Verdinelli, I. and Wasserman, L. Feature Importance: A Closer Look at Shapley Values and LOCO. Statistical Science, 39, 2023. doi: 10. 20.1214/24-ST5337
Molnar, C. Interpretable machine learning. Lulu. com, 2020.
Watson, D. S. and Wright, M. N. Testing conditional in- dependence in supervised learning algorithms. Machine Learning, 110(8)2117-211)
Nadeau, C. and Bengio, Y. Inference for the generaliza- tion error. Advances in Neural Information Processing Systems, 12, 1999.
2021.
Williamson, B. D., Gilbert, B., Carone, M., and Simon, N. Nonparametric variable importance assessment using machine learning techniques. Biometrics, 77:9-22, 2021.
Nguyen, B. T., Thirion, B., and Arlot, S. A conditional randomization test for sparse logistic regression in high- dimension. Advances in Neural Information Processing Systems, 15:13691-, 2022.
Williamson, D., Gilhert, P. B., Simon, N. R., and Carone, M. General Framework for Inference Algorithm- Agnostic Variable Importance. Journal of the American Statistical Association, 2023.
Pang, M., Gabelle, A., Baha-Chaudhung P., Huijbers, W., Gafson, A., Matthews, P. M., Tian, L.,, Rubino, I., Hughes, R., Moor, C. d., Belachew, S., and Shen, C. Preci- sion medicine analysis of heterogeneity in individual- level treatment response amyloid beta removal in early Alzheimer's ' disease. Alzheimer's & Dementia, 2023. doi: 10002/lz..143.11
Wu, W., Zhang, Y., Jiang, J., Lucas, M. V., Fonzo, G. A., Rolle, C. E., Cooper, C., Chin-Fatt, C., Krepel, N., Cor- nelssen, C. A., et al. An electroencephalograpiie signa- ture predicts antidepressant response in major depression. Nature biotechnology, 38(4))439-447) 2020.
Parrish, R. L., Buchman, A. S., Tasaki, S., Wang. Y., Avey, D., Xu, J., De Jager, P. L., Bennett, D. A., Epstein, M. P., and Yang, J. Sr-twas: leveraging multiple reference panels to improve transcriptome-win wide association study power by ensemble machine learning. Nature Communi- cations, 15(1)66646, 2024.
Ratitch, B., Bell, J., Mallinckrodt, C., Bartlett, J. W., Goel, N., Molenberghs, G., O'Kelly, M., Singh, P., and Lip- kovich, I. Choosing estimands in clinical trials: putting the ich e9 (r1) into practice. Therapeutic innovation & regulatory science, 54:324 341, 2020.
Reyero Lobo, A., Neuvial, P., and Thirion, B. Sobol-cpi: doubly robust conditional permutation importance statis- tic. 2025. doi: 10448550/ar2i..2511.775
Sanchez, P., Voisey, J. P., Xia, T., Watson, H. I., O'Neil, A. Q., and Tsaftaris, S. A. Causal machine learning for healthcare and precision medicine. Royal Society Open Science, 9, 2022.

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
A. Supplementary Materials
A.I. Notations
X
Covariates. The matrix is of dimension x d with n the sample size and d the number of covariates
A
Treatment assignment
Y
Outcome
Ao
Control response
P1 Treatment response
Conditional Average Treatment Effect
Propensity score
Finite sample estimate of a function f
X- Set of covariates excluding covariate j
D3 E,, |X Finite sample estimate of covariate j from other
2 Xj Dy. Residual covariate j that cannot be predicted from the others
f shuffle(r,). Permuted version of the residual
AB_ d dimensional vector containing as entries, for (k) ر
= By 3e where are the coefficients of the model fitted on the full set of covariates and are the coefficients of the model fitted on the subset of covariates X-j
V3 E [[ Total Sobol index, quantity of interest estimated by both LOCO and PermuCATE.
Table I.
A.2. R- and pseudo-outcomee-rsst decomposition
In a scenario where the outcome is not deterministic: m=x)= + (a (())))((()))+++++ + x), given a CATE estimate f, the expectation of the pscudo-outcome risk can be formulated as,
E ElRpo X, E [(((Z) =(X))²)
((AAX))(A-A(X) ((1,
(X) - (x())
(0Y
m(A, x(X))
+ +((x) - =(X))
27
x(X(1 - T(X))
)
[e(X, A) T xX)) (X) - (X))
2
(T(X) +(X)) +E
+ 2
t T (X)(1-
risk
resculed noise term
=0, since B/c(A,
[5(X, A) xx)) (((X,A)(A~xXX)) +(X(1 - -
E -(x))² +E
(SI)
Similarly, as shown in Doutreligne & Varoquaux 2025, the R-risk RR(F X, A. can be re-written as,
B[RR(F, X, A,Y) = R ((V m(X)) - (4 xx)))()
EE ((A (x) (A (X))F((X))
= E ((A ((x))((=+(x))) + ((X) + [=(A,))]]
reweighted 7 risk
=0, since 3|E(A.X))=0
noise variance
= ==((1-=(x0)) ((11==(M)))
-(x((()(((X))-- ( + + + [((AA
(S2)
12

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
Both risks are therefore comparable in the sense that they are feasible when estimated values used for the quantities it Aa and that they involve combination of r-MSE and noise term (termed Bayes error in Doutreligne & Varoquaux 2025). However, the presence of the term π (1 x (XX) at the denominator in SI makes the quantity less stable and thus desirable for CATE estimator selection. For instance, in a scenario where the propensity score is poorly fitted and can take extreme values (close to or 1) a DR-learner might still provide reliable estimation of the CATE. However, the term x(X]j ) will likely take extreme values and lead to inconsistencies between the Rpo-risk and the oracle 7-MSE. While this consideration holds for model selection of CATE estimators and is consistent with the findings from Doutreligne & Varoquaux 2025, as shown in S1 we did not notice any difference for the estimation of variable importance. It is likely that noise terms cancel when taking the difference in the variable importance for either LOCO or CP I. The following figure depicts the comparable experimental results up to scaling factor for the scenario LD
EI( 7[X)( E(X,A)
PO-risk
R-risk
LOCO
LOCO
4
1.0
2
0.5
0
0.0
True positive
-2
0.5
False negative
10a
10a
10
102
10a
10
X2
PermuCATE
PermuCATE
X3
1.0
2
0.5
0
0.0
-2
-0.
102
10a
10
102
10a
10a
N
Figure S1. The R- and PO-risk provide, up to a scaling factor, the same measures of variable importance. This plot reproduces the results from Figure 1 (left panel) which depicts the importance measured for the truly important variables (X1, X2, X3) in the simulation LD from Hines et al. 2022. In addition, the same experiment was performed using the R-risk (right panel).
A.3. Proof of Proposition 3.2
Proof. The causal risk considered in this proof is the pstudo-outoome risk described in subsection A.2. Here again, XP.. denotes the features matrix with the jth feature sampled from the conditional distribution [Xr.sl ~ X x²xx(-3)
2
Tiast Σ
[( 8 A,Y) R(F, A,
(F(X)-((8P.))
<(X,A)
27
(()(( - π x(X))
T4tsst Nreet
Titest Meantt
E(X,A)
(r(x)- - -A(X))
(T))(1 - +(XI)
Noess Tisess
Durst
Renes
(r(x) (((x)=F(RE) ((x))=(8Pp))) -=(XPe))))
Σ (((x)- ((x)))
Tttst
Then the law of large numbers provides,
2
[[[[(X)- -
CPT
13

---

Mcasuring Variable Importance in Heterogencous Treatment Efffett with Confidence
Then using the consistency hypothesis for f. ensures that the second error term BI(r(X) - Borain- asymptotically converges to the true CATE.
vanishes
Finally, thanks to the consistency hypothesis on 2/ and Assumption 3.1 we can apply Proposition 3.3 from Reyero Lobo et al. (2025) which implies that the conditional distribution p(X+IXx(-)) is approximated without error, meaning that [Xpsls - 2i (x((5)) ind ej with Ej TL x(=)). The tower property then provides
E (x(X) ( - (XP..)) = E ( (7(x) ((xx)))-
Using that Xrj id x|x(-),, we finally obtain
E E (T(X) ((Xp)))) [x E 2E (E ((xy)x(()))) [x(--]]]]]]
2E v [(X))x([=1]
A.4. Proof of Proposition 3.4
Proof. Expectations are taken over the test set, conditionally to the training set Dtrain. The causal risk considered in this proof is the pseudo-outcome risk described in subsection A.2. The conditional notation is omitted for readability. Here again, Xr. denotes the features matrix with the jth feature sampled from the conditional distribution [Xr,yls ~ p(X*|x Additionally, we denote X P. it's estimate using the estimator ún for the conditional sampling step.
Similarly to Appendix A.3,
1 (B[R(T, X, A, A, B[R(T, xr.: A, (,
El PPPelIPPran
(T(X) (F(x))(((xr))
E
Perain
E ((xx)=+(x]] Pirain]
Developing the first term gives
E (T(x) 1000--(00 {(Xp)))
D Dtrain
= E
(T(X) ((Pps))
+E (x(Xp)) -((p))-(((Rpp)) -(Kpp)))
Derain
(A)
(1)
+ 2E
- (xx)0
(C)
The crossed-term (C) involves the estimation bias of r and v-j- This errors can be easily decomposed as
E (r(( - ((X,)) -(())) train] = [[[((()) - -XP:)) P(XP))))
- -(X.)) -(XP,))) D Perain
so we first have the eтrоr of T and then of ny.
Using the proof of Proposition 3.2 (A) is exactly 2E [V(((X))[X((3) = -207,, twice the total Sobol index. Also, (B) can be decomposed in,
E
- -((P))) (F(XP)))(({Ps)
Perain E [TT(PP))) - F(XP:))) +E (F(XP:)) - =((Pps))) Perainn
2E (T(Xps) - (Psaia)
14

---

Mcasuring Variable Importance in Heterogencous Treatment Effeett with Confidence
L.d. Then, using that XP ~ X, we have that E [(T(Xp)) F(XP.)) the error estimating r cancels out for PermuCATE, in contrast to LOCO, which depends on this error. Finally the cross term being a product zero-mean errors coming from independent optimizations can be neglected, leading to,
- Derain = E ( -(X))) Therefore,
E CeellPaanii = + E (f((r.) =((r,))) Derain
(9)
Lastly i being K-Lipschitz, get
(FXxs))=(8Pss)) Durann
SE K
XPs
Dirain
Then, recalling that the jth covariate Xyy is obtained by sampling from the conditional distribution, it corresponds to Bj (x(()))) + (X's - Vy (x were X' and X is its jth covariate. therefore comes that
E (f(Xp))-F(3) -(8P)))
2
Derain
< KE
( -Xx Derain
2K ) Derain
where we have used that X i.d. X'. Finally, denoting 10 1/2 E ||| - -2,11 | obtain,
E
V CPI Perain
v; + Op(lly - -B211²),
where the Op notation denotes bounded in probability.
A.5. Proof of Proposition 3.5 Proof.
E
Th Fooco Demin
= E (F =E(R(F A,Y - R(,, A, [erasn]]
=E((-=_((X-)) - -T-- - Duraia)
On the one hand, we have that
E (T - -x,
Perain]
=E (T(X)- - (X-))) Derain
F [T-r(X-1) ()( Perain
+I [[(((x) - -x-((x ((mmi))
where we have used that
E
-x
= [T (x - [(T(X) -T (X- Dirain] [
0.
We note that IE (((x)---s(x-1) is exactly E [V(((X))(x(--) = MJ, the total Sobol index. Combining both terms, completes the proof,
E
V
Dtrain
= W, + E (T- - Dtrain] E (TT(X) 2 (10)
w + Or([IT 7 -ー 2 ||| -
Were || 1 is the 1.2-norm with the fixed train set Drrain
15

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
A.6. Finite sample variance of LOCO and CPI in linear scenario
This subsection presents the analytical variance of the CPI and LOCO estimators in the linear setting it builds on top of results from Reyero Lobo al. 2025. In such scenario, assume the cate to be linear r(X) X B and consider linear estimators ((X)= X
A.7. LOCO
For LOCO, recalling Equation (10)
E
T
Loco Ptrain
= V + E
- Durain] - [(=(x) - -((x)) Perain] -F-;(x (x-
Considering linear models, we get:
Pirain = ( - -x -x-3((11)
E
(T- - -x
D train
Denoting - 3((-) - -}(-)) get,
(T-)(X - 1)) 2))² Ptrain = || ( 11
with = [X E x(
and | ||::
EA. Using a similar reasoning leads to,
E
T
D Perain]
Vi ||
LOCO
Finally, using assuming that independent samples are used to fit r and т-j (for instance using sample splitting similarly to Williamson et al. 2023) we get
var
T
LoCo
Dermin)) 1) var (144)-))13.. + var (14
A.8. PermuCATE
Regarding CPI, we recall from Equation (9),
D Vd T
CPIl
Derainl
(f(X)-r(Xx:)) + (f(Xr))- F(Zr,))) Dersin))
Concerning A, recall that the jth covariate of Xp, corresponds to 2y (x((-)) + - -,,(x(((0))) with X' 1l X. That all other covariates remain unchanged and that we assume 3x.. Consequently.
E (F(X..) =(KPp))) Dtrain
=E EE 33 (1)(x((1))
+ (X" (x - (x((-)) - -D,(x'
Ptrain)
= 33 (1s(x(-))) - - -2((x(-))) (x(((-1)) 2 E ((XPg) =(8P:))) ( train| = -33 ( x((-))) 2
Derain
Then considering Gaussian covariates, the model v is linear, giving 1j (x(x(-1) = x(-1) + +
ID (Derain
= 23 A
where A9 = 7, is the coefficient difference when estimating the linear function vy and Σ., EX To obtain the last line, we used that X" 11 Therefore when using the rescaled version T CPT we obtain,
var Cr
var(3>1177133)
16

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
A.9. Variance comparison in a linear scenario
An empirical comparison of the variance of LOCO and PermuCATE shown in Fig S2.
a
PermuCATE
LOCO
Vn.va . var
102
Vn
10
b 100
||1i/1 -O-
114B-1 -nl13
10-1
10-2
10-3
102
103
104
sample size
Figure S2. The lower variance of PermuCATE explained by finitc-sample estimation error. (a) /n times the variance of the importance estimate of X2 (2(2) on the LD dataset using LOCO and rmuCATT different sample sizes. (b) Dependence of different error terms on sample size. For PermuČATE: this noise can be attributed to the the error in the estimation of i IE,, [X,|X* -J] and captured by its || For Loco, it captured by the fluctuations the estimation the model coefficients of model introduced by the refitting procedure and is captured by 1A33,1133 The correlation structure from the LD dataset has been removed for clarity.
A.10. Datasets
Low Dimensional (LD) dataset
Taken from the work of Hines et al. 2022:
N(O, 1
0.5
(Xs, X6) (X3,
0.5
T(X) X, + 2X, + X3
π 7X)) = expit (-0.4X, + X2
~ Bernoulli(z(X))
(x( X33 X
Y px() AT(X)+
17

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
The High Dimensional Linear (HI) dataset Create the , S of important variables, respectively for the functions n, /o; 7 by randomly drawing without replacement dimp 10 important variables from the set (1.... dd of possible variables. We denote the vector of covariates X = (* za]. Covariates follow random normal distribution and with a correlation coefficient that can be modified.
d Σ
(X) expit
i $ S B ~ Rademacher(0.5)
if E S
A Rademacher variable takes value in {-1,1} with equal probability.
MO(X)=
Σ T
if ş S
3,10
g
(ademacher(0.5)
if E S
T(X)
1=1
if $ ST
87
BT
(ademacher((.55)
if E ST
We then construct,
A ~ (er)ouli(T(X)))
Y = (1 e) (X( + e Ar()) + ( (0,
Where E [0, 1] denotes the effect size which controls the strength of the treatment effect.
The High Dimensional Polynomial (HLL dataset features. Consequently,
This dataset has similar construction as HLL but uses polynomial
d
x (X) expit
Σ
if ST VJ if ST A AJ S
B
8 j.j.k
(.demacher(0.5))
Mo and 5 follow similar construction as in LD using the same polynomial features.
In addition, we also control the proportion of treated units by computing the pdh quantile of the propensity scores that we denote Qp in order to avoid extreme imbalances between the number units in the treated and control group.
A ~ Bernouli - Qp)
Y = (1 - e) (X + Ar(X) +N
Infant Health and Development Program (IHDP) This benchmark consists in data from a clinical trial studuying premature the effect of intervention population premature infants. The Intervention Group received home visits, attendance special child development center, and pediatric follow-up. It contains 747 subjects and 25 covariates (6
18

---

Mcasuring Variable Importance in Heterogencous Treatment Effeett with Confidence
continuous, 19 binary) describing infants and their mothers. It is imbalanced, with 139 treated and 608 controls, making the CATE estimation more challenging than in matched clinical trial. While covariates are real, the outcomes are simulated using.
Y ~ AN(exp w. (es (X + 1)
where is samples from described in the publication from Shalit et al. 2017. Similar to Shalit et al. 2017; Curth & Schaar 2021, used 100 repetitions of the simulations.
0,,,0.22 (0,0.1,0.2,03,, 0.4) with probabilities {0.6, 0.1,0.,, W and are offset terms as
A.11. Shyper-parameter search
optimized using nested cross-validation loop. For linear models, used the scikit- learn implementation RidgeCV for regression and Logist itsiiiiio ceegressersiicccooo onCV with range of penalization strength from 10-3 to 10 with 10 logarithnically spaced values. For the gradient boosting tree we used the implementations Hi Bos Ingllassiier and Btt.rrrieeettttttineeeeeeeeeeeee
All typer-parameters
respectively for regression and classification. After using randomized search for hyper- prperpparameters we used a learning rate of 0.1 (range explored: [10-3 3, 10 and a maximum number of leaves for each tree of 10 (range explored: [10, 100]).
A.12. Variable importance using grouping
Although not emphasized in the main text, the PermuCATE approach compatible with grouping strategies, which can be useful when dealing with highly correlated variables that lead to vanishing conditional importance. To address this, Chamma al. 2024 proposed extending the conditional permutation importance framework to groups of variables. To empirically demonstrate this, we increased the dimensionality of the IHDP dataset by adding noisy copies of the variables- Gaussian noise was added to continuous covariates and binary labels were flipped with 0.1 probability. Groups were formed by combining original variables with their noisy copies, and conditional importance was measured for these groups The previous metrics were adapted accordingly: type-1 error involved identifying group of null variables, and a true positive involved a group with at least one predictive variable. As shown in Figure S3, which is analogous to Figure 4 but with dimensionality on the y-axis, increasing dimensionality negatively impacts performance since no additional information provided. Interestingly, doubling the dimensionality did not significantly harm power AUC for predicting importance but did reduce the PEHE. As before, CATENet models exhibited higher type-1 error rates, especially with increasing dimensionality, likely due to their highly stochastic optimization process.
A.13. Non-Gaussian variable distributions
In this section, propose an additional experiment where revisit the simulation study presented in Figure 2 also to vary the complexity of the conditional distributions. To do so we use simulations inspired by (Hollmann et al., 2025), we sample the covariates using a latent variable model. We first sample latent variables from mixtures of Gaussians and then generate observed covariates through a non-linear transformation with interaction terms of the latent variables. This results in non-Gaussian covariates with complex conditional distributions, illustrated in Figure S5, a) presents the marginal distributions, b) the kurtosis of generated variables (blue) compared with the kurtosis for Gaussians (orange) for comparison, and e) their correlation matrix. Similarly to Figure 2, we generate CATE functions, one linear and one non-linear and compared the variance of Pe rmuCATE and LOCO for three different models varying sample sizes. The results shown in Figure 56 reveal lower variance for PermuCATE compared to LOCO similar to Figure 2. Our interpretation is that the complexity of covariates distributions also affects LOCO, by making the CATE harder to estimate.
19

---

Mcasuring Variable Importance in Heterogencous Treatment Effects with Confidence
a
C
26
KH
52
78
104
130
-0.8
-0.6
-0.4
0.4
0.6
negative PEHE
power (a = 0.05)
b
d
26
model
52
CF
78
H0111
CATENet
TabPFN
104
HHI
130
0.6
0.8
1.0
0.0
0.1
0.2
AUC
type-l error (a = 0.05)
Figure S3. Comparison of variable importance methods the IHDP benchmark using groups of variables. The CATE was estimated using Causal Forest (CF, Athey & Wager 2019), deep neural network (CATENet, Curth & Schaar 2021). For each learner, variable importance was estimated with PermUCATE. The dimensionality was artificially increased by creating (a) Displays the negative Precision in Estimating Heterogeneous Effects (PEHE) actoss the models. (b) presents the AUC for estimating variable importance, which be framed as a classification task. The predicted importance represented by 1 - p-value, which is subsequently utilized to compute statistical power (c) and the type-1 error (d) at the significance level alpha=0.05.
20

---

Mcasuring Variable Importance in Heterogencous Treatment Effeets with Confidence
a
C
26
52
78
KH
104
130
-0.8
-0.6
-0.4
0.2
0.4
negative PEHE
power (a 0.05)
b
d
26
model
52
CF
78
HH
CATENet
TabPFN
104
HH
130
0.6
0.8
1.0
0.0
0.2
0.4
AUC
type-I error (a 0.05)
Figure S4, Same Figure S3 but using Leave One Group Out (LOGO)
Veriable 20
Variable
Variable x2
Variable x3
Variable x4
Kurtosis test
70
53
132
60
Mattoo
xp-4
20-
30
11-1
muttmdde
normall
Variable 25
Variable x6
Variable x7
Variable x8
Variable 29
50-
Correlation Matrix
1.00
3.72
6.50
012
025
222
Figure S5. Simulation with non-Gaussian variable distributions a) presents the marginal distributions, b) the kurtosis of generated variables (blue) compared with the kurtosis for Gaussians (orange) for comparison, and c) their correlation matrix.
21

---

Mcasuring Variable Importance in Heterogencous Treatment Effects with Confidence
linear
b
non-linear
Ridge
variance
PermuCATE
LOCO
10-1
104
d
10
CateNet
variance
10a
10-4
10-1
f
100
TabPFN
105
variance
10
102
103
104
103
104
sample size
sample size
Figure S6. Same Figure 2 but using the variable distributions described in Figure $5
22


Reviews:
[{"review_id": "qXjRuJrS5R", "reviewer": "Reviewer_jT2T", "summary": "The paper proposes a variable importance measure (VIM) to understand the variables driving the conditional average treatment effect (CATE) function. The measure is based on the principle of conditional permutation importance where the variable of interest is permuted while keeping matching its conditional distribution and then checking the impact on CATE through a quantity known as total Sobol index. The work studies the estimation error, finite-sample variance, and type-I and II error of the proposed VIM to an existing measure named LOCO from Hines et al. 2022. Main contributions are to develop a natural permutation-based VIM and to show that it performs favorably in extensive experiments.\n\n### update after rebuttal\n\nI have read the response and other reviews. The response clarifies my concerns on interpretation of the Sobol index and performance of the method for more complex covariate distributions. Thanks for including an additional experiment.\n\nHowever, the theoretical claims on estimation errors/variance for LOCO vs proposed method require substantive changes based on the discussion with Reviewer wUZa. Asymptotically, LOCO seems to be better whereas in finite data experiments proposed methods is consistently better. The authors provide an hypothesis for why finite sample performance could be better, namely, MSE error in tau is worse. I would suggest the authors to carefully check the hypothesis in experiments and thus provide an explanation grounded in their theoretical analysis.\n\nMore importantly, choosing between LOCO and proposed method is not straight-forward, a concern that has not been addressed. Authors should discuss tests for the assumption or procedures to choose between the two methods. More examples of settings where the data might favor LOCO will help.\n\nGiven the experimental results on the proposed method, I am positive that the work is significant. The theoretical justification and guidance for using the method could be improved. Hence, I raise my score to 3 weak accept.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thorough review. The comments and questions are relevant and will help us improve the paper. \n# Dependence of PermuCATE on the estimate of $\\tau$\n- We refer the reviewer to the response provided to Reviewer wUZa, in the section \"Dependence of PermuCATE on the estimate of $\\tau$ \". We explain why second-order error terms cancel out in the estimation of PermuCATE, and propose to rephrase proposition 3.4 to include first-order error terms.\n - We agree with the comment that eq 7 does not explicitly point to a worse dependence on the dimension. We plan on rephrasing the L342 to clarify that this is an empirical observation.\n# Complexity of the conditional distribution estimation\n - The IHDP benchmark uses real-world data, with non-Guassian continuous features and imbalanced categorical (up to 4 categories) features. Because the ground truth for the conditional distribution is not available, the validity of Assumption 3.1 cannot be assessed. However, Figure 4 shows that PermuCATE outperforms LOCO with better statistical power and AUC for detecting the important features while controlling the type-1 error on this dataset. \n - In addition to this benchmark, we propose an additional experiment where we revisit the simulation study presented in Figure 2 also to vary the complexity of the conditional distributions. To do so we use simulations inspired by [1], we sample the covariates using a latent variable model. We first sample latent variables from mixtures of Gaussians and then generate observed covariates $X$ through a non-linear transformation with interaction terms of the latent variables. This results in non-Gaussian covariates with complex conditional distributions, illustrated in this figure: https://imgur.com/a/JqT04in In a), presents the marginal distributions, b) the kurtosis of generated variables (blue) compared with the kurtosis for Gaussians (orange) for comparison, and c) their correlation matrix. Similarly to Figure 2, we generate 2 CATE functions, one linear and one non-linear and compared the variance of PermuCATE and LOCO for three different models at varying sample sizes. The results shown in this figure: https://imgur.com/a/t0qhb61 reveal a lower variance for PermuCATE compared to LOCO similar to Figure 2. Our interpretation is that the complexity of covariates distributions also affects LOCO, by making the CATE harder to estimate. \n\n[1] Hollmann et al. Accurate predictions on small data with a tabular foundation model. Nature 2025\n# Description of the total Sobol index\n - We refer the reviewer to the response provided to Reviewer VGy6, in the section \"Description of the total Sobol index\". We provide an intuitive presentation of the total Sobol index. We will integrate this presentation into the final version of the paper.\n# Discussion of model-x knockoff\nThe proposed approach shares similarities with model-X knockoffs (KO) as both model the conditional distribution of covariates, assuming this task is easier than estimating the quantity of interest (e.g., CATE). This point is discussed in Remark 3.3 of the main paper. A key difference is that model-X KO is designed for variable selection, whereas PermuCATE and LOCO focus on variable importance estimation. The latter provides richer information by measuring importance (total Sobol index) rather than just binary selection. Additionally, KO handles multiple testing, while our method ensures type-I error control. Also, constructing KO variables requires complete exchangeability, a stricter condition than the conditional independence needed for CPI. The conditional randomization test (CRT) from Candès et al. 2018 is more comparable to our approaches for individual conditional independence testing but is much more computationally expensive.\n\nBesides the work from Hines et al. 2020, we are unaware of other methods that provide model agnostic variable importance for CATE. \n# Extension to groups of variables\n - As mentioned in L235 and demonstrated in Appendix A.12, CPI-based approaches can be extended to groups of variables. CPI with grouping in the context of prediction problems was formally treated in https://doi.org/10.1609/aaai.v38i10.28997. There is no particular complexity when performing multi-output regression. For each variable of the group, the conditional distribution can be estimated independently and in parallel: $\\forall j\\in G$ estimate $\\nu_j = \\mathbb{E}[X^j|X^{-G}]$, with $X^{-G}$ the complementary of group G. \n# Ranking of important variables\n - The ranking of important variables is indeed richer information than the AUC for binary classification of important variables. However, except for simple scenarios (Gaussian covariates and linear CATE) such as those presented in Figure 1, the ground truth importance (or ranking) is not known. For IHDP, two problems prevent the analytical computation of the total Sobol index: the true conditional distribution is unknown, and the CATE function is non-linear.", "claims_and_evidence": "- Claims on inferential properties of the VIM are demonstrated convincingly by experiments.\n\n- Claims on benefits of the proposed method are not very rigorous. The comparison of bias between LOCO and PermuCATE in eps (5) and (7) is a bit misleading. Proposition 3.4 hides the fact that PermuCATE also depends on how good is the CATE estimate tau since it is assumed to be consistent and Lipschitz. Therefore, both methods are susceptible to estimation errors in tau. The eps (5) and (7) hide dependence on sample sizes and assume consistency for tau. Therefore, it is misleading to say they are finite sample analyses. Relatedly, eq (7) does not clearly point to a worse dependence on dimension of covariates. The claim in line 342 is not supported by the theoretical results. Both (5) and (7) will depend on dimension.", "methods_and_evaluation_criteria": "- The permutation-based variable importance is a natural concept and the method skillfully develops it for CATE.\n\n- Experiment setup including the metrics and baselines are relevant.", "theoretical_claims": "I read proof for Proposition 3.2 carefully and skimmed the proofs for Proposition 3.4 and 3.5. For 3.2, the cited result from Reyero Lobo et al. 2025 shows convergence in Wasserstein which is written as an exact equality in the proof. I think the main result is ok, but the statement should be written carefully.", "experimental_designs_or_analyses": "- I checked the design for synthetic experiments which control for complexity of conditional outcome function. \n\n- The complexity of conditional covariate density function is not tested adequately. Since Assumption 3.1 on covariate density is important to the results, I think the experiments need to vary complexity of density function as a way to check sensitivity to the assumption, including settings when the assumption is not true.", "relation_to_broader_scientific_literature": "- The method is new and experimental comparison of finite-sample behavior against LOCO are impressive.\n\n- Please discuss the relation to model-x knock off literature in more detail since the methods similarly model the covariate density. Does there exist CATE importance measure from the literature that are natural baselines?", "essential_references_not_discussed": "The paper can discuss other methods that are based on the principle of leave one covariate out such as Zhang and Janson 2020 and more papers on VIMs based on model-x knock off if any.\n\nLu Zhang, Lucas Janson. Floodgate: inference for model-free variable importance. arXiv 2020 https://arxiv.org/abs/2007.01283", "other_strengths_and_weaknesses": "# Strengths\n\n- Method is conceptually simple and inherits the statistical inference guarantees on type-I error from prior results.\n\n- Proposed method provides a useful alternative to existing variable importance measures by leaning more towards modeling the conditional covariate distribution instead of the conditional outcome distribution. Therefore, the method might be more suitable in some applications.\n\n- Experimental validation is thorough. Authors rigorously estimate the importance measures by using flexible CATE-learners wrapped in ensemble methods like super learner and use cross-fitting to reduce variance. Throughout the results report confidence intervals or p-values over sufficiently many repeated samples of the data.\n\n\n# Weaknesses\n\n- It is unclear when to use LOCO vs PermuCATE since they have similar statistical inference properties but differ in assumptions. Consider discussing a test for checking Remark 3.3 that states that the conditional covariates are easier to model than conditional outcome function. Was it the case for IHDP dataset? The claim that modeling conditional covariate distributions are easier for some applications should be more carefully demonstrated or discussed through citing literature.\n\n- The presentation can be improved. I felt that CPI could be explained in more detail before the methods to give the context for comparisons to LOCO. Please present the statistical guarantees for CPI. Please introduce total Sobol index and discuss why it is a good measure for CATE.\n\n- Compared to existing method LOCO, the proposed method is not readily extensible to computing variable importance for variable subsets. The nuisance parameter \\nu will be a challenging multivariate regression, whereas LOCO still requires a regression on univariate response. Computing importance for subsets is important to handle highly-correlated variables and to limit the computation since often variable subsets can be grouped together and treated as one variable.\n\n- Experiments on IHDP data did not test whether variables are ranked in the true order. AUC only tests for ranking between important and not important variables. A more prevalent use of variable importance measures is to identify which of the important variables are the most important. Please consider reporting metrics like precision at k or Kendall tau.", "other_comments_or_suggestions": "Minor comments, no response is requested for the following\n\n- Some terms in the introduction were not introduced. Please describe what is the purpose of variable importance measures, in what way they help in interpretation, feasible risk, and what is the outcome or effect in CATE in the biomedical applications discussed in introduction. Explain how showing variable importance for CATE will help in biomedical applications. \n\n- Please discuss CPI in some more detail in related work or introduction to make the reader appreciate its importance.\n\n- Please explain the remark on total Sobol index in line 135.\n\n- Consider labeling Eq. (3) as R(…)=…\n\n- I really like the Figure 1. It conveys information on effect of sample size on VIM, comparison to true value, and hypothesis test results quite clearly.\n\n- CPI and PermuCATE are used interchangeably, CPI in notation and PermuCATE in text. Please consider using the same name everywhere.\n\n- Please define support in line 370. Does it mean variables used in the true outcome function?\n\n- Please complete the sentences at the end of Remark 3.3. It is evident they support the claim, however, can be written more explicitly.\n\n- Motivation for the method is from high-dimensional settings like SNPs data whereas the IHDP datasets in evaluation study has moderate number of covariates.\n\n\n\nTypos\nLine 066 an\nLine 162 estimated\nLine 333 pseudo", "questions_for_authors": "- Please discuss validity of the claims on variance of the two methods.\n\n- Please intuitively describe total Sobol index.\n\n- Please provide more evidence either from literature or IHDP data that covariate density is easier to model or suggest when to apply the proposed method.\n\n- Please clarify whether the experiment setup checks Assumption 3.1 systematically.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "VPySi5WzKh", "reviewer": "Reviewer_wUZa", "summary": "The paper proposes a new explainability method for the conditional average treatment effect (CATE) estimators, namely, PermuCATE. PermuCATE measures global features importances and is based on a conditional permutation importance (CPI) methods. The authors demonstrated that PermuCATE aims to estimate the expected conditional variance of CATE, given the subset of covariates (total Sobol index). Also, they showed, that under additional assumptions,  PermuCATE achieves lower variance than the existing leave-one-covariate-out (LOCO) method. Ultimately, the authors empirically demonstrate the superiority of their method in terms of statistical power and precision.\n\n\n**Post-rebuttal response**: I thank the authors for addressing the majority of my concerns. Still, I have two remaining concerns: (1) the presentation and (2) the efficiency of the method. \n1. The theoretical results (namely, Propositions  3.4 and 3.5) need to be carefully revisited, especially, considering the large number of corrections and edits during the rebuttal. Here, I also mean streamlining the notation. \n2. I read the latest authors' response and now it looks like the proposed method is asymptotically strictly worse than the existing LOCO: the proposed method contains the error of the first order wrt. the estimated $\\hat{\\tau}$ whereas the  LOCO has the same error squared (= hence, less sensitive to the estimation of $\\tau$). Therefore, I don't see how the method with the apriori asymptotically worse performance can be preferred in practice given that we cannot reliably choose or validate Sobol index predictions simply from observational data.  \n\nHence, I keep my current score and encourage authors to further adjust their method so that it achieves the lowest variance (= semi-parametric efficiency bound).", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thorough review. The comments and questions are relevant and will help us improve the paper. \n# Dependence of PermuCATE on the estimate of $\\tau$\n - To clarify the dependence of PermuCATE on the estimation of $\\tau$, we first provide some intuition: LOCO and PermuCATE quantities are estimated by computing a difference between two risks (eq 1 and 4). Contrary to LOCO, PermuCATE uses the same estimator ($\\hat{\\tau}$) in both terms of the difference, leading to estimation terms that cancel each other out. This intuition is demonstrated in Supplement A.4. Specifically, in L715, the second-order $\\tau$ estimation term cancels out using that $\\\\mathbb{E}[(\\tau (X_{P,j})-\\hat{\\tau}(X_{P,j}))^2|\\mathcal{D_{train}}]=\\mathbb{E}[(\\tau (X)-\\hat{\\tau}(X))^2|\\mathcal{D_{train}}]$ because by construction $X_{P,j}\\overset{i.i.d.}{\\sim}X$.\n - The only term involving the estimation error is term C (L697), a first-order error term: $O(\\tau-\\hat\\tau)$. We initially argued that it is negligible compared to second-order terms and proposed to ignore it. However, we agree that clarifying the dependence on the estimation error of $\\tau$ is important and will rephrase proposition 3.4 accordingly. \n - The Lipschitz assumption aimed to show that in equation 9, the bias term carries the error in estimating the conditional distribution $\\nu_j$. As confirmed by experiments, this assumption is not restrictive.\n\nTo clarify these last two points, we propose to rephrase the proposition using the intermediate step (eq 9) of the proof in A.4.\n\n## Proposition 3.4:\nUnder Assumption 3.1, for a consistent CATE estimator $\\hat{\\tau}$, the finite sample estimation of the importance for the $j^{\\text{th}}$ covariate is: \n$$\\frac{1}{2}\\mathbb{E}[\\widehat{\\Psi_{\\mathrm{CPI}}^j}|\\mathcal{D_{train}}] - \\Psi_j^* = \\mathbb{E}[(\\hat{\\tau}(X_{P,j}) - \\hat{\\tau}(\\widehat{X_{P,j}}))^2| \\mathcal{D_{train}}] + O_P(\\tau - \\hat{\\tau})$$\n\n**Remark**: While the first term involves the CATE estimate $\\hat{\\tau}$, it does not involve the estimation error.\n\n## Corollary:\nUnder the additional assumption that $\\hat{\\tau}$ is Lipschitz continuous and $\\hat\\nu$ is consistent, the right term becomes $$O_P(||\\nu_{j} - \\widehat{\\nu}_{j}||^2) + O_P(\\tau - \\hat{\\tau})$$\n\n# Variance comparison between LOCO and PermuCATE\n - Indeed, LOCO is asymptotically efficient. However, please note that this convergence rate explains its asymptotic behavior and relies on the convergence rates of the models used. In finite samples, the estimation error of these models drives the variance of PermuCATE and LOCO. As shown in eq. 5 and 7, PermuCATE has a smaller dependence on the complex model $\\tau$ and, therefore, a smaller finite sample variance. This was experimentally confirmed in Figure 2 and the additional simulation study mentioned below.\n\n# Subset used for conditional perturbations\n - The goal of the conditional perturbation is to sample from the conditional distribution $X_{P,j} \\sim p(X^j|X^{-j})$. Under assumption 3.1, each covariate can be decomposed as a function of the other covariates plus an additive noise term $X^j =\\nu_j(X^{-j}) + \\epsilon_j$. As presented in algorithm 1, the function $\\nu_j$ is estimated on the training data. This avoids overfitting, which could lead to perfectly predicting $X^j$, not perturbing the data and leading to vanishing importance. Then, to sample the additive noise term, we use permutations of the residuals: $\\epsilon_j = \\mathrm{shuffle(X^j - \\nu_j(X^{-j}))}$. We would like to insist that no parameter estimation is involved in the permutation of the residuals and, therefore, no information leakage from the test set. This also maintains scikit-learn API compatibility by using the `fit` method (with training data) to estimate $\\nu_j$, and the `score` method (with test data) to sample perturbations.\n# Benchmarks that do and do not satisfy Assumption 3.1\n -  The simulation scenarios LD, HL, and HP (see datasets in the main paper) used multivariate Gaussians, which satisfied this assumption. \n - We refer the reviewer to the response provided to Reviewer jT2T, section \"Complexity of the conditional distribution estimation.\" We clarify that the IHDP benchmark uses real-world data, where this assumption cannot be verified. Furthermore, we provide an additional experiment where the dependency structure of the covariates is more complex than that of multivariate Gaussians. In both cases, PermuCATE outperforms LOCO, supporting the validity of assumption 3.1.\n - We also refer to the section \"Discussion of model-x knockoff\" in response to Reviewer jT2T, where we discuss that this assumption is common in the literature and review references. This assumption ensures valid conditional sampling from $p(X^j|X^{-j})$. However, any other sampling method could be used as a drop-in replacement. \n\n# Factor 2 in risk difference\n - As demonstrated in Proposition 3.2, PermuCATE estimates the total Sobol index up to a factor of 2.", "claims_and_evidence": "The major claims are supported by proofs and the authors provided an extensive empirical evaluation of their method.\n\nHowever, I have several major concerns regarding the claimed properties of PermuCATE:\n1. I find it hard to believe that the method achieves lower variance than LOCO, given a semi-parametrically efficient variant of LOCO  [1] offers the asymptotically efficient estimator of the target total Sobol index. I encourage the authors to provide a wider discussion on this (also, take into consideration issue 1 in Theoretical Claims, as it seems that Prop. 3.4 is missing an error term).\n2. The residuals for conditional perturbations of covariates are estimated on the same subset (D_test) that is being used to estimate the global variable importance. Doesn’t this introduce estimation bias or compromise statistical tests for variable significance? \n\nIf all of the issues from above and below could be resolved, I would be happy to raise my score.\n\nReferences:\n- [1] Hines, O., Diaz-Ordaz, K., and Vansteelandt, S. Variable importance measures for heterogeneous causal effects. arXiv, 2022. doi: 10.48550/arxiv.2204.06030.", "methods_and_evaluation_criteria": "The chosen baseline methods are relevant and make sense. However, the authors did not discuss whether the chosen benchmarks satisfy Assumption 3.1 needed for the method’s consistency. Thus, I would like to see a clear separation between the benchmarks that do and do not satisfy Assumption 3.1.", "theoretical_claims": "I found several inaccuracies in the theoretical statements:\n1. In Prop. 3.4, it is hard to believe that the bias of estimating the total Sobol index does not depend on the error of estimating the ground truth CATE, namely, $O_P(\\lVert \\tau - \\hat{\\tau} \\rVert^2)$. Upon checking the proof in the Appendix, it seems that indeed, the authors missed this term while transitioning between the proofs of Prop. 3.2 and 3.4.\n2. I don’t think the additive noise assumption is necessary for the outcome, as introduced in Appendix A.2 (line 575). The general DR-/R-learners do not make such an assumption [1].\n\nReferences:\n- [1] Morzywolek, Pawel, Johan Decruyenaere, and Stijn Vansteelandt. \"On a general class of orthogonal learners for the estimation of heterogeneous treatment effects.\" arXiv preprint arXiv:2303.12687 (2023).", "experimental_designs_or_analyses": "See Methods And Evaluation Criteria.", "relation_to_broader_scientific_literature": "The paper suggests a lower variance alternative to the existing CATE variable importance method, LOCO [1]. Yet, this comes at a cost of the additional assumption (Assumption 3.1), namely the additivity of noise for conditional distributions of covariates. In my opinion, this is a pretty strong assumption (although the authors argued differently in Remark 3.3), which limits the application of PermuCATE in practice.\n\nReferences:\n- [1] Hines, O., Diaz-Ordaz, K., and Vansteelandt, S. Variable importance measures for heterogeneous causal effects. arXiv, 2022. doi: 10.48550/arxiv.2204.06030.", "essential_references_not_discussed": "To the best of my knowledge, all the most important works were discussed here.", "other_strengths_and_weaknesses": "Nothing new here.", "other_comments_or_suggestions": "- Did you mean “pseudo-outcomes” instead of “potential outcomes” in line 152?\n- I would be more precise in lines 114-115 (2nd column) by defining what “comparable risks” mean. It is important to mention, that both the DR-learner risk (=PO-risk) and R-learner risk yield the same CATE estimators **only** in population and when the ground-truth nuisance functions are used [1].\n- I had a hard time distinguishing the setting and related work. I suggest the authors split Sec. 2 into two parts for better readability. \n- Some notation is not properly defined or introduced (e.g., $\\beta_j$ in Eq. 6 or $\\Psi_{\\text{LOCO}}$ in Eq. 8. \n\nReferences:\n- [1] Morzywolek, Pawel, Johan Decruyenaere, and Stijn Vansteelandt. \"On a general class of orthogonal learners for the estimation of heterogeneous treatment effects.\" arXiv preprint arXiv:2303.12687 (2023).", "questions_for_authors": "1. Is there a reason why the risk difference (line 13 of Alg. 1) is divided by 2? It seems like this only complicates the derivations.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "sOshbxYH2f", "reviewer": "Reviewer_GH6j", "summary": "The submission describe a variable importance method (PermuCATE) generalizing CPI (Chamma et al. 2023), \na theoretical analysis of PermuCATE is performed showing the behaviour of the estimator in finite sample. \nExtensive experiments are implemented over a variety of datasets comparing PermuCATE against LOCO.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer. The comments and questions are relevant and will help us improve the paper. \n\n# Previously published references \n - We would like to highlight that important references for this work have been published. Specifically, Chamma et al. 2023 in NeurIPS and Verdinelli & Wasserman, 2023 in Statistical Science. We will update the bibliography to include the final published version of the Verdinelli & Wasserman paper instead of the preprint.\n\n # Connection between CATE and variable importance\n  - Following the distinction made by Verdinelli & Wasserman, 2023, our work focuses on population-level variable importance. Estimating how important a particular feature is in the underlying data-generating process $\\mu$ is done by estimating the total Sobol index using the LOCO procedure. For a variable $j$, $\\Psi_j^* = \\mathcal{L}(y, \\mu{-j}(X^{-j})) - \\mathcal{L}(y, \\mu(X))$, where $\\mu_{-j}$ is a sub-model estimated without the variable $j$.\n\n\n - In the present work, the focus on causality is motivated by the type of problem we are interested in: interventional data where patients are assigned to treatment or control groups. In Verdinelli & Wasserman's vocabulary, the underlying data-generating process we wish to study is the individual treatment effect, also known as the conditional average treatment effect (CATE), $\\tau(x) = \\mathbb{E}[Y(1) - Y(0) | X=x]$.\n\n\n - Extending the approach from Verdinelli & Wasserman to CATE estimation is not straightforward because the ground truth CATE is not observable — patients are either in the treatment or control group, not both (the fundamental problem of causal inference). In our main paper, we show in Equations 1 and 4 that instead of the loss difference used in Verdinelli & Wasserman, a difference of feasible causal risks can be used to consistently estimate the total Sobol index for the CATE.\n\n # Motivation and problem setting\n - The motivation for this approach is to understand, at the population level, what drives heterogeneity in the treatment effect. For example, in a clinical setting, a given treatment might lead to adverse events in a subset of patients (e.g., those with a specific genetic profile or lifestyle ...). The proposed approach would then allow us to identify the risk factors (e.g., genetic profiles or lifestyle ...) associated with a higher risk of adverse events when treated.\n\n\nWe appreciate the reviewer's comment and will include this clarification and discussion in the terms of Verdinelli & Wasserman in the final version of the paper.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "yes they make sense", "theoretical_claims": "no i did not check proofs in the appendix", "experimental_designs_or_analyses": "yes the experimental design is appropriate and sound", "relation_to_broader_scientific_literature": "the method is an extension of previous work (Chamma et al. 2023,  Verdinelli & Wasserman, 2023 and  Hines et al., 2022). \nThe basic references of this work are all preprints,", "essential_references_not_discussed": "I am not aware of previous references not discussed.", "other_strengths_and_weaknesses": "while the paper is well written in general I fell it should provide a better description and connection between CATE and variable importance. \n\nFor instance Verdinelli & Wasserman, 2023 make a distinction between 3 types of variables importance: population  (how important is age in a regression function?), algorithmic (how much did age affect the estimated value of the regression) and causal (ow would the outcome have\nchanged if Mary had been 5 years younger). \nThe present submission deals with the causal intepretation, being it the CATE. While  Verdinelli & Wasserman, (2023) focus on the population case. This is justified by the authors of this submission citing Hines et al., 2022 (which is a preprint not published yet?), I skimmed Hines et al., 2022 but I could not pinpoint where the connection can be made. \nMeaning, how feature importance of the CATE tell us something about effect modification or differences in the causal effect? maybe this is naive but I have intuitive feeling that is not the same.", "other_comments_or_suggestions": "no", "questions_for_authors": "For instance Verdinelli & Wasserman, 2023 make a distinction between 3 types of variables importance: population  (how important is age in a regression function?), algorithmic (how much did age affect the estimated value of the regression) and causal (ow would the outcome have\nchanged if Mary had been 5 years younger). \nThe present submission deals with the causal intepretation, being it the CATE. While  Verdinelli & Wasserman, (2023) focus on the population case. This is justified by the authors of this submission citing Hines et al., 2022 (which is a preprint not published yet?), I skimmed Hines et al., 2022 but I could not pinpoint where the connection can be made. \nMeaning, how feature importance of the CATE tell us something about effect modification or differences in the causal effect? maybe this is naive but I have intuitive feeling that is not the same.\n\nMoreover it is not completely clear to me the problem setting, especially the goal of this? the objective is to obtain causal estimates on effect modification?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nUC04FwKyM", "reviewer": "Reviewer_VGy6", "summary": "This paper tackles measuring variable importance in conditional average treatment effect (CATE) functions. One of the few current approaches consists in applying the LOCO method to CATE estimation ; as the CPI method is an alternative to LOCO, authors propose applying CPI to CATE estimation. They prove the consistency of both CPI and LOCO to the total Sobol index for a specific risk function, where the variance of CPI is driven by its less error-prone missing covariate estimation procedure. They empirically show that CPI leads to faster convergence and lower variance than LOCO on CATE variable importance measure problems.\n\n(EDIT : updated my score following the rebuttal)", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer. The comments and questions are relevant and will help us improve the paper.\n\n# Discussion of assumption 3.1\n - Assumption 3.1, that a given covariate can be decomposed as a function of the other covariates plus an additive independent noise term is quite common in the literature. It was especially discussed in (Candlès et al., Panning for Gold, 2017), introducing the model-x knockoff framework. In section 1.3, the authors write, \"assuming we know everything about the covariate distribution but nothing about the conditional distribution $Y |X1, . . . , Xp$\". We discuss the relevance of this assumption in remark 3.3 of the main paper. It is motivated by plausible scenarios where, for instance, the relationship between the covariates is much simpler than the relationship between the covariates and the outcome. For example, in genetics, the relationship between single nucleotide polymorphisms (SNPs) is easier to model than the relationship between SNPs and a disease. Also, when a large number of unsupervised data is available, learning the conditional distribution of the covariates can be easier than learning the outcome, for which only a few labeled examples are available. Finally, a temporal aspect may also be considered, where the covariates are observed simultaneously while the outcome is observed later.\n\n\n - This assumption was also made in previous published work on CPI (Chamma et al, 2023, NeurIPS), which also revealed that in the context of prediction tasks, CPI was a robust method for variable importance estimation.\n\n\n - We also refer the reviewer to the response provided to Reviewer jT2T, section \"Complexity of the conditional distribution estimation.\" We clarify that the IHDP benchmark uses real-world data, where this assumption cannot be verified. And present an additional simulation experiment where the dependency structure of the covariates is made more complex, using a latent variable model. PermuCATE still shows a lower variance than LOCO in this scenario. \n\n\n# Description of the total Sobol index \n - The total Sobol index is a well-studied quantity (see, for instance, [1, 2, 3, 4]) that measures the influence of a variable or group of variables on the output of a (possibly complex) model. For a model $\\mu$ predicting an outcome $y$ given covariates $X$, the total Sobol index of a variable $j$ (or group) is given by: $$\\Psi^* = \\mathbb{E}[\\mathbb{V}[\\mu(X) | X^{-j}]] = \\mathbb{E}[(y-\\mu_{-j}(X^{-j}))^2]-\\mathbb{E}[(y-\\mu(X))^2],$$\nwhere $\\mu_{-j}$ is the model where the variable $j$ is removed. \n\n - By writing the total Sobol index as a loss difference, we can see that it measures the loss increase when removing the variable $j$ from the model, which provides a measure of the importance of the variable. It can also be seen as an unnormalized generalized ANOVA (difference of $R^2$) [5]:\n\n$$\\Psi^* = \\mathbb{V}(y)\\left[\\left[1-\\frac{\\mathbb{E}[(y-\\mu(X))^2]}{\\mathbb{V}(y)}\\right]-\\left[1-\\frac{\\mathbb{E}[(y-\\mu_{-j}(X^{-j}))^2]}{\\mathbb{V}(y)}\\right]\\right] .$$\n\n - While this is true for predictive models, extending this definition to CATE estimation is not straightforward due to the fundamental problem of causal inference, the ground truth CATE is not observable (in clinical trials, patients are either in the treatment or control group, not both). \n\n - It can be seen from the above equation that LOCO is a plug-in estimator of the total Sobol index, where the model $\\mu$ is replaced its finite sample estimate $\\hat{\\mu}$. We show in the paper why this plug-in estimator is not the best choice for estimating the total Sobol index, especially in finite sample.\n\nWe will integrate this intuitive presentation of the total Sobol index in the final version of the paper.\n\n- [1] Sobol, IM. \"Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates.\" Mathematics and computers in simulation, 2001\n- [2] Bénard, C. et al. Mean decrease accuracy for random forests: inconsistency, and a practical solution via the Sobol-MDA, Biometrika, 2022\n- [3] Hooker, G. et al. Unrestricted permutation forces extrapolation: variable importance requires at least one more model, or there is no free variable importance. Stat Comput, 2021\n- [4] Williamson, B. et al. A General Framework for Inference on Algorithm-Agnostic Variable Importance. Journal of the American Statistical Association, 2022\n- [5] Williamson, B et al. Nonparametric variable importance assessment using machine learning techniques. Biometrics, 2021", "claims_and_evidence": "Overall, the method seems well-justified ; I'll just have a slight concern for Assumption 3.1 that looks like the backbone of the method and its advantages and that I find lacks a bit of justification. Does it come from the past literature on CPI? How is it generally justified?", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "I did quickly and did not find any issues", "experimental_designs_or_analyses": "Yes, and everything looks good", "relation_to_broader_scientific_literature": "As the LOCO method was applied to CATE estimation, and the CPI method was not, the submission fills this gap. It does justify theoretically and experimentally the use of CPI compared to LOCO.", "essential_references_not_discussed": "Not to my knowledge", "other_strengths_and_weaknesses": "The paper is original in the sense that it is a combination of existing ideas ; it might be seen as a relatively low-hanging fruit but the theoretical and experimental evaluations done by authors increase the significance of their work. The paper is generally clear, except on the point with Assumption 3.1 above and the total Sobol index, which seems to a critical notion of the paper (see questions below)", "other_comments_or_suggestions": "I strongly suggest you explain in greater depth the total Sobol index, eg where it comes from and why it is very appropriate for variable importance. It seems like a critical quantity in your work, as the ground-truth measure to which either CPI and LOCO is expected to converge. You might explain why it is to be taken as a ground-truth.\n\nI do not see any typos.", "questions_for_authors": "As said previously, can you please give details (literature + relevance) for :\n\n1) Assumption 3.1\n\n2) The total Sobol index?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper focuses on measuring variable importance in conditional average treatment effect (CATE) functions. Existing methods such as LOCO are applied to CATE estimators, and the authors propose the application of CPI as an alternative. A novel method, PermuCATE, based on CPI, is introduced to measure global feature importance in CATE estimators. Empirically, PermuCATE demonstrates superiority over LOCO in terms of faster convergence, lower variance, statistical power, and precision.

Overall, the reviewers' feedback was positive. The method was commended for being conceptually simple and inheriting the statistical inference guarantees on type - I error from previous results. Most of the major concerns raised by reviewers, such as enhancing theoretical justification and providing explanations for the connections with related work, have been addressed. 

After reviewing the paper, the reviewers' comments, and the authors' rebuttal, I recommend the acceptance of this paper at this time.
Decision: Accept (poster)