We would like to thank the reviewer. The comments and questions are relevant and will help us improve the paper.

Discussion of assumption 3.1

• 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.

• 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.

• 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.

Description of the total Sobol index

• 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],$ where $\mu_{-j}$ is the model where the variable $j$ is removed.

• 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]:

$\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] .$

• 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).

• 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.

We will integrate this intuitive presentation of the total Sobol index in the final version of the paper.

• [1] Sobol, IM. "Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates." Mathematics and computers in simulation, 2001
• [2] Bénard, C. et al. Mean decrease accuracy for random forests: inconsistency, and a practical solution via the Sobol-MDA, Biometrika, 2022
• [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
• [4] Williamson, B. et al. A General Framework for Inference on Algorithm-Agnostic Variable Importance. Journal of the American Statistical Association, 2022
• [5] Williamson, B et al. Nonparametric variable importance assessment using machine learning techniques. Biometrics, 2021

We would like to thank the reviewer. The comments and questions are relevant and will help us improve the paper.

Previously published references

• 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.

Connection between CATE and variable importance

• 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$.

• 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]$.

• 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.

Motivation and problem setting

• 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.

We 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.

We thank the reviewer for the thorough review. The comments and questions are relevant and will help us improve the paper.

Dependence of PermuCATE on the estimate of $\tau$

• 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$.
• 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.
• 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.

To clarify these last two points, we propose to rephrase the proposition using the intermediate step (eq 9) of the proof in A.4.

Proposition 3.4:

Under Assumption 3.1, for a consistent CATE estimator $\hat{\tau}$, the finite sample estimation of the importance for the $j^{\text{th}}$ covariate is: $\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})$

Remark: While the first term involves the CATE estimate $\hat{\tau}$, it does not involve the estimation error.

Corollary:

Under 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})$

Variance comparison between LOCO and PermuCATE

• 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.

Subset used for conditional perturbations

• 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.

Benchmarks that do and do not satisfy Assumption 3.1

• The simulation scenarios LD, HL, and HP (see datasets in the main paper) used multivariate Gaussians, which satisfied this assumption.
• 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.
• 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.

Factor 2 in risk difference

• As demonstrated in Proposition 3.2, PermuCATE estimates the total Sobol index up to a factor of 2.

We thank the reviewer for the thorough review. The comments and questions are relevant and will help us improve the paper.

Dependence of PermuCATE on the estimate of $\tau$

• 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.
• 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.

Complexity of the conditional distribution estimation

• 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.
• 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.

[1] Hollmann et al. Accurate predictions on small data with a tabular foundation model. Nature 2025

Description of the total Sobol index

• 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.

Discussion of model-x knockoff

The 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.

Besides the work from Hines et al. 2020, we are unaware of other methods that provide model agnostic variable importance for CATE.

Extension to groups of variables

• 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.

Ranking of important variables

• 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.