Semiparametric Efficiency of Residual Correlation Testing under Gaussian Additive Noise Models

Consider the Gaussian additive noise model (GANM), where $X,Y,{\bf Z}$ satisfy the ANM in (2), where the regression errors $(\epsilon_{x},\epsilon_{y})^{\top}\sim N({\bf 0},{\bm{\Sigma}})$ where

$\displaystyle{\bm{\Sigma}}=\left(\begin{matrix}\sigma_{x}^{2}&\rho\sigma_{x}\sigma_{y}\\ \rho\sigma_{x}\sigma_{y}&\sigma_{y}^{2}\end{matrix}\right).$

Here, $\sigma_{x}^{2}$ and $\sigma_{y}^{2}$ are the variances of $\epsilon_{x}$ and $\epsilon_{y}$, respectively, and $\rho$ is the correlation coefficient between $\epsilon_{x}$ and $\epsilon_{y}$, i.e., $\rho=\mbox{corr}(\epsilon_{x},\epsilon_{y})$.

Clearly, under GANM, since $X$ and $Y$ depend on ${\bf Z}$ only through the conditional mean function $m_{x}({\bf Z})$ and $m_{y}({\bf Z})$, we know that $X\;\,\rule[0.0pt]{0.29999pt}{6.69998pt}\hskip-2.5pt\rule[0.0pt]{6.49994pt}{0.29999pt}\hskip-2.5pt\rule[0.0pt]{0.29999pt}{6.69998pt}\;\,Y|{\bf Z}$ if and only if the regression errors are independent, i.e., $\epsilon_{x}\;\,\rule[0.0pt]{0.29999pt}{6.69998pt}\hskip-2.5pt\rule[0.0pt]{6.49994pt}{0.29999pt}\hskip-2.5pt\rule[0.0pt]{0.29999pt}{6.69998pt}\;\,\epsilon_{y}$, which, under the joint Gaussian assumption for $(\epsilon_{x},\epsilon_{y})$, is further equivalent to $\mbox{corr}(\epsilon_{x},\epsilon_{y})=0$, i.e., $\rho=0$. Therefore, to test whether (1) holds, we can equivalently test

$\displaystyle H_{0}:\rho=0\quad\text{vs.}\quad H_{1}:\rho\neq 0.$

(3)

In this way, under GANM, testing conditional independence of $X$ and $Y$ given ${\bf Z}$ is equivalent to testing the uncorrelatedness of the regression errors $\epsilon_{x}$ and $\epsilon_{y}$.

We first consider the oracle case when the regression functions $m_{x}$ and $m_{y}$ are known. In this case, we could directly calculate the true regression errors

$\displaystyle\epsilon_{xi}=X_{i}-m_{x}({\bf Z}_{i}),\quad\epsilon_{yi}=Y_{i}-m_{y}({\bf Z}_{i}),\qquad i=1,\dots,n,$

Based on $(\epsilon_{x1},\epsilon_{y1}),\dots,(\epsilon_{xn},\epsilon_{yn})$, we then estimate $\rho$ by the Pearson correlation coefficient

$\displaystyle\widehat{\rho}_{\rm{orc}}=\frac{\widehat{\sigma}_{xy,\rm{orc}}}{\widehat{\sigma}_{x,\rm{orc}}\widehat{\sigma}_{y,\rm{orc}}}=\frac{\sum_{i=1}^{n}\epsilon_{xi}\epsilon_{yi}}{(\sum_{i=1}^{n}\epsilon_{xi}^{2})^{1/2}(\sum_{i=1}^{n}\epsilon_{yi}^{2})^{1/2}}.$

Here,

$\displaystyle\widehat{\sigma}_{xy,\rm{orc}}=n^{-1}\sum_{i=1}^{n}\epsilon_{xi}\epsilon_{yi},\quad\widehat{\sigma}_{x,\rm{orc}}^{2}=n^{-1}\sum_{i=1}^{n}\epsilon_{xi}^{2},\quad\widehat{\sigma}_{y,\rm{orc}}^{2}=n^{-1}\sum_{i=1}^{n}\epsilon_{yi}^{2}$

