A Distribution Free Conditional Independence Test with Applications to Causal Discovery

In this section, we propose to characterize the conditional dependence of $X$ and $Y$ given $Z$ through quantifying the mutual dependence among $U$, $V$ and $W$. Although our proposed test is based on the distance between characteristics functions, our proposed test is much simpler and has different asymptotic distribution as well as different convergence rate from the conditional distance correlation proposed by Wang et al. (2015). Let $\omega(\cdot)$ be an arbitrary positive weight function and $\varphi_{U,V,W}(\cdot)$, $\varphi_{U}(\cdot)$, $\varphi_{V}(\cdot)$, and $\varphi_{W}(\cdot)$ be the characteristic functions of $(U,V,W)$, $U$, $V$ and $W$, respectively. Then

			$U$, $V$ and $W$ are mutually independent
		$\displaystyle\Longleftrightarrow$	$\displaystyle\varphi_{U,V,W}(t_{1},t_{2},t_{3})=\varphi_{U}(t_{1})\varphi_{V}(t_{2})\varphi_{W}(t_{3})\textrm{ for all $t_{1},t_{2},t_{3}\in\mathbb{R}$}$
		$\displaystyle\Longleftrightarrow$	$\displaystyle\iiint\big{\|}\varphi_{U,V,W}(t_{1},t_{2},t_{3})-\varphi_{U}(t_{1})\varphi_{V}(t_{2})\varphi_{W}(t_{3})\big{\|}^{2}\omega(t_{1},t_{2},t_{3})dt_{1}dt_{2}dt_{3}=0,$

where $\|\psi\|^{2}=\psi^{\mbox{\tiny{T}}}\overline{\psi}$ for a complex-valued function $\psi$ and $\overline{\psi}$ is the conjugate of $\psi$. By choosing $\omega(t_{1},t_{2},t_{3})$ to be the joint probability density function of three independent and identically distributed standard Cauchy random variables, the integration in the above equation has a closed form,

			$\displaystyle Ee^{-|U_{1}-U_{2}|-|V_{1}-V_{2}|-|W_{1}-W_{2}|}-2Ee^{-|U_{1}-U_{3}|-|V_{1}-V_{4}|-|W_{1}-W_{2}|}$
			$\displaystyle\hskip 85.35826pt+Ee^{-|U_{1}-U_{2}|}Ee^{-|V_{1}-V_{2}|}Ee^{-|W_{1}-W_{2}|},$		(1)

where $(U_{k},V_{k},W_{k})$, $k=1,\ldots,4$, are four independent copies of $(U,V,W)$. Here the choice of the weight function $\omega(t_{1},t_{2},t_{3})$ is mainly for the convenient analytic form of the integration. Different from the distance correlation (Székely et al., 2007), our integration exists without any moment conditions on the data, which is more widely applicable. Furthermore, with the fact that $U\bot\!\!\!\bot W$ and $V\bot\!\!\!\bot W$, (1) boils down to

$\displaystyle E\left\{S_{U}(U_{1},U_{2})S_{V}(V_{1},V_{2})e^{-|W_{1}-W_{2}|}\right\},$

(2)

where $S_{U}(U_{1},U_{2})$ and $S_{V}(V_{1},V_{2})$ are defined as

	$\displaystyle S_{U}(U_{1},U_{2})$	$\displaystyle=$	$\displaystyle E\left\{e^{-|U_{1}-U_{2}|}+e^{-|U_{3}-U_{4}|}-e^{-|U_{1}-U_{3}|}-e^{-|U_{2}-U_{3}|}\mid(U_{1},U_{2})\right\},$
	$\displaystyle S_{V}(V_{1},V_{2})$	$\displaystyle=$	$\displaystyle E\left\{e^{-|V_{1}-V_{2}|}+e^{-|V_{3}-V_{4}|}-e^{-|V_{1}-V_{3}|}-e^{-|V_{2}-V_{3}|}\mid(V_{1},V_{2})\right\}.$

