Machine-Learning-Assisted Comparison of Regression Functions

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

53,062 characters · 9 sections · 40 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

\affil[1]{Department of Statistics and Data Science, Cornell University} \affil[2]{Department of ISOM, Hong Kong University of Science and Technology} \affil[3]{Department of Biostatistics, Epidemiology and Informatics, University of Pennsylvania}

abstractWe revisit the classical problem of comparing regression functions, a fundamental question in statistical inference with broad relevance to modern applications such as data integration, transfer learning, and causal inference. Existing approaches typically rely on smoothing techniques and are thus hindered by the curse of dimensionality. We propose a generalized notion of kernel-based conditional mean dependence that provides a new characterization of the null hypothesis of equal regression functions. Building on this reformulation, we develop two novel tests that leverage modern machine learning methods for flexible estimation. We establish the asymptotic properties of the test statistics, which hold under both fixed- and high-dimensional regimes. Unlike existing methods that often require restrictive distributional assumptions, our framework only imposes mild moment conditions. The efficacy of the proposed tests is demonstrated through extensive numerical studies.

{\smallKeywords: Comparison of regression functions; High dimensionality; Kernel-based conditional mean dependence; Machine learning.}

Introduction

The comparison of two regression functions is a fundamental problem in regression analysis munk1998nonparametric. Consider two independent data sets $\mathcal{D}^{(l)}=\{(Y_{i}^{(l)},X_{i}^{(l)})\}_{i=1}^{N_{l}}$ for $l=1,2$, where $Y_{i}^{(l)}\in\mathbb{R}$ is the response and $X_{i}^{(l)}\in\mathbb{R}^{p}$ is a set of covariates with dimension $p$ allowed to diverge. For $l=1,2$, we assume that $\{(Y_{i}^{(l)},X_{i}^{(l)})\}_{i=1}^{N_{l}}$ are independent and identically distributed samples of $(Y^{(l)},X^{(l)})\sim P^{(l)}\equiv P_{Y\mid X}^{(l)}\otimes P_{X}^{(l)}$, and define the regression function $m^{(l)}(x)=\mathbb{E}(Y^{(l)}\mid X^{(l)}=x)$ and the error $\varepsilon^{(l)}=Y^{(l)}-m^{(l)}(X^{(l)})$. We aim to test

equation[equation omitted — 174 chars of source]

To ensure that the testing problem is nontrivial, we assume that $P_{X}^{(1)}\ll P_{X}^{(2)}$ and $P_{X}^{(2)}\ll P_{X}^{(1)}$, where the symbol $\ll$ stands for absolute continuity. Thus, the hypotheses ((ref)) can be equivalently stated by replacing $P^{(1)}_{X}$ with $P^{(2)}_{X}$.

The problem ((ref)) also naturally arises in a variety of related areas. For instance, in data integration, it is necessary to assess whether two data sets share a common regression function before proceeding with modeling, estimation, and inference. When the null hypothesis $H_{0}$ holds, pooling information across sources can enhance statistical efficiency and lead to more reliable conclusions. In contrast, when heterogeneity is present, naive aggregation of data sets may obscure meaningful differences and result in misleading scientific findings. A similar issue appears in transfer learning, where it is commonly assumed that the conditional mean is identical in the source and target populations wang2025phase—a weaker version of the standard covariate shift assumption. Testing hypotheses ((ref)) is thus crucial for validating methods built upon this assumption. Another example comes from causal inference. In the binary treatment setting, testing the null hypothesis of a zero conditional average treatment effect crump2008nonparametric can be formulated as problem ((ref)) under the standard assumptions of consistency and no unmeasured confounding. In this context, the two independent samples are naturally defined by the treatment indicator.

Much effort has been devoted to the problem ((ref)) in the literature; see Section 7 of gonzalez2013updated for a comprehensive review. Classical methods employ parametric models for the regression functions and assess equality through comparisons of model parameters. However, this approach necessitates correct model specification, which is often unrealistic in practice. The problem of testing the equality of two regression functions in the nonparametric setting was first considered by hall1990bootstrap,king1991testing, and subsequently explored in a series of studies including delgado1993testing,kulasekera1995comparison,young1995non,dette2001nonparametric,lavergne2001equality,neumeyer2003nonparametric,pardo2007testing,srihera2010nonparametric,pardo2015non. The existing literature primarily focuses on the univariate case, i.e., $p=1$. Most available methods rely on smoothing-based estimators of the regression functions, and thus suffer from the curse of dimensionality even if they can be extended to multivariate covariates. Furthermore, restrictive distributional assumptions are often imposed to facilitate theoretical analysis. hall1990bootstrap,king1991testing,delgado1993testing concentrate on identical design points. king1991testing,young1995non assume Gaussian errors. kulasekera1995comparison assumes independence between $\varepsilon^{(l)}$ and $X^{(l)}$, and dette2001nonparametric,neumeyer2003nonparametric,pardo2007testing,pardo2015non posit a heteroscedastic model of the form $\varepsilon^{(l)}=\sigma^{(l)}(X^{(l)})e^{(l)}$, where $e^{(l)}$ is independent of $X^{(l)}$ and $\sigma^{(l)}(\cdot)$ denotes the variance function. As noted by racine2020smooth, these independence assumptions may be unduly stringent and should themselves be subjected to empirical testing. Moreover, several studies dette2001nonparametric,neumeyer2003nonparametric,pardo2007testing,pardo2015non assume that $X^{(l)}$ has compact support (e.g., $[0,1]$) with density bounded away from zero, thereby excluding frequently encountered Gaussian distributions.

In contrast to traditional smoothing-based methods, modern machine learning techniques offer flexible and powerful tools for modeling regression functions while alleviating the curse of dimensionality. In this paper, we propose two tests that leverage machine learning methods to address problem ((ref)). To this end, we introduce a generalized notion of kernel-based conditional mean dependence which yields a novel characterization of the null hypothesis. Our proposed tests possess several attractive features relative to existing approaches:

itemize• They remain valid irrespective of whether the covariate dimension $p$ is fixed or diverging. • They do not rely on any specific functional form of the regression functions and require only mild moment conditions on the covariates and errors. • They allow for heterogeneous covariate and error distributions across the two samples. • Under $H_{0}$, the test statistics admit simple Gaussian limiting distributions, allowing straightforward implementation without resorting to resampling procedures. Additionally, users retain full flexibility in the choice of regression methods.

