Test and Measure for Partial Mean Dependence Based on Machine Learning Methods

In this section, we consider the partial mean independence test problem (1.1). That is, it aims to check whether $Y$ and $W$ are partially mean independent after controlling for the nonlinear effect of $Z$. We first make a new reformulation of the null hypothesis $H_{0}$ as follows

$\displaystyle{H}_{0}^{\prime}:\,\,\,\mbox{E}(\tilde{Y}|X)=\mbox{E}(\tilde{Y}),$

(2.1)

where $\tilde{Y}=Y-\mbox{E}(Y|Z)$. Note that $\mbox{E}(\tilde{Y})=\mbox{E}[Y-\mbox{E}(Y|Z)]=0$, and thus the right-hand side of $H^{\prime}_{0}$ is actually reduced to 0. Based on this reformulation $H^{\prime}_{0}$, a natural idea is to compare the estimators of $\mbox{E}(\tilde{Y}|X)$ and $\mbox{E}(\tilde{Y})$ via measuring whether $\mbox{E}\left[\{\mbox{E}(\tilde{Y}|X)-\mbox{E}(\tilde{Y})\}^{2}\right]$ equals 0. By introducing the sample analog estimator of the zero term $\mbox{E}(\tilde{Y})$, one could avoid the fundamental difficulty that a naive test statistic has a degenerate distribution. Let ${\mathcal{D}}=\{X_{i},Y_{i}\}_{i=1}^{N}$ be a random sample from the population $\{X,Y\}$. We propose a new test statistic based on machine learning methods and data splitting technique. To be specific, we split the data randomly into two independent parts, denoted by ${\mathcal{D}}_{1}$ and ${\mathcal{D}}_{2}$. Let $|{\mathcal{D}}_{2}|=n_{2}=:n$ and $|{\mathcal{D}}_{1}|=n_{1}=:n^{\gamma}$ with $\gamma>1$, where $n_{1}+n_{2}=n^{\gamma}+n=N$. Without loss of generality, we assume that $n^{\gamma}$ is an integer. We first estimate the conditional mean $h(Z_{i})=\mbox{E}(Y_{i}|Z_{i})$ using machine learning methods based on the first subset of data ${\mathcal{D}}_{1}$, denoted by $\widehat{h}_{{\mathcal{D}}_{1}}(Z_{i})$, and then estimate $g(X_{i})=\mbox{E}(\tilde{Y}_{i}|X_{i})$ by regressing $Y_{i}-\widehat{h}_{{\mathcal{D}}_{1}}(Z_{i})$ on $X_{i}$ using machine learning methods based on ${\mathcal{D}}_{1}$ to get its final estimator $\widehat{g}_{{\mathcal{D}}_{1}}(X_{i})$. Then, we construct the test statistic via comparing the estimators of $\mbox{E}(\tilde{Y}|X)$ and $\mbox{E}(\tilde{Y})$ based on ${\mathcal{D}}_{2}$,

$$T_{n}=\frac{1}{n}\sum_{i\in{\mathcal{D}}_{2}}\left[\widehat{g}_{{\mathcal{D}}_{1}}(X_{i})-\frac{1}{n}\sum_{j\in{\mathcal{D}}_{2}}\left\{Y_{j}-\widehat{h}_{{\mathcal{D}}_{1}}(Z_{j})\right\}\right]^{2}.$$

(2.2)

To simplify notations, here we use $i\in{\mathcal{D}}_{2}$ to denote that $i$ belongs to the index set corresponding to ${\mathcal{D}}_{2}$. Intuitively, the larger $T_{n}$ is, the stronger evidence we have to reject $H_{0}$. We call this test the partial mean independence test (pMIT).

