Bias correction for Chatterjee's graph-based correlation coefficient

Bias correction for Chatterjee's graph-based correlation coefficient

{5pt} {5pt} {5pt} {5pt} \hypersetup{colorlinks,breaklinks,urlcolor=blue,linkcolor=blue}

{\bf Keywords}: measure of dependence, nearest neighbor graph, graph-based statistics, regression adjustment, nonparametric ridge regression.

Introduction

Sourav Chatterjee, in his groundbreaking paper chatterjee2020new, introduced a novel and elegant approach to estimating the dependence measure of MR3024030. This measure captures both independence and perfect functional dependence: it equals 0 if and only if the two random scalars are independent, and 1 if and only if one is a measurable function of the other almost surely. Chatterjee’s work led to a new rank-based statistic, now widely known as Chatterjee’s rank correlation.

Among the notable extensions of this work is the contribution by Azadkia21simple, who generalized Chatterjee’s rank correlation to multivariate settings through a sophisticated integration of rank statistics and nearest neighbor (NN) graphs. However, unlike Chatterjee’s original rank correlation—which has been shown to be root-$n$ consistent and asymptotically normal under broad conditions lin2022limit,kroll2024asymptotic—the NN-based correlation coefficient of Azadkia21simple suffers from a bias term that may decay at a rate slower than root-$n$, rendering the estimator not root-$n$ consistent in general.

This article examines the bias of the NN-based correlation coefficient in greater depth and contributes in two main directions. First, building on the thought-provoking work of viel2025convergenceratenearestneighbour, which analyzed the bias of another NN-based estimator in causal inference, we show that the bias of the NN-based correlation coefficient could become asymptotically negligible when the dimension $d \le 3$. This finding advances the current understanding, which establishes analogous asymptotic negligibility only in the univariate case $d=1$.

Second, for more general scenarios in which the bias may no longer be ignorable, we propose a regression-adjustment technique and derive the large-sample properties of the resulting bias-corrected estimator of Chatterjee’s NN graph-based correlation coefficient. Our approach is inspired, once again, by related developments in causal inference, including the influential ideas of rubin1973use and the theoretical foundations laid out in abadie2011bias, lin2021estimation, lin2022regression, and cattaneo2024rosenbaum. Leveraging these insights, we develop a new theory for nonparametric ridge regression and demonstrate that the proposed correction can asymptotically eliminate the bias of the original NN estimator, thereby restoring root-$n$ consistency without affecting the limiting variance.

In both directions, the original NN-based estimator of Azadkia21simple (when $d \le 3$) and the newly proposed bias-corrected estimator (for more general settings) share the same attractive property: each achieves root-$n$ consistency and asymptotic normality with an identical limiting variance, regardless of the dimensionality of the covariates. This robustness allows seamless integration with both analytical variance estimators lin2022limit and bootstrap-based procedures dette2025simple, thereby facilitating valid inference.

{\bf Notation.} In the following, for any two real numbers $a, b$, we write $a \wedge b := \min\{a,b\}$ and $a \vee b := \max\{a,b\}$. For any two real sequences $\{a_n\}_{n\geq 1}$ and $\{b_n\}_{n\geq 1}$, we denote $a_n= O(b_n)$ (or equivalently, $a_n\lesssim b_n$) if there exists a constant $C>0$ such that $|a_n|\le C|b_n|$ for all sufficiently large $n$, and $a_n= o(b_n)$ if $|a_n|/|b_n|\to 0$ as $n\to \infty$. For any two real sequences $\{a_n\}, \{b_n\}$, we write $a_n \asymp b_n$ if both $a_n\lesssim b_n$ and $b_n\lesssim a_n$ hold true. We write $O_{\mathrm{P}}$ and $o_{\mathrm{P}}$ if $O(\cdot)$ and $o(\cdot)$ hold in probability, respectively.

Chatterjee's NN graph-based correlation coefficient