Notation. Denote by $z_{\alpha}$ the $\alpha$ quantile of the standard normal distribution, for $\alpha\in(0,1)$. For two sequences of real numbers $\{a_{n}\}_{n=1}^{\infty}$ and $\{b_{n}\}_{n=1}^{\infty}$, we write $a_{n}=O(b_{n})$ if there exists a finite $n_{0}\ge 1$ and a finite $C>0$ such that $|a_{n}|\le C|b_{n}|$ for all $n\ge n_{0}$; $a_{n}=o(b_{n})$ if $\lim_{n\rightarrow\infty}a_{n}/b_{n}=0$; $a_{n}=\Omega(b_{n})$ if $b_{n}=O(a_{n})$; $a_{n}=\omega(b_{n})$ if $b_{n}=o(a_{n})$; $a_{n}\asymp b_{n}$ if $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. The symbols $\overset{p}{\rightarrow}$ and $\overset{d}{\rightarrow}$ stand for convergence in probability and in distribution, respectively. Denote by $\|\cdot\|$ the Euclidean norm in $\mathbb{R}^{p}$. For any probability measure $P$ on $\mathbb{R}^{p}$, let $L^{2}(P)$ denote the Hilbert space of all square integrable functions on $\mathbb{R}^{p}$. For $f\in L^{2}(P)$, denote by $\|f\|_{L^{2}(P)}=(\int f^{2}dP)^{1/2}$ the $L^{2}(P)$ norm of $f$.

Preliminaries

This section reviews reproducing kernel Hilbert spaces and introduces a generalized notion of the kernel-based conditional mean dependence.

Let $\mathcal{H}$ be a Hilbert space of real-valued functions defined on a topological space $\mathcal{U}$. A function $k:\mathcal{U}\times\mathcal{U}\rightarrow\mathbb{R}$ is called a reproducing kernel of $\mathcal{H}$ if (i) $\forall u\in\mathcal{U}$, $k(\cdot,u)\in\mathcal{H}$, and (ii) $\forall u\in\mathcal{U},\forall f\in\mathcal{H}$, $\langle f,k(\cdot,u)\rangle_{\mathcal{H}}=f(u)$, where $\langle\cdot,\cdot\rangle_{\mathcal{H}}$ is the inner product associated with $\mathcal{H}$. If $\mathcal{H}$ has a reproducing kernel, it is said to be a reproducing kernel Hilbert space (RKHS). According to the Moore-Aronszajn theorem, for every symmetric, positive definite function (henceforth kernel) $k:\mathcal{U}\times\mathcal{U}\rightarrow\mathbb{R}$, there is an associated RKHS $\mathcal{H}_{k}$ with reproducing kernel $k$. For $\theta>0$, define $\mathcal{M}_{k}^{\theta}(\mathcal{U})=\{\mu\in\mathcal{M}(\mathcal{U}):\int k^{\theta}(u,u)d|\mu|(u)<\infty\}$, where $\mathcal{M}(\mathcal{U})$ denotes the set of all finite signed Borel measures on $\mathcal{U}$. The kernel mean embedding of $P\in\mathcal{M}_{k}^{1/2}(\mathcal{U})$ into the RKHS $\mathcal{H}_{k}$ is defined by the Bochner integral $\Pi_{k}(P)=\int k(\cdot,u)dP(u)\in\mathcal{H}_{k}$. The kernel $k$ is said to be characteristic if the mapping $\Pi_{k}$ is injective. Conditions under which kernels are characteristic have been studied by sriperumbudur2008injective,sriperumbudur2010hilbert. In particular, the Gaussian kernel and the Laplace kernel on $\mathbb{R}^{p}$ are characteristic.

Let $V\in\mathbb{R}$ be an arbitrary random variable and $U\in\mathcal{U}$ be an arbitrary random element. For a given kernel $k:\mathcal{U}\times\mathcal{U}\rightarrow\mathbb{R}$, under the assumptions that $\mathbb{E}(V^{2})<\infty$ and $U\sim P\in\mathcal{M}_{k}^{1}(\mathcal{U})$, the kernel-based conditional mean dependence (KCMD) of $V$ on $U$ is defined as lai2021kernel \[ \text{KCMD}(V\mid U)=\mathbb{E}[\{V-\mathbb{E}(V)\}\{V'-\mathbb{E}(V)\}k(U,U')], \] where $(V',U')$ denotes an independent copy of $(V,U)$. The martingale difference divergence proposed in shao2014martingale is a special case of the KCMD associated with the kernel induced by the Euclidean distance lai2021kernel. When the kernel $k$ is characteristic, we have $\text{KCMD}(V\mid U)\ge 0$, and $\text{KCMD}(V\mid U)=0$ if and only if $\mathbb{E}(V\mid U)=\mathbb{E}(V)$ almost surely.

The definition of the KCMD can be generalized as follows. For a given $v_{0}\in\mathbb{R}$, define \[ \text{KCMD}^{*}(V\mid U)=\mathbb{E}\{(V-v_{0})(V'-v_{0})k(U,U')\}, \] which replaces $\mathbb{E}(V)$ with $v_{0}$.

lemmaAssume that $\mathbb{E}(V^{2})<\infty$ and $U\sim P\in\mathcal{M}_{k}^{1}(\mathcal{U})$. When the kernel $k$ is characteristic, for any $v_{0}\in\mathbb{R}$, we have $\textup{KCMD}^{*}(V\mid U)\ge0$, and $\textup{KCMD}^{*}(V\mid U)=0$ if and only if $\mathbb{E}(V\mid U)=v_{0}$ almost surely.

This generalization of the KCMD, which is not considered in lai2021kernel, serves as the foundation of our proposed test statistics. The proof of Lemma (ref) is provided in the appendix.

Test statistic and asymptotic properties

We first make a new reformulation of the null hypothesis $H_{0}$. Define

align*[align* omitted — 99 chars of source]

which become the true errors $\varepsilon^{(1)}$ and $\varepsilon^{(2)}$ under $H_{0}$. Then $H_{0}$ holds if and only if