Recall that $U$, $V$ and $W$ are uniformly distributed on $(0,1)$. With further calculations based on (2), we obtain a normalized index and define it as $\rho$ to measure the mutual dependence:

	$\displaystyle\rho(X,Y\mid Z)$	$\displaystyle=$	$\displaystyle c_{0}E\big{\{}\left(e^{-|U_{1}-U_{2}|}+e^{-U_{1}}+e^{U_{1}-1}+e^{-U_{2}}+e^{U_{2}-1}+2e^{-1}-4\right)$		(3)
			$\displaystyle\left(e^{-|V_{1}-V_{2}|}+e^{-V_{1}}+e^{V_{1}-1}+e^{-V_{2}}+e^{V_{2}-1}+2e^{-1}-4\right)e^{-|W_{1}-W_{2}|}\big{\}},$

where $c_{0}=(13e^{-3}-40e^{-2}+13e^{-1})^{-1}$. Several appealing properties of the proposed index $\rho(X,Y\mid Z)$ are summarized in Theorem 1.

The step-by-step derivation of $\rho(X,Y\mid Z)$ and proof of Theorem 1 are presented in the appendix. Property (1) indicates that the index $\rho$ ranges from zero to one, equals zero when the conditional independence holds, and is equal to one if $Y$ is a strictly monotone transformation of $X$ conditional on $Z$. Property (2) shows that the index $\rho$ is a symmetric measure of conditional dependence. Property (3) illustrates that the index $\rho$ is invariant to any strictly monotone transformation. In fact, $\rho$ is not only invariant to marginal strictly monotone transformations, but also invariant to strictly monotone transformations conditional on $Z$. For example, it can be verified that $\rho(X,Y\mid Z)=\rho[m_{1}\{X-E(X\mid Z)\},Y\mid Z]$.

In this section, we establish the asymptotic properties of the sample version of the proposed index under the null and alternative hypothesis. Consider independent and identically distributed samples $\{X_{i},Y_{i},Z_{i}\}$, $i=1,\dots,n$. To estimate the proposed index $\rho(X,Y\mid Z)$, we apply kernel estimator for the conditional cumulative distribution function. Specifically, define

	$\displaystyle\widehat{f}_{Z}(z)=n^{-1}\sum_{i=1}^{n}K_{h}(z-Z_{i}),$
	$\displaystyle\widehat{U}=\widehat{F}_{X\mid Z}(x\mid z)=n^{-1}\sum_{i=1}^{n}K_{h}(z-Z_{i})\mathbbm{1}(X_{i}\leq x)/\widehat{f}_{Z}(z),$
	$\displaystyle\widehat{V}=\widehat{F}_{Y\mid Z}(y\mid z)=n^{-1}\sum_{i=1}^{n}K_{h}(z-Z_{i})\mathbbm{1}(Y_{i}\leq y)/\widehat{f}_{Z}(z),$

where $K_{h}(\cdot)=K(\cdot/h)/h$, $K(\cdot)$ is a kernel function, and $h$ is the bandwidth. Besides, we use empirical distribution function to estimate the cumulative distribution function, i.e., $\widehat{W}=\widehat{F}_{Z}(z)=n^{-1}\sum_{i=1}^{n}\mathbbm{1}(Z_{i}\leq z)$. The sample version of the index, denoted by $\widehat{\rho}(X,Y\mid Z)$, is thus given by

	$\displaystyle\widehat{\rho}(X,Y\mid Z)$	$\displaystyle=$	$\displaystyle c_{0}n^{-2}\sum_{i,j}\left\{\left(e^{-|\widehat{U}_{i}-\widehat{U}_{j}|}+e^{-\widehat{U}_{i}}+e^{\widehat{U}_{i}-1}+e^{-\widehat{U}_{j}}+e^{\widehat{U}_{j}-1}+2e^{-1}-4\right)\right.$
			$\displaystyle\left.\left(e^{-|\widehat{V}_{i}-\widehat{V}_{j}|}+e^{-\widehat{V}_{i}}+e^{\widehat{V}_{i}-1}+e^{-\widehat{V}_{j}}+e^{\widehat{V}_{j}-1}+2e^{-1}-4\right)e^{-|\widehat{W}_{i}-\widehat{W}_{j}|}\right\}.$