are estimators of $\sigma_{x}^{2}=\mbox{var}(\epsilon_{x})$, $\sigma_{y}^{2}=\mbox{var}(\epsilon_{y})$, and $\sigma_{xy}=\mbox{cov}(\epsilon_{x},\epsilon_{y})$, respectively. Since $(\epsilon_{xi},\epsilon_{yi})^{\top}$ are i.i.d. bivariate Gaussian random vectors, by Theorem 5.1.6 of 23, the asymptotic distribution of $\widehat{\rho}_{\rm{orc}}$ is

$\displaystyle n^{1/2}(\widehat{\rho}_{\rm{orc}}-\rho)\overset{d}{\to}N\{0,(1-\rho^{2})^{2}\}.$

(4)

In particular, under the null hypothesis $H_{0}:\rho=0$, the asymptotic null distribution of $\widehat{\rho}$ is

$\displaystyle n^{1/2}\widehat{\rho}_{\rm{orc}}\overset{d}{\to}N(0,1).$

In practice, the regression functions $m_{x}$ and $m_{y}$ are unknown, and we define $\widehat{m}_{x}$ and $\widehat{m}_{y}$ as nonparametric estimators of $m_{x}$ and $m_{y}$, respectively. Based on $\widehat{m}_{x}$ and $\widehat{m}_{y}$, We calculate the fitted residuals as

$\displaystyle\widehat{\epsilon}_{xi}=X_{i}-\widehat{m}_{x}({\bf Z}_{i}),\quad\widehat{\epsilon}_{yi}=Y_{i}-\widehat{m}_{y}({\bf Z}_{i}),\qquad i=1,\dots,n,$

(5)

and the residual correlation estimator is given by the Pearson correlation coefficient based on $(\widehat{\epsilon}_{x1},\widehat{\epsilon}_{y1}),\dots,(\widehat{\epsilon}_{xn},\widehat{\epsilon}_{yn})$ is

$\displaystyle\widehat{\rho}=\frac{\widehat{\sigma}_{xy}}{\widehat{\sigma}_{x}\widehat{\sigma}_{y}}=\frac{\sum_{i=1}^{n}\widehat{\epsilon}_{xi}\widehat{\epsilon}_{yi}}{(\sum_{i=1}^{n}\widehat{\epsilon}_{xi}^{2})^{1/2}(\sum_{i=1}^{n}\widehat{\epsilon}_{yi}^{2})^{1/2}}.$

(6)

Here,

$\displaystyle\widehat{\sigma}_{xy}=n^{-1}\sum_{i=1}^{n}\widehat{\epsilon}_{xi}\widehat{\epsilon}_{yi},\quad\widehat{\sigma}_{x}^{2}=n^{-1}\sum_{i=1}^{n}\widehat{\epsilon}_{xi}^{2},\quad\widehat{\sigma}_{y}^{2}=n^{-1}\sum_{i=1}^{n}\widehat{\epsilon}_{yi}^{2}$

(7)

are corresponding estimators of $\sigma_{xy}$, $\sigma_{x}^{2}$ and $\sigma_{y}^{2}$ when $m_{x}$ and $m_{y}$ are estimated.

To make the asymptotic theory applicable to a broad class of machine learning methods for estimating $m_{x}$ and $m_{y}$, without requiring explicit asymptotic expansions of the nuisance regression estimators, we employ sample splitting and cross-fitting. Specifically, one part of the data is used to estimate the regression functions, while the other part is used to construct the residual correlation estimator $\widehat{\rho}$. The two resulting estimators are then averaged to obtain the final estimator. More details are given in Section 3.5.

Since $(\epsilon_{x},\epsilon_{y})\;\,\rule[0.0pt]{0.29999pt}{6.69998pt}\hskip-2.5pt\rule[0.0pt]{6.49994pt}{0.29999pt}\hskip-2.5pt\rule[0.0pt]{0.29999pt}{6.69998pt}\;\,{\bf Z}$, the conditional density of $X,Y|{\bf Z}$ can be written in terms of the joint density of the regression errors $(\epsilon_{x},\epsilon_{y})$. That is,

$\displaystyle f_{X,Y|{\bf Z}}(x,y,{\bf z})=f_{\epsilon_{x},\epsilon_{y}}\{x-m_{x}({\bf z}),y-m_{y}({\bf z})\}=f_{\epsilon_{x},\epsilon_{y}}(\epsilon_{x},\epsilon_{y}),$