equation[equation omitted — 110 chars of source]

Based on this reformulation, for a given kernel $k$ on $\mathbb{R}^{p}$, define \[ \Delta^{(l)}=\text{KCMD}^{*}(\eta^{(l)}\mid X^{(l)})=\mathbb{E}\left\{\eta^{(l)}\eta^{(l)\prime}k(X^{(l)},X^{(l)\prime})\right\}\quad(l=1,2), \] where $(\eta^{(l)\prime},X^{(l)\prime})$ denotes an independent copy of $(\eta^{(l)},X^{(l)})$. In this paper, we make the following assumption on our kernel. Examples include the Gaussian kernel and the Laplace kernel.

assumptionThe kernel $k$ is characteristic. Besides, there exists a constant $K>0$ such that $\sup_{x,x'}|k(x,x')|\le K$.

theoremAssume that $\mathbb{E}(|\eta^{(l)}|^{2})<\infty$ for $l=1,2$. Under Assumption (ref), we have $\Delta^{(l)}\ge 0$, and $\Delta^{(l)}=0$ if and only if the null hypothesis $H_{0}$ in ((ref)) holds.

We now introduce our test statistic:

enumerate[(i)] • For $l=1,2$, estimate the regression function $m^{(l)}(\cdot)$ using machine learning methods based on the data set $\mathcal{D}^{(l)}$, denoted by $\widehat{m}^{(l)}(\cdot)$. • Obtain \begin{align*} \widehat{\eta}^{(1)}_{i}&=Y^{(1)}_{i}-\widehat{m}^{(2)}(X^{(1)}_{i})\quad(i=1,\ldots,N_{1}),\\ \widehat{\eta}^{(2)}_{i}&=Y^{(2)}_{i}-\widehat{m}^{(1)}(X^{(2)}_{i})\quad(i=1,\ldots,N_{2}). \end{align*} Note that $\widehat{\eta}^{(l)}_{i}$ is different from the residual $\widehat{\varepsilon}^{(l)}_{i}=Y^{(l)}_{i}-\widehat{m}^{(l)}(X^{(l)}_{i})$. • Without loss of generality, we assume that $N_{l}=2n_{l}\ (l=1,2)$. Define \[ \widehat{\Delta}^{(l)}=\frac{1}{n_{l}}\sum_{i=1}^{n_{l}}\widehat{\eta}^{(l)}_{i}\widehat{\eta}^{(l)}_{i+n_{l}}k(X^{(l)}_{i},X^{(l)}_{i+n_{l}}). \] • Construct the test statistic \[ T=\widehat{\Delta}^{(1)}+\widehat{\Delta}^{(2)}. \]

We make the following assumptions throughout the analysis.

assumptionFor $l=1,2$, assume that $\mathbb{E}(|\varepsilon^{(l)}|^{4})\le C_{1}<\infty$ for some $C_{1}>0$, and \[ \mathbb{E}\left\{\left|\varepsilon^{(l)}\varepsilon^{(l)\prime}k(X^{(l)},X^{(l)\prime})\right|^{2}\right\}\ge c_{1}, \] for some $c_{1}>0$, where $(\varepsilon^{(l)\prime},X^{(l)\prime})$ denotes an independent copy of $(\varepsilon^{(l)},X^{(l)})$.

assumptionSuppose \begin{gather*} \mathbb{E}\left[\left\{\widehat{m}^{(1)}(X^{(2)})-m^{(1)}(X^{(2)})\right\}^{2}\mid\mathcal{D}^{(1)}\right]=o_{p}(n_{1}^{-1/2}),\\ \mathbb{E}\left[\left\{\widehat{m}^{(2)}(X^{(1)})-m^{(2)}(X^{(1)})\right\}^{2}\mid\mathcal{D}^{(2)}\right]=o_{p}(n_{2}^{-1/2}). \end{gather*}

We emphasize that the condition $n_{1}\asymp n_{2}$ is not required. Hence, our procedure remains valid under imbalanced sample sizes. Assumption (ref) guarantees the applicability of the Lyapunov central limit theorem (CLT). By comparison, zhang2018conditional,li2023testing,he2025goodness impose a stronger condition: for some constants $c_{2}$ and $C_{2}$, \[ 0<c_{2}\le\mathbb{E}(|\varepsilon^{(l)}|^{2}\mid X^{(l)})\le\mathbb{E}^{1/2}(|\varepsilon^{(l)}|^{4}\mid X^{(l)})\le C_{2}<\infty\quad\text{almost surely}, \] which directly implies Assumption (ref). Assumption (ref) accommodates potential covariate shift. In the special case where $P^{(1)}_{X}=P^{(2)}_{X}$, Assumption (ref) reduces to \[ \mathbb{E}\left[\left\{\widehat{m}^{(l)}(X^{(l)})-m^{(l)}(X^{(l)})\right\}^{2}\mid\mathcal{D}^{(l)}\right]=o_{p}(n_{l}^{-1/2})\quad(l=1,2), \] a condition commonly assumed in the literature on double/debiased machine learning and non- or semiparametric estimation chernozhukov2018double,kennedy2024semiparametric. Flexible machine learning methods can be employed to estimate the regression functions $m^{(l)}(\cdot)$. Importantly, Assumption (ref) is substantially weaker than conditions requiring uniform convergence in the supremum norm, as in zhang2023classification,he2025goodness. In fact, a sufficient condition for Assumption (ref) is \[ \sup_{x}\left|\widehat{m}^{(l)}(x)-m^{(l)}(x)\right|=o_{p}(n_{l}^{-1/4})\quad(l=1,2). \] However, achieving such a uniform convergence rate typically requires additional conditions, such as compact covariate support or covariate densities bounded away from zero.

theoremUnder $H_{0}$ in ((ref)) and Assumptions (ref)-(ref), we have \[ \left(\frac{\sigma_{1}^{2}}{n_{1}}+\frac{\sigma_{2}^{2}}{n_{2}}\right)^{-1/2}T\overset{d}{\longrightarrow}\mathcal{N}(0,1),\quad\text{as }n_{1},n_{2}\rightarrow\infty, \] where $\sigma_{l}^{2}=\mathbb{E}[\{\varepsilon^{(l)}\varepsilon^{(l)\prime}k(X^{(l)},X^{(l)\prime})\}^{2}]$ and $(\varepsilon^{(l)\prime},X^{(l)\prime})$ is an independent copy of $(\varepsilon^{(l)},X^{(l)})$.