One can also obtain a normalized index $\rho_{0}$, which is a direct normalization based on (1) without considering $U\bot\!\!\!\bot W$ and $V\bot\!\!\!\bot W$:

	$\displaystyle{\rho}_{0}(X,Y\mid Z)$	$\displaystyle=$	$\displaystyle c_{0}\left\{Ee^{-|U_{1}-U_{2}|-|V_{1}-V_{2}|-|W_{1}-W_{2}|}+8e^{-3}\right.$
			$\displaystyle\left.-2E\left(2-e^{-U_{1}}-e^{U_{1}-1}\right)\left(2-e^{-V_{1}}-e^{V_{1}-1}\right)\left(2-e^{-W_{1}}-e^{W_{1}-1}\right)\right\}.$

The corresponding moment estimator is

	$\displaystyle\widehat{\rho}_{0}(X,Y\mid Z)$	$\displaystyle=$	$\displaystyle c_{0}\left\{n^{-2}\sum_{i,j}e^{-|\widehat{U}_{i}-\widehat{U}_{j}|-|\widehat{V}_{i}-\widehat{V}_{j}|-|\widehat{W}_{i}-\widehat{W}_{j}|}+8e^{-3}\right.$
			$\displaystyle\left.-2n^{-1}\sum_{i=1}^{n}\left(2-e^{-\widehat{U}_{i}}-e^{\widehat{U}_{i}-1}\right)\left(2-e^{-\widehat{V}_{i}}-e^{\widehat{V}_{i}-1}\right)\left(2-e^{-\widehat{W}_{i}}-e^{\widehat{W}_{i}-1}\right)\right\}.$

Although $\rho(X,Y\mid Z)=\rho_{0}(X,Y\mid Z)$ at the population level, those two statistics $\widehat{\rho}(X,Y\mid Z)$ and $\widehat{\rho}_{0}(X,Y\mid Z)$ exhibit different properties at the sample level. This is because $\widehat{\rho}(X,Y\mid Z)$ considers the fact that $U\bot\!\!\!\bot W$ and $V\bot\!\!\!\bot W$. But on the other hand, $\widehat{\rho}_{0}(X,Y\mid Z)$ is only a regular mutual independence test statistic, where $\widehat{U}_{i}$, $\widehat{V}_{i}$, and $\widehat{W}_{i}$ are exchangeable. When $X\bot\!\!\!\bot Y\mid Z$, under Conditions 1-4 listed below, $\widehat{\rho}(X,Y\mid Z)$ is of order $(n^{-1}+h^{4m})$, while $\widehat{\rho}_{0}(X,Y\mid Z)$ is of order $(n^{-1}+h^{2m})$ because of the bias caused by nonparametric estimation. Note that $m$ is the order of kernel functions and equal to 2 when using regular kernel functions such as Gaussian and epanechnikov kernels. This indicates that $\widehat{\rho}(X,Y\mid Z)$ is essentially $n$ consistent without under-smoothing while $\widehat{\rho}_{0}(X,Y\mid Z)$ typically requires under-smoothing. In addition, our statistic $\widehat{\rho}(X,Y\mid Z)$ has the same asymptotic properties as if $U,V$ and $W$ are observed, but $\widehat{\rho}_{0}(X,Y\mid Z)$ does not. See Figure 1 for a numerical comparison between the empirical null distributions of the two statistics.

We next study the asymptotical behaviors of the estimated index, $\widehat{\rho}(X,Y\mid Z)$, under both the null and the alternative hypotheses. The following regularity conditions are imposed to facilitate our subsequent theoretical analyses. In what follows, we derive the limiting distribution of $\widehat{\rho}(X,Y\mid Z)$ under the null hypothesis in Theorem 2.

Condition 1. The univariate kernel function $K(\cdot)$ is symmetric about zero and Lipschitz continuous. In addition, it satisfies

$\displaystyle\int K(\upsilon)d\upsilon=1,\quad\int\upsilon^{i}K(\upsilon)d\upsilon=0,1\leq i\leq m-1,\quad 0\neq\int\upsilon^{m}K(\upsilon)d\upsilon<\infty.$