where

$\displaystyle\epsilon_{x}=x-m_{x}({\bf z}),\quad\epsilon_{y}=y-m_{y}({\bf z}).$

(8)

Then, for a realization $(x,y,{\bf z})$, the likelihood function can be written as

$\displaystyle f(x,y,{\bf z})=f_{X,Y|{\bf Z}}(x,y,{\bf z})f_{\bf Z}({\bf z})=f_{\epsilon_{x},\epsilon_{y}}(\epsilon_{x},\epsilon_{y})f_{\bf Z}({\bf z}),$

where $f_{\bf Z}$ is the marginal distribution of ${\bf Z}$. Furthermore, since $(\epsilon_{x},\epsilon_{y})^{\top}\sim N({\bf 0},{\bm{\Sigma}})$, we plug in the Gaussian density function to get

$\displaystyle f(x,y,{\bf z})=\frac{1}{2\pi\sigma_{x}\sigma_{y}\sqrt{1-\rho^{2}}}\exp\left\{-\frac{1}{2(1-\rho^{2})}\left(\frac{\epsilon_{x}^{2}}{\sigma_{x}^{2}}-2\rho\frac{\epsilon_{x}\epsilon_{y}}{\sigma_{x}\sigma_{y}}+\frac{\epsilon_{y}^{2}}{\sigma_{y}^{2}}\right)\right\}f_{\bf Z}({\bf z}),$

(9)

where $\epsilon_{x}$ and $\epsilon_{y}$ are given by (8).

Note that, in the likelihood, our parameter of interest is $\rho$, and all other parameters, including $\eta=\{\sigma_{x},\sigma_{y},m_{x}(\cdot),m_{y}(\cdot),f_{\bf Z}(\cdot)\}$, can be viewed as nuisance parameters. As we see, in the nuisance parameters, $\sigma_{x}$ and $\sigma_{y}$ are parametric, and the rest ones $m_{x}(\cdot)$, $m_{y}(\cdot)$ and $f_{\bf Z}(\cdot)$ are nonparametric.

The log likelihood is given by

	$\displaystyle\mbox{log}f(x,y,{\bf z})$	$\displaystyle=$	$\displaystyle\mbox{log}(2\pi)-\mbox{log}(\sigma_{x})-\mbox{log}(\sigma_{y})-\frac{1}{2}\mbox{log}(1-\rho^{2})$
			$\displaystyle-\frac{1}{2(1-\rho^{2})}\left(\frac{\epsilon_{x}^{2}}{\sigma_{x}^{2}}-2\rho\frac{\epsilon_{x}\epsilon_{y}}{\sigma_{x}\sigma_{y}}+\frac{\epsilon_{y}^{2}}{\sigma_{y}^{2}}\right)+\mbox{log}\{f_{\bf Z}({\bf z})\}.$

To derive the efficient score for $\rho$, we first calculate the nuisance tangent space associated with each of the the nuisance parameter in $\eta=\{\sigma_{x},\sigma_{y},m_{x}(\cdot),m_{y}(\cdot),f_{\bf Z}(\cdot)\}$. We then derive its orthogonal decomposition and its orthogonal complement, on which we will need to project the score function with respect to $\rho$ to find the efficient score (3; 37).

Denote ${\cal H}=\{f(x,y,{\bf z}):E\{f(X,Y,{\bf Z})\}=0,\mbox{var}\{f(X,Y,{\bf Z})\}<\infty\}$ the Hilbert space of all possible influence functions. We first give the nuisance tangent space as the next proposition.

Based on the nuisance tangent space in Proposition 1, we derive the orthogonal complement of the nuisance tangent space as the following proposition.

In the next theorem, we give the efficient score function, which is defined as the projection of the score function with respect to $\rho$ onto $\Lambda^{\perp}$, which is given in Proposition 2. The efficient score is written as $S_{\rm eff}=\Pi(S_{\rho}|\Lambda^{\perp})$, where $\Pi$ is the projection operator. Based on the efficient score, we can further derive the efficient Fisher information and the efficient influence function as the next theorem.