To proceed, we make some discussions on the reformulation (2.1). Different from existing approaches to compare the two conditional expectations $\mbox{E}(Y|X)$ and $\mbox{E}(Y|Z)$, we compare one conditional expectation $\mbox{E}(\tilde{Y}|X)$ and one unconditional expectation $\mbox{E}(\tilde{Y})$ of transformed $\tilde{Y}$. With the unbalanced sample splitting strategy, under the null hypothesis, the estimation errors from $\sum_{i\in{\mathcal{D}}_{2}}\widehat{g}_{{\mathcal{D}}_{1}}^{2}(X_{i})$ are controlled and $\sum_{i\in{\mathcal{D}}_{2}}\widehat{g}_{{\mathcal{D}}_{1}}^{2}(X_{i})$ becomes negligible. While with the formulation $H^{\prime}_{0}$, $\sum_{i\in{\mathcal{D}}_{2}}(n^{-1}\sum_{j\in{\mathcal{D}}_{2}}\tilde{Y}_{j})^{2}$ converges to a chi-squared distribution and determines the asymptotic distribution of $T_{n}$, making it valid for inference. If instead, we do not estimate $\mbox{E}(\tilde{Y})$ and consider $n^{-1}\sum_{i\in{\mathcal{D}}_{2}}\widehat{g}^{2}_{{\mathcal{D}}_{1}}(X_{i})$ directly, this natural statistic still has the degeneracy issue. Thus, introducing the estimation of ${\rm E}(\tilde{Y})$ is crucial and is a new complement to existing methods. Essentially, we tackle the issue of degeneracy by introducing additional randomness in the error term of $Y$. However there are two differences between our procedure and Dai et al. (2024)’s procedure. Firstly, our procedure does not involve data perturbation and does not require to select perturbation size. Secondly and more subtlely, we introduce the randomness in $\tilde{Y}$ directly in the comparison with $g(X)$. While Dai et al. (2024)’s procedure involves data perturbation with the difference of prediction errors of $m(X)$ and $h(Z)$. This seemingly minor change has a major beneficial effect: theoretically our procedure can detect local alternative hypothesis which converges to the null at faster rate. Intuitively, our procedure directly detects shift from $g(X)$, while their procedure focuses on $\mbox{E}[g^{2}(X)]$. Furthermore, our reformulation is general and can be adopted to deal with other conditional mean restriction testing problems. For instance, we may extend our procedure to test the conditional independence. Actually the conditional independence of $Y$ and $W$ given $Z$ is equivalent to $\mbox{E}[I(Y\leq y)|X]=\mbox{E}[I(Y\leq y)|Z]$ for any $y$. We can then apply our method for each $y$ and integrate all the information at different $y$’s to construct the final test statistic.

Then, we derive the asymptotic distributions of the statistic $T_{n}$. Let $\sigma^{2}_{Y|Z}=\mbox{Var}\left\{Y-\mbox{E}(Y|Z)\right\}=\mbox{E}\left[\left\{Y-h(Z)\right\}^{2}\right]$ and it can be estimated by $\widehat{\sigma}^{2}_{Y|Z}={n}^{-1}\sum_{i\in{\mathcal{D}}_{2}}\left\{Y_{i}-\widehat{h}_{{\mathcal{D}}_{1}}(Z_{i})\right\}^{2}.$ We impose the following assumption on the estimation errors of the conditional means using machine learning methods. For some $a\in(0,1)$ and a generic positive constant $c$, we suppose

• (C1)

  $\mbox{E}\left[\left\{\widehat{h}_{{\mathcal{D}}_{1}}(Z)-h(Z)\right\}^{2}\right]\leq cn_{1}^{-a}=cn^{-\gamma a}$ and $\mbox{E}\left[\left\{\widehat{g}_{{\mathcal{D}}_{1}}(X)-g(X)\right\}^{2}\right]\leq cn_{1}^{-a}=cn^{-\gamma a}$.