Condition 2. The bandwidth $h$ satisfies $nh^{2}/\log^{2}(n)\to\infty$, and $nh^{4m}\to 0$.

Condition 3. The probability density function of $Z$, denoted by $f_{Z}(z)$ is bounded away from $0$ to infinity.

Condition 4. The $(m-1)$th derivatives of $F_{X\mid Z}(x\mid z)f(z)$, $F_{Y\mid Z}(y\mid z)f(z)$ and $f_{Z}(z)$ with respect to $z$ are locally Lipschitz-continuous.

The proof of Theorem 2 is given in the appendix. To understand the asymptotic distributions intuitively, we showed in the proof that $n\widehat{\rho}(X,Y\mid Z)$ can be approximated by degenerate V-statistics, i.e.,

$$n\widehat{\rho}(X,Y\mid Z)=c_{0}n^{-1}\sum_{i,j}h(U_{i},V_{i},W_{i};U_{j},V_{j},W_{j})+o_{p}(n^{-1}),$$

and $E\{h(U_{i},V_{i},W_{i};U_{j},V_{j},W_{j})\mid(U_{i},V_{i},W_{i})\}=0$. By the spectral decomposition,

$$h(u,v,w;u^{\prime},v^{\prime},w^{\prime})=\sum_{j=1}^{\infty}\lambda_{j}\Phi_{j}(u,v,w)\Phi_{j}(u^{\prime},v^{\prime},w^{\prime}).$$

Therefore, $n\widehat{\rho}(X,Y\mid Z)=c_{0}\sum_{j=1}^{\infty}\lambda_{j}\{n^{-1/2}\sum_{i=1}^{n}\Phi_{j}(U_{i},V_{i},W_{i})\}^{2}+o_{p}(n^{-1}),$ which converges in distribution to the weighted sum of independent chi squared distributions provided in Theorem 2 because $n^{-1/2}\sum_{i=1}^{n}\Phi_{j}(U_{i},V_{i},W_{i})$ is asymptotically standard normal (Korolyuk and Borovskich, 2013). Moreover, the $\lambda_{j}$s, $j=1,2,\dots$ are real numbers associated with the distribution of $U$, $V$ and $W$, all of which follow uniform distributions on $[0,1]$. In addition, $U$, $V$ and $W$ are mutually independent under the null hypothesis. This indicates that the proposed test statistic is essentially distribution free under the null hypothesis. Therefore, we suggest a simulation procedure to approximate the null distribution and decide the critical value. The simulation procedure can be independent of the original data and hence greatly improved the computation efficiency. In what follows, we describe the simulation-based procedure in detail to decide the critical value $c_{\alpha}$.

1. 1.

   Generate $\{U_{i}^{*},V_{i}^{*},W_{i}^{*}\}$, $i=1,\dots,n$ independently from mutually independent standard uniform distributions;

2. 2.

   Compute the statistic $\widehat{\rho^{*}}$ based on $\{U_{i}^{*},V_{i}^{*},W_{i}^{*}\}$, $i=1,\dots,n$, i.e.,

   	$\displaystyle\widehat{\rho^{*}}$	$\displaystyle=$	$\displaystyle c_{0}n^{-2}\sum_{i,j}\left\{\left(e^{-|U_{i}^{*}-U_{j}^{*}|}+e^{-U_{i}^{*}}+e^{U_{i}^{*}-1}+e^{-U_{j}^{*}}+e^{U_{j}^{*}-1}+2e^{-1}-4\right)\right.\hskip 65.44142pt$		(4)
   			$\displaystyle\left.\left(e^{-|V_{i}^{*}-V_{j}^{*}|}+e^{-V_{i}^{*}}+e^{V_{i}^{*}-1}+e^{-V_{j}^{*}}+e^{V_{j}^{*}-1}+2e^{-1}-4\right)e^{-|W_{i}^{*}-W_{j}^{*}|}\right\}.$

3. 3.

   Repeat Steps 1-2 for $B$ times and set $c_{\alpha}$ to be the upper $\alpha$ quantile of the estimated $\widehat{\rho^{*}}$ obtained from the randomly simulated samples.