As shown in (10), the efficient score function has a simple quadratic form involving only the standardized regression errors, and the corresponding efficient estimating equation leads to a simple explicit estimator in Section 3.4.

Based on semiparametric theory, the efficient estimator can be obtained through implementing

$\displaystyle\sum_{i=1}^{n}S_{\rm eff}(X_{i},Y_{i},{\bf Z}_{i})=0.$

Plugging in the efficient score function (10), we have

$\displaystyle\sum_{i=1}^{n}\left(\rho\frac{\epsilon_{xi}^{2}}{\sigma_{x}^{2}}-2\frac{\epsilon_{xi}\epsilon_{yi}}{\sigma_{x}\sigma_{y}}+\rho\frac{\epsilon_{yi}^{2}}{\sigma_{y}^{2}}\right)=0.$

(12)

The solution to the estimating equation (12) can be explicitly written as

$\displaystyle\widehat{\rho}=\frac{2\sigma_{x}\sigma_{y}\sum_{i=1}^{n}\epsilon_{xi}\epsilon_{yi}}{\sigma_{y}^{2}\sum_{i=1}^{n}\epsilon_{xi}^{2}+\sigma_{x}^{2}\sum_{i=1}^{n}\epsilon_{yi}^{2}}.$

(13)

In practice, all nuisance parameters need to be estimated. Firstly, $m_{x}$ and $m_{y}$ are estimated by $\widehat{m}_{x}$ and $\widehat{m}_{y}$. Then, based on $\widehat{m}_{x}$ and $\widehat{m}_{y}$, we need to calculate the residuals $\widehat{\epsilon}_{xi}$ and $\widehat{\epsilon}_{yi}$ in (5) as the estimated versions of regression errors $\epsilon_{xi}$ and $\epsilon_{yi}$. Next, based on the residuals $\widehat{\epsilon}_{xi}$ and $\widehat{\epsilon}_{yi}$, we can estimate the unknown parameters $\sigma_{x}$ and $\sigma_{y}$ by $\widehat{\sigma}_{x}$ and $\widehat{\sigma}_{y}$ as defined in (7).

Surprisingly, after replacing all nuisance parameters in (13) by their empirical estimators, we obtain

$\displaystyle\widehat{\rho}=\frac{2(n^{-1}\sum_{i=1}^{n}\widehat{\epsilon}_{xi}^{2})^{1/2}(n^{-1}\sum_{i=1}^{n}\widehat{\epsilon}_{yi}^{2})^{1/2}\sum_{i=1}^{n}\widehat{\epsilon}_{xi}\widehat{\epsilon}_{yi}}{n^{-1}\sum_{i=1}^{n}\widehat{\epsilon}_{yi}^{2}\sum_{i=1}^{n}\widehat{\epsilon}_{xi}^{2}+n^{-1}\sum_{i=1}^{n}\widehat{\epsilon}_{xi}^{2}\sum_{i=1}^{n}\widehat{\epsilon}_{yi}^{2}}=\frac{\sum_{i=1}^{n}\widehat{\epsilon}_{xi}\widehat{\epsilon}_{yi}}{(\sum_{i=1}^{n}\widehat{\epsilon}_{xi}^{2})^{1/2}(\sum_{i=1}^{n}\widehat{\epsilon}_{yi}^{2})^{1/2}},$

which coincides with residual correlation estimator as in (6). Therefore, the residual correlation estimator in (6) can be equivalently interpreted as the estimator induced by the efficient influence function under GANM.

Note that the nuisance regression functions $m_{x}$ and $m_{y}$ need to be estimated nonparametrically from the data. To facilitate the asymptotic analysis, we employ sample splitting and cross-fitting as introduced in Section 2.2. This construction separates the nuisance estimation step from the estimation step for $\rho$, which avoids requiring explicit asymptotic expansions of the nonparametric regression estimators and allows the asymptotic theory to rely primarily on suitable convergence rate conditions.

Specifically, we use the first $n_{1}$ data to estimate $m_{x}$ and $m_{y}$, denote the corresponding estimators by $\widehat{m}_{x1}$ and $\widehat{m}_{y1}$, and use the remaining $n_{2}$ data to estimate $\sigma_{x}^{2}$, $\sigma_{y}^{2}$, $\sigma_{xy}$ and $\rho$. That is,