The first part of (C1) reasonably assumes that machine learning methods are able to estimate the conditional mean accurately with a certain convergence rate of $n_{1}^{-a}$. For example, Bauer and Kohler (2019), Schmidt-Hieber (2020) and Kohler and Langer (2021) studied the high-level estimation errors of estimators of conditional means via DNNs. Specifically, Bauer and Kohler (2019) showed that ${\rm E}\left[\left\{\widehat{h}_{{\mathcal{D}}_{1}}(Z)-h(Z)\right\}^{2}\right]=O((\log n_{1})^{3}n_{1}^{-2q/(2q+d^{*})})$, where $q$ is the smoothness parameter and $d^{*}$ is the intrinsic dimensionality. From their corollary 1, the above rate can be even in the order of $(\log n_{1})^{3}n_{1}^{-2/3}$. Further, rates of convergence of various boosting procedures have been widely studied. For instance, Bühlmann and Yu (2003) derived that the estimation error achieves rate of convergence $O(n^{-2q/(2q+1)})$ for one-dimensional regression, when adopting smoothing splines as weak learners. Bühlmann (2006) established the consistency of $L_{2}$-Boosting in high-dimensional setting without stating the rate. Kueck et al. (2023) proved that for high-dimensional linear regression model, iterated post-$L_{2}$-Boosting and orthogonal $L_{2}$-Boosting achieve the rate of convergence $O(s\log p\log n/n)$ with $s$ being the sparsity level. More discussions can be seen after Condition (C2). For the second part of (C1), we note that $g(X_{i})={\rm E}(Y_{i}|X_{i})-{\rm E}(Y_{i}|Z_{i})$ can be estimated by $\widehat{g}_{{\mathcal{D}}_{1}}(X_{i})=\widehat{m}_{{\mathcal{D}}_{1}}(X_{i})-\widehat{h}_{{\mathcal{D}}_{1}}(Z_{i})$, where $\widehat{m}_{{\mathcal{D}}_{1}}(X_{i})$ is an estimator of the conditional mean $m(X_{i})={\rm E}(Y_{i}|X_{i})$. Thus, the second part of (C1) is deduced from the first assumption in (C1) and a similar assumption on $\widehat{m}_{{\mathcal{D}}_{1}}(X_{i})$. In the literature, Williamson et al. (2021) assumed that $\mbox{E}\left[\left\{\widetilde{h}_{{\mathcal{D}}_{1}}(Z)-h(Z)\right\}^{2}\right]\leq cn_{1}^{-a}$ with $1/2<a$, where $\widetilde{h}(Z)$ is obtained by regressing $\widehat{m}_{{\mathcal{D}}_{1}}(X)$ on $Z$ based on ${\mathcal{D}}_{1}$. Lundborg et al. (2022) also imposed similar conditions about estimation errors for the relevant estimators in their assumptions 3-4. In a special case, for linear regression model with ordinary least squared estimator, it can be shown that our proposed estimator satisfies $\mbox{E}\left[\left\{\widehat{g}_{{\mathcal{D}}_{1}}(X)-g(X)\right\}^{2}\right]\leq cn_{1}^{-1}$. Further for linear regression model with LASSO estimator, we show that our proposed estimator satisfies $\mbox{E}[\left\{\widehat{g}_{{\mathcal{D}}_{1}}(X)-g(X)\right\}^{2}]\leq cs\log p/n_{1}$ with $s$ being the sparsity level. Thus in both cases, the second assumption in (C1) holds.

Although we are not constrained to any particular estimation method to estimate $g(X)$, we have found that the proposed one works well in practice. In fact under the null hypothesis, it generally has very small overall estimation error $\sum_{i\in\mathcal{D}_{2}}\widehat{g}^{2}_{{\mathcal{D}}_{1}}(X_{i})$ and thus can control empirical sizes very well. While under the alternative hypothesis, as long as the over-fitting is not too severe, our procedure can still detect corresponding alternative hypothesis. As in Lundborg et al. (2022), we do not claim that our choice is always the best, but we find this choice works well in practice and set it to be a sensible default choice.