We begin with a brief review of the current understanding of Azadkia21simple's NN graph-based correlation coefficient. Consider $(X,Y) \in \mathbb{R}^d \times \mathbb{R}$ to be a pair of random variables. To quantify the strength of dependence between $X$ and $Y$, MR3024030 introduced the following population quantity, now referred to as the Dette–Siburg–Stoimenov dependence measure: \[ T = T(Y,X) := \frac{\int \mathrm{Var}\big(\mathrm{E}[\mathbbm{1}(Y \geq t) \mid X]\big) \mathrm{d} \mu(t)}{\int \mathrm{Var}\big(\mathbbm{1}(Y \geq t)\big) \mathrm{d} \mu(t)}, \] where $\mu$ denotes the law of $Y$ and $\mathbbm{1}(\cdot)$ is the indicator function. As shown in MR3024030 and further developed in Azadkia21simple, $T$ satisfies the desirable properties of being

Let $(X_1,Y_1),\ldots,(X_n, Y_n)$ be the sample. For each $i \in [n] := \{1,2,\ldots,n\}$, define \[ R_i := \sum_{j=1}^n \mathbbm{1}(Y_j \leq Y_i), \quad \text{and} \quad N_i := \argmin_{j \ne i} \|X_i - X_j\|, \] where $R_i$ is the rank of $Y_i$ among $\{Y_1, \ldots, Y_n\}$ and $N_i$ indexes the nearest neighbor (NN) of $X_i$ under the Euclidean metric, $\|\cdot\|$. To estimate $T$ based on the sample $\{(X_i, Y_i)\}_{i=1}^n$, Azadkia21simple proposed the following NN-based rank correlation coefficient: \[ \hat T_n = \hat{T}_n(Y,X) := \frac{6}{n^2 - 1} \sum_{i=1}^n \min\big\{R_i, R_{N(i)}\big\} - \frac{2n + 1}{n - 1}. \]

The asymptotic properties of $\hat T_n$ are summarized in the following theorem under the following assumption of the data generating scheme.

Theorem (ref) characterizes the asymptotic behavior of $\hat T_n$. In particular, Part (ref) highlights the main obstacle preventing $\hat T_n$ from achieving root-$n$ consistency: the presence of a theoretically non-negligible bias term. This is not unexpected for researchers familiar with NN-based statistics. Indeed, bias correction has played a central role in related areas such as NN graph-based estimation for causal inference abadie2011bias, lin2021estimation and entropy estimation MR3909934. In the next section, we will present our approach to addressing this final gap for $\hat T_n$.

Beyond the above theoretical results, it is worth acknowledging ongoing efforts in this research direction, including deeper investigations of $T$ bucher2024lack, further analyses of $\hat T_n$ shi2020power, bickel2022measures, lin2024failure, han2024azadkia, auddy2021exact, zhang2024asymptotic, exploration of new dependence measures satisfying the defining properties of $T$ griessenberger2022multivariate, ansari2025directextensionazadkia, strothmann2024rearranged, azadkia2025new, and extensions of $\hat T_n$-type ideas to broader statistical contexts huang2020kernel, azadkia2021fast, gamboa2022global, ansari2023quantifying, gao2024family, fuchs2024quantifying, tran2024rank, zhang2025extensions. More recent developments include zhou2025association and olivares2025powerful. For broader overviews and related perspectives, we refer the reader to the surveys conducted by han2021extensions and chatterjee2022survey.

Bias analysis

We begin with a closer analysis of the bias of $\hat T_n$. To this end, we first recall the closed-form expression for the bias term $\mathrm{E}[\hat T_n]-T$, as derived in lin2022limit.

Define the leading bias term in (ref) by \[ L \coloneqq \mathrm{E}\Big[\mathbbm{1}\!\big(Y_2 \le Y_1 \wedge Y_{N(1)}\big)\Big] - \mathrm{E}\Big[\mathbbm{1}\!\big(Y_2 \le Y_1 \wedge \tilde Y_1\big)\Big]. \] Lemma (ref) then directly yields \[ \mathrm{E}[\hat T_n] - T = 6L + O(n^{-1}). \] Thus, analyzing the bias reduces to understanding the behavior of $L$.