It is worth noting that Theorem (ref) remains valid whether the dimension $p$ is fixed or diverges as $n_{1},n_{2}\rightarrow\infty$. A natural plug-in estimator of $\sigma_{l}^{2}$ is given by \[ \widehat{\sigma}_{l}^{2}=\frac{1}{n_{l}}\sum_{i=1}^{n_{l}}\left\{\widehat{\eta}^{(l)}_{i}\widehat{\eta}^{(l)}_{i+n_{l}}k(X^{(l)}_{i},X^{(l)}_{i+n_{l}})\right\}^{2}\quad(l=1,2), \] The following theorem shows that $\widehat{\sigma}_{l}^{2}$ is ratio-consistent under the null hypothesis. Then we reject $H_{0}$ at a significance level $\alpha$ if $(\widehat{\sigma}_{1}^{2}/n_{1}+\widehat{\sigma}_{2}^{2}/n_{2})^{-1/2}T>z_{1-\alpha}$.

theoremUnder $H_{0}$ in ((ref)) and Assumptions (ref)-(ref), we have $\widehat{\sigma}_{l}^{2}/\sigma_{l}^{2}\overset{p}{\longrightarrow} 1$ for $l=1,2$. Consequently, $(\widehat{\sigma}_{1}^{2}/n_{1}+\widehat{\sigma}_{2}^{2}/n_{2})^{-1/2}T\overset{d}{\longrightarrow}\mathcal{N}(0,1)$ as $n_{1},n_{2}\rightarrow\infty$.

remarkOne might consider constructing a test statistic based on the U-statistic estimator \[ \widehat{\Delta}^{(l)}_{U}=\frac{1}{n_{l}(n_{l}-1)}\sum_{1\le i\neq j\le n_{l}}\widehat{\eta}^{(l)}_{i}\widehat{\eta}^{(l)}_{j}k(X^{(l)}_{i},X^{(l)}_{j})\quad(l=1,2). \] We instead adopt $\widehat{\Delta}^{(l)}$ for two principal reasons. First, employing $\widehat{\Delta}^{(l)}_{U}$ requires strengthening the convergence rate in Assumption (ref) from $o_{p}(n_{l}^{-1/2})$ to $o_{p}(n_{l}^{-1})$—a condition that fails even under simple parametric models. Without this stronger rate, it becomes necessary to use a substantially smaller test set. For example, Assumption 3 and Theorem 1 in he2025goodness require the test set size to be of order $o(\sqrt{n_{l}})$. Second, the asymptotic distribution of $\widehat{\Delta}^{(l)}_{U}$ differs markedly between fixed- and high-dimensional regimes. When $p$ is fixed, $\widehat{\Delta}^{(l)}_{U}$ is a degenerate U-statistic whose null distribution is an infinite weighted sum of chi-squared variables, where the weights depend on the unknown underlying distribution. This necessitates a resampling calibration procedure such as the wild bootstrap lai2021kernel, which incurs considerable computational cost. In contrast, when $p$ diverges, a martingale CLT argument hall2014martingale yields the asymptotic normality of $\widehat{\Delta}^{(l)}_{U}$ under Assumption (ref) and the additional condition \[ \frac{\mathbb{E}\left(\left[\mathbb{E}\left\{h(Z_{1}^{(l)},Z_{2}^{(l)})h(Z_{1}^{(l)},Z_{3}^{(l)})\mid Z_{2}^{(l)},Z_{3}^{(l)}\right\}\right]^{2}\right)}{\left[\mathbb{E}\left\{h^{2}(Z_{1}^{(l)},Z_{2}^{(l)})\right\}\right]^{2}}\longrightarrow0,\quad\text{as }n_{l}\rightarrow\infty, \] where $Z_{i}^{(l)}=(\varepsilon_{i}^{(l)},X_{i}^{(l)})$ and $h(Z_{1}^{(l)},Z_{2}^{(l)})=\varepsilon^{(l)}_{1}\varepsilon^{(l)}_{2}k(X^{(l)}_{1},X^{(l)}_{2})$. Similar conditions appear in zhang2018conditional,li2023testing, but they are generally unverifiable without strong prior knowledge. By contrast, Theorem (ref) indicates that $\widehat{\Delta}^{(l)}$ is asymptotically normal under the null hypothesis in both fixed- and high-dimensional settings, without the need for resampling procedures or reliance on such practically unverifiable technical conditions.

Next, we study the asymptotic behavior of $T$ under the alternative hypothesis.

theoremFor $l=1,2$, assume that $\mathbb{E}[\{m^{(1)}(X^{(l)})-m^{(2)}(X^{(l)})\}^{2}]\le C_{3}<\infty$ for some $C_{3}>0$. Under $H_{a}$ in ((ref)) and Assumptions (ref)-(ref), we have \begin{align*} T&=\Delta^{(1)}+\Delta^{(2)}+O_{p}(n_{1}^{-1/2}+n_{2}^{-1/2})\\ &\quad+\|m^{(1)}-m^{(2)}\|_{L^{2}(P_{X}^{(2)})}o_{p}(n_{1}^{-1/4})+\|m^{(1)}-m^{(2)}\|_{L^{2}(P_{X}^{(1)})}o_{p}(n_{2}^{-1/4}). \end{align*} Thus, sufficient conditions for consistency, i.e., $(\widehat{\sigma}_{1}^{2}/n_{1}+\widehat{\sigma}_{2}^{2}/n_{2})^{-1/2}T\overset{p}{\longrightarrow}\infty$, are $\Delta^{(1)}+\Delta^{(2)}=\omega(n_{1}^{-1/2}+n_{2}^{-1/2})$ and $\Delta^{(1)}+\Delta^{(2)}=\Omega(\|m^{(1)}-m^{(2)}\|_{L^{2}(P_{X}^{(2)})}n_{1}^{-1/4}+\|m^{(1)}-m^{(2)}\|_{L^{2}(P_{X}^{(1)})}n_{2}^{-1/4})$.