Next, we derive the asymptotic distributions of $T_{n}$ in the following theorem.

Theorem 2.1 demonstrates that $T_{n}$ is $n$-consistent under the null hypothesis and the asymptotic null distribution is a standard chi-squared distribution, which provides a computationally efficient way to compute the p-values without any resampling procedures. Under the fixed alternative hypothesis, $T_{n}$ is root-$n$-consistent and asymptotically normal. It also reveals that the proposed test has non-trivial power against the local alternative hypothesis which converges to the null at the root-$n$ rate. Although the root-$n$ rate is generally slower than the root-$N$ rate, it is very close to the root-$N$ rate under parametric regression models where we can take $a=1$ and $\gamma$ can be very close to 1. The condition $\gamma>1/a$ reveals the relative roles of the sample size $n_{1}$ used to estimate $g(X)$. If $g(X)$ can be estimated accurately, then $a$ could be close to 1 and the sample size $n_{1}$ can be reduced. However, it is challenging to determine $\gamma$ from a theoretical perspective, since $a$ is generally unknown, often related to the smoothness of the conditional mean function and the dimensionality. Practically, we suggest achieving unbalanced sample splitting by tuning the ratio of the sample size of $\mathcal{D}_{1}$ to the total sample size, i.e., $\xi=n_{1}/N$. We provide an adaptive procedure for selecting proper $\xi$ in section 2.2, through which Type I error can be well controlled without much loss of power.

Now suppose that $\mbox{E}[\{\widehat{g}_{{\mathcal{D}}_{1}}(X)-\tilde{g}(X)\}^{2}]=O(n^{-a}_{1})$. Here $\tilde{g}(X)$ may not be equal to the true function $g(X)$. As long as $\mbox{E}[\tilde{g}^{2}(X)]>c_{g}$ for some positive constant $c_{g}$ under the alternative hypothesis, we have $V_{n}=nT_{n}/\widehat{\sigma}_{Y|Z}^{2}=\sigma_{Y|Z}^{-2}\sum_{i\in{\mathcal{D}}_{2}}\left[\tilde{g}(X_{i})-n^{-1}\sum_{j\in{\mathcal{D}}_{2}}\left\{Y_{j}-h(Z_{j})\right\}\right]^{2}+o_{p}(1)$ and thus $V_{n}\rightarrow\infty$. That is, we still can detect the alternative hypothesis even when $\widehat{g}_{{\mathcal{D}}_{1}}(X)$ may not be a consistent estimator of $g(X)$.

We further note that Theorem 2.1 can be strengthened to uniform results. If we enhance the condition (C1) to be valid for all distributions $P$’s of $(X,Y)$ and assume that $c_{Y}^{\prime}<\mbox{E}|Y-h(Z)|^{2+\tau}<C_{Y}^{\prime}$ for all $P$’s with $\tau>0$ and $c_{Y}^{\prime}$, $C_{Y}^{\prime}$ being positive constants. Then, from the Lemmas 18-20 in Shah and Peters (2020) (see also Lemmas 16-17 in Lundborg et al. (2022)), we can show that under the null hypothesis, $V_{n}$ converges to chi-squared distribution uniformly and the test can uniformly control the nominal level $\alpha$.

From Theorem 2.1, in terms of ${\rm E}[g^{2}(X)]$, our procedure can detect local alternative hypothesis which converges to the null at $n$ rate. With assumption that $1/2<a<1$ as in condition (C2), the convergence rate is faster than root-$N$ and slower than $N$ rate. While from Lundborg et al. (2022), the convergence rates of Williamson et al. (2023) and Dai et al. (2024) are generally root-$N$. Although the procedure in Lundborg et al. (2022) can achieve the fastest possible $N$ rate, the theorem is obtained with fixed dimension $W$ under linear regression models. Furthermore, compared with Lundborg et al. (2022), we establish the asymptotic distribution of the proposed test statistic under alternative hypothesis within the model-free framework.

