Testing Conditional Mean Independence Using Generative Neural Networks
Authors: Yi Zhang
Linjun Huang
Yun Yang
Xiaofeng Shao
Venue: ICML 2025 poster (2025)
Keywords: Conditional Distribution
Maximum Mean Discrepancy
Kernel Method
Double Robustness
Primary area: general_machine_learning
TL;DR: We introduce a novel method for conditional mean independence testing using deep generative neural networks with strong performance in high-dimensional setting.

Abstract: Conditional mean independence (CMI) testing is crucial for statistical tasks including model determination and variable importance evaluation. In this work, we introduce a novel population CMI measure and a bootstrap-based testing procedure that utilizes deep generative neural networks to estimate the conditional mean functions involved in the population measure. The test statistic is thoughtfully constructed to ensure that even slowly decaying nonparametric estimation errors do not affect the asymptotic accuracy of the test. Our approach demonstrates strong empirical performance in scenarios with high-dimensional covariates and response variable, can handle multivariate responses, and maintains nontrivial power against local alternatives outside an $n^{-1/2}$ neighborhood of the null hypothesis. We also use numerical simulations and real-world imaging data applications to highlight the efficacy and versatility of our testing procedure.

Testing Conditional Mean Independence Using Generative Neural Networks
Yi Zhang *1
Linjun Huang"
Yun Yang² Xiaofeng Shao
Abstract
variables. In machine learning, CMI testing has broad ap- plications in areas like interpretable machine learning (Mur- doch et al., 2019), representation learning (Bengio et al., 2013; Huang et al., 2024) and transfer learning (Maqsood et al., 2019; Zhuang et al., 2020).
Conditional mean independence (CMI) testing is crucial for statistical tasks including model de- termination and variable importance evaluation. In this work, we introduce a novel population CMI measure and a bootstrap-based testing proce- dure that utilizes deep generative neural networks to estimate the conditional mean functions in- volved in the population measure. The test statis- tic is thoughtfully constructed to ensure that even slowly decaying nonparametric estimation errors do not affect the asymptotic accuracy of the test. Our approach demonstrates strong empirical per- formance in scenarios with high- high-dimensional co- variates and response variable, can handle multi- variate responses, and maintains nontrivial power against local alternatives outside an n 112 neigh- borhood of the null hypothesis. We also use nu- merical simulations and real-world imaging data applications to highlight the efficacy and versatil- ity of our testing procedure.
In this paper, we address the problem of CMI testing for multivariate response variables and covariates. Specifically, for random vectors X RRdx, Y Rdy, and Z E Rd2, we test the null hypothesis Ho that Y is conditionally mean independent of X given Z, i.e., B[Y]X x, Z =| E[Y 2 = 2]. a.c. (I, E Rdx+d2 For example, consider predicting age Y from a facial image. To test whether Y can be predicted using images with a specific facial region (potentially containing sensitive information) covered in black, X can represent the covered region and Z the re- maining facial image. Similarly, to test whether Y can be predicted using a low-resolution version or extracted fea- tures of a facial image, X can represent the original image and Z a low-dimensional feature vector derived from X using standard extraction methods (e.g.. Autoencoder, Av- erage pooling, or PCA (Berahmand et al., 2024)). Under Hg, removing X from the predictive model f(X, Z), where f : Rd - dy denotes the nonparametric regression function, does not reduce prediction accuracy.
1. Introduction
For machine learning tasks, identifying which (functions of) covariates contribute to predicting the response is par- ticularly important, as deep neural networks (DNNs) often process high-dimensional data, such as images and text, that may include irrelevant features (Li et al., 2017; Van Lan- deghem et al., 2010). Performing dimensionality reduction before DNN training provides several advantages: it en- hances interpretability. improves prediction accuracy, and reduces computational costs (Cai et al., 2018). In evaluating DNN training, covariate significance is typically assessed by comparing performance metrics, such as mean squared error or classification accuracy, for DNNs trained with and without specific covariates (Dai et al., 2022; Williamson et al., 2023). As discussed in Section 2.3 of Williamson et al. (2023), many common performance metrics, such as the coefficient of determination (R²) for regression and empirical risk with cross-entropy or 0-1 loss for binary clas- sification, are functionals of the conditional mean functions of Y. which means Ho will imply that the covariates are not significant. Thus, CMI testing provides valuable tool for
Conditional mean independence (CMI) testing is a funda- mental tool for model simplification and assessing variable importance, and plays a crucial role in statistics and ma- chine learning. In traditional statistical applications, such as nonparametric regression, CMI testing identifies sub- sets or functions of covariates that meaningfully predict the response variable. This is essential for improving model effi- ciency, accuracy, and interpretability by avoiding redundant
'Equal contribution 'Department of Statistics, Univer- sity of Illinois at Vihana-Champpagee Champaign, II, USA. 'Department of Mathematics, University of Maryland, Col- lege Park, MD, USA. 'Department of Statistics and Data Science, and Department of Economics, Washington Uni- versity in St Louis. St Louis, MO, USA. Correspondence to: Yi Zhang syizl9@ilinooi..ee junh2@ illinoiis.ccooco edu>, Yun Yang <yy84@ edu, Xiaofeng Shao < shaox edu>.
Linjun Huang <lin-
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Testing Conditional Mean Independence using GNN
assessing covariate importance across variety of machine learning applications.
under Ho. [(F, m) converges to zero at a rate faster than the n 1/2 rate at which ['r; m) - I converges to a non- degenerate limiting distribution under the alternative (Fan & Li. 1996. Section 1).
1.1. Related literature
To the best of our knowledge, most existing CMI tests focus on univariate Y and face one more of the following three performance deteriorates when some or all of (X, Y, Z) are high-dimensionn tests lack theoretical size guarantees in general nonparamet- ric settings; and (3) they cxhibit weak power in detecting local alternatives. For a recent survey, see Lundborg (2022). We discuss how these challenges arise and how existing CMI tests have partially addressed them.
Second, the nonparametric estimation errors for r(2) and mn(r, x) typically decay slower than the n n=1/2 parametric rate, and the convergence rate of ( m) under Ho depends heavily on how quickly these errors decay. For CMI tests that use kernel smoothing to estimate the conditional mean functions, the estimation error has an explicit form, allow- ing it to be addressed directly, along with the degeneracy issue, when deriving asymptotic results under Ho- As a result, all the aforementioned kemel smoothing-based CMI tests have theoretical size guarantees. In contrast, due the black-box nature of DNNs, the estimation errors for r(z) and m(x, =) cannot be explicitly decomposed or handled in the same way as kernel smoothing estimators. Conse- quently, addressing the degeneracy issue requires additional debiasing procedures to mitigate the impact of these errors and achieve accurate size control. For example, Williamson et al. (2021) and Lundborg et al. (2024) constructed statis- tics based alternative forms of I' to reduce bias. However, the degeneracy issue persists: as shown in Theorem I of Williamson et al. (2021), their estimator is Vn-consistent when r > 0, but no asymptotic results are provided un- der Ho (i.e., 1 = 0). Similarly, the test in Lundborg et al. (2024), which builds on the conditional independence test from Shah & Peters (2020), relies on strong assumptions to bypass the degeneracy issue, outlined in Section 2.1 and Part (a) of Theorem 4 in Lundborg et al. (2024). To address the degeneracy issue, Dai et al. (2022) introduced additional noise of order Oy (nnn21 to the estimator of r, while Williamson et al. (2023) utilized sample splitting to estimate separate components of an equivalent form of T on different subsamples. However, both approaches are ad hoc (Section 6.2 of Verdinelli & Wasserman (2024)) and, as discussed in Appendix S3 of Lundborg et al. (2024), suffer from significant power loss under the alternative hypothesis.
major issues: (I)
(2) the
Performance deterioration in high dimensional setting, This issue primarily arises from the estimation of the condi- tional mean functions r(2) == | / = =] and m m(x z) = E|Y|X = x, Z = =]. Early CMI tests, such as those in Fan & Li (1996); Delgado & Manteiga (2001); Zhu & Zhu (2018); Lavergne& & Vuong (2000); Aut-Sahalia al. (2001), relied on kernel smoothing methods for these estimations. Consequently, these CMI tests suffer from the curse of di- mensionality: their performance declines significantly as the dimensions dz, dx + dz, and/or dy increase (Zhou et al., 2022, Section 1). For instance, Figure 1 in Zhu & Zhu (2018) shows that the empirical power of tests proposed by Fan & Li (1996) and Delgado & Manteiga (2001) decreases sharply with increasing dz and dx, exhibiting trivial power under sample size of 200 with moderate dimension- ality (dz- = 4 and dx 2 8). To address this, recent CMI tests utilize machine learning tools, such as deep neural net- works (DNNs) and the kernel trick, to estimate conditional mean functions. These tools are effective in approximating complex, digeeniimesionnl functions with underlying low- dimensional structures. For example, tests in Williamson et al. (2021); Dai al. (2022); Lundborg et al. (2024); Cai et al. (2022); Williamson et al. (2023); Cai et al. (2024) leverage DNNs for conditional mean estimation, achieving better performance in high-dimensional settings.
Weak power against local alternatives. On one hand, the CMI tests proposed in Fan & Li (1996); Zhu & Zhu (2018); Lavergne Vuong (2000); Ast-Sahalia et al. (2001); Williamson et al. (2021); Dai et al. (2022); Williamson et al. (2023) fail to detect local alternatives with signal strength An v [[[r(Z) m(X, 2) of order 1/2 primarily due to their reliance on the population measure r. Specif- ically, the tests in Fan & Li (1996); Zhu & Zhu (2018); Lavergne & Vuong (2000); Aat-Sahalia et al. (2001) take the form nh./2f, where -> a bandwidth parame- ter used in kernel smoothing. and = dz or dx + dz. Consequently, these tests cannot detect local alternatives converging to the null faster than n -1/2 in Williamson et al. (2021); Dai et al. (2022); Williamson et al. (2023) use the population CMI measure To = r -T2,
Theoretical size guarantee. The main challenge in es- tablishing a theoretical size guarantee stems from two key issues. First, most of the existing CMI tests (ex- cept the one proposed by Delgado & Manteiga (2001)) rely on the sample estimation of the population CMI mea- sure I := [[r(()) - m(X, w 2)] or its equiv- alent forms, where w is a positive weight function. A common plug-in estimator of is given by T(P,m) = mX, Z.)) 2 Zi), where f and m are estimators of the conditional mean functions. For population CMI measure to be valid, it must uniquely characterize CMI, meaning that the measure equals zero if and only if Ho holds. While I' satisfies this requirement, its estimator [(F m) suffers from a degeneracy problem:
n 1
The tests
2

---

Testing Conditional Mean Independence using GNN
train a generative neural network (GNN) to sample from the (approximated) conditional distribution of X given Z. The CME is then estimated using the Monte Carlo method with samples generated from the trained GNN.
Table Summary of existing CMI tests. High-Dim: whether the test has good empirical performance when some or all the dimensions of X, Y and Z are large; Size Guarantee: whether the test has theoretical results under Hu that guarantec accurate asymptotic size control in general nonparametric setting; Local Alt: whether the test can detect local alternatives with signal strength A. converging to zero at the parametric rate n 11/2.
1.3. Organization and notations
The paper is organized as follows: proposed population CMI measure, the test statistic, and the bootstrap calibration procedure, along with its asymptotic properties and consistency results. Section 3 evaluates the test using finite sample simulations, comparing it with other methods. Section 4 presents two real data examples. Section 5 concludes with final remarks. All proofs and additional details are deferred to the Appendix. A Python implementation of our proposed test procedure is available at
Section 2 introduces the
Tests (1996)
High-Dim
Size Guarantee
Local Alt
No
Yes
No
Delgudo & Manteiga (2001)
No
Yes
Yes
Zhu & Zhu (2018)
No
Yes
No
Lavergne & Vuong (2000)
Yes
No
Sat-Suhalia et al. (2001)
No
Yes
No
Williamson ot al. (2021)
Yes
No
No
al. (2022)
Yes
No
Lundborg et (2024)
Yes
No
Yos
Williamson et al. (2023)
Yes
Yes
No
Cai et al. (2024)
Yes
Yes
No
wtt/://innnnnnnnninnoonnnnnnnnna 6 6749/
Our method
Yes
Yes
Yes
The following notations will be used throughout the pa- per. For any positive integer d and random vectors (x' "x?.....,X"" X²,, Z) defined on the same probability space, Pxixa...x and IP x x212 denote the joint distribution of X1 xx, ......, Xd and its conditional distribution given Z, respectively. Let Ez represent the expectation with re- spect to Izz. and let I'm denote the Lebesgue measure on Rd. For positive integer n, define [n]] {1, 2, ................................ .n.. For probability measure sa(-) on Rd and p 2 1, let (IR( = {f Rd + IR ʃ (u(()) < oo For any Hilbert spaces A and B, let A @ B denote their tensor product, and use (.. and || lla to denote the inner product and induced norm on A, respectively. For random vectors a, b E Rd, the Gaussian kernel is defined as K[aa)) exp( l|a bll3/(2o²) and the Laplacian kernel as b) = exp( - ||a - /ll1/o), where a >> 0 is the band- width parameter, || ||2 is the Euclidean fg-norm, and ||- llx is the nnorm.
2))%, which is equivalent to r. Since the quadratic terms l11 and 22 can only be estimated at the 112 rate, these tests are limited detecting local alternatives with An of order -1/4, On the other hand, the test proposed in Cai et al. (2024) employs sample splitting and requires the size of the training subsample used to estimate the conditional mean functions to be substantially larger than that of the testing subsample. As a result, their test is limited to detecting local alternatives with An converging to zero slower than n 11/2, which can still result in significant power loss in practice.
where 11 = B( and =
1.2. Our contributions
Table summarizes limitations of existing CMI tests. To overcome these challenges, we propose a novel CMI testing procedure with the following advantages:
1. The test demonstrates strong empirical performance even when some or all dimensions of X, Y, and Z are high. Notably, it is well-suited for scenarios where imaging data serve as covariates, responses, or both.
2. Conditional Mean Independence Testing
In this section, we introduce a novel population measure for evaluating CMI and propose corresponding sample-based statistic for conducting the CMI test. We then establish theoretical guarantees for the proposed procedure, including size control and power against local alternatives.
2. The test achieves precise asymptotic size control under Ho-
3. The test exhibits nontrivial power against local alterna- tives outside an n 1 htt.///niihhboroood of Ho.
2.1. Population conditional mean independence measure Recall that the goal is to test the null hypothesis of Ho E|Y X, Z] = B a.s.-Pxz xx against the alternative hy- pothesis of H₁ PEEY X, Z] # 2 > 0. To moti- vate our population measure for evaluating CMI, we begin with the following result, which provides equivalent charac- terizations of CMI (Daudin, 1980).
To achieve these features, we propose a new population CMI measure closely related to the conditional independence measure introduced in Daudin (1980). Additionally, we develop sample version of this population measure in a multiplicative form, which is key to mitigating the impact of estimation errors in nonparametric nuisance parameters (i.e., the conditional mean functions). Our test not only requires estimating r(z) but also the conditional mean embedding (CME, Song et al., 2009, Definition 3) of X given Z into a reproducing kernel Hilbert space (RKHS) on the space of X. Instead of directly estimating the CME using DNNs, we
Proposition 1. If ( < co, then the following prop- erties are equivalent to each other:
(a) E|Y 2) = || axx-
(b) [[[[(XX)) 7) F[S((, Z) V] = 0 for any f E
3

---

Testing Conditional Mean Independence using GNN
Ly Lg((Rx Pxz))
Due to the multiplicative form inside the expectation in equation (1), when constructing a sample version of T as the test statistic, the estimation errors of the two nuisance parameters (gx and gr) do not affect the asymptotic prop- erties of the statistic, as long as the product of these errors decays faster than n -1/2; see Lemma 17 in Appendix E. This desirable property is commonly referred to as double robustness (Bang & Robins, 2005; Zhang et al., 2024).
(c) 2) BI 2I-1IXX any I E L2( +22 ,xz).
(Y-E)) = 0 for
Remark 2. It is straightforward to see that (a) + (b) (c), and (c) implies (a) by laking f(X,7) B[yTe| c| X 2] over all e E Rdv The only difference between (c) and (b) in Proposition 1 is that Y is centered at B| Z in (c). This additional centering is crucial for reducing biases from the estimation of the conditional mean functions; our proposed population CMI measure will be derived from (c) by consid- ering all f within a dense subset of L2 (xxRRxxd Pxz)).
Our I bears similarity to the maximum mean discrepancy- based conditional independence measure (MMDCI) pro- posed in Zhang et al. (2024), which was designed for testing the stronger null hypothesis of conditional independence (CI). While CI assumes that the conditional distribution of the response variable Y is independent of an additional covariate X given an existing covariate Z, the CMI assump- tion only requires that the conditional mean function of Y does not depend on X given Z. This means that includ- ing X in the regression function of Y does not improve its predictive ability in the mean squared error sense. In con- trast, rejecting CI does not provide a direct interpretation in terms of predictive ability (Lundborg et al., 2024, Scction 1). Furthermore, the MMDCI statistic compares the joint distribution of (X, Y, Z) under the null with its true distri- bution using maximum mean discrepancy (MMD) (Gretton et al., 201 2), requiring stronger assumptions to fully char- acterize CI (see Assumption 1 in Zhang et al. (2024)). In contrast, our population CMI measure avoids using MMD to capture the full distributional properties of (X, Y, Z). mak- ing it more suitable for CMI testing compared to existing population CMI measures. Computationally, the test pro- posed in Zhang et al. (2024) requires training two generative neural networks (GNNs) to sample from the conditional dis- tributions of X and Y given Z, respectively. In contrast, our proposed test statistic only requires sampling from the conditional distribution of X given Z, resulting in lower computational cost and reduced memory usage.
Let Kx : Rdx x Rdxx and Kzz : Rdz x Rd denote two symmetric positive-definite kernel functions that define two reproducing kernel Hilbert spaces (RKHS) Hy and Hiz over the spaces of X and Z, respectively. Additionally, let Ku- - =Kx x X2 represent the product kernel, with H&g being the corresponding RKHS induced by Ko over the product space RRdx Rdz, Motivated by part (c) of Proposition 1, HHlo
we define linear operator D R
>== (<o 2), - ([Ko 2), )12]
Ec
[X - -[ c). anr c E R
From the reproducing property, we see that for any ʃ E Ho and any c C Rdy,
().Ec)n = = - 2) - E|y EY "d)
Under the assumption that the RKHS Ho is dense in Ta (Rdx+dz, Pxz), which holds if Ixx and Kx are L-- or co-universal kernels (Szabó& & Sriperumbudur, 2018, The- orem 5), such as the Gaussian and Laplacian kernels consid- ered in this paper, the preceding display along with Propo- sition implies that the null hypothesis Ho holds if and only if Σ is the zero operator (i.e., Ee = 0 E IHl for any e E E ). The following proposition formalizes this intu- ition and serves as the foundation for our proposed popula- tion CMI measure,
2.2. Testing procedure
To construct a sample version of r it is necessary to estimate gy and gx. To address the curse of dimension- ality, we use deep neural networks (DNNs) to estimate gy. For gx, which is an RKHS- function, we esti- mate gx gx(2) for fixed E RR by sampling M copies
m x")
2)),
(I)
| 1 (z)) and where V(KYY' = (r)) T Y
= -x((X, =xxxx))=xx(xxxxx)
(xX', ,)ux - -(xx KX, 11x + (ex (Z), gx max
(Xa, J.M. from the conditional distribution of X given Z = z. The estimate of (2) is then given by the sample average M -- (x the noise-outsourcing lemma (Theorem 6.10 of Kallenberg (2002); see also Lemma 2.1 of Zhou et al. (2022)), for any integer m 2 1, there exists a measurable function G such that, for any 17 ~ Im,, independent of Z, we have G(1, 2) Z ~ Pxiz. Based on this result, we train a generative moment matching network (GMMN, Dziugaite et al., 2015; Li al., 2015) to approximate G. For any fixed 2 E dd we generate i.i.d. copies na from N (0,
Here, (X', (','' is an independent copy of (X, (X,Y, 2), 2). gy and gx are defined as gy ()-B(y]-= E dy and (x (X -== E Hx. respectively. Proposition 3. If Assumption 4(a) holds, then (see Defi- nition 11 in Appendix C the definitions of the Hilbert- Schmidt operator and norm)
M
Hxx Specifically, using
(a) Σ is a Silbert-Schmidt operator, and its Hilbert- Hiibert-Schmidt norm, denoted as [[E][us, satisfies |Eiliss = -
(b) The null Ho holds if and only if = O.
4

---

Testing Conditional Mean Independence using GNN
input them into the trained GMMN & along with 2, and treat the outputs G()::)) =) as approximately sampled from the conditional distribution of X given Z = We refer to the trained GMMN as the (conditional) generator, and de- tailed training procedures are provided in Appendix A. Note that other generative models, such as generative adversarial networks (GAN) (Goodfellow et al., 2014) and diffusion models (Yang et al.. 2023), can also be applied to approx- imate G, as long the approximation errors satisfy mild decaying conditions. The field of generative modeling has recently seen significant advances, with notable examples including Athey et al. (2024); Baptista et al. (2020); Shi et al. (2022); Nguyen et al. (2024). Incorporating these newly de- veloped models into our testing procedure is straightforward and may potentially improve both the computational cost and the performance of the proposed test.
Pxyz. Consequently, we cannot directly determine the rejection threshold of nT,, without knowing PXY z. Instead, we employ a wild bootstrap method to approximate the distribution of nT,,. Following Section 2.4 of Zhang et al. (2018), for a positive integer B and b E [B]. we generate {<m))== from a Rademacher distribution, and define the
bootstrap version of Tr as
M 2
1
Σ
Ü(X,, (Y), Y2)
-1)
Jtk J.KE(O
. KZz Zz )esse }
(2)
Since ( Ma.... M can be viewed as samples from the distribution of Tr (c.f. Theorem 9), we reject Hg level Y E (0, 1) if B Σ>= 1 [[nT)>nf. >nf.) <7- As a default choice, we set B 1, 000 throughout our numerical experiments.
In order to eliminate dependence on conditional gener- ator estimation and improve test size accuracy, we fol- low Shi et al. (2021) and adopt a sample-splittng and cross-fitting framework to train the GMMN and DNN. For easy presentation, we consider two-fold splitting in this paper, although the test can be readily generalized to the multiple-fold splitting setting. Speciffically, for posi- tive integers M and m E [M]. let {(x,, Y, }:))== and {n," in Ti=1 be i.d. copies of (X, Y, 2) and such that ([(X,, Yú, 2.) (n. 1a li (n(" T 12 } are mutually
2.3. Theoretical properties
In this part, we evaluate the asymptotic performance of the test as the sample size increases, focusing on its empiri- cal size (Type-I error) control and power against (local) alternatives. Specifically, we aim to determine whether the proposed test satisfies two desirable theoretical properties: 1. the probability of incorrectly rejecting the null hypothesis converges to the specified significance level yy as the sample size grows; 2. the test maintains the capability to detect (local) alternative hypotheses with a deviation from null that diminishes at the parametric rate of n 11/2.
independent. Let X X((m) = G(7,"" ; 21) denote the data sam- pled from Px.zz. We divide [n] into two equal folds ((1) and (2) where J(1)- = J(-2) = {1,2,.... [n/2]] and = = [m// For E [2]. we train a GMMN generator G, using data (X, Z,) for i E J(-)), Similarly, we train DNN Øj using data (Y,, 2:) for i E J(-5) as an estimator of gy (four neural networks are trained in total).
From our procedure, we define the event Rej,, ==
[BI 1 (n.7. (n..):.ff.) { as rejecting Ho. To demon- strate that our proposed test has an accurate asymptotic size, it suffices to prove that the probability of rejecting Ho. given that Ho is true, converges to the specified level 7, i.c., lim, P(Rej, Ho) = 7. To achieve this, we first derive the limiting null distribution of nef,, in Theorem 5. Subsequently, we show the consistency of our bootstrap procedure in Thcorem 9, that is, conditional on the sample, the rescaled bootstrap statistic converges to the same limit- ing null distribution. To study the asymptotic power of our proposed test, define a local alternative hypothesis in equation (4) whose deviation from null scales as n n- for some a 2 0. The asymptotic behavior of T,, and the boot- strap counterpart under different values of a is analyzed in Theorem 8 and Theorem 9, respectively. Based on these re- sults, we derive the asymptotic power of the test in Theorem 9. Let T, denote the oracle test statistic, defined similarly to Th except that the estimated nuisance parameters (gy:9x) are replaced with their true values; see equation (14) in Appendix D for the precise definition.
Let X (m))) = Gn,", Z,) and (Z.) = (Z,) if e J(), then we define the test statistic as
T.
Σ
(XX,,
K2) 2.)
(.
where
Xa) 1
Σ Xx(XX,, 8(m))
M
M M Σ
M
8(m)) ) Kx (Xa,
(x( {(()))
M
and (Y Given the trained neural networks and assuming Gaussian or Laplacian kernel functions are used, T, resembles the average of two U-statistics of degree two, with its value depending on pairwise distances between samples. As result, the computational complexity scales linearly with the dimensions (dx + dy + d2) and quadratically with the sample size n. In comparison, the computational complexity of the statistics proposed in Cai et al. (2024) and Dai et al. (2022), both designed for univariate Y, scales linearly with n. As will be shown in Theorem 5 below, the limiting null distribution of nf,, depends on the unknown distribution
[YL-ōv (Zw)
The following assumption is needed to derive the asymptotic properties of our statistic.
Assumption 4. Assume M + Do as n2 + 00. For
Co > O. 01, 02 E (0, })) such that an a1 > 2. D. e
5

---

Testing Conditional Mean Independence using GNN
{x,, X,", (,), Fg E gr (Z.).0v (Z,) and E (2) where is E g(s) for each s E [2], we have:
Assumption 4 by the total variation distance and the mean squared error, i.e., B B[d (Px x Px,1z.) O(n 201) and B[lloy (2+) gy (Z.) [|]
O(n-202), which com- mon assumptions made in existing works for characterizing qualities of conditional generators and nonparametric regres- sion functions. However, the total variation metric may not be a suitable metric for characterizing the closeness between nearly mutually singular distributions, which happens when data are complex objects such as images or texts exhibiting low-dimensional manifold structures (Tang & Yang, 2023). Moreover, due to the double robustness property of our test statistic, we do not impose separate constraints on the respective estimation errors of gy and the mean embedding of Px|z. Instead, we only require that their product decays faster than 11/2, Notably, when gy is sufficiently smooth, Og can approach 1 1/2, allowing to remain very small to accommodate complex and high-dimensional distributions of Px|z (c.g., when both X and Z are images).
(a) E { (<x(DD; D. D.II 2,) <Co:
(b) F (11E))x((...))))2=EX 8 I (K2(Z2, Z,)
+ +1E.1BK2(Z, 2,)1 -001
=0(n ::);
[ (1-00(22-00022)) (4) Kz (.)7, =00 =0(=
(c)
It it worth noting that the required technical assumptions are relatively mild and are commonly used in the literature of CMI testing. Part (a) of Assumption 4 ensures that the (conditional) mean embeddings into the RKHSs, as well as the operator Σ, are well-defined and holds for bounded kernels such as the Gaussian and Laplacian kernels. Part (b) of Assumption requires the estimation errors of the neural networks to decay to zero at rates -0- and n n for 02,02 E (0, oo) such that o1 + 02 > 1/2. Similar rate requirements appear in Cai et al. (2024) and Lundborg et al. (2024).
Theorem 5. SSuppose Assumption 4 holds, then under Ho. 02
we have that T,, -T, = (n- 1/2 +n -a +n + 1=1/2-(a, +a2)) as -> bo; moreover, nf, B. Td. where {T(, T. } are i.i.d. random variables following a mixture of centered chi-squares distribution (see Appendix D.I).
When estimating gY with DNNs, Bauer& & Kohler (2019) demonstrated that the decay rate for TE [[ggr (Z) - Ty [ can be bounded from above by n 22=/(2+''"),, ) where 8 is the smoothness parameter of gy, and d' denotes its intrinsic dimensionality. Regarding the estimation error of the con- ditional mean embedding of Pxiz z. our requirement is less restrictive than the assumptions based on the total varia- tion distance between Pxiz and its estimator, which are commonly adopted in the literature.
From this theorem, we observe that the "plugged-in" test statistic T, becomes asymptotically equivalent to the oracle statistic T, as long as a1 + 02 1/2. This demonstrates the property of double robustness: the slowly decaying non- parametric estimation errors in (gy, (gy,9x) ) do not compromise the asymptotic accuracy of the test, provided the product of these errors decays faster than n -1/2
Now let us switch to the asymptotic power of the test under local alternatives. We use the triple (X", y z²) de- note (X, Y, Z) under Ho and consider sequence of triples (X", VVA. 7 under the alternative hypothesis:
Specifically, if Kx. K2 are bounded kernels and ||УII is bounded (without loss of generality, assume they are bounded by 1), then Assumption 4 reduces to [[[[E[K(x( X,) Z,] - F[Kx((, x(2) Zi] Hx
O(n 221) and E[llgy (.)) -x (Z,) [ = O(n 202)). Note that
Hin: YY, =
+ " 2 + Rn.
(4)
Here, the -valued function g and the random vector RS satisfy E[G(X", | = 0 and | 2"| O, respectively. so that the exponent 2 0 determines the decay rate of the deviation from null under Hi.; for instance, setting a = 0 and R., = R corresponds to a fixed alternative. The following assumption is needed to derive the asymptotic results under Hnn when a > O.
E [[Kx(
X.) Z1 - E[ 271 3 ]
B[f(X(")) - f(x)) |
E
sup
SEHx: |15 nx $1
E [(((8(())) - f(x.) z }
<E
sup
f:R²x -R: lf $1
Assumption 6. For D DE (Kx(X") -). 9..x 9x(2))) there ex-
-2 di (P 12. Px)
(3) ists random variable ( E Rdd such that B |IDI R|
<lEKZ 2n) (2)) }
where || . ala denotes the function supreme norm, and drv(..) denotes the total variation distance.
Remark 7. Assumption 6 iimples that when > 0, the demeaned random vector Rn = Y,A - EE |Xx" 1 converges to some random vector Ç. Instead of fixing 3P-EEY' - |Z| as in nonparametric regression mod- els (see equation (1.I) in Zhu & Thu (2018)), we allow R...
Here, the inequality in the second line is implied by the fact that ll = sup [f(x)]
supzeR*x (x, < If vxx(2,2) V
R2 = = Y
|xx- Therefore, we can also replace the error metries in
6

---

Testing Conditional Mean Independence using GNN
to change with and ( can be different from Y - Theorem 8. Suppose Assumptions 4 holds, then under Hin (see Appendix D.2 for precise definitions of relevant con- stants and random variables),
2. Under Hin with a = O, n.MaN L T - E T,
=
where {Ti, T2) are i.i.d. random variables following a mixture of centered chi-squares distribution.
With Assumption 6 further satisfied, we have
1. If a = O, then vā - co) D 1 1E=19
900 as n + 0o for some positive constant CO and i..dd mean
V2
3. Under Hin with a > 0, n.Ma Furthermore, if we let Mass.y denote the (1 71 quantile of ffM conditioning on the data, then the power (probabil- ity of detecting the alternative) of our testing procedure sat- isfies: ifa/1/2 then P(nIn 2 -inMa)) - 1; if a=fff114 1/2.
MM D TT-
zero normal random variables ( -2 ):
With Assumption 6 further satisfied, we have
2. 1f0-a<<1/2 then Tn 2n cfoor some c > 0;
3. Ifa = 1/22 then nf, Ds C+TA+
v2 Ej=1 J3. where = Ei Tıj and (AA::TA! are i.i.d. random variables following a mixture of centered chi chi-squares distribution.
then P(nIn 2. >MAMa) + +
N
To where To., denotes the (1 - 7)th quantile of T fa> > 1/22 then P(nTT 2 MMM,)) -
7-
Since P(nTT, 2 MMM..) = P(Rej,,) item I of Theorem 9 demonstrates that the proposed test achieves asymptotically correct size; while items 2 and 3, along with the second part of the theorem, establish that the proposed test is consistent (the probability of rejecting Hg converges to 1) when the alternative lies outside an n -1/2 neighbioiooo of Ho (i.e.. when a < 1/2). Similar to the optimal test based on nT,, discussed after Theorem 8, the proposed test will also ex- hibit trivial power in detecting alternatives that too close to the null (i.e., when a > 1/2).
1/2, 4. 1fa> 1/2, then nf, T
Theorem 5 establishes that rejecting Ho if nf, exceeds the rejection threshold as the (1 -- 7)t quantile of the limiting distribution T constitutes a valid test procedure for Ho with an asymptotic size oy for any y E (0, 1). Addition ally, items 1 and 2 of Theorem 8 imply that n. nT diverges to infinity under in if the deviation from the null hypothesis decays slower than n 1/2 ensuring the power of this test procedure to approach one as n + oc. Together. these prop- erties imply the minimax-optimaliyy of the test procedure based on the test statistic T.. provided the rejection thresh- old can be computed. On the other hand, items 3 and 4 of Theorem 8 show that the test has trivial power in detecting alternatives that are too close to the null. This is expected, as even for parametric linear model, a meaningful CMI test cannot detect local alternatives with deviations from the null decaying faster than the n- 1/2 parametric rate (e.g., see Section 1.1 of Lundborg et al. (2024)).
3. Simulation Results
To evaluate the finite-sample performance of our test, we adopt the two examples from Cai et al. (2024), where the response variable Y is univariate. Each experiment is re- peated 500 times with sample sizes n E (200, 400, 8003. The nominal significance level is set at 5%.
Example A1: Consider the linear regression model Y;
87 2, + 3xX, + f i E [n]. where E i.s.d. ~ (.,0 are independent of {Xi, Di) dx = dz = 25 and
Finally, we demonstrate that our proposed test in Section 2.2 based on bootstrap climinates the need to compute the rejec- tion threshold, and the resulting test procedure is asymptoti- cally indistinguishable from the optimal test based on nff, For a generic bootstrapped statistic M, defined according to equation (2) with )im replaced by independent sam- ple {ei from standard normal distribution, we say that MM converges in distribution in probability to random variable B if, for any subsequence MM*x there exists -M;; further subsequence Mnny such that the conditional dis-
xT) iid. NOO, Σ) with the (i, element of being Eng - 0.3 for E [50]. We set the first two compo- nents of Bz one and the rest as zero. For this example, the null Ho corresponds to Bx = =0. Under the sparse alter- native, the first two components of 3x are 0.2/ and the rest as zero. Under the dense alternative, every component of 3x is fixed at 0.2/v25- = 0.04.
Example A2: Consider the nonlinear regression model Y 8 Z; + (B]xx,² + E, for i E [m], under the same B2 and (X,, Zi, ] as in Example Al. The null Ho also corresponds to 3x= = O. Under the sparse alternative, we set the first five components 3x to be 10 1/2 and the rest as zero. Under the dense alternative, we set the first twelve components of Bx as 142 and the rest zero.
tribution of Milyyy given the data {Xi, Yú, Zi, 1/,"" 1o0 Ji.m=1 converges in distribution to B* almost surely (c.g.. see Defi-
nition 2.1 of Zhang et al. (2018)). We use the notation + to denote convergence in distribution in probability.
For T,., we opt to use the Laplacian kernel and the band- width parameter for each kernel is selected according to the median heuristic (Gretton et al., 2012, Section 8). For comparison, we also include the simulation results for the CMI test proposed in (Williamson et al., 2023, Algorithm
Theorem 9. Suppose Assumptions 4 holds, then we have (see Appendix D.3 for precise definitions of relevant con- stants and random variables),
1. Under Ho, nMM D. Tt.
7

---

Testing Conditional Mean Independence using GNN
Table 2. Empirical size and size adjusted power for Examples A1 and A2.
PMIT
pMrr
pMIT M
PMIT XGB
PCM
PCM
VIM
DSP
DSP M
To
T,
XGB
DNN
XGB
DNN
XGB
DNN
DNN
Oracie
5.6
6.4
5.8
7.0
6.2
8.4
8.0
9.6
2.0
0.0
3.2
0.0
7.0
Ho
400
6.0
4.2
8.0
6.8
11.8
9.66
14.2
13.2
3.4
0.0
4.8
0.0
0.0
6.8
7.0
800
7.0
11.2
9.8
12.8
14.2
17.2
20.2
21.0
1.2
0.0
4.8
0.0
0.0
7.2
5.4
200
3.8
9.4
5.6
9.1
3.2
5.6
4.8
7.6
11.0
25.4
12.6
15.8
15.8
88.4
36.0
Example A1 H1 spanse
400
18.6
34.6
30.4
35.8
11.8
35.2
24.8
56.8
27.0
81.8
14.8
14.8
29.8
100
98.6
800
91.6
66.0
97.8
66.4
96.0
96.2
99.8
95.8
87.4
99.8
31.6
52.8
91.5
100
100
200
6.0
11.6
7.0
11:2
3.0
15.8
4.6
20.0
15.6
30.0
13.2
22.4
43.8
95.4
56.0
deuse
17.6
63.6
32.2
65.2
11.0
83.6
20.6
86.3
29.2
17.0
47.2
76.8
100
100
800
86.6
94.8
93.6
95.2
85.6
100
99.0
100
91.2
100
38.6
80.2
98..8
100
100
200
5.6
6.4
5.8
7.0
6.2
8.4
8.0
9.6
2.0
0.0
3.2
0.0
0.0
7.0
8.8
Ho
400
6.0
4.2
8.0
5.8
11.8
9.6
14.2
13.2
3.4
0.0
4.8
0.0
0.0
6.8
7.0
800
7.0
11.2
9.8
12.8
14.2
17.2
20.2
21.0
1.2
0.0
4.8
0.0
0.0
7.2
5.4
200
7.4
7.6
13.0
12.2
17.6
10.2
27.0
21.0
10.6
17.2
10.6
46.4
74.6
36.8
19.2
Example A2 H2 sparse
400
27.2
24.8
17.0
29.8
28.4
50.6
57.4
00.4
23.2
82.0
26.0
73.0
90.8
80.0
76.2
92.2
42.5
98.6
45.6
97.6
90.0
91.6
94.2
100
71.8
99.6
100
200
1.6
8.2
12.0
12.8
14.4
9.4
25.0
24.8
6.6
11.2
10.6
59.4
79.4
23.8
17.6
H1 dense
400
14.0
32.0
27.2
36.8
15.0
67.4
33.0
75.8
13.4
43.0
12.6
84.0
96.4
57.4
57.8
800
44.2
60.0
67.6
61.6
39.8
98.4
79.4
99.2
58.0
97.2
312
99.6
100
99.88
99.2
3) (denoted as VIM), the single and multiple split statistics proposed in (Lundborg et al., 2024, Algorithm and PB) (denoted as PCM and PCM respectively), the single and multiple split statistics proposed in (Dai et al., 2022, equa- tions (3) and (7)) (denoted as DSP and D respectively), the single/multiple split CMI tests proposed in Cai et al. (2024) (denoted as PMIT and pMITM M respectively) as well as their power enhanced versions (denoted as PMIT. and PMITeM respectively). For the four tests proposed in Cai et al. (2024), we include the simulation results when the conditional mcan functions are learned using eXtreme Gra- dient Boosting (XGB) and DNN. In addition, we also show the simulation results for an oracle version of our statistic Tu, where the true CME and conditional mean function of Y are used instead of their estimators.
proposed in Cai et al. (2024) depends heavily on the sparsity of the alternative and the estimation method used. The XGB estimation method has better performance under sparse al- ternative while the DNN estimation method outperforms under the dense alternative.
4. Imaging Data Applications
We apply our proposed CMI method to identify important facial regions for two computer vision tasks: recognizing facial expressions and predicting age.
4.1. Facial expression recognition
We examine whether covering specific facial regions affects facial expression prediction accuracy using the FER2013 dataset (Goodfellow al., 2013), which contains 48 48 grayscale images labeled with one of seven expressions: angry, disgust, fear, happy, sad, surprise, neutral. Following Dai et al. (2022), we analyze seven hypothesized regions (HRs): top left corner (TL), nose, right eye, mouth, left eye, eyes, and face (see Figure 2). After preprocessing as in Section 6.D of Dai al. (2022), we use 11,700 image-label pairs, denoted {(X,, 11700 and generate images Z by masking HRs.
For Example A1, as shown in Table 2, VIM and T, have relatively accurate size under Ho, while PCM and DSP (as well as their multiple split versions) are severely undersized. The empirical size for PMIT with XGB estimation method is close to the nominal level, but it is oversized with DNN es- timation method. The PMIT tests with multiple split and/or power enhancement all have large size distortions and the size distortion gets larger as increases. For the size ad- justed power, our test Tn outperforms all other tests for all values of under both sparse and dense alternatives, and VIM has the largest power loss although it also utilizes sam- ple splitting and cross fitting. Note that PMIT with DNN estimation method actually has larger power than PMIT which is supposed to be the power enhanced version of PMIT.
The test statistic T, is computed on ten subsamples (size n = 2000), with the resulting p-vauues plotted in Figure 1, alongside p-values for DSP statistics (Dai et al., 2022) under 0-1 loss and CE loss. To assess testing results, we compare test accuracies from a VGG network (Khaireddin & Chen, 2021) trained (Y,, Z4) for each HR against the baseline accuracy using (X,, Y.). As shown in Figure I, T, correctly identifies the nose and TL as non discrimmiiniaiee regions (p-values above the 5% level), while rejecting Ho for other HRs, consistent with their lower test accuracies. DSP p-values vary by loss function, with CE loss ex- hibiting stronger detecting power but inflated type-I error for TL and nose. Additionally, DSP M results do not align with accuracy trends, highlighting inconsistencies in region
For Example A2, as shown in Table 2, the empirical size results are the same in Example A1. For the size adjusted power, DSP and DSPM have the best overall performances (especially when the sample size is small) and our test T, have similar power performance as DSP and DSPM As when n = 800. Note that the power performances of the tests
8

---

Testing Conditional Mean Independence using GNN
100
baselin acc: 70.89
100
Face MAE: 8.42
78
p-value
accuraoy
D
p-value
PMIT
00
USP GE loss
80
5.7
pMi
acc
bato MAE
MAE
25
5.8
1808 5%
Face 8DC. 37.80
level 556
6.5
T1 Nove Right Fye Morrh
Left Fyn
Fyex FACE
T1
Nown
Minut
Fyns
Face
region
region
Figure 1. Box plot of the p-values (left y-axis) and the test accura- cies (red line, right y-axis) for different HRs. The blue dashed line represents the baseline accuracy. The red dashed line represents the 5% nominal level. The test accuracy for the face-coverey case (face acc: 37.60) shown the bottom right corner.
Figure 3. Box plot of the p-values (left y-axis) and the test MAE (red line, right y-axis) for different HRs. The blue dashed line represents the baseline MAE, The red dashed line represents the 5% nominal level. The test MAE for the face-covered case (face MAE: 8.42) shown the bottom right comer.
Figure 4. Original facial images in UTKFace (first column) and the covered images with HRs: TL, nose, mouth, eyes, face (Columns 2-6).
Figure 2. Original facial images in FER2013 (first column) and the covered images with HRs: TL, nose, right eyc, mouth, left eye, eyes, face (Columns 2-8).
rejected for eyes and face due to significantly higher test MAEs. pMIT as consistently produces low p-values, even when TL is covered (test MAE close to baseline), indicating oversizing. Meanwhile, pMIT fails to reject Hg in some cases (e.g., eyes, mouth) despite higher test MAEs than the nose-covered case, highlighting inconsistencies. Further details and more discussions of the results can be found in Appendix B.2
importance detection. Further details and more in-depth discussion of the results are provided in Appendix B.I.
4.2. Facial age estimation
We investigate the impact of covering specific fa- cial regions age prediction accuracy using the cropped and aligned UTKFace dataset (Zhang et al., com/
2017), available at http::://ww..kggee datasets/ a hatjjaaatk-ae farck-face--ao cropped. Five hypothesized regions (HRs) are analyzed: top left corner (TL), nose, mouth, eyes, and face (see Figure 4 and Ap- pendix A.3). After converting images to grayscale and selecting age labels between 20-59 years, we obtain 16,425 image-label pairs, {(X,, 16425 and the masked images {2,)]6425 The statistic T,, is computed on ten subsamples (n 2000 each), with p-values and mean absolute errors (MAE) from EfficientNet BO model (Tan & Le, 2019) plotted in Figure 3. Baseline MAES from (X,, Y;) are also included for comparison.
5. Conclusion
In this paper, we propose a novel conditional mean inde- pendence test that addresses limitations of existing methods. Using RKHS embedding, sample splitting, mooss-itting, and deep generative neural networks, our fully non-parametric test handles multivariate responses, performs well in high- dimensional settings, and maintains power against local alternatives converging at the n parametric rate. Simu- lations and data applications demonstrate its effectiveness. Some potential future directions include exploring condi- tional quantile independence testing, developing variable selection methods with false selection rate control, and per- forming diagnostic checks for high-dimensional regression models.
Figure 3 shows that T, p-values are above the 5% level when TL, nose, or mouth covered, indicating these re- gions are not critical for age estimation. However, Ho is

---

Testing Conditional Mean Independence using GNN
Impact Statement
Dai, B., Shen, X., and Pan, W. Significance tests of feature relevance for black-box learner. IEEE transactions neural networks and learning systems, 35(2):1898-1911 2022.
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which feel must be specifically highlighted here.
Daudin, J. Partial association measures and an application to qualitative regression. Biometrika, 67(3):581-55) 1980.
References
Att-Sahalia, Y., Bickel, P. J., and Stoker, T. M. Goodness- of-fit tests for kernel regression with an application to option implied volatilities. Journal of Econometrics, 105 (2)3363 412, 2001.
Dehling, and Mikosch, Random quadratic forms and the bootstrap for u-statistics. Journal of Multivariate Analysis, 51(2):392-413, 1994.
Delgado, M. A. and Manteiga, W. G. Significance testing in nonparametrie regression based on the bootstrap. The Annals of Statistics,
Athey, S., Imbens, G. W., Metzger, J., and Munro, E. Using wasserstein generative adversarial networks for the design of monte carlo simulations. Journal of Econometrics, 240 (2):105076, 2024.
2001.
Dziugaite, G. K., Roy, D. M., and Ghahramani, Z. Training generative neural networks via maximum mean discrep- ancy optimization. arXiv preprint arXiv: 505.03906, 2015.
Aubin, J.-P. Applied functional analysis. John Wiley & Sons, second edition, 2000.
Bang, H. and Robins, J. M. Doubly robust estimation in missing data and causal inference models. Biometrics, 61 (4):962-973, 2005.
Fan, Y. and Li, Q. Consistent model specification tests: omitted variables and semiparametrie functional forms. Econometrica: Journal of the econometric society, pp. 865-890, 1996.
Baptista, R., Hosseini, B., Kovachki, N. B., and Marzouk, Y. Conditional sampling with monotone gans: from genera- tive models to likelihood-fie inference. arXiv preprint Stti:22006.06555. 2020.
Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. Generative adversarial nets. In Advances in Neural Information Processing Systems, volume 27, 2014.
Bauer, B. and Kohler, M. deep learning as remedy for the curse of dimensionality in nonparametric regression. The Annals of Stattistics, 47(4))2261 2285, 2019.
Goodfellow, I. J., Erhan, D., Carrier, P. L., Courville, A., Mirza, M., Hamner, B., Cukierski, W., Tang, Y., Thale,, D., Lee, D.-H., et al. Challenges in representation learn- ing: A report on three machine learning contests. In Neural information processing: 20th international con- ference, ICONIP 2013, daegu, korea, november 3-7, 2013. Proceedings, Part III 20, pp. 117-124. Springer, 2013.
Bengio, Y., Courville, A., and Vincent, P. Representation learning: A review and new perspectives. IEEE transac- tions on pattern analysis and machine intelligence, 35(8): 1798-1828, 2013.
Berahmand, K., Daneshfar, F., Salchi, E. S., Li, Y., and Xu, Y. Autoencoders and their applications in machine learning: a survey, Artificial Intelligence Review, 57(2): 28, 2024.
Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A. A kernel two-sample test. Journal of Machine Learning Research, 13(25):723-75) -7333, 2012.
Hsing, T. and Eubank, R. Theoretical foundations of func- tional data analysis, with an introduction to linear opera- tors, volume 997. John Wiley Sons, 2015.
Cai, J., Luo, J., Wang, S., and Yang, Feature selection in machine learning: A new perspective. Neurocomputing, 300:70-79. 2018.
Huang, J., Jiao, Y., Liao, X., Liu, and Yu, Z. Deep di- mension reduction for supervised representation learning. IEEE Transactions Information Theory, 2024.
Cai, I.,, Guo, X., and Zhong, W. Test and measure for partial mean dependence based machine learning methods. Journal of the American Statistical Association, O(ja): 1-32, 2024.
Kallenberg, O.
Foundations of Modern Probability. Springer, second edition, 2002.
Cai, Z., Lei. J., and Roeder, K. Model-free prediction test with application to genomics data. Proceedings of the National Academy of Sciences, 119(3)):2205511819, 2022.
Khaireddin, Y. and Chen, Z. Facial emotion recognition: State of the art performance on fer2013, arXiv preprint arXiv:2 00.808888 2021.
10

---

Testing Conditional Mean Independence using GNN
Lavergne, P. and Vuong, Q Nonparametric significance testing. Econometric Theory, 16(4):576-601, 2000.
Shi, C., Xu, T., Bergsma, W., and Li, L. Double generative adversarial networks for conditional independence testing. The Journal of Machine Learning Research, 22(1):13029-
Lee, A. J. U-statistics: Theory and Practice. Routledge, 2019.
13060, 2021.
Shi, Y., De Bortoli, V., Deligiannidis, G., and Doucet, A. Conditional simulation using diffusion schrödinger bridges. In Uncertainty in Artificial Intelligence, pp. 1792-1802. PMLR, 2022.
Lee, C. F., Zhang, X., and Shao, X. Testing conditional mean independence for functional data. Biometrika, 107 (2):331-344 346, 2020.
Lehmann, E. L. and Romano, J. P. Testing Statistical Hy- potheses. Springer Science & Business Media, 2006.
Song, L.., Huang. J., Smola, A., and Fukumizu, K. Hilbert space embeddings of conditional distributions with appli- cations to dynamical systems. In Proceedings of the 26th Annual International Conference on Machine Learning. pp. 961-968, 2009. ISBN 9781605585161.
Li, B., Lai, Y.K.., and Rosin, P. L. Example-based image colorization via automatic feature selection and fusion. Neurocomputing. 266:687-698, 2017.
Szabó, Z. and Sriperumbudur, B. K. Characteristic and universal tensor product kernels. Journal of Machine Learning Research. 18(233):1-29 2018.
Li, Y., Swersky, K., and Zemel, R. Generative moment matching networks. In International Conference on Ma- chine Learning, pp. 1718-1727, 2015.
Tan, M. and Le, Q. Efficientnet: Rethinking model seal- ing for convolutional neural networks. In International conference on machine learning, pp. 6105-6114. PMLR, 2019.
Lundborg, A. R. Modern methods for variable significance lesting. PhD thesis, Apollo University of Cambridge Repository, 2022.
Lundborg, A. R., Kim, I., Shah, D., and Samworth, R. The projected covariance measure for assumption-leat variable significance testing. The Annals of Statistics, 52 (6):2851-2878, 2024.
Tang, R. and Yang. Y. Minimax rate of distribution estima- tion on unknown submanifolds under adversarial losses. The Annals of Statistics, 51(3):1282-1308 2023.
Van Landeghem, S., Abeel, T., Saeys, Y., and Van de Peer, Y. Discrininative and informative features for biomolecular text mining with ensemble feature selection. Bioinformat- ics, 26(18 :::55 i560, 2010.
Maqsood, M., Nazir, F., Khan, U., Aadil, F., Jamal, H., Mehmood, I., and Song, O.-y. Transfer learning assisted classification and detection of alzheimer's disease stages using 3d mri scans. Sensors, 19(1 1):2645, 2019.
Verdinelli, I. and Wasserman, L., Feature importance: A closer look at shapley values and loco. Statistical Science, 39(4):623-6)) 03633-6366661 2024.
Murdoch, W.J. Singh, C., Kumbier, K., Abbasi-Asl, R., and Yu, B. Interpretable machine learning: definitions, meth- ods, and applications. arXiv preprint arXiv:1901.04592. 2019.
Weidmann, J. Linear operators in Hilbert spaces. Springer Science & Business Media, 1980.
Nguyen, B., Nguyen, B.. Nguyen, H. T., and Nguyen, V. A. Generative conditional distributions by neural (en- tropic) optimal transport. In Proceedings of the 41st In- ternational Conference on Machine Learning, ICML'24. JMLR.org, 2024.
Williamson, B. D., Gilbert, B., Carone, M., and Simon, N. Nonparametric variable importance assessment using machine learning techniques. Biometrics, 77(1):9-22, 2021.
Park, J. and Muandet, K. A theouric approach to kernel conditional mean embeddings. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), Advances in Neural Information Processing Sys- tems, volume 33, pp. 21247-21259. Curran Associates, Inc., 2020.
Williamson, B. D., Gilbert, P. B., Simon, N. R., and Carone, M. A general framework for inference on algorithm- agnostic variable importance. Journal of the American Statistical Association, (64544)):665-16655 2023.
Yang, L., Zhang, Z., Song, Y., Hong, S., Xu, R., Zhao, Y., Zhang, W., Cui, B., and Yang, M.-H. Diffusion models: A comprehensive survey of methods and applications. ACM Computing Surveys, 56(4):1-3) 2023.
Rudin, W. Functional analysis. McGraw-Hill, New York, second edition. edition, 1991.
Shah, R. D. and Peters, J. The hardness of conditional inde- pendence testing and the generalised covariance measure. The Annals of Statistics, 48(3):1514-1533) 2020.
Zhang, X., Yao, S., and Shao, X. Conditional mean and quantile dependence testing in high dimension. The An- nals of Statistics, 46(1):219-244 246. 2018.

---

Testing Conditional Mean Independence using GNN
Zhang, Y., Huang, L.,, Yang, Y., and Shao, X. Doubly robust conditional independence testing with generative neural networks. arXiv preprint arXi::2407.1 17694, 2024.
Zhang, Z., Song, Y., and Qi, H. Age progression regression by conditional adversarial autoencoder. In Proceedings of the IEEE conference computer vision and paltern recognition, pp. 5810-5818, 2017.
Zhou, X., Jiao, Y., Liu, J., and Huang, deep gener- ative approach to conditional sampling. Journal of the American Statistical Association, 118(543):1-12, 2022.
Zhu, X. and Zhu, L. Dimension reduction- based signif- icance testing in nonparametric regression. Electronic Journal of Statistics, 12(1):1468-1506 2018.
Zhuang, F., Qi, Z., Duan, K., Xi, D., Zhu, Y., Zhu, H., Xiong, H., and He, Q. A comprehensive survey on transfer learning. Proceedings of the IEEE, 109(1):43-76, 2020.
12

---

Testing Conditional Mean Independence using GNN
Appendix to "Testing Conditional Mean Independence Using Generative Neural Networks"
The Appendix is organized as follows: Appendix details the training procedure and hyperparameter selection for the GMMN generator and DNN used in Sections 3 and Appendix B provides an in-depth analysis and discussion of the results from Sections 4.1 and 4.2. Appendix C contains the proof of Proposition 3. Appendix D formalizes and proves Theorems 5-9. Finally, Appendix E includes an essential auxiliary lemma and its proof.
A. Implementation details In this section, we offer additional details about the training procedure and network structure of the networks used in the paper.
In practical applications, the conditional distribution Px|z is often unknown, and we use a generative neural network (GNN) to train a conditional generator G R x 2 - Rxx for approximate sampling from Pxiz- By sampling latent variable from simple distribution (c.g., standard normal) over IR",, the conditional distribution of X G(n 2) given Z (denoted as 2z)z) serves as a good approximation of Pxiz-
We adopt the GMMN framework for its strong performance and alignment with our test. The conditional generator G is trained by minimizing the sample-based squared MMD between Pxz and Px using training data (X,, Z,) 2 and Mnt latent variables {nt 1, ny. 1, I.....M):
M):
t = arg min GES
Σ Kef k.EE(TT]
U(X,, (Z Zt),
(5)
ur (nT
M
with - Kx
1
Σ
Kx (Xes (x 2e)
M
m=1
Σ [xe Gx Zx)) + M
1
1
Σ
(xx( Zι), Gx(r Zt)).
M
where G represents an approximation family, such as (deep) neural networks. In the original GMMN framework by Dziugaite et al. (2015), the Monte Carlo sample size M fixed at one. However, our results indicate that increasing M significantly enhances the empirical performance and stability of the training process.
Let na denote the batch size, Nr is the total training epoch, 7x is the learning rate of Adam optimizer, Md, is the Monte Carlo sample size mentioned above, dm is the dimension of the input noise 7, and p is the distribution of the input noise. The table below shows how the above mentioned hyperparameters are selected in each section. Other hyperparameters of the GMMN used only in some specific section is discussed in Appendix A. 1-A3.. The training procedure is outlined in Algorithm I.
Section
726
N/r
dm
p (0, 1/3 1s0)
3
128
1000
4.9 1.90100 3
10
50
4.1
64
300
1 1-103
10
2304
N(O,, 1/10- 12304)
4.2
16
300
1. 10-3
10
50
N(0,, 1/10- Iso)
Table Hyperparameter selection for G
A.1. Numerical simulation in Section
Model structures and hyperparameters
To train the DNN πy we use the model structure and hyperparameters outlined in the following table with the MSE loss. To train the GMMN G, we follow Algorithm with loss Lx = L'x + L%x where L'x and L% are defined in equation (5) using the Laplace Kernel and Gaussian Kernel, respectively. Let y be the learning rate of Adam optimizer, BN represents whether batch normalization is used, and I "RLL is the activation coefficient of the leaky ReLU activation function. The hyperparameter values used in Section summarized in the following table.
13

---

Testing Conditional Mean Independence using GNN
Algorithm 1 Training Conditional Generator G
Input: batch size n,, data {(X,, Y, Z.) Hejn), total training epoch Nr, learning rate of Adam optimizer 7x, number of synthetic samples My. distribution of the input noise p Initialize the conditional generator G Compute the batch number b, [n/no]
for i 1 to Nr do Randomly split the data {XX,, Y, ki(m) into b, batches for to n, do
Select (.,, Z,) from the b-th batch
Sample 11," ~ p for [mo]. m [Mr]
Compute = Z.) and loss Lx based on equation (5)
Update G based on loss Lx using Adam optimizer with learning rate XX
end for end for Output: G
77b
Nr
Y
# of hidden layer
# of nodes
BN
T RELL
drop out
128
1000
4.9 - 10
128
No
0.8
0.05
QY
7
1000
2.2
10-3
256
No
0.7
0.15
Table Network structure and hyperparameter selection for section 3
Baseline algorithms
We implemented the PMIT, pMIT., PMIT M. PMIT PCM, PCMM. VIM, DSP, DSP methods using the code provided by Cai et al. (Section 5.1, 2024) with the default settings.
A.2. Facial expression recognition application in Section 4.1
Sampling procedure
To compute the test accuracy, we adopt the following train-test split procedure: select half of the images from each emotion label (recall that there are 7 emotions) as the testing set (with a sample size of 5850), and use the remaining images as the training set (with a sample size of 5850). We train the VGG network following the settings in Khaireddin & Chen (2021). The trained network evaluated on the testing set to get the test accuracy.
To obtain a subsample with a size of 2000, randomly sample 2000/1 700 of the triples from each emotion. To compute Tn, then split this subsample equally into two folds, with each fold containing approximately 1000/11700 of the triples from each emotion.
Model structures and hyperparameters We use the VGG network from Khaireddin & Chen (2021) to estimate gy and follow their default hyperparameters to train the network.
The GMMN generator is designed to process input images and noise vectors and is trained by Algorithm I with loss Lx -'x, where L'x is defined in equation (5) using the Laplace Kernel. It starts with fully connected layer that transforms the noise vector into 48 x 48 feature map, followed by ReLU activation. This feature map is then concatenated with the input image, resulting in a tensor with dimensions 2 x 48 x 48. The concatenated tensor is passed through an intermediate block consisting of four convolutional layers. Fach convolutional layer has 128 output channels, a kernel size of 3, and padding of 1, with ReLU activations applied after each layer. This output of the intermediate block is a feature map of dimensions 128 x 48 x 48, which is passed through a final convolutional layer to reduce the number of channels from 128 (the final convolutional layer has kernel size of 3 and padding 1). A Sigmoid activation function is applied to produce the final output, which has dimensions 1 x 48 x 48.
For each HR, the synthetic data ((m) used to calculat the statistic î,, is constructed by replacing the covered region in Z, with the corresponding region from the output image of the trained GMMN.
14

---

Testing Conditional Mean Independence using GNN
Figure 5. Original facial images in FER2013 (first column) and the covered images with HRs: TL, nose, right eyc, mouth, left eye, eycs, face (Columns 2-8). From row to 7, expressions are "angry', 'disgust,, 'fear', 'happy'. 'sad", 'surprise". 'neutral'. The pixel locations (height range. width range) for each HR are: TL: (1:6, 1:6); nose: (25:33, 21:30); right eye:(11:27, 26:46); Case mouth: (35:46, 19:31); Case left eye: (11:27, 6:26); Case eyes: (11:27, 6:46); Case face: (5:45, 5:45).
Baseline algorithms
We implemented the DSP M (with 0-1 loss and Cross Entropy (CE) loss) using the code provided by Dai et al. (Section 6.D, 2022) with the default settings.
A.3. Facial age estimation application in Section 4.2
Sampling procedure
To compute the test MAE, we adopt the following train-test split procedure: select 50 images from each age label (recall that age labels ranging from 20 to 59) as the testing set (with a sample size of 2000), and use the remaining images as the training (with sample size of 14,425). We train the EfficientNet BO network (Tan & Le, 2019) (with the final layer modified to output single value) the training by minimizing the MAE loss. The trained network is evaluated the testing set get the test MAE.
To obtain a subsample with a size of 2000, we randomly sample 50 triples from each age label. To compute F,,, we then split this subsample equally into two folds, with each fold containing approximately 25 triples from each age label.
Model structures and hy perparameterss
We the EfficientNet BO network from (Tan & Le, 2019) estimate gy, with the following hyperparameters to train the network: batch size nu, = 128, total epochs Nr 300, learning rate 7 = 1 10-3, and the loss metric is MAE loss. We employ the ReduceLROnPP Patnnlaateau learning rate scheduler, which reduces the learning rate by a factor of 0.2 if the MAE of the testing fold does not improve for 10 consecutive epochs. The minimum learning rate is set to 1 x 10-6, To save computation time and avoid overfitting, we adopt the early stopping technique based on the MAE the testing fold, with an early stopping patience 25 epochs.
The GMMN generator trained by Algorithm with loss Lx = L'x. where L'x is defined in equation (5) using the Laplace Kernel and it adopts the same ReduceLROnPIaa Pateau learning rate scheduler and early stopping technique mentioned above. It starts with fully connected layer that transforms the noise vector into a 224 x 224 feature map, followed by a ReLU activation. This feature map is then concatenated with the input image, resulting in a tensor with dimensions 2 x 224 x 224. The concatenated tensor is passed through an intermediate block consisting of four convolutional layers. Each convolutional layer has 128 output channels, a kernel size 3, and padding of 1, with ReLU activations applicd after each layer. This output of the intermediate block is feature map of dimensions 128 224 224, which passed through
15

---

Testing Conditional Mean Independence using GNN
final convolutional layer to reduce the number of channels from 128 (the final convolutional layer has kernel size of 3 and padding of 1). A Sigmoid activation function is applied to produce the final output, which has dimensions x 224 224.
Figure 6. Original facial images in UTKFace (first column) and the covered images with HRs: TL, nose, mouth, eyes, face (Columns 2-6). The pixel locations (height range, width range) for each HR are: TL: (1:21, 1:21); nose: (101:141, 86:141); mouth: (131:211, 21:201); eyes: (31:101, 21: 21:201); face: (1:201,1:2244
For each HR, the synthetic data x/()) used to calculate the statistic Tr constructed by replacing the covered region in Z with the corresponding region from the output image of the trained GMMN.
Baseline algorithms
We implemented the PMIT, PMIT., e. PMITM- M. and PMIT methods using the code provided by Cai et al. (Section 6.1, 2024). We modified their model to use the same EfficientNet BO architecture that we employed. Additionally, we changed the loss function to MAE loss and did not scale the age label Y.
B. Detailed analysis of imaging data applications B.1. Facial expression recognition
In this application, we examine whether covering some region of facial image will influence the prediction accuracy of facial expression. We use the Facial Expression Recognition 2013 Dataset (FER2013, Goodfellow al., 2013) consisting of 48 x 48 pixel grayscale facial images, each attached with label from one of seven facial expressions: angry, disgust, fear, happy. sad, surprise, neutral (denoted as Expression 1-7). As in (Dai et al., 2022), we consider seven cases where different hypothesized regions (HR) are covered: top left corner (TL), nose, right eye, mouth, left eye, eyes, face; see Appendix A.2 and Figure 2 for locations of these HRs.
After applying the same preprocessing procedure as outlined in Section 6.D of Dai et al. (2022), we obtain 11700 image-label pairs which will be the samples used in this application. Let ((x,, 117000 denote the 48 x 48 pixel facial images and their corresponding labels. Note that for any i [11700] and E [7]. Y, E RT is one-hot vector with the jth component being one and the rest being zero if Expression j is associated with X,. For each HR, we use {Z, 2,111100 11700 to denote the facial images with the HR covered in black.
The statistic Tn, is evaluated ten times on different subsamples, with sample size n = 2000, from the 11700 triples {(X,, Y, ) 11700 and the box plot of the ten p-values are plotted in Figure 1. As comparison, we also include the p-values of the DSP statistics from (Dai et al., 2022) with different loss functions: 0-1 loss and Cross Entropy (CH) loss. To evaluate the testing results, calculated the test accuracy of VGG network (Khaireddin & Chen, 2021) trained/evaluated on the whole sample {(Y,, |1)00 since lower accuracy indicates alternative hypothesis with stronger signal. The resulting test accuracies for different HRs, as well as the accuracy for the VGG network trained/evaluated {(X,, Y)) 11700 (denoted as baseline acc), are also plotted in Figure 1; see Appendix A.2 for the sampling procedure and other implementation details.
As shown in Figure 1, the lower quantiles of the p-values from T, are above the 5% nominal level when nose or top left corner of the facial image is covered, indicating that these regions are not discriminative to facial expressions. For all the other HRs, the Ho is rejected since all the p-values from T,, are smaller than the nominal level. The test results from F, is consistent with the test accuracies for different HRs, since the test accuracies when nose or TL is covered are very close to the baseline acc, while for other HRs the test accuracies are noticeably lower than the baseline acc. The test results from M varies drastically when different loss functions are used. For DSP with 0-1 loss, the lower quantiles of the p palues are above (close to) the nominal level for the cases when left eye or both eyes are covered, whereas the test accuracies
16

---

Testing Conditional Mean Independence using GNN
for these cases significantly lower than baseline acc, indicating these regions are indeed important in detecting facial expressions. The DSPM with CE loss has stronger detecting power than DSP with 0-1 loss since the p-values for the former is in general smaller than the p-values for the later, which means stronger However, when TL or nose are covered (which correspond to the cases under Hy), the lower quantiles of the p-values from DSP with CE loss are smaller than the nominal level, which may result in inflated type-I error. In addition, the median p-values for different HRs from DSP and MMPM M do not decrease monotonously as the test accuracies decrease.
B.2. Facial age estimation
In this application, we investigate the impact of covering specific regions of facial images on the accuracy of age prediction using a well-cropped and aligned version of the UTKFace dataset (Zhang et al., 2017), which available at https: //www wwwgleg com/coo/ /atasetssssssseee s aaaaaaaiiiaaaaaaaaaa different HRs are covered: the top left corner (TL), nose, mouth, eyes, face; see Appendix A.3 and Figure for the locations of these HRs.
face-aacee cropped. We examine five scenarios where
After converting the images to grayscale and selecting the age labels ranging from 20 to 59 years old, we obtain 16425 image-label pairs, which will be the samples used in this application. Let {X,, ] denote the 224 x 224 pixel facial images and their corresponding age label Y; E R. For each HR, we use [2:],11 16425 to denote the facial images with the HR covered black.
The statistic Tú is evaluated ten times different subsamples, each with a sample size of 2000, from the triples {(X,, Y,, Z,)] fi=l 16425 The box plot of the ten p- p-values is shown in Figure 3. For comparison, also include the p-values of the PMIT and pMIT As statistics from Cai et al. (2024). To evaluate the testing results, we calculated the MAE (mean absolute error) of EfficientNet BO network (Tan & Le, 2019) trained/evaluated on the whole sample {(Y,, ] 16425, The resulting MAE for different HRs, as well as the MAE for the EfficientNet BO network trained/evaluuted using 16425 (denoted as baseline MAE), are also plotted Figure I; see Appendix A.3 for the sampling procedure and other implementation details.
{(X,, )
As shown in Figure 3, the median of the p-values from T, for different HRs decrease as the MAEs increase, and they above the 5% nominal level when TI., nose, or mouth (three HRs with the smallest MAEs) of the facial image is covered, indicating that these regions are not discriminative for age estimation. For the eyes- and face-covered cases, the null hypothesis Ho is rejected since the median of p-values from T, are smaller than the nominal level, which makes sense since the test MAE in these two cases are significantly larger than the bascline MAE.
The p-values from PMIT M are all lower than the nominal level, even for the TL-covered case where the test MAE is close to the baseline MAE This is consistent with the simulation result in Section where pMIT As is severely oversized. The p-values from PMIT does not change monotonously with the test MAE. Based on the median p- p-values, pMIT rejects Ho in the nose-covered case but failed to reject Hu in the month- and eyes-covered cases, even though the later two have larger test MAE than the nose-covered case.
C. Proof of Proposition 3
The following theorem comes from Theorem 2.2 in Park & Muandet (2020). Note that H l and TH6, are RKHSs defined in Section 2.1. Denote H.x & Hz as the tensor product space of H and Hz (Weidmann, 1980, page 47-48).
Theorem 10. HIx & Hz is generated by the functions 8 g Rdx R, with f E Hix and g E li defined by @ 9(I 3) = (x) Moreover, HIx @ Hz and H. are identical.
The following definition of Hilbert-Schmidr operator and Hilbert-Schmidt norm from Definition 12.1.1 and Theorem 12.1.1 in Aubin (2000).
Definition 11. Let A and B be two separable Hilbert spaces and {r: 1 is orthonormal base of A. A continuous linear
operator T A + B is Hilbert-Schmidt operator if the series
20
1/2
|| h
||Tr; B
(6)
converges, and
I| is called the Hilbert-Schmidt norm.
17

---

Testing Conditional Mean Independence using GNN
The following definition of the mean of Hilbert space valued ranodm element derived from Theorem 3.2.1 and Definition 7.2.1 in Hsing& & Eubank (2015).
Definition 12. Let X be random element taking values in separable Hilbert space A. If Bl]xII < 0, then x is Bochner integrable (Hsing & Eubank, 2015, Definition 2.6.3) and its mean Ex is a unique element in A such that (S, Ax) = F.J,x), x) for any / A.
C.1. Proof of part (a) in Proposition 3
We prove E is a filleet-Schmiie operator by showing that Me E Ho for any E Rdy and || < 0.
To show Ee E Ho, denote h1,,, V, 2, c) [Ko ((x,, =), -) [K ((X, =), ) [V ]|] c. then we have Ec Eh(X, Y, Z, c). By Definition 12, it suffice show: (1), h(zz, V, 2, c) E Hil for any H and almost every (2,
Rdx+dy+da +d2; (2), Y,J,(X,Y, CII,, < co. For part (1), it is clear that Ko 2), [ E[V|=]] E Ho. and
E[Ko 2). )) [] HHl follows from Lemma 15. So part (1) is proved. For part (2), note that
E | Y, Z, 11 -B{1 (o), (X, -) - B[Ko
)
||
[[
B c
H.
<41cBB || ((x, -11 )11 + ||B[[Ko .) 1 [. 1
| 12 + 11 (2)1
<
[||Ko >ll +B[I)Ko((XX,)) 8), )lE. 2] | + 11 (2)
(7)
=4 [CBE
[[X, 2), >Il,, + +|]Ko((( Z), 111 [[ | + 1l 11,r (2)1
(8)
where equation (7) follows from a generalized Jensen's inequality (Park& & Muandet, 2020, Appendix A) and X is random vector such that 2) d (X, Z) and X is conditionally independent of Y given Z. By Assumption 4(a), the right hand side of equation (8) is finite and part (2) is proved since [X|,(X, Y, Z, e)E S [B]]A(XX, Y, Z, )
To show |> < O, let (r: 1r be an orthonormal base of Rd where ry E Rे is vector with äth element being and the rest being zero for any E [dy]. Then we have
dv
Er: 11
1=1
dy
((X, ( -EE .11] [ T r | Ho
=
(9)
2E
|| (1X, 1 +IIE(Ko((X,2) -)) )|]
|| Y Bly 1
Eo
< <4 { [|1Ko [[X, 2), > +B[IKo(((X,2) II. z | + 11 11,y (2)1]
<00,
(10)
where equation (10) follows from equation (8). So Σ is a Hilbert-Schmidt operator.
To prove | its = 1 note that from equation (9) we have
|>| hs -Be
B[Ko 2), [o - 2). )121) f(,,a)
To
-B{(
- (17] [xx(X',, - E[[Xx(X' )12 Kz Kz xxxx)x)))) vorra}
E[xxx
(11)
(Xxx(Xx,))22)) ()11) Hx V(,,Y"))2(222) -E(x [[xxxxx (x, -XXXX) -E(UIX.X"VV,
E[xx(X',
(12)
where equation (11) follows from Lemma 15 and equation (12) follows from the reproducing property of RKHS. So part (a)
of Proposition is proved.
18

---

Testing Conditional Mean Independence using GNN
C.2. Proof of part (b) in Proposition 3
The operator norm (Rudin, 1991, Theorem 4.4) of Σ is defined as
|E sup Ec)a |lell s 1, ll s <1].
(13)
According to Proposition 12.1.2 in (Hsing & Eubank, 2015), |Ellrr S |[ If = |Ellis 0, then we have |E 0 which is equivalent to E being zero operator. If of O, then we have Σr; = 0 IHlo for any i E [dy], which implies 1* | = O. So we have that I* = 0 if and only if > is a zero operator, which is equivalent to Ho holds according to part (c) of Proposition 1.
D. Proofs of main results in Section 2.3
Define two oracle statistics
T
E XVX,
(14)
1)
2
itk j.k3
1 2 Tn
Σ (1X,
2
2
8=1
j.
where
X(, Xx) - En Kxx (X,, x(()) - Eng Kx (x(, .Xx) + Engin ne Kx (X(m)) (R(
(15)
Denote S; = =X;,, Z,), under Hin, it is easy to that Assumption 4 imply the following assumption:
Assumption 13. For i, on and D, defined in Assumption 4, we assume there is large positive number C3 such that
(a) Rxx0; Kx((;
1005)
Kz Zi) <Cs.
(b)
R. [=
2.) O(n-11) and
Rmn ",,, (2 =0(n
D.I. Proof of Theorem 5
First, we give formalized version of Theorem 5: Theorem 4 (formalized). Suppose Assumptions 4 holds, then under Ho.
T. Tr, (0p(n 'M 1/2 + 01 + n ann + n 1/2-
and nfT, Dg T = Ti + Tb, where Tn is defined in the same way as T, with Ú(X,,X)) X.), V(Y,, Yx) replaced by U(X,, Xx), V(Y,, Yx) (see equation (14) in Appendix D) and (T), T are i.id. random variables with Ti Ee -1). Here, X2 are i.i.d. chi chi-square random variables with one degree freedom and X2 's are eigenvalues of compact self-adjoint operator on L2(Rdx didx+dy +dy dy +dz Pxyz) induced by the kernel function 1(((1, YI, Zi), (X2, Y2, Z2)) V (Y1, K2) (Z,, Zz); that is, there exists orthonormal basis {f((X1 Y1, 11) of L2 (dd tdy PXYZ)
such that
E [[ Y, 71), (2, (2, V1, {,f((r,))))
Proof.
Since P, -T, = [M + 02 + 1/2 1/2- (a) +o2)) (p by Lemma 17, it suffice to show
19

---

Testing Conditional Mean Independence using GNN
nTn 2 Tt. Note that
nTn
12
Σ
UXX,, (,)
==1
jek j.kej
2 n T(n). 2
For any E [2]. T! is degenerate U statistic of order two and 2 TT!")) 2 T follows directly follows from asymptotic
theory for degenerate U statistic, see Theorem in Section 3.2.2 of Lee (2019) (also stated as Lemma 16 in Appendix E). Since (T(m), Th") are mutually independent, we have nT, 2 7t.
D.2. Proof of Theorem 8
First, we give formalized version of Theorem 8: Theorem (formalized). Suppose Assumptions 4 holds, then under Hin
1. If a = O. then vn(Ta co) 0, and . are i.d.
V2 -2 where = Xo) (2),,,,,, 22)} >
1
mean zero normal random variables with variance equal to X2, Y2, 21))
4 Var (E(U(XI (2,, Zz)
With Assumption 6 further satisfied, we have
2. 1f0<a<1/2 then T, Bn where c
(Xx, 22) >0.
3. Ifa = 1/2, then nfn & c 1 Ei= where = E}= TÀj and {T), T22t are i.i.d. random variables
with Th= the eigenvalues corresponding to the kernel function h((X1, Sı, 21), (X2, E2)) = ((X1, X2)SI (2K2('1, Zz) Here Gj are independent mean zero normal random variables, possibly correlated with Thy, with variance equal to 4--9
1), Xi are i.i.d. chi-square random variables with one degree of freedom and XAS are
4 Var (x) Zi) (2 G(X2,X2, Zz) Kx Z2) (2. X2, X2,
Zz))
4. Ifa > 1/2, n,, 3T
Proof. If a 0, by Lemma 17 we have Vnfa = VnI, + op(1), so we only need to show that Vn(T, - co) 4
J/2 -1=1 9/ Note that
2
Σ UX,, co) }
(T - co)
2
==1
jck j.kE(())
1
(n) co).
For any E [2], 7(1) is a non-degenerate U statistic with E7(") = cp and ( - -n)) D, ( follows directly follows
from asymptotic theory for non-degenerate U statistic, see Theorem in Section 3.2.1 of Lee (2019) (also stated as Lemma
16 Appendix E). Since (T(()) f(((0) are mutually independent, have Vn(T,, co) 4 +/2 3 9
20

---

Testing Conditional Mean Independence using GNN
If a 0, by Lemma 17 and Assumption 6 we have
nfn =nf, Op +11))+++))))))
2a
Σ
2
2
U(X,, (o)((s))) (w+O(()) (1)
2
1=1
2
fike)())
n1
Σ
X()(())) (Rax (a)+(())) (Ras - -si])
J4L
E
(o) (i) +R Ga) K2) Zn)
2 1)
J.kE
Σ
Xx) - - 22) 2x)
2
i4k
J.kE
Σ
Xilc) (K22 Zn) +o
2
-1)-1)
izk
=n 2-200 Ji n
J2 + n
nJ4 + nJs + n to (n2 Ja
(16)
Note that J1 4 degenerate U statistics. By Assumption 6 and part (a) of Assumption 13, we have J3 (n22) and Js op(n If 0<<a<1/2, nf.,- n 10 (1 + (1)) + cc. If > 1/2, nf. = nJa + (1) & TÁ follows in the same way as in the proof of Theorem 5. a 1/2, nf,, J1 + vnJ2 nJ6 + o1).). Following similar approach as in the proof for Theorem of Lee al. (2020), have vn.J2 nJa 4 1 Σ) 9j T and the proof is finished.
((),) -2)} = 0. 12 is mean zero non-degenerate U statistics and J3 is
v2
D.3. Proof of Theorem 9
First, we give formalized version of Theorem 9: Theorem (formalized). Suppose Assumptions 4 holds, then we have,
1. Under Ho, nMM 2 Tt.
O. nMA Σ = where {Tr,Te} are i.id random variables with Ti Σi- 1x(x - 1), X2 are Lid chi-square rundom variables with one degree of freedom and 7is are eigenvalues of h( Y,,2,)) Z.), (X2, Y2, Z2)) = X2V (2) Zz).
2. Under Hin with
With Assumption further satisfied, we have
3. Under Hin with > 0, nMM 22 T Furthermore, if we let MAM, denote the (1 - 7)th quantile of nM M conditioning on the data, then the power (probability
of detecting the alternative) of our testing procedure satisfies: if a<1122, then P(nf,, 2 MMM,)) + 1: if a = 1/2, then
P(nf, 2 MMM.)) + P(C+TA++ v2 Ez G 2 To where To,, denotes the (1 (1-9) 7)th quantile of if >> 1/2, then
P(nf, 2 M.M.) + 7.
Proof. Define
M, J =
Σ
Zweezet
3 1)
(.KE(()
21

---

Testing Conditional Mean Independence using GNN
and [I1, [,, Iz} according the same equations as (41 I2 [3] in Appendix E with Kz Zạ) replaced by Kzz Zegege. Denote E. and Pa the conditional expectation and probability of random variable conditioning on {X,, Y, Zí, 1/" N.
Ji=1,==1 First, we show that the difference between M, and M, is asymptotically negligible in the sense
that
P(P'(InM - nM, 2 c) 2 6) + 0 for any c <> 0,8 1-8000 0.
(17)
By Chebyshev's inequality, it suffice to show IE" Mn nMM1' 2 (1), which is implied by EEI +I3 op(n (2) Note that E [" +Ig] = op(n 2) is already shown in the proof of Lemma 17. For example,
E'I
Y2) ,
[a 1)]²
2 2
j jKEJ
and the right hand side is op(n=) according to part (a) of the proof of Lemma 17 for 71-
E 7
Σ
( 1]]
jzh j.EJ
and the right hand side op(n according to part (b) of the proof of Lemma 17 for 7
Note that nM, 2 I' under Ho and n.47 2 Ti under Hin with a = 0 is a standard result in bootstrapping for U-statistics (see Theorem in Lee al. (2020) or Theorem 3.1 Dehling & Mikosch (1994)).
remains to show n.M, 2 TA under Hin with a O. For [6]. define D, according to the same formulas as in equation (16) with K, (Zj. Za) replaced by Kzz ( Z ee Then we have
n.M = 20 + n - D2 + aD3 + nD4 + nDg + nDo.
It is easy to see that D* - 0,(n for i E [3] and E = (n for i 4,5 (by part (a) of Assumption 13), So the difference between n.M, and nD2 is asymptotically negligible and nD; D T follows from Theorem 4 in Lee al. (2020) or Theorem 3.1 in Dehling& & Mikosch (1994).
To show the asymptotic power of our test, let Müsa, denote the (1 7)th quantile of nMM conditioning on n.M
(X,, Y,, Zi, nith "/ m==1.m=1. By theorem 9 and part (ii) of Lemma 11.2.1 in Lehmann & Romano (2006), M M RMMr
2 To,
M M, 3. To, MAM, Fs Ti,, under Ho. HI, with and HI, with 0 respectively, where Tu,,, To, and Ti, are the (1 - 7hth quantiles of Ti, I and Ti respectively. Then the asymptotic power results in Theorem 9 follows.
E. Auxiliary Lemmas and Its Proofs
In this section, we give three lemmas that used in the proofs of Proposition 3 and Theorems 5-9.
The following lemma is used in the proof of Proposition 3. Lemma 15. If Assumption 4(a) holds, then E [[ ((X, )
E Hy is well defined and B((oo(XX,)) .)]²] =
(xx] Kz(z, (2.-) for almost every E 122..
Proof. By Assumption 4(a),
B(B[IIK((XX,)), )11oo11] -Ef =E(VIx(XX (E[XxxXXXXXx21 1/2 <0,
(18)
which implies [[]IKo((XXx)) 1|1 < oo for almost every z Rdz, So from Definition 12 we have that
E [[Xo((X, =), -)(2] E Ho is well defined for almost every E RRd2, Similarly, we can show [[[xx (X, []]] E Hi
22

---

Testing Conditional Mean Independence using GNN
is well defined for almost every 2 Rdz, For any E Hlx and y Hiz, by Assumption 4(a) we have [B)(()I 5 E/ (x,,xx and similarly Blf(X)g(Z) co. By Theorem 10 have
(B(Ku((XX, (),(x) n). =E ((x))(2)
EE
[(x) ((2))
(Xx (X, f, | (K2( 9)Ez
= [E[xx (x, )) ) K2 (z, J @ g)m
so we have E |Co ((X, =), ) = [xxxx (x, |]== for almost every 2 Rd and the lemma is proved.
The following lemma summarizes Theorem in Section 3.2.1 and Theorem in Section 3.2.2 of Lee (2019), and is used in the proofs of Theorems 5 and Lemma 16. Let Un be a U-statistic based on kernel 1(x1,, with of Var(B))((X,, and 0 = E[I((X, where {Xi]== are i...d. copies from distribution Px. We have
1. If o > then vn(UUn 0) N00
2. If o 0 and E[ X2)], then nUn, 4 Ee Ao(23 1), where {2,) are independent standard normal random variables, and the X, are the eigenvalues of the integral equation
J (tx))))))
The following lemma gives the rate at which the difference between the statistic T, and the oracle statistic Tn decays to zero, which is used in the proofs of Theorems 5-9. Lemma 17. Suppose Assumption 4 holds, then we have
+ n 11112- under Ho- 1. T, T, O (n [M 1/2 n 01 n 02]
2. Tn Tn (1 mwmma M 1/2
+ n n + n (it/a..//2.1100) under Hin.
Proof. The proof is divided into two parts:
(a) T, - = O0(M-1/22 nnn 1/2 (ax under Ho and Tr -T, = 00 1/22 +n 1/2 (/2 1/2 1/2-a])
under H Hin.
(b) - -p Hın-
under Note that in the following proof, we use C to denote generic positive constant which may vary each time it appears.
1-m (o0 under Ho and Tr Tw = ( 1/2,1) /
Proof of part (a):
suffice to prove the asymptotic order in part (a) for T(L) Ti(1) where
T())
Σ
(X),,,
2
2
j4k (.kEJ()
70
Σ (xX,,
(
j j.kE(()
Since Xx) = Xx) - H2 Xt). where
Wi(j m) =Kx(X, ((m))
(Xx, 8/ml) - Kx(PP)) {(()))
Kx EX,xxxX,... x((1) E, (xf Xxa) + Eygng (P, (()))
1 M
(X,,, M
Σ W.).kmm).
23

---

Testing Conditional Mean Independence using GNN
we have ) - /, where
11
Σ 37k J.KET
2
For any & E (0, 1). by Assumptions 4, 13 and Chebyshev's inequality, there exist large positive number C₁ such that for any positive integer n, P(An.) > where
An
44B.x:/h....4444.0....0001111111 I Can
|.Neraaa00.ccia00a4.../.22.20140 SCm
(D1, (10000/..//|1SS////55///11111111 <Ca:
B |3 (, s Cin
( (2) S Can
(19)
In the following proof, we fix {(2, Y, Z))ee(n) An and treat G and 91 as fixed functions. Note that Kz) XI,, (z) Zx)] = 0,
Hi (X, X2) V (,),,
S n2
-1(1) + (2)
First, we show that 1(1) I n²M under Ho or Hin- For integers m,, m2 E [M] such that m, * m2, we have B(W1(1,2 2, (1, mY S1, Y2, S2 = O, which implies
C
1( 11, 12.
C - (I) (G) (8,",, (P(0)
+xxx X S., S2. M
It can also be shown that ()),,))3 &
(Zi (2))
So from part (a) of Assumption (conditioning on ((x,, Y,, 7,) 7,) )ea((1) E A.). we know 7.")) < W.M For (22) we have that under
24

---

Testing Conditional Mean Independence using GNN
Ho
10 www.iiaeeawaiia.a............... Za)
nM
())xx
nM
x (Kx((2))))
Kx (X3 ,))xx
M 0r (2x)
x I or (2:)] (z (2)
C
B
Kx ",
nM
x
|| Kix
Ex |a (22,
x T I -)BBn Xx ( 22 (Z3, (z)
By part (a) of Assumption 4, E [|xx(X - Kxx Saa.:aaafo....ia.aa----//a/aaa/. < CCs. By part (b) of Assumption 4 and Cauchy-Schward inequality,
R I Ax (X2) -) Kxx / 72)
SCn Snn-a100
[ Xxx (Xg ,-) -)3 Kx Kx (3
2z)
SCn 011-2
So we have 119) S 2 M
under Ho. Under Hin, since E| 0.
nM (1
C
Kx (211) ,))xx
x (Ki(( ,) - Xxx Byg Kxx (8."), :))xx
Yi ãr (Z))] (y) (Z)) (g(())
x (4 22)222 20]
(20)
By Assumptions and 13, we have
(xxx((Xss -B, (88", ,))xx
x (Vi-ov (2.) (1 (2) 20)
(8 xx [|
x E - Kx h1lp::/...4.....././//ee/...//..
xE
SOn -0-00-00
Kix "|/p::/hhhh.aa.a.............r.r (1. M)x
(r (21)] |
9r (9r |, (2)) I
SE || ) E Xx 11
E E ) Kx
(Zz) (l) (,)
Kz (Zz Zz)
| Kx (x . By Xx (xx (Zs) gr (2s) 2/ Kz (Zs. Zs)
SOn 2((2)az))
25

---

Testing Conditional Mean Independence using GNN
and
Kx
( [(t) [1] ( 22) Z3))
< |[xxx (1 :) - / Kx (K
R - Xx 11xx |1k
K2
x [Xxx ) X Xx 1 2z)
<Cn
which implies 1,2) s CM- (n 2(0 (otor))) + n -1) I CM (n-2 + +n -1- under Hin- Combining with the result for I(1), we can conclude that conditioning on {X,, Y, Z;) bey-n E An. E[I] <
CM + n 2(0) under Ho and E[4] < (n 2a + n 2 + n 1 2(at)))) under Hin- Let Ah
denote the complement set of An. Under Hg, since
1
1
P( n ++n 1/2-01-02 C) <P(
+
n -n 1/2 022-
1
P(
- 1/2-11-2-0000 > C[((X, Ys, Zi)kett))) LAn LA. (e)tt-))) + +P(A.)
1 C2
- +
(21)
and e can be arbitrarily small, we can conclude that, without conditioning on {(x,, Y, ki) 11 (M 111=-112(0o+00)) ) under Hg. Similarly, we can prove 1 =
(11202 -n 1/2--o)
under Hin- So part (a) is proved.
Since the same conditioning argument as in equation (21) can be used to derive the asymptotic order without conditioning on the training data, from now on, we treat all the neural networks as fixed functions.
Proof of part (b):
suffice to prove the asymptotic order in part (b) for T())) T,()) where
((1)=
Σ
Xx,,,, (z)K2 (3 7k)-
2
J.REJ(
Since Ü(X,,X*) Xx) = U(X,, - F,, (X,, Xx) and V Vb) = V(S, - F22 Yx), where
(X,, Xx) =R (x,, (((1) + Fyg Kx F,g (x)) , xx) - Bagre Kx (2)0), ((D) , (X, j. x(1) K. (871) .Xs) (X,) + Bugtot Xx (X, (x,)) (
HrYY, Yx) =Y -Y) gr (Zx)
+ gy (2,) Tr (Z,) (Zk)
- YjJ gr (xx) gr (2,) + (Z,) Tgy
gr (Zx),
then we have 7,(1) T-(1) = - - ī + Tsg where
i i= 2 - 1)
Σ XXX,
jfk
j.kE.(")
7,
Σ U(X, (,X,,,,)))B,,,,3))333
2
j&
j.kE)(()
1g
Σ
(
jtk
J.kEg()
26

---

Testing Conditional Mean Independence using GNN
Consider I2. Under Ho.
http://kw....nnoo0....///////311 Spl 0,
(V),,,
so we have
Var(T) (1)
(22)
2(2 1)
Note that V(V,, 12)2 = 1 (3
(2)) (22)18))
Xx), ",
(x, 210)
Cx ( (x())
- By Kx [11))] + [ Kx xx Bot Ex (2)1) (2)))
(23)
and -x
(x,, Xx
X(") Kx (x M2x so by part (a)
and (b) of Assumption 4, we have I = Op 1-m)). Under Hin
i
Σ
X (Rn))) Zn)
2
j
J())
Σ
Xxxx
RRw) (ennjK2(244 KZ((Z
7(2 1)
jtk
),ke
Σ
- 1)
2
1
Σ
Xx), [Ros] [Rai] Kz(Z,, Zz)
- -)
jfk
(1)
= 71777,2 7100+
We now show 7,)) 0p(n-112-00 (same result holds for 7,20 and the proof is omitted). Note that
B [x,, Xx) [ [[S;)] [RRA]] K2 ] O, by equation (23), we have
Var(7 ((11)
E [B Kx (X1, X60) - B (X1, "","":"""""""""""""""""""""""""" "1,.
xx [ N2, x1)) - Kx (x,, "1G.1h/h006111133331333333330010 70)
F
x()) Butty
+ [F (x1) x R (x (1sa://F.A......1o.o3333//3333.1 72) }
=0(n -1-20-211 )
where the last equality follows from Assumption 13. So we have 7. -1/2-a0)) under Hin. For 7(40) the conditional expectations of (X), [[R..][[ Kz (2j, Zx) given (Rnj, ST) or (Rni, Sk) are zero so 7, -1-ai) by Assumption 13. Following similar approach as for 7,1). we can show Var (7(3)) O(n -1-da-a--20 20s ), which implies 7(0) = Op(n 112-22--u under H... So we have 1 Op (n 12122---0--00 +n 11aa) under Hm..
27

---

Testing Conditional Mean Independence using GNN
Now consider Tz. Under Ho, Tg Op(n=1-az) can be shown similar the proof for T1. Under H let
1-22)
Kx (x), x(0), Xx - Kx (x)", xx
-,,,,,,)-,,, gr (Zk);
(1) (2- (2,) Yk:
(2)) ( (Zx),,
then we have
T22
E U(X,, (,, V2) (,) (2,)))2222,
701) +12
7,
7000
1
Σ
U(X, Xx) (r)
(Ze) gx
(Zx) (
Kz( Zx)
2 -
*
J(KE
+
1
Σ
XIX,,xX)(rr -0(s))
Jor (2w) (20) K2 Zn)
2 1)
2
jzk
(1)))(())
XW),,XN(RRn) (0 (2nn))K2222 K2(2,,210 2n)
( 1)
j#k (1))
((11) 2111
31)
For 7211 since the conditional expectations of the summand given Sj or S5 are zero, using Assumption 4, we can show 7,10) = =Op(n (11-a-1) 1=-a)) (see the proof for I, under the null). For 7,301) since the conditional expectations of the summand given (R.j, Sz) or S& are zero, 7(31) Op(n-1-0) follows in the same way T111) (using part (b) of Assumption 4 and part (a) of Assumption 13). For f(21), we only have
(000 s.]
E
E [ [2] (z))
0,
which implies
Var(Z221)) S (1-"s)))((((]]))"">>>>>>>>>>>>>>
(x Kz [3]]
n C
(z))
n 122
By part (b) of Assumption 4 and part (a) of Assumption 13, 7,221 = (n -112-00000) and we can conclude that T.
O (n 1/2-a--02 + n --)). The proof for 7,2) Op(n 1/2- 02 + n 11 is similar 7()) and the proof for
7(3) -10 is similar to I under the null. So we have ī₂ = Op(n 1/2 1/2- a- + n -1-22) under H
Now consider the term T3. Denote k) Ef (X, Xz) Zz } and W2(j, E BY, Yz) Zz}. then we
28

---

Testing Conditional Mean Independence using GNN
have
Is
Σ
3 (
iyk
1
+
Σ
[i (G
2
J.KE((1)
Σ
W
2(2 1)
(,NE(())
= 11 + 1/
7(3
+ 7501
(24)
We now show that 2 130 Op(n 1/2- (0140))))))) and 0
under both Hg and Hin- Note that
E E < WE (1, 2) W2 (1, 2) (21,, 21, Z2) and
(1(1,) = | x(11 Kx(, x, Kx
- x()) ) Z1, 2z. (,, x(11 )x
S E Kxx (,, (1)) Ex (s, x |1 |E Kxx (s, x(11) BE (,, 82 ||ax
Similarly, [W((1,)) S ([gr - (z 12llgr (21) - 3r (Z))l2,.. so by part (b) of Assumption 4, we have E 1
Onn +oo))) For T have ETT") 0 and
Var
n E (1, 1W
W1 |4h.::/0066100000333440011111113
By part (b) of assumption 4 and equation (25), E W1 ",, O(n -((0)
Since
-)) = = [2]] (2)) (z))))9rz
[| 0r -0))--((224 (Z1)
(25)
Y2) |,,-::/|)|11,b00000--0000,,,,1,1 34123//4444433111a1/////33111113 (9))))r())])
/r (2))
by part (a) and (b) of Assumption 4, we have
X)) 33 |1,.://.0.__..11.....0.2222..,.0 O(n 00n00(nnnn)))) -
which implies T:2) Op(n 1o:2n-1/2-aan) (1/2-(a and 7:3) 01 follows similarly. For 7,0) since E7
29

---

Testing Conditional Mean Independence using GNN
under Ho, we can show / = Op(n - 1/2-(cett))) in the same way as 7 Under Hin:
(
2 2
JHk J(FE)(1)
[,a 2x)
2
((kE)()
-1g11) 2(21)
7
(x,))=)) 7x)
( 1)
(1kE
-.. "(0n)))))) 9x K2 Zx)
2(7
-
(111)
((12))
3
(26)
We now show 11)
Op(n 1/2 1/2- (o (same result holds for T/21) and the proof is omitted). Note that
B17 (
n
1E ;) B (x2 Ilux V (12(21,2)K22(2,
1 (2)) (2z) 1(10(5))1b100
Kix (XP) lllax
+ no
|[, - (x ,llx (122(1,,)))2(2222 22)
=0nn =0nn (ataataa))
by Jensen's inequality and part (b) of Assumptions 4 and 13, so we have 71111) Op(n +a+aat()t)) Since
E 0. following similar approach as for F5²) we can show 7(112) = On -/211// and conclude that
P(11) = (n +n-(ata+tar) So part(b) is proved.
30


Reviews:
[{"review_id": "raeuxBez2T", "reviewer": "Reviewer_KpXH", "summary": "This develops a novel method to test for conditional mean independence that works well in high-dimensions, gives asymptotic size control, and has nontrivial power against local alternatives. This depends on using deep learning to learn g_y and g_x using a bootstrap sample. They then use the test statistic in the unnumbered equation before equation 2. They then generate data under a null distribution and use that to approximate the p-value. They show theoretically  They then test the performance on synthetic data with regression examples along with real examples testing whether masking affects prediction accuracy.", "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 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.\n\n**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).\n\nThe results were consistent with those in Table 2 of the manuscript:\nempirical sizes: 5.8\\% and 6\\% (close to the original);\nsize-adjusted power (sparse alternative): 98.6\\% and 99.8\\%;\nsize-adjusted power (dense alternative): 100\\% in both cases.\nThis suggests that the test’s performance is robust to moderate changes in key network parameters.\n\nFor 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.\n    \nRegarding 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. \n\n- **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.\n- **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.\n\n**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.\n\n\n**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.\n|                                | M=5 B=200 | M=5 B=500 | M=20 B=200 | M=20 B=500 | M=60 B=200 | M=60 B=500 |\n|-----------------------|-----------|-----------|-------------|------------|-------------|---------------------|\n| Empirical Size                 | 6.2       | 6.6       | 5.8        | 6.2        | 6.8        | 6.4        |\n| Power under Sparse Alternative | 99.2      | 99        | 99.6       | 99.6       | 98.6       | 99         |\n| Power under Dense Alternative  | 100       | 100       | 100        | 100        | 100        | 100        |\n\n**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.", "claims_and_evidence": "The authors justify the claims that this method works well in testing CMI in high dimensions, both theoretically and empirically. In my opinion their evidence is sufficient but some of the comparisons are disappointing. The authors don't really explore the quality of neural network approximations, but given that proving theoretical guarantees even for more rudimentary tasks is difficult I don't fault them for that. But I would have preferred some material in the appendix seeing how some of the neural network parameters changing affects results. Also would be nice to have a sense of how long these methods take to run.", "methods_and_evaluation_criteria": "Hypothesis testing is straightforward for evaluation criteria. Given that this is designed for high dimensions, high dimensional regression and image questions make sense. One of the problems I do have with this paper is that none of the competitor methods are used on the image data from what I can tell. I would have assumed that it would be a straightforward improvement.", "theoretical_claims": "I checked the theoretical claims to the best of my ability.", "experimental_designs_or_analyses": "I checked the analyses. I suppose this is the most appropriate category to point out that from what I can tell, the paper doesn't analyze how the Monte Carlo method performs along with the number of bootstrap samples. This is important, because one possibility is that this is just a very high-powered poorly sized test but the errors in the approximations introduced address that.", "relation_to_broader_scientific_literature": "This paper improves on extensive literature on conditional mean independence. Many kernel methods struggle with high dimensional data. While most tests have size guarantees, at least of the most pertinent methods also struggle at maintaining power at a parametric rate, or at least being backed with theoretical guarantees.", "essential_references_not_discussed": "The references cited seem quite extensive.", "other_strengths_and_weaknesses": "None that are not previously described", "other_comments_or_suggestions": "Please don't put assumptions in the appendix. It's annoying to have to look there. Also it's weird to have assumption 7 and 9 but not 1-6 and 8.", "questions_for_authors": "Unless I missed it why is the effect of choice of B not evaluated?\n\nAre there direct comparisons on the image data?", "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": "T800PmP2KY", "reviewer": "Reviewer_tmTe", "summary": "This paper proposes a novel statistical method to conditional mean independence (CMI) testing. First, the authors introduce a new population-level CMI measure and develop a bootstrap-based hypothesis testing framework that employs generative neural networks to approximate conditional mean functions. Its test statistic is constructed to reduce the influence of nonparametric estimation errors, ensuring asymptotic precision. The proposed method performs well in high-dimensional settings, and supports multivariate responses. Finally, the experiments on simulated and real-world data demonstrate the effectiveness of the proposed method.\n\n##  update after rebuttal\n After reviewing the rebuttal addressed to me and those for other reviewers, I am willing to maintain my score.", "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 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. \n\n**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.  \nThe 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.\n\n**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.\n\nPart (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.\nRegarding 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).\nFor 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.\nAs suggested, we will include a brief discussion of these assumptions, particularly the more technical ones, in the revised version of our paper.\n\n**Reference**\n\nStone, C. J. (1982): \"Optimal global rates of convergence for nonparametric regression.\" Ann. Statist. \n\nBauer and Kohler (2019): \"On deep learning as a remedy for the curse of dimensionality in nonparametric regression.\" Ann. Statist. \n\nBerrett et al. (2019): \"The Conditional Permutation Test for Independence While Controlling for Confounders.\" Journal of the Royal Statistical Society Series B: Statistical Methodology", "claims_and_evidence": "Yes, the claims made in the submission appear to be supported by clear and convincing evidence.", "methods_and_evaluation_criteria": "Yes, the proposed hypothesis testing framework makes sense for the CMI testing  task.", "theoretical_claims": "Yes. I check some proofs,  including Theorems 4，5，6.", "experimental_designs_or_analyses": "Yes, the experimental designs and analyses are sound.", "relation_to_broader_scientific_literature": "This paper proposes a novel framework for CMI testing problem with strong theoretical guarantee.", "essential_references_not_discussed": "No, the paper includes all essential and relevant references.", "other_strengths_and_weaknesses": "Strengths:\n- The method is shown to control Type I error while maintaining nontrivial  power.\n\n- The proofs are clear and sound.\n\n- Comparisons against existing CMI tests highlight superior empirical  performance of the proposed method.\n\nWeaknesses:\n- The method requires to train multiple deep neural networks, which increases computational cost. \n\n- The theoretical results depends on some strong assumptions, such as assumption 7. The rationality of the assumptions in this paper should be discussed.", "other_comments_or_suggestions": "See Weaknesses.", "questions_for_authors": "No.", "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": "ZMvk74RJ8I", "reviewer": "Reviewer_WEw2", "summary": "The paper introduces a new nonparametric test for conditional mean independence (CMI) that leverages deep generative neural networks to estimate conditional mean embeddings. The proposed method uses a novel population measure based on RKHS embeddings and constructs a test statistic in a multiplicative form that is robust to the slower convergence rates of nonparametric nuisance parameter estimators. To mitigate estimation errors, the authors combine sample splitting and cross-fitting with a generative moment matching network (GMMN) for sampling from the conditional distribution of covariates. The paper provides comprehensive theoretical guarantees (including asymptotic size control and power properties under local alternatives) and supports its claims via extensive simulation studies and real-world imaging applications (facial expression recognition and age estimation).", "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 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. \n\n**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. \n\nFollowing 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.\n    \n**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.\n\n**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. \n \nIn light of this comment, we will include a discussion on computational complexity in the revised version.\n\n**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.\n\n**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.\n\nFortunately, 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.", "claims_and_evidence": "The central claims are well-supported by theoretical analysis and empirical results. In particular:\n\n- The theoretical claims about precise asymptotic size control are supported by rigorous proofs and verified in simulation studies showing empirical sizes close to nominal levels. The claim of detecting local alternatives is validated through theoretical analysis in Theorem 5 and empirical power evaluations in simulations. However, the theoretical results rely on certain technical assumptions (e.g., on the decay rates of estimation errors) might limit the generality of the results in practice.\n\n- The claim of strong empirical performance in high-dimensional settings is demonstrated through comprehensive simulations against multiple baseline methods and by experiment on real-world imaging applications.", "methods_and_evaluation_criteria": "This work proposes to use generative models to approximate conditional distributions for RKHS-based testing, which is novel to overcome challenges in high-dimensional nonparametric estimation. The evaluations on both synthetic experiments (with clear sparse and dense alternatives) and applications on real imaging datasets are appropriate and provide convincing evidence of the method’s effectiveness.", "theoretical_claims": "The paper contains several theoretical contributions, with proofs detailed in the supplementary material. I reviewed the main steps in the proofs of Theorems 4, 5, and 6. Under the given assumptions (e.g., Assumptions 7 and 9), the arguments appear sound. However, some of those technical assumptions seem to be quite strong. Clarification on the practical implications of these assumptions would be beneficial.", "experimental_designs_or_analyses": "The simulation studies (Examples A1 and A2) are well-designed and cover a range of scenarios (both linear and nonlinear models, and sparse versus dense alternatives). The experimental analyses also include a comparison with multiple state-of-the-art methods.\n\nHowever, while the results demonstrate clear benefits of the proposed test in terms of both size control and power, I believe a more detailed ablation study—particularly regarding the sensitivity to hyperparameter choices and kernel bandwidth selection—could strengthen the empirical section further.", "relation_to_broader_scientific_literature": "I believe the authors have done good work on literature review within the CMI testing literature, clearly identifying limitations of existing methods.", "essential_references_not_discussed": "While the paper cites a wide range of relevant literature, it would benefit from discussing recent advances in generative modeling (GANs and diffusion models) for conditional distribution estimation (i.e. trying different design of the generative network $\\hat G$), for example:\n\n1. Athey, S., Imbens, G. W., Metzger, J., and Munro, E. Using Wasserstein generative adversarial networks for the design of Monte Carlo simulations. Journal of Econometrics, 2021.\n\n2. Baptista, R., Hosseini, B., Kovachki, N. B., and Marzouk, Y. Conditional sampling with monotone GANs: From generative models to likelihood-free inference. arXiv preprint arXiv:2006.06755, 2020.\n\n3. Shi, Y., De Bortoli, V., Deligiannidis, G., and Doucet, A. Conditional simulation using diffusion Schrödinger bridges. In Uncertainty in Artificial Intelligence, pp. 1792–1802. PMLR, 2022.\n\n4. Nguyen, B., Nguyen, B., Nguyen, H. T., & Nguyen, V. A. Generative conditional distributions by neural (entropic) optimal transport. In Proceedings of the 41st International Conference on Machine Learning (pp. 37761-37775), 2024.", "other_strengths_and_weaknesses": "Good: \n\n- The paper is in general clear to follow, though can be dense with notation in the first few pages.\n\n\nNeed to address: \n\n- There is limited discussion of hyperparameter sensitivity (kernel bandwidths, network architectures)", "other_comments_or_suggestions": "n/a", "questions_for_authors": "1. How does computational complexity scale with the dimensionality of X, Y, Z, and sample size n compared to existing methods?\n2. Have you explored alternative conditional generative models beyond GMMNs (e.g., conditional GANs, score-based diffusion model, see Essential References Not Discussed)? How might they affect test performance?\n3. How robust is your method to model misspecification when estimating conditional mean functions in highly nonlinear relationships?", "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 introduces a new nonparametric test for conditional mean independence (CMI), that works well in high-dimensional settings, by leveraging deep generative neural networks to estimate the conditional mean.  A new population-level CMI measure is proposed and its sample version serves as the test statistic. To mitigate the estimation errors a combination of sample splitting and cross-fitting with a generative moment matching network (GMMN) for sampling from the conditional distribution of covariates. The theory provides theoretical guarantees for both the asymptotic size control and power under local alternatives. Numerical experiments support the theory.

One concern is the lack of comparisons on image data with competing methods. The authors addressed this in their response.  

Overall, this work makes a significant contribution to the literature of testing CMI.
Decision: Accept (poster)