remarkThe slower rates $o_{p}(n_{1}^{-1/4})$ and $o_{p}(n_{2}^{-1/4})$ under $H_{a}$ arise from the first-order errors of the machine learning estimators $\widehat{m}^{(l)}(\cdot)$. By contrast, under $H_{0}$, only second-order errors contribute to the asymptotics. Under parametric models, these rates improve to $O_{p}(n_{1}^{-1/2})$ and $O_{p}(n_{2}^{-1/2})$.

Although the test statistic $T$ is valid across all dimensional regimes, its power exhibits distinct behaviors in fixed- and high-dimensional settings. When $p$ is fixed and $\Delta^{(1)}+\Delta^{(2)}$ is a constant under $H_{a}$, the test is consistent, i.e., $(\widehat{\sigma}_{1}^{2}/n_{1}+\widehat{\sigma}_{2}^{2}/n_{2})^{-1/2}T\overset{p}{\longrightarrow}\infty$. In contrast, when $p\rightarrow\infty$, the quantity $\Delta^{(1)}+\Delta^{(2)}$ depends on $p$, and the power is determined jointly by the dimension and the sample sizes. To illustrate, consider the Gaussian kernel $k(x,x')=\exp(-(2\gamma^{2})^{-1}\|x-x'\|^{2})$, where $\gamma$ is the bandwidth parameter. Following yan2023kernel,han2024generalized, we assume $\gamma^{2}\asymp p$ and impose the following condition.

assumptionLet $X^{(1)}=\mu_{1}+\Sigma_{1}^{1/2}U$ and $X^{(2)}=\mu_{2}+\Sigma_{2}^{1/2}V$, where $\mu_{l}$ and $\Sigma_{l}$ denote the mean vector and covariance matrix of $X^{(l)}$, respectively. The random vectors $U,V\in\mathbb{R}^{p}$ have independent components with mean zero, variance one and finite fourth moment. Moreover, the spectrum of $\Sigma_{l}$ lies within $[M^{-1},M]$ for some $M>1$.

Assumption (ref) is standard in high-dimensional statistics and random matrix theory. We emphasize that this assumption is introduced solely for theoretical analysis and is not required for the implementation of our test. The bandwidth order $\gamma^{2}\asymp p$ is consistent with the widely used median heuristic ramdas2015adaptivity.

theoremFor $l=1,2$, let $\tau_{l}=\operatorname{tr}(\Sigma_{l})/\gamma^{2}$, where $\operatorname{tr}(\cdot)$ denotes the trace. Under $H_{a}$ in ((ref)) and Assumption (ref), we have \[ \Delta^{(1)}+\Delta^{(2)}=e^{-\tau_{1}}\left\{\mathbb{E}(\eta^{(1)})\right\}^{2}+e^{-\tau_{2}}\left\{\mathbb{E}(\eta^{(2)})\right\}^{2}+\left\{\mathbb{E}(|\eta^{(1)}|^{2})+\mathbb{E}(|\eta^{(2)}|^{2})\right\}O(p^{-1/2}). \] When $\mathbb{E}(\eta^{(l)})=0$ for $l=1,2$, we further have \begin{align*} \Delta^{(1)}+\Delta^{(2)}=\frac{e^{-\tau_{1}}}{\gamma^{2}}\left\|\operatorname{cov}(\eta^{(1)},X^{(1)})\right\|^{2}+\frac{e^{-\tau_{2}}}{\gamma^{2}}\left\|\operatorname{cov}(\eta^{(2)},X^{(2)})\right\|^{2}+\left\{\mathbb{E}(|\eta^{(1)}|^{2})+\mathbb{E}(|\eta^{(2)}|^{2})\right\}O(p^{-1}). \end{align*}

Combining Theorems (ref) and (ref), we obtain the following conclusion. If \[ \max\left\{\left|\mathbb{E}(\eta^{(1)})\right|,\left|\mathbb{E}(\eta^{(2)})\right|\right\}\ge c_{3}, \] for some $c_{3}>0$, the proposed test is consistent regardless of the growth rate of the dimension $p$. By contrast, when $\mathbb{E}(\eta^{(l)})=0$ for both $l=1,2$, the test statistic $T$ primarily detects alternatives satisfying $\operatorname{cov}(\eta^{(l)},X^{(l)})\neq 0$ for some $l=1,2$, as $p\rightarrow\infty$. To capture the nonlinear dependence of $\eta^{(l)}$ on $X^{(l)}$ in high-dimensional settings, we introduce an alternative test statistic in Section (ref).

remarkBy the definition of $\eta^{(l)}\ (l=1,2)$, we have \begin{align*} \mathbb{E}(\eta^{(1)})&=\mathbb{E}\left\{Y^{(1)}-m^{(2)}(X^{(1)})\right\}=\mathbb{E}\left\{m^{(1)}(X^{(1)})-m^{(2)}(X^{(1)})\right\},\\ \mathbb{E}(\eta^{(2)})&=\mathbb{E}\left\{Y^{(2)}-m^{(1)}(X^{(2)})\right\}=\mathbb{E}\left\{m^{(2)}(X^{(2)})-m^{(1)}(X^{(2)})\right\}. \end{align*} In the simplest case where both $m^{(1)}(\cdot)$ and $m^{(2)}(\cdot)$ are linear, the equalities $\mathbb{E}(\eta^{(1)})=\mathbb{E}(\eta^{(2)})=0$ hold under $H_{a}$ provided that $\mu_{1}=\mu_{2}=0$ and the two models share the same intercept. Generally, however, the regression functions are often highly nonlinear in practice. Therefore, under the alternative hypothesis, it is uncommon for both $\mathbb{E}(\eta^{(1)})$ and $\mathbb{E}(\eta^{(2)})$ to vanish simultaneously, especially in the presence of covariate shift (i.e., $P_{X}^{(1)}\neq P_{X}^{(2)}$).

Alternative test statistic in high dimensions