In addition, we comment that a larger training sample could have a negative impact on the statistical power for our testing procedure. This is the price we have to pay for dealing with the estimation errors of machine learning methods and the critical degenerate limiting null distribution issue. Similarly, Dai et al. (2024) also required larger training sample to make their bias-sd-ratio $T_{2}$ negligible. When the alternatives are believed to be complex, the use of machine learning methods including DNN can be helpful to estimate the conditional means accurately. If the relationship between the response and the covariates is believed to be parametric, parametric methods should be used and then it is not necessary to use the larger training sample. The loss of power associated with performing a test on a small testing sample could be offset by efficient estimation of regression functions. In the next session, we will also discuss some procedures for power enhancement.

We first present our pMIT procedure in Algorithm 1 based on single data splitting.

In the following, we make some improvements to enhance the empirical testing power. Inspired by the seminal idea of Fan et al. (2015), we first introduce a power-enhanced version of $V_{n}$, $V_{n}^{*}=V_{n}+\tau\sum_{i\in{\mathcal{D}}_{2}}\widehat{g}^{2}_{{\mathcal{D}}_{1}}(X_{i})$ with $\tau>0$. As discussed before under $H_{0}$, $\sum_{i\in{\mathcal{D}}_{2}}\widehat{g}^{2}_{{\mathcal{D}}_{1}}(X_{i})$ is asymptotically negligible and thus $V_{n}^{*}$ still converges to $\chi^{2}(1)$ under $H_{0}$, which implies that the size can be controlled asymptotically. While under the local alternative hypothesis $H_{1n}:g(X)=\delta(X)/\sqrt{n}$, $\sum_{i\in{\mathcal{D}}_{2}}\widehat{g}^{2}_{{\mathcal{D}}_{1}}(X_{i})\rightarrow\mbox{E}[\delta^{2}(X)]$ and thus $V_{n}^{*}\overset{d}{\rightarrow}\chi^{2}(1)+\sigma^{-2}_{Y|Z}{\rm E}[\delta^{2}(X)]+\tau{\rm E}[\delta^{2}(X)].$ Since $V_{n}^{*}\geq V_{n}$, the power is always enhanced. It reveals that the degeneracy issue can be a blessing in terms of power. Also the same idea can be applied to other recently developed test statistics. How to adaptively choose $\tau$ warrants further investigation. Selecting an appropriate value of $\tau$ involves a trade-off between size control and power enhancement. With larger $\tau$, the power is larger but the empirical size faces the risk of being distorted. In practice, we simply take $\tau=1$.

Next, we discuss how to select the sample splitting ratio $\xi$ to balance the estimation bias and the test efficiency in Algorithm 1, which can be viewed as the trade-off between type I and type II errors. For the selection of $\xi$, Dai et al. (2024) provided an easy-to-use formula and also a novel data-adaptive procedure. Inspired by Dai et al. (2024), we also consider a data-adaptive splitting procedure and select optimal $\xi$ by controlling the estimated Type I error on permutation datasets. To be specific, we first randomize the covariates $W$ and obtain the corresponding $p$-values. After that, we estimate the Type I error by using the proportions of rejections over $M$ permutations, and select the value of $\xi$ that can control the estimated Type I error under some predetermined significance level. This procedure is presented in Algorithm 2. In the process of searching optimal $\xi$ in the candidate set $\Xi$, it stops once the criterion $\widehat{{\rm Err}}(\xi)\leq\alpha$ is met in practice to enjoy the computational efficiency.

Finally, the testing performance of the pMIT procedure outlined in Algorithm 1 may be affected by the additional randomness of the data splitting. To improve the algorithm stability, we introduce an ensemble testing procedure $\rm pMIT_{M}$ based on multiple data splittings in Algorithm 3. Actually we aggregate $p$-values from different sample splits to obtain an adjusted $p$-value $Q^{*}$. The $p$-value aggregation is an important problem in statistics. For recent theoretical investigations, see for instance DiCiccio et al. (2020), Vovk and Wang (2020), and Choi and Kim (2023). Instead of combining p-values, Guo and Shah (2023) studied the combination of test statistics from data splitting and made a deep theoretical investigation.