$\displaystyle\widehat{\sigma}_{xy2}=n_{2}^{-1}\sum_{i=n_{1}+1}^{n}\widehat{\epsilon}_{xi1}\widehat{\epsilon}_{yi1},\quad\widehat{\sigma}_{x2}^{2}=n_{2}^{-1}\sum_{i=n_{1}+1}^{n}\widehat{\epsilon}_{xi1}^{2},\quad\widehat{\sigma}_{y2}^{2}=n_{2}^{-1}\sum_{i=n_{1}+1}^{n}\widehat{\epsilon}_{yi1}^{2},$

where

$\displaystyle\widehat{\epsilon}_{xi1}=X_{i}-\widehat{m}_{x1}({\bf Z}_{i}),\quad\widehat{\epsilon}_{yi1}=Y_{i}-\widehat{m}_{y1}({\bf Z}_{i}),\qquad i=n_{1}+1,\dots,n.$

and the final estimator of $\rho$ is

$\displaystyle\widehat{\rho}_{2}=\frac{\widehat{\sigma}_{xy2}}{\widehat{\sigma}_{x2}\widehat{\sigma}_{y2}}=\frac{\sum_{i=n_{1}+1}^{n}\widehat{\epsilon}_{xi1}\widehat{\epsilon}_{yi1}}{(\sum_{i=n_{1}+1}^{n}\widehat{\epsilon}_{xi1}^{2})^{1/2}(\sum_{i=n_{1}+1}^{n}\widehat{\epsilon}_{yi1}^{2})^{1/2}}.$

(14)

We can also switch the roles of the two parts of data, using the later $n_{2}$ data to estimate $m_{x}$ and $m_{y}$, which are denoted similarly by $\widehat{m}_{x2}$ and $\widehat{m}_{y2}$. We then use the first $n_{1}$ data to estimate $\sigma_{x}^{2}$, $\sigma_{y}^{2}$ and $\sigma_{xy}$, which are denoted similarly by $\widehat{\sigma}_{x1}^{2}$, $\widehat{\sigma}_{y1}^{2}$ and $\widehat{\sigma}_{xy1}$. Then we construct the estimator $\widehat{\rho}_{1}$ as

$\displaystyle\widehat{\rho}_{1}=\frac{\widehat{\sigma}_{xy1}}{\widehat{\sigma}_{x1}\widehat{\sigma}_{y1}}=\frac{\sum_{i=1}^{n_{1}}\widehat{\epsilon}_{xi2}\widehat{\epsilon}_{yi2}}{(\sum_{i=1}^{n_{1}}\widehat{\epsilon}_{xi2}^{2})^{1/2}(\sum_{i=1}^{n_{1}}\widehat{\epsilon}_{yi2}^{2})^{1/2}},$

where

$\displaystyle\widehat{\sigma}_{xy1}=n_{1}^{-1}\sum_{i=1}^{n_{1}}\widehat{\epsilon}_{xi2}\widehat{\epsilon}_{yi2},\quad\widehat{\sigma}_{x1}^{2}=n_{1}^{-1}\sum_{i=1}^{n_{1}}\widehat{\epsilon}_{xi2}^{2},\quad\widehat{\sigma}_{y1}^{2}=n_{1}^{-1}\sum_{i=1}^{n_{1}}\widehat{\epsilon}_{yi2}^{2},$

with

$\displaystyle\widehat{\epsilon}_{xi2}=X_{i}-\widehat{m}_{x2}({\bf Z}_{i}),\quad\widehat{\epsilon}_{yi2}=Y_{i}-\widehat{m}_{y2}({\bf Z}_{i}),\qquad i=1,\dots,n_{1}.$

Taking $n_{1}=n_{2}=n/2$ for simplicity, we can average the two estimators $\widehat{\rho}_{1}$ and $\widehat{\rho}_{2}$ to construct the final estimator

$\displaystyle\widehat{\rho}=(\widehat{\rho}_{1}+\widehat{\rho}_{2})/2.$

(15)

In the next theorem, we give the asymptotic properties of $\widehat{\rho}_{2}$. The analogous properties for $\widehat{\rho}_{1}$ can be similarly derived.

Note that (17) indicates that the estimator $\widehat{\rho}$ is indeed efficient under the convergence rate assumptions in (16).

