Testing Conditional Mean Independence Using Generative Neural Networks

We greatly appreciate your valuable comments, which have helped lead to a much-improved manuscript. In the following, we present our point-by-point responses to your questions and will take into account all your suggestions in a revised version of our manuscript.

Generality and practical implications of the assumptions. For a detailed discussion on Assumptions 7 and 9, please refer to our reply to Reviewer tmTe. In summary, our proposed CMI test is fully nonparametric, and Assumptions 7 and 9 do not impose explicit restrictions on the data distribution (e.g., boundedness, continuity, or sub-Gaussianity), enhancing its practical applicability. Furthermore, the double robustness property of the test statistic allows for mild assumptions on the error decay rates of nuisance parameters. For example, if we assume that $(Y, Z)$ follows a linear regression model, then $g_Y$ can be estimated at the $n^{-1/2}$ rate, which implies that the conditional distribution $P_{X|Z}$ only needs to be consistently estimated without strict rate requirements.

Following your comments, we will include a brief discussion in the revised version of our paper, highlighting the generality and practical implications of these assumptions, particularly the more technical ones.

Hyperparameter sensitivity. For the sensitivity to network hyperparameters, please refer to our reply to Reviewer tmTe. Regarding the choice of kernel bandwidths, we followed the median heuristic in our manuscript, as it is widely used in kernel-based tests. To further evaluate the sensitivity to bandwidth selection, we conducted additional simulations for Example A1 with a fixed sample size $n=400$ using bandwidths determined by either the mean pairwise distance or the “$\gamma$th quantile heuristic” for $\gamma$ $\in$ {25%, 75%}. Specifically, the bandwidth for $\mathcal{K}_X$ was set as the mean or $\gamma$th quantile of {$|X_j - X_k|_1 : j,k \in [n]$}, with similar choices for $\mathcal{K}_Z$. The empirical sizes for the test using the mean, 25% quantile, and 75% quantile bandwidths are 7%, 6.8%, and 7%, respectively. The size-adjusted power under the sparse alternative are 98.6%, 98.4%, and 98.6%, while under the dense alternative, it remains 100% in all cases. These results indicate that the test’s empirical performance remains stable across different bandwidth selection methods.

Computational complexity. Given the trained neural networks and assuming Gaussian or Laplacian kernel functions, our statistic resembles the average of two U-statistics of degree two, with its value depending on pairwise distances between samples. As a result, the computational complexity scales linearly with the dimensions $(d_X + d_Y + d_Z)$ and quadratically with the sample size $n$. For comparison, the computational complexity of pMIT$_M$ (Cai et al., 2024) and DSP$_M$ (Dai et al., 2022), both DNN-based CMI tests focused on univariate $Y$, scales linearly with $n$. The network training complexity is $O(E \cdot n \cdot P)$, where $E$ is the number of epochs and $P$ is the total number of trainable parameters.

In light of this comment, we will include a discussion on computational complexity in the revised version.

Alternative GNN structure and discussion on recent advances in generative modeling. At the early stages of this paper, we experimented with conditional GANs similar to those in Shi et al. (2021) to approximate $P_{X|Z}$. The empirical performance of the test was comparable to the current approach using GMMN. However, a key drawback of GANs is their longer training time, as both the generator and discriminator must be trained simultaneously. In contrast, GMMN has a more efficient training process and yields a test statistic whose empirical performance closely matches the oracle statistic; see Table 2 in Section 3. We will incorporate a more detailed discussion of recent advancements in generative modeling, particularly those highlighted in your review, in the revised version of our paper.

Robustness to model misspecification. A key strength of our proposed test is its fully nonparametric nature, as it does not assume a specific parametric form for the mean functions. Thanks to the universal approximation properties of neural networks, model misspecification is not a concern asymptotically if the network width increases with the sample size. However, in practice, network structures are fixed, which may introduce approximation errors.

Fortunately, the double robustness property of our statistic mitigates sensitivity to these approximation errors, making our method more reliable than approaches that lack this property. As demonstrated by the simulation results in Section 3 for both linear and nonlinear models, our test's performance is comparable to the oracle test.

We greatly appreciate your valuable comments, which have helped lead to a much-improved manuscript. In the following, we present our point-by-point responses to your questions and will take into account all your suggestions in a revised version of our manuscript.

Computational cost. Due to the sample splitting and cross-fitting framework, we need to train four neural networks (two GNNs and two DNNs) to construct our statistic, which is analogous to the two-split pMIT$_M$ test proposed in Cai et al. (2024). Regarding computation time, it takes 41.0 seconds for our method to complete one Monte Carlo simulation for Example A1 with a sample size of $n=400$. This is longer than competing methods but of the same order as pMIT$_M$; please refer to our reply to Reviewer KpXH for more details on computation time.
The computation time for training the neural networks depends on the complexity of the network, which is primarily determined by the data structures. For our numerical results in Section 3, the network structures used are relatively simple (with only one hidden layer). These simple structures are easy to train and still yield satisfactory empirical performance. Importantly, our method does not rely on a specific machine learning algorithm or network structure. As long as the estimation error meets the requirements in Assumption 7, any new or different machine learning techniques and network architectures can be used to reduce the training cost.

Rationality of the assumptions. Part (a) of Assumption 7 ensures that the (conditional) mean embeddings into the RKHSs, as well as the operator $\Sigma$, are well-defined. This assumption is commonly used in the literature of kernel-based conditional (mean) independence testing and holds for bounded kernels such as the Gaussian and Laplacian kernels.