The rejection of $H_{0}$ implies that the sub-vector of covariates $W$ is still able to add some contributions to the conditional mean of $Y$ after controlling for the nonlinear effect of $Z$. Then, it is of great interest to quantify the partial dependence between $W$ and $Y$ given $Z$. In this section, we propose a new partial dependence measure, called partial Generalized Measure of Correlation (pGMC), and derive its theoretical properties and confidence interval.

To quantify the partial dependence between $W$ and $Y$ given $Z$, a commonly used method is partial correlation coefficient, defined as

$$\rho(Y,W|Z)=\frac{\mbox{E}[\{Y-\mbox{E}(Y|Z)\}\{W-\mbox{E}(W|Z)\}]}{\sqrt{\mbox{E}[\{Y-\mbox{E}(Y|Z)\}^{2}]}\sqrt{\mbox{E}[\{W-\mbox{E}(W|Z)\}^{2}]}},$$

which equals the Pearson correlation coefficient between the errors $Y-\mbox{E}(Y|Z)$ and $W-\mbox{E}(W|Z)$. It only measures the linear correlation between $Y-\mbox{E}(Y|Z)$ and $W-\mbox{E}(W|Z)$. It can be zero even there exists explicit partial dependence between $Y$ and $W$ given $Z$. For example, let $W$ and $Z$ be two independent standard normal random variables and set $Y=Z+W^{2}$. In this case, $Y$ and $W$ are perfectly dependent given $Z$ but $\rho(Y,W|Z)=0$. To measure the nonlinear dependence for the conditional mean of $Y$ given $W$ adjusting for $Z$, we consider the following conditional variance decomposition

$\displaystyle\mbox{Var}(Y|Z)=\mbox{E}\left\{\mbox{Var}(Y|X)|Z\right\}+\mbox{Var}\left\{\mbox{E}(Y|X)|Z\right\}.$

It motivates us to consider the following partial measure

$$r^{2}(Y,W|Z)=\frac{\mbox{E}[\mbox{Var}\{\mbox{E}(Y|X)|Z\}]}{\mbox{E}[\mbox{Var}(Y|Z)]}=\frac{\mbox{E}[\{\mbox{E}(Y|X)-\mbox{E}(Y|Z)\}^{2}]}{\mbox{E}[\{Y-\mbox{E}(Y|Z)\}^{2}]}=\frac{\mbox{E}[\{m(X)-h(Z)\}^{2}]}{\mbox{E}[\{Y-h(Z)\}^{2}]},$$

(3.1)

which can be interpreted as the explained variance of $Y$ by $W$ given $Z$ in a model-free framework. To make $r^{2}(Y,W|Z)$ well defined, it is necessary to assume that the denominator of (3.1) is nonzero. That is, $Y$ is not almost surely equal to a measurable function of $Z$. Otherwise, given $Z$, $Y$ is a constant almost surely and it makes no sense to discuss the contribution of $W$ to the conditional mean of $Y$. When there is no conditional covariates set $Z$, $r^{2}(Y,W|Z)$ reduces to $\mbox{GMC}(Y|X)=\mbox{Var}\left\{\mbox{E}(Y|X)\right\}/\mbox{Var}(Y)$, the Generalized Measure of Correlation (GMC) proposed by Zheng et al. (2012). Thus, the pGMC can be considered as an extension of the GMC to measure the partial dependence between $W$ and $Y$ given $Z$.

Next, we discuss more properties of $r^{2}(Y,Z|X)$ in the following proposition.