Note that the variance of $\phi_{\rm eff}(X,Y,{\bf Z})$ is

$\displaystyle E\left\{\phi_{\rm eff}^{2}(X,Y,{\bf Z})\right\}=E(S_{\rm eff}^{2})^{-1}=(1-\rho^{2})^{2}.$

The next theorem gives the the asymptotic distribution of $\widehat{\rho}$.

The asymptotic variance in (18) coincides with that in the oracle result (4), indicating that the convergence rate conditions in (16) are sufficient to make the nonparametric regression estimation errors negligible in terms of first-order asymptotics.

Under the null hypothesis $\rho=0$, the asymptotic distribution of $\widehat{\rho}$ becomes

$\displaystyle n^{1/2}\widehat{\rho}\xrightarrow{d}N(0,1).$

Therefore, we can construct a Wald test for (3) (19). At significance level $\alpha$, our decision rule is to reject $H_{0}$ if $n^{1/2}|\widehat{\rho}|>z_{1-\alpha/2}$, where $z_{1-\alpha/2}$ denotes the $(1-\alpha/2)$-quantile of the standard normal distribution. Equivalently, the corresponding p-value is given by $2\Phi(-n^{1/2}|\widehat{\rho}|)$, where $\Phi$ is the cumulative distribution function of the standard normal distribution.

In general, to construct a confidence interval for $\rho$, since the asymptotic variance $(1-\rho^{2})^{2}$ in (18) involves the unknown parameter $\rho$, we use the plug-in estimator as $(1-\widehat{\rho}^{2})^{2}$. Therefore, by Slutsky’s theorem, an asymptotic $100(1-\alpha)\%$ confidence interval for $\rho$ is given by

$\displaystyle\widehat{\rho}\pm n^{-1/2}z_{1-\alpha/2}(1-\widehat{\rho}^{2}).$

In the simulation studies, we consider two models:

	Model 1:		$\displaystyle{\bf Z}=(Z_{1},\dots,Z_{6})^{\top},\quad\text{where}\ Z_{1},\dots,Z_{6}\overset{\mathrm{iid}}{\sim}\mathrm{Uniform}(0,1),$
			$\displaystyle X=0.4Z_{1}+0.6\sqrt{Z_{3}}+\epsilon_{x},\quad Y=0.6Z_{4}+0.4\sqrt{Z_{5}}+\epsilon_{y},$
	Model 2:		$\displaystyle{\bf Z}=(Z_{1},\dots,Z_{10})^{\top},\quad\text{where}\ Z_{1},\dots,Z_{10}\overset{\mathrm{iid}}{\sim}\mathrm{Uniform}(0,1),$
			$\displaystyle X=(Z_{1}+\dots+Z_{10})/10+\epsilon_{x},\quad Y=\mbox{log}(1+Z_{3}+Z_{6})+\epsilon_{y},$

where in both models, the regression errors $(\epsilon_{x},\epsilon_{y})^{\top}\;\,\rule[0.0pt]{0.29999pt}{6.69998pt}\hskip-2.5pt\rule[0.0pt]{6.49994pt}{0.29999pt}\hskip-2.5pt\rule[0.0pt]{0.29999pt}{6.69998pt}\;\,{\bf Z}$ and $(\epsilon_{x},\epsilon_{y})^{\top}\sim N({\bf 0},{\bm{\Sigma}})$ with

$\displaystyle{\bm{\Sigma}}=\left(\begin{matrix}\sigma^{2}&\rho\sigma^{2}\\ \rho\sigma^{2}&\sigma^{2}\end{matrix}\right).$

We set $\sigma^{2}=0.2$. We take sample sizes to be $n=100,200,500,1000,2000,5000$. Under the null hypothesis, we simply set $\rho=0$. Under the alternative hypothesis, we consider two sets of $\rho$: (1) $\rho=\pm 0.25,\pm 0.5,\pm 0.75,\pm 1$; (2) $\rho=\pm 0.025,\pm 0.05,\pm 0.075,\pm 0.1$. The former set of $\rho$ indicates a strong dependence signal, which represents deviation from the null hypothesis; the later indicates a weak dependence signal which is hard to detect. For each model, we conduct 500 independent experiments, and the significance level is set as $\alpha=0.05$.