Because $(U^{*},V^{*},W^{*})$ has the same distribution as that of $(U,V,W)$ under the null hypothesis, it is straightforward that this simulation-based procedure can provide a valid approximation of the asymptotic null distribution of $\widehat{\rho}(X,Y\mid Z)$ when $B$ is large. The consistency of this procedure is guaranteed by Theorem 3.

Next, we study the power performance of the proposed test under two kinds of alternative hypotheses, under which the conditional independence no longer holds. We first consider the global alternative, denoted by $H_{1g}$, we have

$\displaystyle H_{1g}:\quad X\not\!\perp\!\!\!\perp Y\mid Z.$

We then consider a sequence of local alternatives, denoted by $H_{1l}$,

$\displaystyle H_{1l}:\quad F_{X\mid Z}(x\mid Z=z)-F_{X\mid(Y,Z)}\left\{x\mid(Y=y,Z=z)\right\}=n^{-1/2}\ell(x,y,z).$

The asymptotical properties of the test statistics $\widehat{\rho}(X,Y\mid Z)$ under the global alternative and local alternatives are given in Theorem 4, whose proof is in the appendix. Theorem 4 shows that the proposed test can consistently detect any fixed alternatives as well as local alternatives at rate $O(n^{-1/2})$.

The methodology developed in Section 2 assumes that all the variables are univariate. In this section, we generalize the proposed index $\rho$ to the multivariate case. Let ${\bf x}=(X_{1},\ldots,X_{p})^{\mbox{\tiny{T}}}\in\mathbb{R}^{p}$, ${\bf y}=(Y_{1},\ldots,Y_{q})^{\mbox{\tiny{T}}}\in\mathbb{R}^{q}$ and ${\bf z}=(Z_{1},\ldots,Z_{r})^{\mbox{\tiny{T}}}\in\mathbb{R}^{r}$ be continuous random vectors. More specifically, all elements of ${\bf x}$, ${\bf y}$ and ${\bf z}$ are continuous random variables. Define $F_{X_{1}\mid{\bf z}}(X_{1}\mid{\bf z})$ for the cumulative distribution function of $X_{1}$ given ${\bf z}$ and $F_{X_{k}\mid{\bf z},X_{1},\ldots,X_{k-1}}(X_{k}\mid{\bf z},X_{1},\ldots,X_{k-1})$ for the cumulative distribution function of $X_{k}$ given ${\bf z},X_{1},\ldots,X_{k-1}$ for $k=2,\cdots,p$. Similar notation apply for ${\bf y}$ and ${\bf z}$. Denote $\widetilde{U}_{1}=F_{X_{1}\mid{\bf z}}(X_{1}\mid{\bf z})$, $\widetilde{V}_{1}=F_{Y_{1}\mid{\bf z}}(Y_{1}\mid{\bf z})$, $\widetilde{W}_{1}=F_{Z_{1}}(Z_{1})$,

	$\displaystyle\widetilde{U}_{k}=F_{X_{k}\mid{\bf z},X_{1},\ldots,X_{k-1}}(X_{k}\mid{\bf z},X_{1},\ldots,X_{k-1}),\ k=2,\ldots,p,$
	$\displaystyle\widetilde{V}_{k}=F_{Y_{k}\mid{\bf z},Y_{1},\ldots,Y_{k-1}}(Y_{k}\mid{\bf z},Y_{1},\ldots,Y_{k-1}),\ k=2,\ldots,q,$
	$\displaystyle\widetilde{W}_{k}=F_{Z_{k}\mid Z_{1},\ldots,Z_{k-1}}(Z_{k}\mid Z_{1},\ldots,Z_{k-1}),\ k=2,\ldots,r.$

Further denote ${\bf u}=(\widetilde{U}_{1},\ldots,\widetilde{U}_{p})^{\mbox{\tiny{T}}}$, ${\bf v}=(\widetilde{V}_{1},\ldots,\widetilde{V}_{q})^{\mbox{\tiny{T}}}$, ${\bf w}=(\widetilde{W}_{1},\ldots,\widetilde{W}_{r})^{\mbox{\tiny{T}}}$. Similar to test of conditional independence for random variables, we first establish an equivalence between the conditional independence ${\bf x}\bot\!\!\!\bot{\bf y}\mid{\bf z}$ and the mutual independence of ${\bf u}$,${\bf v}$ and ${\bf w}$, which is stated in Theorem 5.