For $l=1,2$, write $X^{(l)}=(X^{(l)}(1),\ldots,X^{(l)}(p))^{\top}$ and $X^{(l)}_{i}=(X^{(l)}_{i}(1),\ldots,X^{(l)}_{i}(p))^{\top}\ (i=1,\ldots,N_{l})$. Given a kernel $k$ on $\mathbb{R}$, define \[ \Delta^{(l)}_{d}=\text{KCMD}^{*}(\eta^{(l)}\mid X^{(l)}(d))=\mathbb{E}\left\{\eta^{(l)}\eta^{(l)\prime}k(X^{(l)}(d),X^{(l)\prime}(d))\right\}\quad(l=1,2;\ d=1,\ldots,p), \] where $(\eta^{(l)\prime},X^{(l)\prime})$ is an independent copy of $(\eta^{(l)},X^{(l)})$.

Now we propose an alternative test statistic for detecting nonlinear conditional mean dependence of $\eta^{(l)}$ on $X^{(l)}$ in high dimensions: \[ T_{a}=\sum_{l=1}^{2}\sum_{d=1}^{p}\widehat{\Delta}^{(l)}_{d}, \] where \[ \widehat{\Delta}^{(l)}_{d}=\frac{1}{n_{l}}\sum_{i=1}^{n_{l}}\widehat{\eta}^{(l)}_{i}\widehat{\eta}^{(l)}_{i+n_{l}}k(X^{(l)}_{i}(d),X^{(l)}_{i+n_{l}}(d)), \] and $\widehat{\eta}^{(l)}_{i}$ is defined as in Section (ref).

In contrast to ((ref)), which is equivalent to the null hypothesis $H_{0}$, the test statistic $T_{a}$ targets the weaker null hypothesis

equation[equation omitted — 129 chars of source]

Similar idea was first introduced by zhang2018conditional and later adopted by li2023testing. This relaxation is necessary in high-dimensional regimes, where the alternative space is vast due to both the growing dimension and the possibility of nonlinear dependence.

remarkOur problem differs entirely from that of zhang2018conditional,li2023testing. Furthermore, our goal is to test whether $\mathbb{E}(\eta^{(l)}\mid X^{(l)})=0$, rather than whether $\mathbb{E}(\eta^{(l)}\mid X^{(l)})=\mathbb{E}(\eta^{(l)})$. To illustrate, suppose that $m^{(1)}(\cdot)-m^{(2)}(\cdot)\equiv a$ for some $a\neq 0$; that is, the two regression functions differ only by a constant shift. In this case, we have $\mathbb{E}(\eta^{(l)}\mid X^{(l)})\neq 0$ but $\mathbb{E}(\eta^{(l)}\mid X^{(l)})=\mathbb{E}(\eta^{(l)})$.

We next establish the asymptotic properties of $T_{a}$. To this end, Assumption (ref) is reformulated as follows.

assumptionp{assu2} Define \[ G(X^{(l)},X^{(l)\prime})=\sum_{d=1}^{p}k(X^{(l)}(d),X^{(l)\prime}(d))\quad(l=1,2). \] For $l=1,2$, assume that $\mathbb{E}(|\varepsilon^{(l)}|^{4})\le C_{1}<\infty$ for some $C_{1}>0$, and \begin{equation} \mathbb{E}\left\{\left|\varepsilon^{(l)}\varepsilon^{(l)\prime}G(X^{(l)},X^{(l)\prime})\right|^{2}\right\}\ge c_{4}p^{2}, \end{equation} for some $c_{4}>0$, where $(\varepsilon^{(l)\prime},X^{(l)\prime})$ denotes an independent copy of $(\varepsilon^{(l)},X^{(l)})$.

Condition ((ref)) holds automatically when $p$ is fixed. For diverging $p$, a sufficient condition for ((ref)) is that \[ \min_{1\le d_{1},d_{2}\le p}\mathbb{E}\left\{\left|\varepsilon^{(l)}\varepsilon^{(l)\prime}\right|^{2}k(X^{(l)}(d_{1}),X^{(l)\prime}(d_{1}))k(X^{(l)}(d_{2}),X^{(l)\prime}(d_{2}))\right\}\ge c_{4}, \] for some $c_{4}>0$.

theoremUnder $H_{0}$ in ((ref)) and Assumptions (ref), (ref) and (ref), we have \[ \left(\frac{\xi_{1}^{2}}{n_{1}}+\frac{\xi_{2}^{2}}{n_{2}}\right)^{-1/2}T_{a}\overset{d}{\longrightarrow}\mathcal{N}(0,1),\quad\text{as }n_{1},n_{2}\rightarrow\infty, \] where $\xi_{l}^{2}=\mathbb{E}[\{\varepsilon^{(l)}\varepsilon^{(l)\prime}G(X^{(l)},X^{(l)\prime})\}^{2}]$ and $(\varepsilon^{(l)\prime},X^{(l)\prime})$ is an independent copy of $(\varepsilon^{(l)},X^{(l)})$.

remarkHeuristically, the statistic $T_{a}$ can be regarded as a version of the statistic $T$ with kernel $G$ defined in Assumption (ref). Accordingly, the theoretical analysis proceeds in a manner analogous to that presented in Section (ref).

A natural plug-in estimator of $\xi_{l}^{2}$ is given by \[ \widehat{\xi}_{l}^{2}=\frac{1}{n_{l}}\sum_{i=1}^{n_{l}}\left\{\widehat{\eta}^{(l)}_{i}\widehat{\eta}^{(l)}_{i+n_{l}}\sum_{d=1}^{p}k(X^{(l)}_{i}(d),X^{(l)}_{i+n_{l}}(d))\right\}^{2}\quad(l=1,2). \] The following theorem shows that $\widehat{\xi}_{l}^{2}$ is ratio-consistent under $H_{0}$. Then we reject $H_{0}$ at a significance level $\alpha$ if $(\widehat{\xi}_{1}^{2}/n_{1}+\widehat{\xi}_{2}^{2}/n_{2})^{-1/2}T_{a}>z_{1-\alpha}$.

theoremUnder $H_{0}$ in ((ref)) and Assumptions (ref), (ref) and (ref), we have $\widehat{\xi}_{l}^{2}/\xi_{l}^{2}\overset{p}{\longrightarrow} 1$ for $l=1,2$. As a result, $(\widehat{\xi}_{1}^{2}/n_{1}+\widehat{\xi}_{2}^{2}/n_{2})^{-1/2}T_{a}\overset{d}{\longrightarrow}\mathcal{N}(0,1)$ as $n_{1},n_{2}\rightarrow\infty$.