To this end, introduce the conditional mean function \[ G_{x}(t) \coloneqq \mathrm{E}\big[\mathbbm{1}(Y \ge t)\mid X=x\big]. \] A simple calculation shows that \[ L = \mathrm{E}\!\Big[F_Y\!\big(Y_1 \wedge Y_{N(1)}\big) - F_Y\!\big(Y_1 \wedge \tilde Y_1\big)\Big] = \int \mathrm{E}\!\left[G_{X_1}(t)\,G_{X_{N(1)}}(t) - G_{X_1}^2(t)\right] \mathrm{d}\mu(t), \] highlighting the central role played by the conditional mean function $G_x(\cdot)$ in the analysis of $\mathrm{E}[\hat T_n]-T$. This is a key insight already emphasized in the original paper Azadkia21simple.

Equation (ref) is a classical representation of the bias of NN-based functional estimators. As shown in Azadkia21simple (cf.\ Theorem (ref)(ref)), a direct analysis of (ref) implies that the resulting estimator converges at a rate slower than $n^{-1/2}$ whenever $d \ge 2$. In the broader literature, singh2016finite were the first to suggest that such functional NN estimators may nevertheless achieve parametric convergence rates in multivariate settings, provided the underlying density is sufficiently smooth, even though the NN density estimator itself generally does not adapt automatically to such smoothness zhao2022analysis. Subsequently, and more closely related to our work, viel2025convergenceratenearestneighbour analyzed the bias of an NN-based estimator in causal inference and demonstrated that the phenomenon identified by singh2016finite extends to this class of NN estimators as well.

In this section, we follow the argument of viel2025convergenceratenearestneighbour and show that the bias term in (ref) can also be asymptotically negligible. To prepare for the analysis, we introduce some additional notation. For a twice continuously differentiable function \[ f: U \to \mathbb{R}, \qquad U \subseteq \mathbb{R}^d \ \text{open}, \] we denote its gradient and Hessian at $x$ by $\nabla f(x)$ and $\nabla^{2} f(x)$, respectively. Let $\mathcal{X} \coloneqq {\rm supp}(X)$ denote the support of $X$, and define \[ \delta(x) \coloneqq

\] where \[ \mathrm{dist}(x,\partial \mathcal{X}) = \inf_{z \in \partial \mathcal{X}} \|x - z\| \] denotes the Euclidean distance from $x$ to the boundary $\partial \mathcal{X}$ of $\mathcal{X}$. Finally, for any positive integer $b$, we use $\mathcal{H}^b$ to denote the $b$-dimensional Hausdorff measure.

The following assumption regulates the behavior of $X$ and imposes smoothness conditions on the conditional mean function $G_x(t)$ with respect to $x$. Similar conditions appear in viel2025convergenceratenearestneighbour.

Assumption (ref) yields the following result, which establishes the asymptotic negligibility of the bias even when $d>1$.

Theorem (ref) provides a convergence rate for $\mathrm{E}[\hat{T}_n]-T$ that is typically sufficient for establishing the central limit theorem (ref). In particular, it shows that—under the stated conditions—the bias is $o(n^{-1/2})$ whenever $d \le 3$. Under additional smoothness assumptions, however, one can derive a more refined expansion of the bias and show that the above rate is, in general, unimprovable and therefore sharp. Such refinements follow from a higher-order expansion of the bias term, paralleling the arguments of viel2025convergenceratenearestneighbour.

Bias correction

Method

Section (ref) outlines the landscape of the bias of $\hat T_n$ and provides both positive and negative insights. On the positive side, Theorem (ref) shows that the bias is asymptotically negligible whenever $d \le 3$. On the other hand, Theorem (ref) demonstrates that, in general, the bias cannot converge at a rate faster than that established in Theorem (ref). Consequently, a bias correction becomes necessary whenever $d \ge 4$.

