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.
62,242 characters · 19 sections · 65 citation commands
A Powerful Bootstrap Test of Independence in High Dimensions
\thispagestyle{empty}
{\bf Keywords:} block multiplier bootstrap, Chatterjee's rank correlation, family-wise error rate, high-dimensional data, independence test.
\newgeometry{tmargin=2cm,lmargin=2cm,rmargin=2cm}
This paper is concerned with nonparametric testing of pairwise independence between a random variable $X$ and many other random variables $Y_1,\ldots,Y_p$, $$H_0\colon Y_j\perp X\:\text{ for all }\:j\in\{1,\dots,p\}, $$ against the alternative $H_1$, which is the negation of $H_0$. The goal is to propose a powerful test of $H_0$ allowing for $p$ to be much larger than the sample size while at the same time not restricting the dependence among $Y_1,\ldots,Y_p$ in any way. In a second step, we want to combine the new test with a stepwise procedure for screening out variables from $Y_1,\ldots,Y_p$ that violate independence so as to control the family-wise error rate.
There are many applied examples in which testing $H_0$ and, in particular, screening out variables that violate independence is of interest. For instance, in causal inference, one might want to test whether a treatment indicator has an effect on various outcomes and then select those outcomes on which there is an effect. Such a test could also be applied to “placebo” outcomes, i.e. pre-treatment outcomes that the researcher knows cannot have been affected by the treatment, to validate unconfoundedness assumptions. Another example concerns testing fairness of machine learning predictions, where one might want to test independence of a prediction from a set of protected characteristics. As a final example, one might want to test independence of a measure of environmental exposure (e.g., whether or not a person smokes) from a vector of genetic markers. In our empirical application, we use the proposed test for identifying genes whose transcript levels oscillate during the cell cycle.
As test statistic we consider the maximum of $p$ rank correlation coefficients by chatterjee2021new for testing independence between $Y_j$ and $X$, $j=1,\ldots,p$. Critical values for the test are computed via a block multiplier bootstrap that perturbs an asymptotically linear representation of the rank correlations. We then show that, as the sample size grows, the proposed test controls size uniformly over a large class of data-generating processes. This result allows for the dimension $p$ to grow exponentially in sample size and the dependence among $Y_1,\ldots,Y_p$ is left completely unrestricted.
While it has been shown that the nonparametric bootstrap does not consistently estimate the distribution of a single Chatterjee's rank correlation Lin:2024ui, the proposed block multiplier bootstrap achieves this by accounting for the dependence of the summands in the asymptotically linear representation. The existence of an asymptotic linear representation for Chatterjee's rank correlation together with suitably strong control of the remainder terms allows us to employ recent results for the approximation of maxima of sums of high-dimensional vectors chernozhukov2013gaussian,chernozhukov2015comparison,chernozhukov2019inference to establish the validity of the block multiplier bootstrap also for the maximum of many Chatterjee's rank correlation coefficients.
The block multiplier bootstrap requires a tuning parameter choice, namely the block size that should tend to infinity with sample size. We are able to provide a simple choice that enjoys a certain optimality property: it minimises the distance of the bootstrap variance and the asymptotic variance for each individual test statistic. In simulations, however, we find that the size and power of our test are almost insensitive to the particular value of the block size.
The test is consistent against all fixed alternatives under which $Y_1,\ldots,Y_p,X$ are continuously distributed. Finally, our test is computationally attractive even for very large sample sizes and high dimensions $p$.
The test can be combined with a stepwise multiple testing procedure romano2005exact for screening out variables from $Y_1,\ldots,Y_p$ that violate independence so as to asymptotically control the family-wise error rate uniformly over a large set of data-generating processes. This result allows for the dimension $p$ to grow exponentially in sample size and the dependence among $Y_1,\ldots,Y_p$ is left completely unrestricted. The companion R package \href{https://github.com/mauolivares/hdIndep}{hdIndep} facilitates the implementation of the methodology developed in this paper.
\paragraph{Related Literature}
Our paper is most closely related to the literature on testing independence of two random vectors, $H_0^{joint}\colon Y \perp X$, where $Y\in\mathbb{R}^p$ and $X\in\mathbb{R}^q$. This is a different hypothesis from the one we consider; for $q=1$, it implies, but is not implied by, $H_0$; when $H_0$ holds, but the copula of $Y_1,\ldots,Y_p$ given $X$ is not independent of $X$, then $H_0^{joint}$ is violated. Therefore, tests that control the rejection probability under $H_0^{joint}$ are not guaranteed to control it under our hypothesis $H_0$. Sinha:1977re, Taskinen:2005ty, Bakirov:2006ui, Szekely:2007io, Heller:2012oi, Heller:2012op, Shi:2022oi propose nonparametric tests of $H_0^{joint}$, where $p$ and $q$ are of arbitrary, but fixed (with sample size) dimensions. Szekely:2013re show that the test statistic of Szekely:2007io is biased in high dimensions and an independence test based on it therefore does not control size in high dimensions. They also propose a bias-corrected test statistic and derive its asymptotic distribution under the null of independence when the dimensions of both vectors grow with the sample size. The asymptotic regime under which their test is valid requires the dimensions of both vectors to grow, so it is not clear (at least to us) whether it is also valid when one of the two dimensions remains constant as the sample size grows. In addition, their derivation of the test statistic's limiting distribution requires the elements of the two vectors to be exchangeable, a condition we do not require for $Y_1,\ldots,Y_p$. Ramdas:2015aa show that both independence tests, Szekely:2007io and Szekely:2013re have low power against “fair alternatives” in high dimensions. Our simulations in Section (ref) show that these tests do not necessarily control size under $H_0$.
More recently, zhou2024rank and Wang:2024ui propose other rank-based tests, e.g. based on Hoeffding's D, Blum-Kiefer-Rosenblatt's R and Bergsma-Dassios-Yanagimoto's $\tau$ among others, of $H_0^{joint}$ in high dimensions. The validity of these tests relies on at least one of the dimensions $p$ and/or $q$ diverging so that a central limit theorem across the elements of, say, $Y$ can be invoked. This approach necessarily restricts the dependence of $Y_1,\ldots,Y_p$, while our validity results leave the dependence completely unrestricted.
In simulations, we find that our test is typically more powerful than the alternatives by Szekely:2007io, Szekely:2013re, zhu2020distance, and zhou2024rank (in scenarios in which they are valid) in high dimensions (large $p$) or when there is dependence among the variables $Y_1,\ldots,Y_p$.
Our proposed test has two additional advantages over competitors (in scenarios in which they are valid): (i) unlike tests by zhu2020distance and zhou2024rank ours is computationally inexpensive and (ii) we develop a stepwise procedure for selecting hypotheses that are not rejected so as to control the family-wise error rate.
There is a large literature on nonparametric tests of mutual independence among the elements of a random vector. Some examples are Leung:2018iu, Shun:2018ui, Drton:2020oi, Wang:2024op, Bastian:2024io; see also references therein. Xia:2024aa propose such a test based on Chatterjee's rank correlation. While the hypothesis considered in our paper also involves many nonparametric independence tests, it substantially differs from the null of mutual independence. This is because, in our testing problem, $X$ occurs in every independence hypothesis and the dependence of $Y_1,\ldots,Y_p$ is left unrestricted.
Finally, since our proposed test is based on the rank correlation for two random variables proposed by chatterjee2021new, our paper is also related to a recent and fast-growing literature that examines the rank correlation coefficient's properties. chatterjee2021new shows asymptotic normality of the correlation coefficient under independence of the two random variables. Lin:2022aa and kroll2024asymptotic show that it is also asymptotically normal under dependence. Shi:2021oi and Lin:2022io examine and propose improvements of the power of tests of independence based on the rank correlation coefficient. Based on the earlier work by chatterjee2021new, Azadkia:2021oi introduce a graph-based correlation coefficient that can be seen as a multivariate extension of Chatterjee's correlation coefficient. Han:2024io derive its asymptotic distribution under the null that a vector $Y$ (with fixed dimension) is independent of a random variable $X$, and they find, in simulations, that the test may be more powerful than competitors in higher dimensions. For a recent review of this literature, see Chatterjee:2024op.
This section first introduces the new test, establishes asymptotic size control uniformly over a large class of data-generating processes, and then consistency against all fixed alternatives under which $Y_1,\ldots,Y_p$ and $X$ have continuous distributions. The section concludes with the development of an optimal tuning parameter choice.
Let $\mathbb{D}\coloneqq \{(X_i,Y_{1,i},\ldots,Y_{p,i})\}_{i=1}^n$ be an i.i.d. sample drawn from the distribution of $(X,Y_1,\ldots,Y_p)$. For each individual hypothesis $H_{0,j}\colon Y_j \perp X$ there are many available tests in the literature. In this paper, we focus on the test statistic by chatterjee2021new. The motivation for this choice will become clear later in this section. To define the test statistic let $X_{(k)}$ be the $k$-th order statistic of $X_1,\ldots,X_n$, i.e. $X_{(1)} \leq \dots \leq X_{(n)}$, and $Y_{j,(k)}$ be the concomitant of $X_{(k)}$, i.e. if $X_{(k)}=X_l$, then $Y_{j,(k)}=Y_{j,l}$. Denote by $F_{Y_j}$ the cumulative distribution function (cdf) of $Y_j$ and by $\hat{F}_{Y_j}$ the empirical cdf. Then, Chatterjee's rank correlation for testing an individual hypothesis $H_{0,j}$ is
chatterjee2021new shows that $\hat\xi_j$ is a consistent estimator of $$\xi_j \coloneqq \frac{\int Var(\operatorname*{\mathbb{E}}[\mathds{1}\{Y_j\geq t\}|X])f_{Y_j}(t)dt}{\int Var(\mathds{1}\{Y_j\geq t\})f_{Y_j}(t)dt}, $$ a measure of dependence introduced by Dette:2013yy in the case in which $Y_j$ has a continuous distribution with density $f_{Y_j}$. This measure has several attractive features chatterjee2021new. Two features that are particularly important for the test proposed in this paper are that (i) $\xi_j$ is equal to zero if, and only if, $X$ and $Y_j$ are independent, and (ii) the estimator $\hat\xi_j$ admits an asymptotic linear representation (shown in (ref) below) with a remainder that we can show to be sufficiently small. Our arguments for validity of the proposed test in high dimensions crucially depend on property (ii).
The proposed test statistic for $H_0$ is the maximum of the individual Chatterjee's rank correlations:
We propose to compute critical values for the test statistic via a block multiplier bootstrap. To describe the procedure consider first an individual hypothesis $H_{0,j}$. If the null $H_{0,j}$ holds and both random variables are continuously distributed, then the arguments in angus1995coupling imply that Chatterjee's rank correlation has an asymptotically linear representation of the form
where $r_{j,n}$ is a negligible remainder term and $W_{j,i}$ is defined as $$W_{j,i}\coloneqq 2 - 3\@ifstar{\oldabs}{\oldabs*}{U_{j,i+1}-U_{j,i}} - 6 U_{j,i}(1-U_{j,i}) $$ with $U_{j,i} \coloneqq F_{Y_j}(Y_{j,(i)})$. A naive application of the multiplier bootstrap idea would be to repeatedly draw bootstrap multipliers $\varepsilon_1,\ldots,\varepsilon_n$ as independent standard normal random variables that are independent of the data $\mathbb{D}$ and then compute a critical value from the distribution of $\frac{1}{\sqrt{n}}\sum_{i=1}^{n-1}\varepsilon_i W_{j,i}$ given the data. However, there are two problems with this approach. First, $\{W_{j,i}\}_{i=1}^n$ is not an i.i.d. sequence but a 1-dependent process. Second, $\{W_{j,i}\}_{i=1}^n$ is not directly observed because $F_{Y_j}$ is unknown.
To address the first challenge, we decompose the sum $\sum_{i=1}^{n-1} W_{j,i}$ into the sums over “big” and “small” blocks formed of $\{W_{j,i}\}_{i=1}^n$ with the property that the big blocks are independent of each other. Formally, let $q\ge 1$ denote the size of big blocks. It is a tuning parameter to be chosen by the researcher; in Section (ref) we develop an optimal choice of $q$. Since $\{W_{j,i}\}_{i=1}^n$ is a 1-dependent sequence, we let the small blocks be of length one. Further, let $m\coloneqq \lfloor(n-1)/(q+1)\rfloor$, where $\lfloor\nu\rfloor$ is the integer part of $\nu$, denote the number of big blocks. Lastly, $r \coloneqq n-1 - m(q+1)$, $0\le r<q+1$, is the number of remaining summands that are not allocated to any of the small or big blocks. With this notation, we can thus write $$\sum_{i=1}^{n-1}W_{j,i} = \sum_{k=1}^m A_{j,k} + \sum_{k=1}^m B_{j,k} + R_j,$$ where \[ A_{j,k} \coloneqq \sum^{kq+(k-1)}_{l=(k-1)(q+1)+1} W_{j,l}, \quad B_{j,k} \coloneqq W_{j,k(q+1)}, \quad \text{ and } \quad R_j \coloneqq \sum_{k=1}^r W_{j,m(q+1)+k}. \] The big block sums $A_{j,1},\ldots,A_{j,m}$ are independent of each other, and we show that the terms $B_{j,1},\ldots,B_{j,m}$ and $R_j$ are asymptotically negligible. It then follows that $\sqrt{n}\,\hat{\xi}_j $ can be approximated by $\frac{1}{\sqrt{mq}}\sum_{k=1}^m A_{j,k}$, the sum of independent components.
We now need to address the second challenge, which is that $W_{j,i}$, and thus also $A_{j,k}$, are not observed. $W_{j,i}$ can be estimated by \[ \hat W_{j,i}\coloneqq 2 - 3\@ifstar{\oldabs}{\oldabs*}{ \hat U_{j,i+1}- \hat U_{j,i}} - 6 \hat U_{j,i}(1-\hat U_{j,i}), \] where $\hat U_{j,i}\coloneqq \hat{F}_{Y_j}(Y_{j,(i)})$, and $A_{j,k}$ by \[ \hat{A}_{j,k} \coloneqq \sum^{kq+(k-1)}_{l=(k-1)(q+1)+1} \hat{W}_{j,l}. \] While $A_{j,1},\ldots,A_{j,m}$ are independent, $\hat{A}_{j,1},\ldots,\hat{A}_{j,m}$ are only asymptotically independent, i.e. in the limit as $n,m\to\infty$.
Finally, bootstrap multipliers $\varepsilon_1,\ldots,\varepsilon_m$ are drawn as independent standard normal random variables that are independent of the data $\mathbb{D}$. The bootstrap statistic is then defined as
For a nominal level $\alpha\in(0,1)$, the proposed critical value $\hat c(\alpha)$ for our test is the conditional $(1-\alpha)$-quantile of $\hat{T}^B$ given the data $\mathbb{D}$. The test rejects $H_0$ if, and only if, the test statistic $\hat T$ exceeds the critical value $\hat c(\alpha)$.
In this subsection, we show that our proposed test asymptotically controls size uniformly over a large class of data-generating processes.
Assuming all random variables have a continuous distribution simplifies the presentation, but is not essential. Remark (ref) below discusses extensions to the case with discrete distributions.
This assumption requires $q$ to diverge as the sample size grows, but restricts its rate to be neither too slow nor too fast. The assumption also restricts the rate at which the dimension $p$ is allowed to grow with sample size. However, the upper bound on the rate is very large: $p$ may be an exponential function of sample size and thus is allowed to be much larger than sample size. For instance, there are positive constants $\delta_1,\delta_2$ so that $p = e^{n^{\delta_1}}$ and $q = n^{\delta_2}$ satisfy Assumption (ref).
The result in (ref) implies that the proposed test asymptotically controls size. In fact, the asymptotic size is equal to the nominal level $\alpha$ and, in this sense, the test is not conservative. Furthermore, the probability of rejecting $H_0$ when $H_0$ is satisfied can deviate from the nominal level only by a term that is polynomially small in $n$. Importantly, the constants $c$ and $C$ depend on the data-generating process only through the constants $\gamma$ and $C_1$ from Assumption (ref). Therefore, under $H_0$, the inequality in (ref) holds uniformly over all data-generating processes that satisfy the assumption with the same constants, denoted by $\mathbf{P}_{\gamma,C_1}$: $$\limsup_{n\to\infty}\sup_{P\in\mathbf{P}_{\gamma,C_1}} \left|{\mathrm{P}}\left( \hat{T} > \hat{c}(\alpha)\right) - \alpha\right| =0, $$ i.e. our test has asymptotic size equal to $\alpha$ uniformly over those data-generating processes.
It is worth noting that the validity of our test is guaranteed in high dimensions without restricting the dependence of $Y_1,\ldots,Y_p$ in any way.
To establish this result we need to show that the distribution of the bootstrap statistic $\hat{T}^B$ given the data is close to the distribution of the test statistic $\hat{T}$. This is achieved in several steps: (i) show that $\hat{T}$ is close to $$T_0\coloneqq \max_{1\le j\le p}\,\frac{1}{\sqrt{n}}\sum_{i=1}^{n-1}W_{j,i}~,$$ (ii) show that $\hat{T}^B$ is close to $$T_0^{B}\coloneqq\max_{1\leq j \leq p} \frac{1}{\sqrt{mq}}\sum_{k=1}^m \varepsilon_k A_{j,k},$$ and then (iii) show that both $T_0$ and $T_0^B$ are close to the maximum of Gaussian random variables, $$Z_0\coloneqq \max_{1\le j\le p} V_j,$$ where $V\coloneqq (V_1,\dots,V_p)'$ is a Gaussian vector such that $\operatorname*{\mathbb{E}}(V)=0$ and $\operatorname*{\mathbb{E}}(V\,V')=\frac{1}{mq}\sum_{k=1}^m\operatorname*{\mathbb{E}}(A_k A_k')$, $A_k\coloneqq (A_{1,k},\dots,A_{\,p,k})'$. Steps (i) and (ii) are developed in the proof of the theorem and step (iii) is delegated to Lemmas (ref) and (ref). These three steps establish that the distribution of the test statistic is close to that of the bootstrap statistic. However, this result does not yet imply the statement in (ref) because $\hat{c}(\alpha)$ is random and may be correlated with $\hat{T}$. The final step of the proof therefore consists of passing from a deterministic to a random critical value.
Step (iii) of the derivation combines several results on high-dimensional Gaussian approximations from chernozhukov2013gaussian,chernozhukov2015comparison,chernozhukov2019inference. The main challenge of the proof is contained in steps (i) and (ii), where we need to establish that all remainder terms $r_{1,n},\ldots,r_{p,n}$ from the representation (ref) vanish at a suitably fast rate. angus1995coupling shows that, for each $j$, $r_{j,n}=o_P(n^{-1/2})$, but in our high-dimensional setting, we need stronger control of these remainder terms, and these results are developed in Lemma (ref).
The following theorem shows that the proposed test is consistent against any (fixed) violation of the null:
In light of recent discussions related to the local asymptotic power of Chatterjee's test of independence (e.g., Shi:2021oi and Lin:2022io), it would be interesting to study the power of our proposed test to local alternatives. We view this analysis as beyond the scope of the paper and relegate it to future research.
Theorem (ref) shows that our proposed block multiplier bootstrap test asymptotically controls size under rate conditions on $q$, the size of blocks used to construct the critical value $\hat c(\alpha)$. These rate conditions specify that $q\to \infty$ as $n \to \infty$ at a rate that is neither too fast nor too slow, but they do not provide any guidance on how to choose $q$ in finite samples. In this section, we develop a choice of $q$ that is optimal in a certain finite-sample sense.
To describe the optimality criterion, recall from the discussion following Theorem (ref) that the distribution of the bootstrap test statistic $\hat T^B$ is, to first order, determined by its infeasible counterpart $T_0^B=\max_{1\leq j \leq p}T^B_{0,j}$, where $$T^B_{0,j} \coloneqq \frac{1}{\sqrt{mq}} \sum_{k=1}^m \varepsilon_k A_{j,k} $$ with bootstrap multipliers $\varepsilon_1,\ldots,\varepsilon_m$ that are i.i.d. standard normal and independent of the data. By construction, $T^B_{0,j}$ follows a normal distribution with mean zero and variance that depends on the data and, implicitly, on $q$, \[ T^B_{0,j} | \mathbb{D} \sim \mathcal{N}\left(0,V^B_j\right), \quad V^B_j \coloneqq {\rm Var}\left( T^B_{0,j}| \mathbb{D} \right). \] In our bootstrap procedure, $T^B_{0,j}$ mimics the behaviour of the individual test statistic $\sqrt{n} \hat \xi_j$, which is asymptotically normally distributed and, under the null, has variance $v_n$ given in Remark (ref). Since $q$ affects the conditional distribution of $T^B_{0,j}$ only through the bootstrap variance $V^B_j$ and it does not affect the test statistic itself, we aim to choose $q$ such that $V^B_j$ is close to $v_n$.
The following lemma characterizes the expectation and the variance of the bootstrap variance $V^B_j$ under the null $H_{0,j}$.
Lemma (ref) shows that the expectation of $V_j^B$ is above 2/5 for any fixed $q$ and it monotonically converges to 2/5 as $q$ grows large. Since the target variance $v_n$ approaches 2/5 from below (see Remark (ref)), this means that $V_j^B$ is biased upwards with a bias that vanishes asymptotically when $q,n \to \infty$. The variance of $V_j^B$, in turn, is generally increasing in $q$. Figure (ref) illustrates these relationships for samples of size $n=500$ and $n=1000$. We resolve this bias-variance trade-off by considering the mean squared error of the bootstrap variance $V_j^B$: $$MSE_{j,n}(q)\coloneqq \operatorname*{\mathbb{E}}\left[\left(V_j^B - v_n\right)^2\right] = \left\lfloor\frac{n-1}{q+1}\right\rfloor^{-1}\left.
\right\} + \left(\frac{2}{5} + \frac{1}{10q} - v_n\right)^2. $$ Minimising $MSE_{j,n}(q)$ over $q\in\mathbb{N}^+$ yields our proposed optimal choice of $q$: \[ q^*(n) \coloneqq \operatorname*{arg\,min}_{q\in\mathbb{N}^+} MSE_{j,n}(q). \] Since $MSE_{j,n}(q)$ is a known function of the sample size, the optimal choice $q^*(n)$ does not depend on the data (beyond $n$). The reason for that is that the individual bootstrap statistic $T^B_{0,j}$ depends on the data only through the ranks of the concomitants, $U_{j,i}\coloneqq F_{Y_j}(Y_{j,(i)})$, and these are independent and uniformly distributed under the null. In consequence, the optimal choice $q^*(n)$ is independent of $j$ and thus the same for each independence hypothesis $H_{0,j}$.
The minimiser $q^*(n)$ can be computed for each $n$ by evaluating the function $MSE_{j,n}(q)$ over a grid of values $q \in \{1,2,\ldots, \lfloor (n-1)/2 \rfloor\}$. On most of its domain, $q^*(n)$ is a step function with large flat regions, but it can oscillate in small transition regions. For example,
$$q^*(n) =
$$ where $q^*(225)=3$, $q^*(n)=2$ for $n \in \{226,\ldots,232\}$, $q^*(233)=3$, etc. The non-monotone ranges make it impractical to tabularise $q^*(n)$. Instead, we approximate $q^*(n)$ using a smooth function. To motivate this approximation, note that in a setting where $n$ and $q$ are large but $q$ is much smaller than $n$, which is consistent with the asymptotic regime of Assumption~\ref{ass:rates}, $MSE_{j,n}(q)$ is close to $ 0.32 q/n + 0.01/q^2$, which is minimised at $$ \tilde q(n) \coloneqq (n/16)^{1/3}. $$ This expression turns out to provide a good approximation to $q^*(n)$ even for small $n$ in the sense that $|\tilde{q}(n) - q^*(n)| < 1$ for any $n \in \mathbb{N}^+$. This property implies that \[ q^*(n) =
\] The above approximating property of $\tilde q(n)$ is illustrated in Figure (ref).
We note that the goal of this paper is to test the hypothesis $H_0$ that all individual hypotheses $H_{0,j}$, $j=1,\ldots,p$, hold simultaneously, while we derived the optimal $q^*(n)$ for an individual test statistic. Therefore, this choice does not necessarily minimise the distance between the distributions of the max-test statistic $\hat{T}\coloneqq \max_{1\leq j\leq p} \sqrt{n}\hat\xi_j$ and the bootstrap statistic $\hat{T}^B \coloneqq \max_{1\leq j \leq p} \hat{T}^B_j$ in any sense. However, minimising the distance between these two distributions is considerably more difficult because the individual statistics may be arbitrarily dependent and the optimal $q$ would then depend on their (unknown) dependence structure. Developing a feasible version of this approach would require estimation of the copula of a high-dimensional random vector, which our proposal above avoids.
In simulations in Section (ref), we find that our proposed choice of $q$ not only optimises the bootstrap approximation of the individual statistics, but also yields a good bootstrap approximation for the max-statistic.
In the previous section, we considered testing the null $H_0$ that all variables $Y_1,\ldots,Y_p$ are independent of $X$. Now, consider the problem of selecting individual hypotheses $H_{0,j}\colon Y_j\perp X$ that are violated. The previous section already yields a (“single-step”) method of selection: simply select all hypotheses $H_{0,j}$ for which $\sqrt{n}\hat \xi_j > \hat c(\alpha)$. This section introduces a stepdown procedure that improves upon the single-step procedure by possibly rejecting more hypotheses in finite samples. In addition, we show that the stepdown procedure (and, thus, by extension also the single-step procedure) guarantees asymptotic control of the family-wise error rate.
Let $J(P)\subseteq \{1,\ldots,p\}$ denote the set of hypotheses $H_{0,j}$ that are true under $P$. The family-wise error rate is defined as the probability of rejecting at least one true hypothesis, $$FWER_P \coloneqq {\mathrm{P}}(\text{reject at least one } H_{0,j}\colon j\in J(P)). $$ For any $I\subseteq \{1,\ldots,p\}$, let $$\hat T(I) \coloneqq \sqrt{n} \max_{j\in I} \hat \xi_j\qquad\text{and}\qquad \hat T^B(I) \coloneqq \max_{j\in I} \frac{1}{\sqrt{mq}} \sum_{k=1}^m \varepsilon_k \hat A_{j,k}, $$ where $\varepsilon_k$ and $\hat A_{j,k}$ are the multipliers and estimated blocks as introduced in the previous section. Finally, define $\hat c(\alpha;I)$ as the ($1-\alpha$)-quantile of $\hat T(I)$ given the data $\mathbb{D}$.
The following algorithm introduces the stepdown procedure.
The theorem shows that the stepdown procedure in Algorithm (ref) asymptotically controls the family-wise error rate uniformly over data-generating processes in $\mathbf{P}_{\gamma,C_1}$.
In this section, we report results from a series of extensive simulation experiments in which we studied the finite-sample performance of our new test and compared it to existing tests.
First, we examine the influence of the choice of block size $q$ on the finite-sample performance of our test. Having established that the test's rejection frequency is fairly insensitive to $q$, in subsequent experiments, we only consider our test with the optimal choice derived in Section (ref). In the second set of experiments, we extensively study size control and power of our test. We consider a variety of data-generating processes. Key parameters that we vary in the simulations are the degree of dependence among $Y_1,\ldots,Y_p$, the dimension $p$, and whether alternatives are sparse or dense. Having confirmed the theoretical findings that our test controls size and has power against all alternatives, we then show that existing tests based on distance covariance do not control size under our null hypothesis. They are valid only under the stronger null $H_0^{joint}$ mentioned in the introduction, i.e. when the copula of $Y_1,\ldots,Y_p$ given $X$ does not depend on $X$. Finally, we compare our test to existing tests in the special case in which the copula of $Y_1,\ldots,Y_p$ given $X$ does not depend on $X$, because these other tests are valid in that case.
Let $X$ be distributed uniformly on $[-1,1]$, and let $(\epsilon_1,\ldots,\epsilon_p) \sim \mathcal N(0, \Sigma_\tau)$ be a random vector independent of $X$, where $\Sigma_\tau \in \mathbb{R}^{p \times p} $ has diagonal elements equal to one and off-diagonal elements equal to $\tau$. We consider a range of deviations from the null of pairwise independence that differ in the functional form of the association between $X$ and elements of $Y_1, \ldots, Y_p$, as well as in the number and strength of violations of individual hypotheses. In all considered models, the magnitude of the violations is parameterised by $\rho \in \mathbb{R}$, with the null hypothesis corresponding to $\rho=0$, while $\tau$ is the correlation between $Y_1, \ldots, Y_p$ under the null.
In Models 1 and 2, $\rho$ determines the magnitude of the violation of the first hypothesis $H_{0,1}\colon Y_1\perp X$, while all other hypotheses, $H_{0,j}\colon Y_j\perp X$ for $j \in \{2,\ldots, p\}$, hold. In Models 5 and 6, all individual hypotheses are violated to the same degree, while Models 3 and 4 offer an intermediate scenario where the magnitude of the violations is strongest for the first hypothesis and decays exponentially with the index of the hypothesis. The individual violations are reparameterised versions of data-generating processes considered in the simulation study by chatterjee2021new. Note that under the null, when $\rho = 0$, all six models are identical.
We include cosine alternatives in these experiments partly because, in our empirical application, we conjecture such alternatives to be reasonable departures from independence.
We consider three variants of our proposed test that employ different types of studentisation. Implementations can be found in our accompanying R package \href{https://github.com/mauolivares/hdIndep}{hdIndep}.
For the special case in which the copula of $Y_1,\ldots,Y_p$ given $X$ does not depend on $X$, we also compare our tests to the following alternatives:
All simulations are based on $B=499$ bootstrap replications, and $S=5000$ Monte Carlo draws, except in Experiment 4.2, where $S$ is reduced to 1000 due to mitigate computational burden. The significance level is set to $\alpha=0.05$ throughout. Additional results are reported in Appendix (ref) and are omitted from the main text for brevity.
The first set of simulations concerns the rejection rates of our proposed test for different choices of the block size $q$. Figure (ref) presents the results under the null in Panel A, and under the alternatives specified by Models 1--6 in Panels B--G. We focus on the case in which there is no correlation between $Y_1,\ldots,Y_p$ under the null ($\tau=0$). The results are very similar when $\tau=0.5$; see Figure (ref) in Appendix (ref).
First, we note that our baseline procedure BMB0 controls the size across different scenarios for all considered values of $q$. It is, however, conservative in higher dimensions. The fact that the rejection rate is particularly low for $q=1$ is consistent with the upward bias in the bootstrap variance characterised in Lemma (ref). The optimal rule from Section (ref) yields $q^*(n)=3$ for $n=500$, marked by the vertical dashed lines in the graphs. This choice proves reasonable in all the considered scenarios. The issue of conservativeness is alleviated by the studentisation. Both BMB1 and BMB2 effectively shrink the bootstrap distribution and yield rejection probabilities very close to the nominal level of 5% for small values of $q$. Since the correction in BMB1 is negligible for large $q$, the blue and black lines approach each other as $q$ increases. BMB2 maintains rejection rates closer to 5% as $q$ grows, but it generally slightly overrejects.
Panels B--G of Figure (ref) indicate that all tests have high power in low dimensions. Under sparse and decaying alternatives, the power decreases as $p$ grows, while the power increases for dense alternatives.
Since, in Experiment 1, our test was seen to be fairly insensitive to the particular choice of block size $q$, we now analyse size control and power only for the optimal choice, which for $n=500$ is $q^*(n)=3$.
Figures (ref) and (ref) present the power curves of our test under linear and cosine alternatives, respectively. For each model, we consider the cases in which $Y_1,\ldots,Y_p$ are mutually independent ($\tau=0$) or dependent ($\tau=0.5$). In each graph, we show rejection frequencies of our test as we vary $\rho$. The null hypothesis corresponds to $\rho=0$, and as $\rho$ increases the violation of the null becomes larger.
First, the three variants of our test, BMB0, BMB1, and BMB2, control size across scenarios, including high-dimensional settings (large $p$) and those with dependence among $Y_1,\ldots,Y_p$ ($\tau=0.5$).
Second, our test has power against all considered alternatives. Interestingly, comparing the power curves for $\tau=0$ and $\tau=0.5$, we see that the power of our test is not or only slightly affected by dependence among $Y_1,\ldots,Y_p$.
In Appendix (ref), we present analogous results for sample sizes $n=200$ and $n=1000$. While the power of the test increases with sample sizes, the qualitative conclusions remain the same.
As indicated above, the validity of existing tests dcorT.test, dcov.test, ZXZL, and ZZYS-agg.dcov have been established under the stronger null $H_0^{joint}$ that $(Y_1,\ldots,Y_p)\perp X$. We now investigate the size control properties of these tests when our null hypothesis $H_0\colon Y_j\perp X$ for $j=1,\ldots,p$ holds, but $H_0^{joint}$ is violated.
To this end we consider a modified version of the data-generating processes Models 1--6, which are identical under the null, in which the dependence parameter $\tau$ depends on $X$. Specifically, we consider Model 1 with $\rho=0$ and instead of $(\epsilon_1,\ldots,\epsilon_p) \sim \mathcal N(0, \Sigma_\tau)$, we now consider $(\epsilon_1,\ldots,\epsilon_p) | X \sim \mathcal N(0, \Sigma_{\tau}(X))$, where $\Sigma_{\tau}(X)$ has diagonal elements equal to one and off-diagonal elements equal to $\tau(X) := \gamma (1+X)/2$. Since $\rho=0$, the marginal distributions of $Y_j$ are all independent of $X$, but the correlation parameter among $Y_1,\ldots,Y_p$ is a function of $X$. The parameter $\gamma$ governs the strength of that dependence.
Figure (ref) presents the rejection rates of all tests under the null hypothesis $H_0$ when $p=50$. Our test controls size across all values of $\gamma$, while the distance covariance tests dcorT.test and dcov.test do not and their rejection rates increase with $\gamma$. Interestingly, the ZXZL and ZZYS-agg.dcov tests also control size. Since they are based on marginal distance covariance and rank-based statistics, we conjecture that these tests are not only valid under $H_0^{joint}$, but also under $H_0$.
In this section, we compare our test to all other tests in the special case in which the copula of $Y_1,\ldots,Y_p$ given $X$ does not depend on $X$, both under the null and the alternative hypotheses. We divide this experiment into two sub-experiments because the ZXZL and ZZYS-agg.dcov tests are computationally so demanding that we were not able to run experiments in high dimensional settings ($p$ large). Therefore, we first consider an experiment in which we compare our test to the distance covariance tests dcorT.test and dcov.test for a range of values of $p$ up to very high dimensions. In a second experiment, we then compare our test to the ZXZL and ZZYS-agg.dcov tests in more moderate dimensions up to $p=500$.
Figures (ref) and (ref) present the power curves of our test under linear and cosine alternatives, respectively. For each model, we consider the cases in which $Y_1,\ldots,Y_p$ are mutually independent ($\tau=0$) or dependent ($\tau=0.5$). In each graph, we show rejection frequencies of our test as we vary $\rho$. The null hypothesis corresponds to $\rho=0$, and as $\rho$ increases the violation of the null becomes larger.
Several interesting findings emerge from these figures. First, while our test is almost insensitive to dependence among $Y_1,\ldots,Y_p$, the distance covariance tests dcorT.test and dcov.test are not. For instance, comparing Panel A and B in Figure (ref) for sparse alternatives, we see that their power may suffer substantially when $Y_1,\ldots,Y_p$ are dependent. A similar but less extreme pattern is observed under decaying linear alternatives (Panels C and D in Figure (ref)). Under dense alternatives, under which these tests are known to perform particularly well, their power is less affected by dependence among $Y_1,\ldots,Y_p$.
Second, as expected, the distance covariance tests perform particularly well under dense alternatives, dominating the power curves of our tests. Our test performs particularly well under sparse and decaying alternatives, where it tends to be more powerful than the distance covariance tests, particularly so in high dimensions.
Third, considering the cosine alternatives (Figure (ref)), we see that our test is more powerful than the distance covariance tests in all scenarios, including even the dense alternatives in which one would expect the distance covariance tests to perform well.
In Appendix (ref), we present analogous results for sample sizes $n=200$ and $n=1000$. While the power of all tests increases with sample sizes, the qualitative conclusions remain the same.
In this section, we compare our test to the ZXZL and ZZYS-agg.dcov tests.\footnote{The tests ZXZL-D, ZXZL-R, ZXZL-$\tau^*$, and ZZYS-agg.dcov were run using the implementation of zhou2024rank provided at \url{https://github.com/Yeqing-TJ/Rank-based-test-in-high-dimension} [Accessed on March 10, 2025].} Due to the high computational complexity of these additional tests, we limit the simulation setup to samples of size $n=200$ with the maximum of $p=200$ individual hypotheses.\footnote{Estimation of the variance of the rank-based indices of zhou2024rank in one sample with $p=500$ and $n=500$ takes over 40 minutes using an Intel i7-1185G7 @ 3.00 GHz processor, which renders simulations for such settings infeasible. For comparison, our bootstrap test with $B=499$ takes less than a second to compute in this setting.}
Figures (ref) and (ref) show power comparisons analogous to Figures (ref) and (ref), respectively. The broad patterns are similar to those in Figures (ref) and (ref) in that the ZXZL and ZZYS-agg.dcov tests exhibit similar performance to the distance covariance tests studied in Experiment 4.1.
In the low-dimensional settings with linear alternatives (Figure (ref)), the ZXZL and ZZYS-agg.dcov tests dominate our test when $p$ is small, there is no dependence among $Y_1,\ldots,Y_p$, or the alternatives are dense. Our test performs relatively better the larger $p$ and dominates the other tests when $p$ is sufficiently large (except for the dense alternatives). As in Experiment 4.1, our test's power curves are barely affected by dependence among $Y_1,\ldots,Y_p$, but the alternative tests' power curves decrease substantially when $Y_1,\ldots,Y_p$ are dependent; compare for instance Panels A and B in Figure (ref).
Under cosine alternatives (Figure (ref)), the ZXZL and ZZYS-agg.dcov tests have virtually no power, and our test is more powerful than the competitors in all scenarios.
The simulation experiments reveal several interesting findings:
We illustrate the practical usefulness of our proposed test by revisiting the study conducted by hughes2009harmonics, which investigates transcriptional oscillations from mouse liver, NIH3T3, and U2OS cells. The goal is to identify genes whose transcript levels oscillate during the cell cycle. These cycles are crucial as they play a key role in regulating metabolism and liver function, e.g., many liver genes operate on daily cycles, influencing vital processes such as detoxification, energy metabolism, and hormone regulation.
hughes2009harmonics collected liver tissue samples from mice at hourly intervals over a 48-hour period, pooling samples from 3--5 mice at each point. Their gene-level expression data set is available from Gene Expression Omnibus (GEO). In this section, we focus on the liver data set (accession GSE11923). We extracted the data using the GEOquery and BiocManager R packages. The final data set contains $p=45,101$ genes. For each gene $j$, we observe $n=48$ transcript level measurements $Y_{j,1},\ldots,Y_{j,n}$ at different points in time, recorded in $X_1,\ldots,X_n$.
The goal is to test whether gene transcriptions $Y_j$ are independent of the time of measurement $X$ and, in particular, to identify genes that violate independence.
We apply our proposed bootstrap test based on the studentised test statistic BMB1 described in Remark (ref), combined with the stepwise procedure developed in Section (ref), both implemented via the R package \href{https://github.com/mauolivares/hdIndep}{hdIndep}. We set $\alpha=0.05$, the nominal level at which the family-wise error rate is to be controlled, and the number of bootstrap samples to $B=1,000$. We choose the optimal block size developed in Section (ref), which in this application is $q^*(n)=1$.
The stepdown procedure identifies $4,554$ genes that violate the independence hypothesis. These correspond to approximately $10.1\%$ of all gene transcripts. The stepdown procedure rejects $4,052$ hypotheses in the first and $502$ in the second step, demonstrating the power gains from employing multiple steps.
The original study by hughes2009harmonics identified $3,667$ gene transcripts showing oscillatory behaviour, but based on a different methodology for testing a different hypothesis than ours.\footnote{The authors identify oscillatory patterns by testing for hidden periodicities in gene transcriptions using Fisher's G and straume2004dna tests. These tests are designed to test the null that the data are generated by a Gaussian white noise process against the alternative that the data is generated by a Gaussian white noise with a deterministic sinusoidal component. The authors then obtain so-called q-values following storey2003statistical so as to select gene transcripts whose q-values are less than $\alpha=0.05$.} Their procedure is guaranteed to control the false discovery rate at $\alpha=0.05$. Interestingly, our procedure finds more genes than the original study while at the same time providing stronger guarantees in the form of family-wise error rate control.
Among the genes identified by our method, $1,919$ were also identified in the original study, thus demonstrating substantial overlap. On the other hand, it singled out $2,635$ genes not identified in the original study (about $5.8\%$ of all gene transcripts), thus highlighting its ability to detect novel rhythmic gene activity. Figure (ref) shows transcript level measurements of a random sample of six genes identified by our method, but not in the original study. A complete list of identified gene transcripts is available in the replication code (see replication vignette in accompanying R package \href{https://github.com/mauolivares/hdIndep}{hdIndep}).
\ifthenelse{\boolean{arxiv}}{