Although the statistic $T_{a}$ is initially designed for high-dimensional settings, Theorems (ref) and (ref) remain valid whether the dimension $p$ is fixed or diverging. In the fixed-dimensional case, however, the test based on $T$ introduced in Section (ref) is preferable, as it is consistent against all alternatives, whereas $T_{a}$ can merely detect alternatives for which the weaker null hypothesis ((ref)) is violated.

theoremFor $l=1,2$, assume that $\mathbb{E}[\{m^{(1)}(X^{(l)})-m^{(2)}(X^{(l)})\}^{2}]\le C_{3}<\infty$ for some $C_{3}>0$. Under $H_{a}$ in ((ref)) and Assumptions (ref), (ref) and (ref), we have \begin{align*} T_{a}&=\sum_{d=1}^{p}\Delta^{(1)}_{d}+\sum_{d=1}^{p}\Delta^{(2)}_{d}+O_{p}(n_{1}^{-1/2}p+n_{2}^{-1/2}p)\\ &\quad+\|m^{(1)}-m^{(2)}\|_{L^{2}(P_{X}^{(2)})}o_{p}(n_{1}^{-1/4}p)+\|m^{(1)}-m^{(2)}\|_{L^{2}(P_{X}^{(1)})}o_{p}(n_{2}^{-1/4}p). \end{align*} Hence, sufficient conditions for consistency, i.e., $(\widehat{\xi}_{1}^{2}/n_{1}+\widehat{\xi}_{2}^{2}/n_{2})^{-1/2}T_{a}\overset{p}{\longrightarrow}\infty$, are $\sum_{d=1}^{p}\Delta^{(1)}_{d}+\sum_{d=1}^{p}\Delta^{(2)}_{d}=\omega(n_{1}^{-1/2}p+n_{2}^{-1/2}p)$ and $\sum_{d=1}^{p}\Delta^{(1)}_{d}+\sum_{d=1}^{p}\Delta^{(2)}_{d}=\Omega(\|m^{(1)}-m^{(2)}\|_{L^{2}(P_{X}^{(2)})}n_{1}^{-1/4}p+\|m^{(1)}-m^{(2)}\|_{L^{2}(P_{X}^{(1)})}n_{2}^{-1/4}p)$.

When $\max\{|\mathbb{E}(\eta^{(1)})|,|\mathbb{E}(\eta^{(2)})|\}\ge c_{3}$ for some $c_{3}>0$, then $\sum_{d=1}^{p}\Delta^{(1)}_{d}+\sum_{d=1}^{p}\Delta^{(2)}_{d}\asymp p$. Consequently, similar to the test based on $T$, the test constructed from $T_{a}$ remains consistent regardless of the growth rate of the dimension $p$. Moreover, when the dimension $p$ increases and $\mathbb{E}(\eta^{(l)})=0$ for both $l=1,2$, the statistic $T_{a}$, through the quantity $\sum_{d=1}^{p}\Delta^{(1)}_{d}+\sum_{d=1}^{p}\Delta^{(2)}_{d}$, is able to capture componentwise nonlinear dependence between $\eta^{(l)}$ and $X^{(l)}$, in contrast to the test based on $T$.

remarkIn practice, without prior knowledge of the underlying dependence structure, both tests can be applied irrespective of the dimensionality. The extensive numerical studies demonstrate the reliable performance of both test statistics.

Numerical studies

Monte Carlo simulations

In this section, we conduct several Monte Carlo simulations to evaluate the finite-sample performance of our proposed tests, denoted by $T$, the statistic from Section (ref), and $T_a$, the alternative statistic from Section (ref). We assess their empirical size and power across various settings, including low-dimensional and high-dimensional covariates, and different forms of nonlinear regression functions. All results are based on 1000 replications for each setting, and the nominal significance level is set to $\alpha = 0.05$.

For our proposed methods, we employ XGBoost chen2016xgboost to estimate the regression functions $m^{(1)}(\cdot)$ and $m^{(2)}(\cdot)$. We opt for XGBoost owing to its strong predictive performance across diverse data structures, its computational efficiency, and its relative simplicity in tuning compared with more complex methods like neural networks, which typically require larger sample sizes and more extensive hyperparameter optimization. We set a learning rate of $0.15$ and a maximum tree depth of $3$. To prevent both underfitting and overfitting, the optimal number of boosting iterations is determined using 5-fold cross-validation. At the same time, we incorporate an early stopping mechanism, which terminates the training process if the validation performance does not improve for 20 consecutive rounds, with an overall limit of $1000$ iterations. For both test statistics, we use a Gaussian kernel, $k(x, x') = \exp(-(2\gamma^2)^{-1} \| x - x' \|^2)$. The bandwidth parameter $\gamma$ is selected via the median heuristic gretton2012kernel, and is computed separately for each sample. Specifically, for the statistic $T$, a scalar bandwidth is calculated for each sample based on the median of pairwise Euclidean distances within that sample. For the statistic $T_a$, a dimension-specific bandwidth vector is computed for each sample, where each component is derived from the median of pairwise distances along the corresponding covariate dimension within that sample.

Across all examples, we set the sample sizes to be equal, $N_1 = N_2 = N$, and consider different values of $N$ in the simulations. For each group $l \in \{1, 2\}$, we consider the regression model

equation*[equation* omitted — 68 chars of source]

The covariates for the first group are drawn from a standard normal distribution, $X^{(1)} \sim \mathcal{N}(0, I_p)$. To introduce a covariate shift, the covariates for the second group are drawn from a normal distribution with an autoregressive covariance structure, $X^{(2)} \sim \mathcal{N}(0, \Sigma)$, where $\Sigma_{ij} = 0.3^{|i - j|}$. The errors are generated as $\varepsilon^{(1)} \sim \mathcal{N}(0, 0.5^2)$ and $\varepsilon^{(2)}$ from a t-distribution with $5$ degrees of freedom, scaled to have a variance of $0.5^2$. The difference between the two regression functions is controlled by a vector $\beta\in\mathbb{R}^{p}$, where the null hypothesis $H_0$ corresponds to $\|\beta\| = 0$.