As a matter of fact, Lemma (ref) naturally motivates a bias-corrected version of $\hat T_n$. Recall that \[ \mathrm{E}[\hat T_n] - T = 6L + O(n^{-1}). \] To obtain a root-$n$ consistent estimator of $T$, it suffices to construct an estimator $\hat L_n$ such that \[ \hat L_n - L = o_{\mathrm{P}}(n^{-1/2}), \] in which case \[ n^{1/2}(\hat T_n - 6\hat L_n - T) = n^{1/2}(\hat T_n - \mathrm{E}[\hat T_n]) + o_{\mathrm{P}}(1). \]

Toward this goal, recall that \[ L = \mathrm{E}\!\left[F_Y\!\big(Y_1 \wedge Y_{N(1)}\big) - F_Y\!\big(Y_1 \wedge \tilde Y_1\big)\right] = \int \mathrm{E}\!\left[G_{X_1}(t)\, G_{X_{N(1)}}(t) - G_{X_1}^2(t)\right] \mathrm{d} \mu(t), \] suggesting that a natural estimator of $L$ may be obtained by replacing the integral, expectations, and the conditional mean function $G_x(\cdot)$ with their empirical counterparts.

The procedure then proceeds as follows.

The resulting {\it bias-corrected} NN graph-based correlation coefficient is then given by \[ \hat T_n^{\mathrm{bc}} := \hat T_n - 6 \hat L_n. \]

Theory

Before stating the main result of this section, we introduce additional notation. For any function $f: \mathbb{R}^d \to \mathbb{R}$, define its $L^{\infty}$ norm over the support of $X$ as \[ \norm{f}_{\infty} := \sup_{x \in {\rm supp}(F_X)} |f(x)|, \] where $F_X$ denotes the CDF of $X$ and ${\rm supp}(\cdot)$ denotes the support of a distribution or law. For any positive integer $q$, let $\Lambda_q$ denote the set of multi-indices of total degree $q$: \[ \Lambda_q := \left\{ \alpha = (\alpha_1, \ldots, \alpha_d) : \sum_{i=1}^d \alpha_i = q,\ \alpha_i \in \mathbb{Z}^{\geq 0}~~\text{for all }i \in [d] \right\}, \] where $\mathbb{Z}^{\geq 0}$ is the set of nonnegative integers. For each $\alpha \in \Lambda_q$, we adopt the following multi-index notation for derivatives: \[ D^\alpha f = \frac{\partial^{|\alpha|} f}{\partial x^\alpha} = \frac{\partial^{\alpha_1 + \cdots + \alpha_d} f}{\partial x_1^{\alpha_1} \cdots \partial x_d^{\alpha_d}},\quad |\alpha| = \sum_{i=1}^d \alpha_i,\quad \alpha! = \prod_{i=1}^d \alpha_i!,\quad x^\alpha = \prod_{i=1}^d x_i^{\alpha_i}. \] For any positive integer $K$, let $I_K$ denote the identity matrix of dimension $K$. Lastly, for any nonnegative random variable $Z$ and any $s>0$, we denote \[ \mathrm{E}^{s}[Z] := (\mathrm{E}[Z])^{s}. \]

A general theory for bias correction

To facilitate presentation, we begin by introducing general conditions for the bias-corrected estimator $\hat T_n^{\rm bc}$ under certain high-level assumptions on the estimator of $G_x(\cdot)$. These assumptions will be verified in Section (ref) for a specific class of ridge least squares estimators.

We first impose regularity conditions on $F_X$.

Assumption (ref)(ref) is standard in the NN literature. On the other hand, we refer the readers of interest to evans2002asymptotic for distributions that satisfy Assumption (ref)(ref).

Next, we assume sufficient smoothness of the function $G_{\cdot}(t)$, requiring a higher degree than Azadkia21simple, but similar to what was posed in lin2021estimation and cattaneo2024rosenbaum.

Finally, we assume regularity of the estimator $\hat G_{\cdot}(t)$, in line with lin2021estimation and cattaneo2024rosenbaum.

Ridge least squares