Part (b) of Assumption 7 requires the estimation errors of the neural networks to decay to zero at rates $n^{-\alpha_1}$ and $n^{-\alpha_2}$ for $\alpha_1, \alpha_2 \in (0, \infty)$ such that $\alpha_1 + \alpha_2 > 1/2$. Similar rate requirements appear in Cai et al. (2024) and Lundborg et al. (2024), where "black-box" estimators such as DNNs are used. As shown in Stone (1982), the minimax nonparametric regression decay rate for $\mathbb{E} [ |g_Y(Z) - \hat g_Y(Z)|^2]$ is $n^{2p/(2p+d_Z)}$. When estimating $g_Y$ with DNNs, Bauer and Kohler (2019) demonstrated that the decay rate can be $n^{2p/(2p+d^\ast)}$, where $p$ is the smoothness parameter and $d^\ast$ represents the intrinsic dimensionality. Regarding the estimation error of the conditional mean embedding of $P(X|Z)$, our requirement is actually less restrictive than the assumptions on the total variation distance between $P(X|Z)$ and its estimator (see Remark 8 in Appendix C), which were used in Shi et al. (2021). Shi et al. (2021) argue that their assumption holds in a wide range of settings, with examples provided in Berrett et al. (2019). Moreover, due to the double robustness property of our test statistic, we do not impose explicit constraints on the individual estimation errors of $g_Y$ and the mean embedding of $P_{X|Z}$. Instead, we only require their product to decay faster than $n^{-1/2}$. Notably, when $g_Y$ is sufficiently smooth, $\alpha_2$ can approach $1/2$, allowing $\alpha_1$ to remain very small to accommodate complex and high-dimensional distributions of $P_{X|Z}$ (e.g., when both $X$ and $Z$ are images). For Assumption 9, we allow the residual vector $Y - \mathbb{E}[Y|Z]$ to vary under local alternatives. This is a more general setting than in nonparametric regression models, where the residual is assumed to remain the same as under the null hypothesis; see Remark 10 in Appendix C. As suggested, we will include a brief discussion of these assumptions, particularly the more technical ones, in the revised version of our paper.

Reference

Stone, C. J. (1982): "Optimal global rates of convergence for nonparametric regression." Ann. Statist.

Bauer and Kohler (2019): "On deep learning as a remedy for the curse of dimensionality in nonparametric regression." Ann. Statist.

Berrett et al. (2019): "The Conditional Permutation Test for Independence While Controlling for Confounders." Journal of the Royal Statistical Society Series B: Statistical Methodology

We greatly appreciate your valuable comments, which have significantly contributed to improving the quality of our manuscript. Below, we provide point-by-point responses to your questions and will incorporate all of your suggestions in the revised version of our manuscript.

Sensitivity to network parameters and computation time. To evaluate how changes in neural network parameters affect the performance of the proposed test, we repeated the simulation for Example A1 with a fixed sample size of $n = 400$ but varied the network configurations. In the first case, we increased the number of hidden layers in both networks to two, keeping all other parameters unchanged. In the second case, we reduced the number of nodes in $\mathbb{\widehat G}$ and $\widehat g_Y$ by half (to 56 and 128, respectively).

The results were consistent with those in Table 2 of the manuscript: empirical sizes: 5.8% and 6% (close to the original); size-adjusted power (sparse alternative): 98.6% and 99.8%; size-adjusted power (dense alternative): 100% in both cases. This suggests that the test’s performance is robust to moderate changes in key network parameters.

For other hyperparameters (e.g., batch size, learning rate), we used the Optuna automated search package to optimize them by minimizing the loss function defined in Equation (4) of Appendix A.

Regarding computation time, our method takes 41.0 seconds to complete one Monte Carlo simulation for Example A1 with sample size $n=400$ using NVIDIA T4 GPU, which is of the same order as the DNN based test pMIT$_M$. Detailed computation times for our competitors under the same setting are listed below.

• On Intel Xeon CPU. DSP: 3.21 seconds; DSP$_M$: 14.9 seconds; DNN-pMIT: 5.13 seconds; DNN-pMIT$_e$: 5.13 seconds; DNN-pMIT$_M$: 25.0 seconds; DNN-pMIT_e$$_M: 25.0 seconds.
• On 11th gen intel core i7-11800h CPU. XGB-pMIT: 0.167 seconds; XGB-pMIT$_e$: 0.167 seconds; XGB-pMIT$_M$: 1.56 seconds; XGB-pMIT_e$$_M: 1.56 seconds; PCM: 0.20 seconds; PCM$_M$: 1.07 seconds; VIM: 6.667 seconds.

Competitor methods on image data application. We applied some competing methods to the image datasets in our initial submission, but due to space limitations, we have included some details in the appendix. For the image data application in Section 4.1, we compared our method with the DSP$_M$ approach developed by Dai et al. (2022), where a similar application was examined. For the application in Section 4.2, we compared with the pMIT$_M$ method introduced by Cai et al. (2024). The p-values for these comparison methods are included in Figures 1 and 3, and a detailed comparison is provided in Appendix B.

Sensitivity to numbers of Monte Carlo data and bootstrap sample. We have conducted additional simulations to investigate how the number of Monte Carlo synthetic data ($M$) and the bootstrap number ($B$) influence the performance of our test. The results suggest that both the size and power of our testing procedure remain robust against these tuning parameters. Specifically, we repeated the simulation for Example A1 with $n=400$, varying $M$ in {$5,20,60$} and $B$ in {$200,500$}. The results are presented in the following table.

Numbering and location of assumptions. Originally, we chose to place the assumptions in the Appendix to stay within page limits and to keep the focus on explaining the core ideas of our proposed test without introducing excessive technical details. For the next version, we plan to maintain them in the Appendix due to space constraints, but will add a few sentences in the main text to briefly summarize these assumptions. Regarding assumption labeling, we initially used the default formatting from the ICML LaTeX template, but we would be happy to relabel them separately from theorems and remarks if preferred.