exampleWrite $x=(x(1),\ldots,x(p))^{\top}\in\mathbb{R}^{p}$. Let \begin{align*} m^{(1)}(x) &= \sqrt{\{x(1)\}^2 + \{x(2)\}^2 + \sum_{d=1}^{p} \beta_{d} \{x(d)\}^2}, \\ m^{(2)}(x) &= \sqrt{\{x(1)\}^2 + \{x(2)\}^2}. \end{align*}

exampleLet \begin{align*} m^{(1)}(x) &= \exp\{ x(1) \} + \exp\{ x(2) \} + \sum_{d=1}^{p} \beta_d \exp\{ x(d) \}, \\ m^{(2)}(x) &= \exp\{ x(1) \} + \exp\{ x(2) \}. \end{align*}

exampleLet \begin{align*} m^{(1)}(x) &= x(1) + x(2) + \left\{ \sum_{d=1}^{p} \beta_d x(d) \right\}^2, \\ m^{(2)}(x) &= x(1) + x(2). \end{align*}

For each example, we consider a low-dimensional setting with $p = 5$ and a high-dimensional setting with $p = 50$. In the low-dimensional case, all the elements of signal vector $\beta$ are non-zero and have equal magnitude. In the high-dimensional case, we investigate two types of alternatives: a sparse alternative, where only the first two elements of $\beta$ are non-zero with equal magnitude, and a dense alternative, where the first $20$ elements of $\beta$ are non-zero, also with equal magnitude.

For the low-dimensional case ($p = 5$), our tests are compared with the smoothing-based procedure of lavergne2001equality. Consistent with their simulation study, a uniform kernel is employed, and the bandwidth is selected via the rule-of-thumb method.

The simulation results are presented in Tables (ref)--(ref). Under the null hypothesis (i.e., $\|\beta\| = 0$), the empirical sizes of our proposed tests, $T$ and $T_a$, are consistently close to the $5\%$ nominal level across all settings. This indicates that our tests maintain proper size control even in the presence of covariate shift, non-Gaussian errors, and high dimensionality.

As expected, the power of both tests, $T$ and $T_a$, increases with the signal strength $\|\beta\|$ and the sample size $N$, approaching one as either becomes sufficiently large. In high-dimensional settings ($p=50$), both procedures perform well against both sparse and dense alternatives. By contrast, even in low-dimensional settings ($p=5$), the test of lavergne2001equality is overly conservative, exhibiting an empirical size of zero across all scenarios. As a result, it displays no empirical power and fails to detect the nonlinear alternatives in any of the three examples.

table[table omitted — 2,994 chars of source]

Airfoil data example

We apply our proposed tests to the Airfoil Self-Noise dataset airfoil, available from the UCI Machine Learning Repository. This dataset, which has also been used in previous studies such as tibshirani2019conformal,hu2024two, includes $1503$ measurements from wind tunnel experiments on NASA airfoil sections. The response is the scaled sound pressure level, and there are five covariates: frequency, angle of attack, chord length, free-stream velocity, and suction side displacement thickness. Following common practice, we apply a logarithmic transformation to the frequency and thickness variables.

Since the original dataset does not have a predefined two-sample structure, we construct two scenarios to evaluate our tests. First, to assess the empirical size, we simulate the null hypothesis by randomly partitioning the entire dataset into two nearly equal groups with sizes $N_1=752$ and $N_2=751$. Under this setup, the underlying regression functions for the two groups are expected to be identical. We repeat this procedure $500$ times to estimate the Type I error rate. Second, to evaluate the power, we construct an alternative hypothesis by partitioning the data based on the response. Specifically, one group comprises observations with sound pressure levels at or below the median, while the other contains the remaining observations. This procedure intentionally creates two groups with distinct regression functions. To avoid trivial separation of the covariate supports, we randomly swap $5\%$ of the samples between the two groups.

The results of our analysis demonstrate the effectiveness of the proposed tests. Under the null hypothesis, where the data is randomly partitioned, the empirical rejection rates for the statistics $T$ and $T_a$ are $4.6\%$ and $5.8\%$, respectively. Both rates are close to the nominal $5\%$ level, confirming that our tests maintain proper size control on this real-world dataset. Under the alternative hypothesis, constructed by splitting the data based on the response, both tests strongly reject the null hypothesis. The test statistic $T$ yields a $p$-value of $1.574 \times 10^{-19}$, while the statistic $T_a$ gives a $p$-value of $4.726 \times 10^{-17}$.

Communities and Crime data example

We further validate our tests on the Communities and Crime dataset communities_and_crime, also from the UCI Machine Learning Repository. This dataset, which has also been used by joshi2025conformal, combines socio-economic data from the 1990 US Census, law enforcement data from the 1990 US LEMAS survey, and crime data from the 1995 FBI Uniform Crime Report. The dataset contains 1994 communities and 128 attributes. The response is the per capita violent crime rate, and covariates include information like population, median income, and employment rates. After removing non-predictive attributes and any covariates containing missing values, we proceed with the analysis. This preprocessing results in a final covariate dimension of $p = 99$.

To create a two-sample testing problem, we again construct two scenarios following the same logic as in the Airfoil data analysis. To evaluate the empirical size, we create a null scenario by randomly splitting the dataset into two samples of equal size, a procedure repeated 500 times. To assess the power, we generate an alternative hypothesis by dividing the communities into two groups based on their crime rates (at or below vs. above the median) and then swapping $5\%$ of the data points between the groups to ensure overlapping covariate distributions.

Our tests perform well in this setting. Under the null hypothesis, the empirical rejection rates are $5.6\%$ for the statistic $T$ and $5.2\%$ for the statistic $T_a$, both of which are close to the nominal $5\%$ level, demonstrating accurate size control. Under the alternative hypothesis, the tests show decisive power, yielding $p$-values of $2.128 \times 10^{-76}$ for $T$ and $9.260 \times 10^{-83}$ for $T_a$.

Discussion

To conclude, we outline several directions for future research. First, it would be valuable to extend the framework to the comparison of $L$ regression functions for a general $L\ge 2$, or even for $L\rightarrow\infty$. Second, beyond the conditional mean, it would be of interest to employ machine learning methods to test the equality of two conditional variance functions pardo2015tests, or more generally, two conditional distributions yan2022distance. Finally, our framework may find utility in other contexts, such as change point detection nie2022detection.