The remaining task is to verify Assumption (ref), which cannot be directly addressed using standard nonparametric regression results newey1997convergence, chen2018optimal, belloni2015some. In particular, unlike previous work, our setting requires control of estimation error in expectation rather than in probability, necessitating a more stable estimator.

To this end, we employ ridge regularization via ridge least squares tuo2024asymptotic, kurisu2024series. Following the setup of cattaneo2024rosenbaum, let \[ p_K(\cdot) = (p_{1K}(\cdot), \ldots, p_{KK}(\cdot))^\top \in \mathbb{R}^K \] be a $K$-dimensional vector of basis functions capable of approximating \[ \psi_t(\cdot) := G_{\cdot}(t), \quad \text{for any } t \in {\rm supp}(\mu). \] We estimate \[ G_x(t) = \mathrm{E}[\mathbbm{1}(Y \ge t) \mid X = x] \] by projecting $\mathbbm{1}(Y \ge t)$ onto the span of the basis functions. This corresponds to nonparametric linear probability models.

Define \[ \psi_{t,K}(x) := p_K(x)^\top \beta_{t,K}, \quad \hat{\psi}_{t,K}(x) := p_K(x)^\top \hat{\beta}_{t,K}, \] where $\beta_{t,K}$ is the population $L^2$ projection: \[ \beta_{t,K} := \argmin_{b \in \mathbb{R}^K} \mathrm{E}\left[(\psi_t(X_1) - p_K(X_1)^\top b)^2\right], \] and $\hat\beta_{t,K}$ is the ridge estimate: \[ \hat{\beta}_{t,K} := \argmin_{b \in \mathbb{R}^K} \left\{ \frac{1}{n} \sum_{i=1}^n (\mathbbm{1}(Y_i \ge t) - p_K(X_i)^\top b)^2 + \lambda_n \|b\|^2 \right\}, \] with regularization parameter $\lambda_n > 0$. Introduce further

The following quantities play a key role in characterizing the behavior of series estimators and their associated approximation errors: \[ \underline{\lambda}_K := \lambda_{\min}(Q), \quad \zeta_{r,K} := \max_{\alpha \in \Lambda_r} \sup_{x} \|D^\alpha p_K(x)\|, \quad \vartheta_{r,K}^t := \max_{\alpha \in \Lambda_r} \|D^\alpha \psi_t - D^\alpha \psi_{t,K}\|_{\infty}. \] Here, \(\lambda_{\min}(Q)\) denotes the smallest eigenvalue of \(Q\); \(\zeta_{r,K}\) measures the smoothness of the basis functions \(p_K\); and \(\vartheta_{r,K}^t\) captures the best possible order-\(r\) approximation error of the basis \(p_K\) in representing the target function \(\psi_t(\cdot)\).

The above assumption ensures that $Q$ is invertible and is standard in the literature newey1997convergence, cattaneo2024rosenbaum. The expressions for $\beta_{t,K}$ and $\hat{\beta}_{t,K}$ can then be simplified to \[ \beta_{t,K} = Q^{-1} \mathrm{E}[p_K(X) \psi_t(X)], \quad \hat{\beta}_{t,K} = (P^\top P + n \lambda_n I_K)^{-1} P^\top \mathbbm{1}(Y_{[n]} \ge t). \]

Our goal is to derive a bound for \[ \mathrm{E}\left[ \max_{\alpha \in \Lambda_r} \| D^\alpha \hat{\psi}_{t,K} - D^\alpha \psi_t \|_\infty^2 \right]. \] Such a bound can be directly used to verify Assumption (ref). To this end, we introduce the following assumption, which imposes regularity conditions on the smoothness of the basis functions $p_K$.

Assumption (ref) is mild. For instance, $\zeta_{r,K} = O(K^{1+r})$ for power series, and $\zeta_{r,K} = O(K^{1/2 + r})$ for Fourier series, splines, compactly supported wavelets, and piecewise polynomials. Moreover, $\underline{\lambda}_K$ is typically bounded away from zero uniformly in $K$; see newey1997convergence and cattaneo2024rosenbaum. Importantly, we do not require the “noise” term $\mathbbm{1}(Y \geq t) - \psi_t(X)$ to be independent of the covariates $X$.