In Proposition 3.1, properties (1)-(3) demonstrate that the new partial dependence $r^{2}(Y,W|Z)$ characterizes nonlinearly partial dependence. $r^{2}(Y,W|Z)$ is 0 if and only if $Y$ and $W$ are conditionally mean independent given $Z$, and is 1 if and only if $Y$ is equal to a measurable function of $W$ given $Z$. For example, for $Y=Z+W^{2}$ where $W$ and $Z$ be two independent standard normal random variables, $r^{2}(Y,W|Z)=1$. In property (4), $r^{2}(Y,W|Z)\geq\rho^{2}(Y,W|Z)$ further indicates that $r^{2}(Y,W|Z)$ is capable of revealing more conditional association information between $Y$ and $W$ beyond the linear relation. The relationship between the new partial dependence and the GMC in Zheng et al. (2012) is established in property (5). As the GMC is viewed as a generalization of the classical $R^{2}$ measurement to a nonparametric model, the new partial dependence can then be viewed as the relative difference in the population $R^{2}$ obtained using the full set of covariates $X$ as compared to the subset of covariates $Z$ only. Moreover, when the conditional mean of $Y$ does not depend on $Z$, the new partial dependence $r^{2}(Y,W|Z)$ then reduces to $\mbox{GMC}(Y|X)$. Thus $r^{2}(Y,W|Z)$ is a natural extension of $\mbox{GMC}(Y|X)$ for partial mean dependence. From property (6), $r^{2}(Y,W|Z)$ measures the relative reduction of population residual sum of squares with additional $W$ in the model. The last property (7) means that when the conditional information is more, the remaining variables provide less additional information about the conditional mean of the response. On the other hand, the partial dependence is always bounded by the $\mbox{GMC}(Y|X)$.

Next, we estimate the pGMC $r^{2}(Y,W|Z)$ via data splitting. We split the data randomly into two independent parts with equal sample size, denoted by ${\mathcal{D}}_{1}$ and ${\mathcal{D}}_{2}$, and estimate the conditional means $h(Z)$ and $m(X)$ using ${\mathcal{D}}_{1}$. According to the representation of $r^{2}(Y,W|Z)$ in property (6), we estimate its numerator $\mbox{E}\left[\{Y-h(Z)\}^{2}\right]-\mbox{E}\left[\{Y-m(X)\}^{2}\right]$ by the following statistic using ${\mathcal{D}}_{2}$,

$\displaystyle R_{n}=\frac{1}{n}\sum_{i\in{\mathcal{D}}_{2}}\left\{Y_{i}-\widehat{h}_{{\mathcal{D}}_{1}}(Z_{i})\right\}^{2}-\frac{1}{n}\sum_{i\in{\mathcal{D}}_{2}}\left\{Y_{i}-\widehat{m}_{{\mathcal{D}}_{1}}(X_{i})\right\}^{2}.$

To improve the estimation efficiency, we use the cross-fitting method to obtain another estimator $R_{n}^{\prime}$ similarly by swapping the roles of ${\mathcal{D}}_{1}$ and ${\mathcal{D}}_{2}$. Then, we define $R_{n}^{\ast}=(R_{n}+R_{n}^{\prime})/2$ to estimate the numerator of $r^{2}(Y,W|Z)$. We remark that the cross-fitting method has been shown to improve the estimation efficiency asymptotically in Fan et al. (2012), Chernozhukov et al. (2018) and Vansteelandt and Dukes (2022). Similarly, we can estimate the denominator $\mbox{E}\left[\left\{Y-h(Z)\right\}^{2}\right]$ of $r^{2}(Y,W|Z)$ by $\widehat{\sigma}^{2*}_{Y|Z}=\left(\widehat{\sigma}^{2}_{Y|Z}+\widetilde{\sigma}^{2}_{Y|Z}\right)/2$, where $\widehat{\sigma}^{2}_{Y|Z}={n}^{-1}\sum_{i\in{\mathcal{D}}_{2}}\left\{Y_{i}-\widehat{h}_{{\mathcal{D}}_{1}}(X_{i})\right\}^{2},$ and $\widetilde{\sigma}^{2}_{Y|Z}$ is similarly defined by swapping the role of ${\mathcal{D}}_{1}$ and ${\mathcal{D}}_{2}$. Thus, the pGMC $r^{2}(Y,W|Z)$ can be estimated by