In the simulation studies, we implement our method using residual correlation estimator with sample splitting (RPCS). As is mentioned in Section 3.5, the nonparametric regression estimators $\widehat{m}_{x}$ and $\widehat{m}_{y}$ are fitted on the first part of samples, while the residual correlation estimator $\widehat{\rho}$ is computed on the second part. We further switch the roles of two parts and construct another estimate, and finally average over the two estimates to get the final result. We also consider a full-data version, referred to as residual correlation estimator with full data (RPCF), in which both the estimation of $\widehat{m}_{x},\widehat{m}_{y}$ and the computation of $\widehat{\rho}$ are performed using the entire dataset without sample splitting. As an oracle benchmark, we further include a test based on the Pearson correlation of the true regression errors $\epsilon_{x}$ and $\epsilon_{y}$, which we refer to as residual correlation estimator in the oracle setting (RPCO).

On performing the nonparametric regression for $m_{x}$ and $m_{y}$, we use Super Learner (38; 28; 29), which combines a collection of candidate regression algorithms including both parametric and nonparametric methods. In our implementation, we include the mean estimator (SL.mean), generalized linear models (SL.glm), penalized linear regression (SL.glmnet), random forests (SL.ranger), gradient boosting trees (SL.xgboost), and neural networks (SL.nnet). The Super Learner is implemented under the Gaussian loss using nonnegative least squares based on the Lawson–Hanson algorithm (method.NNLS) to estimate the ensemble weights for combining the candidate learners. In addition, 5-fold cross-validation is used to evaluate and combine the individual learners.

We compare our method with several existing approaches. The first is the partial correlation test (PaCo), which can be viewed as a special case of our framework in which the regression functions $m_{x}({\bf z})$ and $m_{y}({\bf z})$ are restricted to be linear in ${\bf z}$. Partial correlation plays a fundamental role in multivariate analysis (1) and graphical models (18). Moreover, its asymptotic distribution coincides with (18); see Section 5.3 of 23. The second method is based on residual Hilbert–Schmidt independence criterion (RHSIC), which measures the dependence between the residuals $\widehat{\epsilon}_{x}$ and $\widehat{\epsilon}_{y}$ using a kernel-based independence criterion; see, for example, 12; 13. The implementation of RHSIC utilizes the R package dHSIC (27). The next three methods are the residual Randomized Independence Test (RRIT), the Randomized Conditional Independence Test (RCIT), and the Randomized conditional Correlation Test (RCoT). Here, RCIT and RCoT were proposed by 33, while RRIT denotes the unconditional version of RCIT (or RCoT) applied to the residuals $\widehat{\epsilon}_{x}$ and $\widehat{\epsilon}_{y}$. As shown by 33, RCIT and RCoT achieve comparable or better empirical performance than the Kernel Conditional Independence Test (KCIT) of 42, including similar power and Type I error control, while being substantially more computationally efficient. Therefore, in our simulation studies, we include only RCIT and RCoT as representatives of nonparametric conditional independence tests. We implement RRIT, RCIT and RCoT by the R package RCIT (34).

Under the null hypothesis ($\rho=0$), the empirical levels for Model 1 are reported in Table 1, while the corresponding boxplots of p-values are shown in Figure 1. The analogous results for Model 2 are provided in Table S.1 and Figure S.1 in the supplement. Overall, most methods achieve empirical levels close to the nominal level of 0.05 when the sample size is sufficiently large. For smaller sample sizes, the residual-based tests still maintain levels reasonably close to 0.05. In contrast, the two nonparametric methods, RCIT and RCoT, appear unable to well control the Type I error rate, particularly when $n=100$ or $200$. A possible explanation is that these methods do not explicitly consider the additive noise structure and may therefore be more sensitive to spurious dependence in small-sample settings.

Under the alternative hypothesis, we report in Table 2 the empirical powers in Model 1 when $\rho$ is positive, and the results when $\rho$ is negative are symmetric and are reported in Table S.2 of the supplement. We also present the boxplots of p-values for Model 1 in Figures S.2–S.5 in the supplement. Analogous results for Model 2 are reported in Tables S.3–S.4 and Figures S.6–S.9 in the supplement.