Simulation studies

In this section, we present a series of illustrative simulations to demonstrate the performance of the bias-corrected estimator and compare it with the original NN graph-based statistic proposed by Azadkia21simple. To this end, we consider the following simulation setup:

In each simulation round, we generate independent copies \( \{(\widetilde{X}_i, \widetilde{Y}_i)\}_{i=1}^n \) from the above multivariate normal distribution. Then we define $X= \big(\Phi(\widetilde{X}^{(1)}), \ldots, \Phi(\widetilde{X}^{(d)})\big)^\top$ and $Y=\Phi(\widetilde{Y})$, where $\Phi$ is the standard normal CDF. In this case, by ansari2025directextensionazadkia, we have the closed-form formula: \[ T(Y,X)\;=\;\frac{3}{\pi}\,\arcsin\!\Big(\frac{1+\rho^{2}}{2}\Big)\;-\;\frac{1}{2}. \]

We consider \( \rho = 0, 0.3, 0.5, 0.7, 0.9 \), dimensions \( d = 2, 4, 6, 8, 10 \), and sample sizes \( n = 300, 600, 900 \). For the bias correction, we use power series basis functions of polynomial degree equal to $2$. We have set the ridge penalty parameter $\lambda_n = n^{-0.85}$. Variance estimation for both \( \hat{T}_n \) and the bias-corrected estimator \( \hat{T}_n^{\text{bc}} \) is performed via the \( m \)-out-of-\( n \) bootstrap with \( m = \lfloor n^{1/2} \rfloor \), as recommended in dette2025simple.

Tables (ref)--(ref) report the root mean square errors (RMSE) and empirical coverage probabilities (ECP) for both estimators over 1,000 independent repetitions. We consider the nominal significance level \( \alpha = 0.05 \) for the empirical coverage. The empirical results show that as the dimension or correlation increases, the performance of \( \hat{T}_n^{\text{bc}} \) clearly dominates that of \( \hat{T}_n \). For small dimensions (especially when $d=2$) or correlation (e.g., when $\rho=0$ when no bias is present), the two (original and bias-corrected) correlation estimators have more comparable performance. All these simulations give necessary supplement to the derived theory.

All reproducible code is available in the associated GitHub repository: \url{https://github.com/chenleihaomars/Simulation-for-ACbc}.

{

}

{

}

{

}

{

}

{

}

Proofs

Proof for (ref)

We introduce some notation. Let $\mathbb{S}^{d-1}$ denote the $(d-1)$-dimensional sphere in $\mathbb{R}^{d}$ and $B(x,r)\subseteq \mathbb{R}^d$ denote the $d$-dimensional open ball in $\mathbb{R}^d$ with center $x$ and radius $r$. Let $v_d$ denote the volume of unit Euclidean ball in $\mathbb{R}^d$. Denote by $\sigma$ the surface measure on $\mathbb{S}^{d-1}$. Let the real number $\mathfrak{D}$ denote the diameter of $\mathcal{X}$, i.e., $\mathfrak{D}=\sup_{x,y\in \mathcal{X}}\|x-y\|$.

Proofs for (ref)

We start by introducing the additional notation used in the following. In the following, let $\|\cdot\|_2$ be the spectral (operator) norm of a matrix, $\mathrm{Tr}(M)$ denote the trace of matrix $M$, and $\mathrm{Diag}(m_1,\ldots,m_n)$ represent an $n$-dimensional diagonal matrix with diagonal elements $m_1,\ldots,m_n$. For symmetric matrices $M_1$ and $M_2$, we write $M_1\preceq M_2$ if $M_2-M_1$ is a positive semi-definite matrix. Let $\|\cdot\|_{L^2}$ represent the $L^2(F_X)$-norm of functions.

Proofs of supporting lemmas