The proof of Theorem 5 is illustrated in the appendix. Theorem 5 established an equivalence between the conditional independence ${\bf x}\bot\!\!\!\bot{\bf y}\mid{\bf z}$ and the mutual independence among ${\bf u}$, ${\bf v}$ and ${\bf w}$. It is notable that when $p,q,r$ are relatively large, $({\bf u},{\bf v},{\bf w})$ may be difficult to estimate because of the curse of dimensionality. In this paper, we mainly focus on the low dimensional case. Next, we develop the mutual independence test among ${\bf u}$,${\bf v}$ and ${\bf w}$. Similar as the univariate case, we set the weight function to be the joint density of $(p+q+r)$ independent and identically distributed standard Cauchy random variables. We may further derive the closed form expression of $\rho({\bf x},{\bf y}\mid{\bf z})$

$\displaystyle\rho({\bf x},{\bf y}\mid{\bf z})=E\left\{S_{{\bf u}}({\bf u}_{1},{\bf u}_{2})S_{{\bf v}}({\bf v}_{1},{\bf v}_{2})e^{-\|{\bf w}_{1}-{\bf w}_{2}\|_{1}}\right\},$

where $({\bf u}_{1},{\bf v}_{1},{\bf w}_{1})$ and $({\bf u}_{2},{\bf v}_{2},{\bf w}_{2})$ are two independent copies of $({\bf u},{\bf v},{\bf w})$. Moreover,

	$\displaystyle S_{{\bf u}}({\bf u}_{1},{\bf u}_{2})$	$\displaystyle=$	$\displaystyle E\left\{e^{-\|{\bf u}_{1}-{\bf u}_{2}\|_{1}}+e^{-\|{\bf u}_{3}-{\bf u}_{4}\|_{1}}-e^{-\|{\bf u}_{1}-{\bf u}_{3}\|_{1}}-e^{-\|{\bf u}_{2}-{\bf u}_{3}\|_{1}}\mid({\bf u}_{1},{\bf u}_{2})\right\},$
	$\displaystyle S_{{\bf v}}({\bf v}_{1},{\bf v}_{2})$	$\displaystyle=$	$\displaystyle E\left\{e^{-\|{\bf v}_{1}-{\bf v}_{2}\|_{1}}+e^{-\|{\bf v}_{3}-{\bf v}_{4}\|_{1}}-e^{-\|{\bf v}_{1}-{\bf v}_{3}\|_{1}}-e^{-\|{\bf v}_{2}-{\bf v}_{3}\|_{1}}\mid({\bf v}_{1},{\bf v}_{2})\right\}.$

and $\|\cdot\|_{1}$ is the $\ell_{1}$ norm. Then $\rho({\bf x},{\bf y}\mid{\bf z})$ is nonnegative and equals zero if and only if ${\bf x}\bot\!\!\!\bot{\bf y}\mid{\bf z}$. By estimating $\rho({\bf x},{\bf y}\mid{\bf z})$ consistently at the sample level, the resulting test is clearly consistent. To implement the test, it is still required to study the asymptotic distributions under the conditional independence using independent and identically distributed samples $\{{\bf x}_{i},{\bf y}_{i},{\bf z}_{i}\}$, $i=1,\dots,n$. We also apply kernel estimator for the conditional cumulative distribution functions when estimating ${\bf u}_{i},{\bf v}_{i}$ and ${\bf w}_{i}$. Specifically, we estimate $F_{A\mid B_{1},\ldots,B_{\ell}}(a\mid b_{1},\ldots,b_{\ell})$ with