$\displaystyle\widehat{r}^{2}(Y,W|Z)=\frac{R_{n}^{\ast}}{\widehat{\sigma}^{2*}_{Y|Z}}.$

It is worth noting that the balanced data splitting and the cross-fitting method is used to estimate the pGMC, which is different from the unbalanced data splitting for the significance testing procedure in the previous section. This is because, when the partial dependence between $W$ and $Y$ given $Z$ truly exists, i.e. $r^{2}(Y,W|Z)\neq 0$, we can use the balanced data splitting to improve the estimation efficiency under the following condition (C2).

• (C2)

  $\mbox{E}\left[\left\{\widehat{m}_{{\mathcal{D}}_{1}}(X)-m(X)\right\}^{2}\right]=o(n^{-1/2}_{1})$ and $\mbox{E}\left[\sigma^{2}(X)\left\{\widehat{m}_{{\mathcal{D}}_{1}}(X)-m(X)\right\}^{2}\right]=o(1)$, where $\epsilon=Y-m(X)$ and $\sigma^{2}(X)=\mbox{E}(\epsilon^{2}|X)$; $\mbox{E}\left\{\widehat{h}_{{\mathcal{D}}_{1}}(Z)-h(Z)\right\}^{2}=o(n^{-1/2}_{1})$ and
  $\mbox{E}\left[\delta^{2}(X)\left\{\widehat{h}_{{\mathcal{D}}_{1}}(Z)-h(Z)\right\}^{2}\right]=o(1)$, where $\eta=Y-h(Z)$ and $\delta^{2}(X)=\mbox{E}\left(\eta^{2}|X\right)$.

Here $n_{1}=n_{2}=n$. The assumptions $\mbox{E}\left[\left\{\widehat{m}_{{\mathcal{D}}_{1}}(X)-m(X)\right\}^{2}\right]=o(n^{-1/2})$ and $\mbox{E}\left[\left\{\widehat{h}_{{\mathcal{D}}_{1}}(Z)-h(Z)\right\}^{2}\right]=o(n^{-1/2})$ in (C2) are also adopted in Chernozhukov et al. (2018) (equation 3.8), Williamson et al. (2021) (assumption A1) and in Vansteelandt and Dukes (2022) (assumptions in Theorem 2). For a full discussion about the estimator of nuisance functions, see Chernozhukov et al. (2018). Chernozhukov et al. (2018) showed that the $l_{1}$-penalized and related methods in various sparse models can satisfy the above condition. For the DNN, Bauer and Kohler (2019) showed that ${\rm E}\left[\left\{\widehat{m}_{{\mathcal{D}}_{1}}(X)-m(X)\right\}^{2}\right]=O((\log n)^{3}n^{-2q/(2q+d^{*})})$, where $q$ is the smoothness parameter and $d^{*}$ is the intrinsic dimensionality. For other recent developments, see also Schmidt-Hieber (2020) and Kohler and Langer (2021). While the other assumptions $\mbox{E}\left[\sigma^{2}(X)\left\{\widehat{m}_{{\mathcal{D}}_{1}}(X)-m(X)\right\}^{2}\right]=o(1)$ and $\mbox{E}\left[\delta^{2}(X)\left\{\widehat{h}_{{\mathcal{D}}_{1}}(Z)-h(Z)\right\}^{2}\right]=o(1)$ follow if we assume that the conditional variance functions $\sigma^{2}(X)$ and $\delta^{2}(X)$ are bounded above.

Next, we establish the asymptotic normality for the $\widehat{r}^{2}(Y,W|Z)$ in Theorem 3.1.