$\displaystyle\widehat{F}_{A\mid B_{1},\ldots,B_{\ell}}(a\mid b_{1},\ldots,b_{\ell})=\frac{\sum_{i=1}^{n}\mathbbm{1}(A_{i}\leq a)\prod_{k=1}^{\ell}K_{h}(B_{ik}-b_{k})}{\sum_{i=1}^{n}\prod_{k=1}^{\ell}K_{h}(B_{ik}-b_{k})}.$

where $(A,B_{1},\ldots,B_{\ell})^{\mbox{\tiny{T}}}$ are $(Z_{k},Z_{1},\ldots,Z_{k-1})^{\mbox{\tiny{T}}}$, $k=2,\ldots,r$, $(X_{\ell},X_{1},\ldots,X_{\ell-1},{\bf z}^{\mbox{\tiny{T}}})^{\mbox{\tiny{T}}}$, $\ell=1,\ldots,p$, or $(Y_{j},Y_{1},\ldots,Y_{j-1},{\bf z}^{\mbox{\tiny{T}}})^{\mbox{\tiny{T}}}$, $j=1,\ldots,q$ when estimating ${\bf w}$, ${\bf u}$ and ${\bf v}$, respectively. The sample version of $\rho({\bf x},{\bf y}\mid{\bf z})$ is given by

	$\displaystyle\widehat{\rho}({\bf x},{\bf y}\mid{\bf z})$	$\displaystyle=$	$\displaystyle n^{-2}\sum_{i,j}E\Big{[}\left\{e^{-\|\widehat{\bf u}_{i}-\widehat{\bf u}_{j}\|_{1}}+e^{-\|{\bf u}-{\bf u}^{\prime}\|_{1}}-e^{-\|\widehat{\bf u}_{i}-{\bf u}\|_{1}}-e^{-\|{\bf u}-\widehat{\bf u}_{j}\|_{1}}\mid(\widehat{\bf u}_{i},\widehat{\bf u}_{j})\right\}$
			$\displaystyle E\left\{e^{-\|\widehat{\bf v}_{i}-\widehat{\bf v}_{j}\|_{1}}+e^{-\|{\bf v}-{\bf v}^{\prime}\|_{1}}-e^{-\|\widehat{\bf v}_{i}-{\bf v}\|_{1}}-e^{-\|{\bf v}-\widehat{\bf v}_{j}\|_{1}}\mid(\widehat{\bf v}_{i},\widehat{\bf v}_{j})\right\}e^{-\|\widehat{\bf w}_{i}-\widehat{\bf w}_{j}\|_{1}}\Big{]},$

where $({\bf u}^{\prime},{\bf v}^{\prime})$ is an independent copy of $({\bf u},{\bf v})$, and further calculations yield that $\widehat{\rho}({\bf x},{\bf y}\mid{\bf z})$ is equal to

	$\displaystyle n^{-2}\sum_{i,j}\left[\left\{e^{-\|\widehat{\bf u}_{i}-\widehat{\bf u}_{j}\|_{1}}+(\frac{2}{e})^{p}-\prod_{k=1}^{p}(2-e^{-\widehat{\widetilde{U}}_{ik}-1}-e^{-\widehat{\widetilde{U}}_{ik}})-\prod_{k=1}^{p}(2-e^{-\widehat{\widetilde{U}}_{jk}-1}-e^{-\widehat{\widetilde{U}}_{jk}})\right\}\right.$
	$\displaystyle\left.\left\{e^{-\|\widehat{\bf v}_{i}-\widehat{\bf v}_{j}\|_{1}}+(\frac{2}{e})^{q}-\prod_{k=1}^{q}(2-e^{-\widehat{\widetilde{V}}_{ik}-1}-e^{-\widehat{\widetilde{V}}_{ik}})-\prod_{k=1}^{q}(2-e^{-\widehat{\widetilde{V}}_{jk}-1}-e^{-\widehat{\widetilde{V}}_{jk}})\right\}e^{-\|\widehat{\bf w}_{i}-\widehat{\bf w}_{j}\|_{1}}\right].$

We next study the asymptotical behaviors of $\widehat{\rho}({\bf x},{\bf y}\mid{\bf z})$ under the null hypothesis in Theorem 6, whose proof is given in the appendix. We begin by providing some regularity conditions for the multivariate data.