A Kernel-based Test of Independence for Cluster-correlated Data
Year: 2022

Abstract:
The Hilbert-Schmidt Independence Criterion (HSIC) is a powerful kernel-based statistic for assessing the generalized dependence between two multivariate variables. However, independence testing based on the HSIC is not directly possible for cluster-correlated data. Such a correlation pattern among the observations arises in many practical situations, e.g., family-based and longitudinal data, and requires proper accommodation. Therefore, we propose a novel HSIC-based independence test to evaluate the dependence between two multivariate variables based on clustercorrelated data. Using the previously proposed empirical HSIC as our test statistic, we derive its asymptotic distribution under the null hypothesis of independence between the two variables but in the presence of sample correlation. Based on both simulation studies and real data analysis, we show that, with clustered data, our approach effectively controls type I error and has a higher statistical power than competing methods.

heading: 1 Introduction
text: We are often interested in studying the dependence between two multivariate variables. For example, in genetic studies, we may want to assess the association between multiple genetic variants within a gene and a group of traits that likely share a common genetic mechanism [1, 2]. In microbiome studies, we may wish to investigate the association between the overall composition of human microbiota, including hundreds of microbial taxa, and multiple host metabolites from a particular metabolic pathway [3, 4]. Such multivariate analyses aggregate information across variables and are often more powerful than univariate analyses. Meanwhile, correlated observations arise in many practical situations. Family-based designs are common in genetic studies [5], where multiple family members are recruited together into the study. Longitudinal data are common in epidemiological studies [6], where variables of interest are measured on the subjects repeatedly over time. Such study designs introduce clustered dependence among the observations and require proper accommodation. In this work, we aim to develop an approach for assessing the dependence between two multivariate variables, based on cluster-correlated data.
A variety of parametric and semi-parametric methods has been proposed to study the association between one or multiple exposure variables and a multivariate longitudinal (or other cluster-correlated) outcome. These methods often extend upon existing tools for univariate longitudinal data. For example, many studies stack the multivariate outcome into a single response vector and then apply the
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
usual methods that account for clustered data, such as generalized estimating equations (GEE) [7, 8] and random effects models [9, 10]. Such approaches generally apply to a low-dimensional setting [11] and are subject to limitations typical of parametric and semi-parametric methods. Randomeffect-based methods require assumptions on the distribution of the multivariate outcome. GEE-based methods rely on good estimation of the correlation structure within clusters as well as across different outcome variables to achieve a high efficiency. Finally, both approaches assume parametric (often linear) relationships between the exposure and the outcome. Therefore, they can only evaluate a limited number of dependence patterns, and are not sufficient as independence tests.
Here we base our approach on the Hilbert-Schmidt Independence Criterion (HSIC), a non-parametric kernel-based measure for assessing the generalized dependence between two multivariate, and potentially high-dimensional, variables [12]. By mapping the two variables into reproducing kernel Hilbert spaces (RKHS’s), the population HSIC can be viewed as a measure of maximized covariance between functions in the two RKHS’s. When the RKHS’s being used are characteristic [13], the population HSIC is zero if and only if the two variables are independent. This measure makes no assumption on the distributions of the variables or the nature of the dependence.
The original HSIC-based independence test [14] applies to independent and identically distributed (i.i.d.) observations. Several extensions have been made to accommodate non-i.i.d. data, but none of the tests, to our knowledge, directly applies to clustered data at an observation level. Zhang et al. (2008) [15] extended the HSIC to certain sequence data, such as XOR sequence and Gaussian process, by specifying the correlation structure of the data as a graphical model and deriving the test statistic based on the maximal cliques. Chwialkowski et al. [16] and Wang et al. [17] developed HSIC-based tests to evaluate the dependence between two time series or random processes in general. Flaxman et al. [18] considered spatial and temporal data; they proposed to first use Gaussian process regression to remove dependence on space and time from the raw variables, and then perform the HSIC test on the resulting de-correlated residuals. However, their approach generally applies to independence testing between two univariate variables. A study with an aim closest to ours is by Rudra et al. [19]: They analyzed the association between multiple genetic variants and a multivariate longitudinal outcome. Rudra et al. concatenated the outcome measurements from different time points at the subject level, and then applied the HSIC test to the subject-level data. Although this is a straightforward approach to deal with clustered correlation, there could be a loss of statistical power by analyzing data at the subject/cluster level rather than observation level.
In this work, we present the first HSIC-based independence test for cluster-correlated data. Using the empirical HSIC [14] as our test statistic, we derive its asymptotic distribution under the null hypothesis of independence between the two variables but in the presence of clustered correlation among observations. We also examine the behavior of the test statistic under the alternative hypothesis and establish the consistency of our test. Furthermore, we provide a way to approximate the asymptotic null distribution of the test statistic and allow for statistical testing in practice. In simulation studies, our proposed test controls type I error rates well and has a much higher statistical power than competing methods across a range of scenarios. In an application to a longitudinal microbiomemetabolite data set, compared to other approaches, our proposed test identifies a larger number of metabolic pathways significantly associated with the overall microbiome composition, highlighting the value of our test in scientific studies.
The remaining sections are organized as following. In Section 2, we provide our background assumption on clustered data and an overview of the HSIC statistic. In Section 3, we study the asymptotic behavior of the HSIC statistic under null and alternative hypotheses, and construct a statistical test of independence for cluster-correlated data. In Section 4, we demonstrate the performance of our proposed test on both simulated and real data. In Section 5, we summarize our work, discuss the limitations of our proposed test and provide a conclusion.

heading: 2 Background
text: In this section, we introduce our assumption on cluster-correlated data and give an overview of the HSIC statistic.

heading: 2.1 General setting
text: Let PXY be a probability measure defined on a sample space X × Y , where both X and Y can be multi-dimensional. Let PX and PY be the marginal distributions on X and Y , respectively. The variables X and Y are statistically independent if PXY = PXPY (equivalently, we can write X ⊥⊥ Y ). We consider a sample of clustered data {(Xj , Yj)}nj=1 drawn from PXY , where the pattern of clustered correlation is balanced and complete: Assumption 2.1. The observations (X1, Y1), · · · , (Xn, Yn) are identically distributed according to PXY , and can be divided into m clusters of fixed size d (i.e., n = md). In particular, the m clusters{[
(Xdi−d+1, Ydi−d+1), · · · , (Xdi, Ydi) ]}m
i=1
are independent from one another while having identical within-cluster correlation structure.
The specific correlation structure among the observations in each cluster can be arbitrary. We are interested in studying the dependence between X and Y based on the sample {(Xj , Yj)}nj=1.

heading: 2.2 Hilbert-Schmidt Independence Criterion
text: We briefly review the Hilbert-Schmidt Independence Criterion (HSIC) proposed by Gretton et al. (2005a) [12]. The HSIC measures the generalized dependence between two variables X and Y , by embeddingX and Y into reproducing kernel Hilbert spaces (RKHS’s) and maximizing the covariance between functions of X and Y in the RKHS’s.
Let HX be an RKHS on X with associated kernel function (i.e., inner product in the RKHS) kX : X × X → R, and letHY be an RKHS on Y with associated kernel function kY : Y × Y → R. Following Gretton et al. (2007) [14], the cross-covariance operator CXY : HY → HX can be defined such that, for any f ∈ HX and g ∈ HY ,
〈f, CXY g〉HX = EXY [( f(X)− EX [f(X)] )( g(Y )− EY [g(Y )] )] = Cov ( f(X), g(Y ) ) .
As shown by Gretton et al. (2005b) [20], the operator norm (i.e., the largest singular value) of CXY , defined by ‖CXY ‖ := supf∈HX ,g∈HY ,‖f‖∞≤1,‖g‖∞≤1 Cov ( f(X), g(Y ) ) , is zero if and only if X ⊥⊥ Y , given that the kernels kX and kY are universal (see Definition 5 of [20]). In this sense, ‖CXY ‖ is a measure of independence between X and Y . The largest singular value of CXY becomes zero when the sum of all squared singular values, denoted as the squared Hilbert-Schmidt norm [12], is zero. Therefore, the squared Hilbert-Schmidt norm of CXY , ‖CXY ‖2HS , is also an independence criterion. This measure is defined as the population HSIC, which can be expressed conveniently in terms of kernel functions:
HSIC(PXY ) := ‖CXY ‖2HS = EXX′Y Y ′ [kX(X,X ′)kY (Y, Y ′)] + EXX′ [kX(X,X ′)]EY Y ′ [kY (Y, Y ′)]− 2EXY [ EX′ [kX(X,X ′)]EY ′ [kY (Y, Y ′)] ] ,
where X ′ is an independent copy of X . It is obvious that, if X is independent from Y , then we have HSIC(PXY ) = 0. Furthermore, for certain characteristic kX and kY [13], HSIC(PXY ) = 0 if and only if X ⊥⊥ Y . Example characteristic kernels include Gaussian kernels and Laplacian kernels [21]. To estimate the population HSIC from a sample {(Xj , Yj)}nj=1, the empirical HSIC can be used:
HSIC(Pn) := 1
n2 n∑ i,j kX(Xi, Xj)kY (Yi, Yj) + 1 n4 n∑ i,j,q,r kX(Xi, Xj)kY (Yq, Yr)
− 2 n3 n∑ i,j,q kX(Xi, Xj)kY (Yi, Yq).
Define the kernel matrices KX and KY such that the (i, j)-th element of KX is kX(Xi, Xj) and the (i, j)-th element ofKY is kY (Yi, Yj). Then the empirical HSIC can also be written in terms of KX andKY :
HSIC(Pn) = 1
n2 tr(HKXHKY ),
whereH = I − 1n11 T is a centering matrix.
Both Gretton et al. (2007) [14] and Zhang et al. (2012) [22] have derived the asymptotic distribution of HSIC(Pn) under the null hypothesis of independence between X and Y as a weighted sum of chi-square variables, when the observations are i.i.d. In Section 3, we examine the asymptotic behavior of HSIC(Pn) based on cluster-correlated observations. It turns out that the null distribution in this case is still a weighted sum of chi-square variables, where the weights are now modified.

heading: 3 HSIC for cluster-correlated data
text: Based on the clustered data setting in Section 2.1, we aim to test the null hypothesis H0 : X ⊥⊥ Y using the empirical HSIC statistic. We first define some useful parameters and statistics.
Assume that the kernel matricesKX andKY defined in Section 2.2 are positive semi-definite. We focus our attention on the centered kernel matrices: K̃X := HKXH and K̃Y := HKYH . Let k̃X and k̃Y be the centered kernel functions 1 derived from kX and kY , with associated RKHS’s H̃X and H̃Y , respectively. Note that the empirical HSIC can be written as HSIC(Pn) = 1n2 tr(K̃XK̃Y ).
Let γX,r be the r-th largest eigenvalue and uX,r = (uX,r(X1), · · · , uX,r(Xn))T be the r-th eigenvector of K̃X . Similarly, we define the eigenvalues γY,r’s and eigenvectors uY,r’s for K̃Y . On the other hand, let λX,r be the r-th largest eigenvalue of the kernel k̃X with respect to PX , with associated eigenfunction φX,r(·), such that ∫ k̃X(x, x ′)φX,r(x ′)dPX(x
′) = λX,rφX,r(x). Similarly, we define the eigenvalues λY,r’s and eigenfunctions φY,r’s for the kernel k̃Y with respect to PY .
For any fixed R ∈ N, let k̃X,R(x, x′) := ∑R r=1 λX,rφX,r(x)φX,r(x
′). For each r where γX,r > 0, let gX,r(x) := √ n
γX,r
∑n j=1 k̃X(x,Xj)uX,r(Xj). Define k̃Y,R and gY,r similarly. The upcoming
theorems will rely on the following assumption:
Assumption 3.1. Suppose that E[k̃2X(X,X ′)] < ∞ and E[k̃2Y (Y, Y ′)] < ∞. Assume that, for each R ∈ N, the classes CX := {x 7→ (k̃X − k̃X,R)2(x, x′) : x′ ∈ X} and CY := {y 7→ (k̃Y − k̃Y,R)2(y, y′) : y′ ∈ Y} are PX -Donsker and PY -Donsker [23], respectively. Further assume that, for each r, the functions x 7→ gX,r(x) and y 7→ gY,r(y) converge uniformly in probability as m→∞, with their limit functions in L2(PX) and L2(PY ), respectively.
In general, Assumption 3.1 ensures that the data-dependent eigenvalues and (elements of) eigenvectors of the kernel matrices K̃X and K̃Y converge in probability to eigenvalues and eigenfunctions of the kernels k̃X and k̃Y . We show in Appendix A.1 that the Donsker class condition in Assumption 3.1 holds for Gaussian kernels. Now we can establish the asymptotic distribution of HSIC(Pn) under H0 : X ⊥⊥ Y based on clustered data. Theorem 3.2. Suppose that, for two multivariate random variables X and Y , we have centered kernels k̃X and k̃Y with discrete eigenvalues. Suppose that Assumption 2.1 and Assumption 3.1 hold. Under the null hypothesis H0 : X ⊥⊥ Y , we have
nHSIC(Pn) = 1
n tr(K̃XK̃Y ) d→ ∞∑ t=1 `tz 2 t as m→∞, (1)
where zt’s are i.i.d. standard normal variables, and `t’s are the solutions to the eigenvalue problem
`tψt,rs
= 1
d ∞∑ p,q=1 E [( d∑ i=1 √ λX,rλY,sφX,r(Xi)φY,s(Yi) )( d∑ i=1 √ λX,pλY,qφX,p(Xi)φY,q(Yi) )] ψt,pq
for some double sequence {ψt,rs}∞r,s=1 ∈ R.
1For the kernel function kX , the corresponding centered kernel function k̃X is: k̃X(x, x′) = kX(x, x′)− EX [kX(X,x′)]− EX′ [kX(x,X ′)] + EXX′ [kX(X,X ′)].
The proof of Theorem 3.2 is provided in Appendix B. To prove the theorem, we first show the convergence of eigenvalues and eigenvectors of K̃X and K̃Y in the presence of clustered data. We then adopt a strategy similar to that of Zhang et al. (2012) [22]: The test statistic nHSIC(Pn) can be expressed as a sum of squared terms, ∑n r,s=1Q 2 rs, where the terms Qrs’s depend on eigenvalues and eigenvectors of the kernel matrices. We could show that Qrs’s are asymptotically jointly normal with mean zero under H0, and the asymptotic variances and covariances of these terms depend on eigenvalues and eigenfunctions of the kernels k̃X and k̃Y .
As a result, the asymptotic distribution of HSIC(Pn) under H0 is a weighted sum of chi-square variables. In particular, we require the number of clusters, m, to be sufficiently large. Knowing the null distribution of the test statistic enables us to construct a statistical test at a given significance level. To examine the power of the proposed test, we further explore the behavior of the test statistic when the null hypothesis is violated. The next theorem states the asymptotic behavior of HSIC(Pn) under the alternative hypothesis H1 : X 6⊥⊥ Y . Theorem 3.3. Suppose that, for two multivariate random variables X and Y , we have centered kernels k̃X and k̃Y with discrete eigenvalues. Suppose that Assumption 2.1 and Assumption 3.1 hold. If there exists some r, s ∈ N such that E[φX,r(X)φY,s(Y )] 6= 0, (2) then
nHSIC(Pn) = 1
n tr(K̃XK̃Y )
p→∞ as m→∞.
When k̃X and k̃Y are characteristic kernels, Condition (2) is equivalent to H1 : X 6⊥⊥ Y .
Here Condition (2) is a sufficient condition for X 6⊥⊥ Y : If X ⊥⊥ Y , then E[φX,r(X)φY,s(Y )] = E[φX,r(X)]E[φY,s(Y )] = 0 for all r, s ∈ N; as a contrapositive, (2) implies X 6⊥⊥ Y . When characteristic kernels are used, based on the definition and property of the population HSIC, we can show that X 6⊥⊥ Y also implies (2). The proof of Theorem 3.3 is provided in Appendix C. To prove the theorem, we show that, under Condition (2), there exists a statistic smaller than dHSIC(Pn) that converges in probability to a positive constant, which results in nHSIC(Pn) = mdHSIC(Pn) going to infinity, as the number of clusters (m) goes to infinity. When the test statistic goes to infinity, the rejection rate of the test would approach one. Hence, based on Theorem 3.3, we have established the consistency of the proposed test.
In practice, the weights `t’s in (1) of Theorem 3.2 are unknown and we need to estimate them with empirical counterparts. In a similar spirit to Theorem 4 of Zhang et al. (2012) [22], the following proposition provides an approximation for the asymptotic null distribution of HSIC(Pn) and allows for independence testing in clustered data. Proposition 3.4. Assume that the conditions in Theorem 3.2 hold. To test the null hypothesis H0 : X ⊥⊥ Y at a significance level α, we can compare the statistic nHSIC(Pn) = 1n tr(K̃XK̃Y ) against the (1− α)-quantile of the distribution of
T̃ = 1
m n2∑ t=1 ˜̀ tz 2 t ,
where zt’s are i.i.d. standard normal variables and ˜̀t’s are eigenvalues of Ṽ Ṽ T , with Ṽ = [ṽ1, · · · , ṽm]. Each vector ṽi is obtained by vectorizing (i.e., stacking the columns of) the n × n matrix M̃i, whose (r, s)-th entry is
M̃i,rs = 1√ d di∑ j=di−d+1 √ γX,rγY,suX,r(Xj)uY,s(Yj).
The proof of Proposition 3.4 is provided in Appendix D, where we show that T̃ has the same asymptotic distribution as nHSIC(Pn) under H0. Note that the eigenvalues of Ṽ Ṽ T , an n2 × n2 matrix, are the same as the eigenvalues of Ṽ T Ṽ , an m×m matrix. In practice, we can calculate the eigenvalues of Ṽ T Ṽ instead to avoid excessive computational burden.
The distribution of T̃ , which is a mixture of chi-square variables, can be efficiently approximated by Davies’ exact method [24]. This method is shown to work well in previous studies [1, 25] that have statistical tests based on a mixture of chi-square distributions.

heading: 4 Experiments
text: In this section, we conduct simulation studies and real data analysis to demonstrate the performance of our proposed test.

heading: 4.1 Simulation studies

heading: 4.1.1 Methods
text: We consider a longitudinal data setting, where a set of exposure variables X ∈ Rp and a set of outcome variables Y ∈ Rq are measured on m subjects at 3 time points. In other words, the observations {(Xj , Yj)}nj=1 are grouped into m clusters with cluster size d = 3. To introduce correlation across different time points as well as across different variables, we use a Kronecker product-based covariance structure [7], which has often been used to model multivariate longitudinal data [11].
The general simulation setting is as following. For each cluster, let xij denote the i-th variable in X measured at the j-th time point, for i = 1, · · · , p and j = 1, 2, 3. Let yij be defined similarly. Within each cluster, we let (x11, x12, x13, · · · , xp1, xp2, xp3)T ∼ N (5×13p,ΣX), where ΣX = RX⊗Rcl, with
RX =  1 ρX · · · ρX ρX 1 · · · ρX
... ...
. . . ...
ρX ρX · · · 1  p×p , Rcl =  1 ρc ρ2cρc 1 ρc ρ2c ρc 1  . Here ⊗ is the Kronecker product. Marginally, we have imposed an exchangeable correlation structure RX across the p variables in X and an AR(1) correlation structureRcl across the three time points. The correlations between distinct variables at different time points are products of the marginal correlations: e.g., Corr(x11, x22) = ρXρc. We simulate a situation where a single exposure (say, the r-th variable inX) affects multiple outcomes, with different effect sizes on different outcomes. Within each cluster, we use the model:
(y11, y12, y13, · · · , yq1, yq2, yq3)T
= (β1f(xr1), β1f(xr2), β1f(xr3), · · · , βqf(xr1), βqf(xr2), βqf(xr3))T + , (3)
where βs, with s = 1, · · · , q, is the effect size of the chosen exposure on the s-th outcome, and ∼ N (0,ΣY ), with ΣY = RY ⊗Rcl. RY is the correlation matrix for the q variables in Y . In simulations, we set p = q = 20, ρX = 0.5 and consider various levels of within-cluster correlation: ρc = 0.3, 0.5 or 0.7. We also letRY = RX .
Type I error simulation To evaluate the type I error rate (rejection rate under H0), we let β1 = · · · = βq = 0, so that the null hypothesis H0 : X ⊥⊥ Y is true. We perform both the proposed HSIC test with proper accommodation for clustered correlation (HSICcl), and the original HSIC test without any adjustment (HSICorig) as in [22] and [1]. These two methods are applied to the data {(Xj , Yj)}nj=1 at the observation level. We consider m = 500, 1000 or 1500 clusters and calculate the empirical type I error rates in each setting based on 1000 simulated data sets.
From Model (3), both X and Y have multivariate normal distributions under H0. While we focus on normal data here, additional Type I error simulations based on non-normal data are considered in Appendix F.1.
Power simulation To evaluate the power (rejection rate under H1), we randomly select one exposure variable from X to be the causal exposure, and make the first η proportion (η = 10%, 20%, 30%, 40%) of outcomes in Y depend on that exposure (with nonzero βs’s). We let the function f(x) take two forms: f(x) = x (Power Scenario 1) and f(x) = log((x− 4)2) (Power
Scenario 2). For s = 1, · · · , ηq, the effect sizes βs’s are generated from a Uniform(0, √ 25/m)
distribution under Power Scenario 1, and from Uniform(0, √ 10/m) under Power Scenario 2.
In the power simulation, we perform HSICcl and two other HSIC-based competing methods. The two competing methods analyze data at the cluster/subject level. In the first method (HSICmean), for each cluster, we take an average of observations at different time points: We consider the new variables X∗ := ( 13 ∑3 j=1 x1j , · · · , 1 3 ∑3 j=1 xpj) T and Y ∗ := ( 13 ∑3 j=1 y1j , · · · , 1 3 ∑3 j=1 yqj) T and then perform the original HSIC test based on {(X∗i , Y ∗i )}mi=1. In the second method (HSICcat), we follow the strategy of Rudra et al. [19] and concatenate the observations at different time points for each cluster: We consider the new variables X∗∗ := (x11, x12, x13, · · · , xp1, xp2, xp3)T and Y ∗∗ := (y11, y12, y13, · · · , yq1, yq2, yq3)T and then perform the original HSIC test based on {(X∗∗i , Y ∗∗i )}mi=1. We consider m = 500 clusters and calculate the empirical power in each setting based on 1000 simulated data sets.
Kernel choices For both X and Y , we consider two different kernels: the Gaussian kernel kX(z1, z2) = kY (z1, z2) = exp(−‖z1 − z2‖22/τ) and the linear kernel kX(z1, z2) = kY (z1, z2) = zT1 z2. For the Gaussian kernel, the shape parameter τ is chosen as the median of the Euclidean distance between each sample pair. While the Gaussian kernel is a characteristic kernel [21], the linear kernel is not characteristic and is designed to detect linear or close-to-linear relationships between X and Y . Nevertheless, linear kernels have been shown to be reasonably powerful in previous association studies [25, 2] and can be computationally efficient (see Appendix G.2).
Additional simulation studies are provided in Appendix F. Additional implementation details including computation time and code availability are provided in Appendix G.

heading: 4.1.2 Results
text: Table 1 shows the empirical type I error rates of HSICorig and HSICcl for normal data. The type I error rate of HSICorig is inflated in each setting, where the inflation becomes greater as the withincluster correlation (ρc) increases. In contrast, HSICcl has a well-controlled type I error rate across all levels of within-cluster correlation. Using the linear kernel, HSICcl has type I error rates close to the nominal α in all situations. Using the Gaussian kernel, HSICcl is conservative when the number of clusters is moderate (m = 500), but its type I error rate gets close to the nominal α at a larger sample size (m = 1500). For non-normal data (Figure F1 in Appendix), the pattern is similar: HSICcl is able to control the type I error rate, either with the Gaussian kernel or with the linear kernel.
HSICcl based on the Gaussian kernel is more conservative, likely because the Gaussian kernel is associated with a larger number of non-zero eigenvalues in finite samples than the linear kernel in our simulation setting. The null distribution for Gaussian-kernel-based HSICcl thus involves more terms in the weighted sum of chi-square variables (as the weights depend on eigenvalues), which might aggregate more finite-sample errors and make the test statistic converge slower to the asymptotic distribution.
Figure 1 shows the empirical power of HSICcl and the two competing methods under Power Scenario 1. In all situations, HSICcl has a higher power than both HSICmean and HSICcat, regardless of the level of within-cluster correlation or the kernel being used. In addition, the power of HSICcl improves quickly as a higher proportion of variables in Y is associated with X . Since X and Y are linearly associated in Scenario 1, it is not surprising that the linear kernel is powerful in detecting this dependence. The Gaussian kernel has a comparable performance as the linear kernel.
Figure 2 shows the empirical power under Power Scenario 2, where X and Y have a non-linear relationship. Similar to Power Scenario 1, for both the Gaussian kernel and the linear kernel, HSICcl achieves a higher power than the competing methods under all levels of within-cluster correlation. When compared between kernels, HSICcl based on the Gaussian kernel is more powerful than HSICcl based on the linear kernel in each setting, showing the advantage of the Gaussian kernel as a characteristic kernel to detect general dependence patterns.
Overall, both Power Scenario 1 and 2 show the considerable power gain of HSICcl over analyzing data at the cluster level. We also note that, the power gain of HSICcl decreases as the within-cluster correlation increases. This is expected since there will be less pronounced information loss in
averaging or concatenating the data at the cluster level if observations within a cluster are highly correlated.
While the above results are based on a fixed cluster size, we have also investigated the effect of cluster size on the performance of HSICcl (Appendix F.2). When the number of clusters (m) and the level of within-cluster correlation (ρc) are fixed, type I error control is similar for different cluster sizes (Figure F2-F3), suggesting that the convergence speed of the test statistic under H0 is likely not affected by cluster size. However, a larger cluster size tends to result in a higher statistical
power (Figure F4), possibly due to an increase in the overall sample size, which allows for additional information gain.

heading: 4.2 Application to real data
text: The vaginal microbiota plays an important role in maintaining vaginal homeostasis. Common vaginal conditions, such as bacterial vaginosis, are characterized by shifts in vaginal microbiome composition and changes in vaginal metabolites [4, 26]. Studying the association between the microbiome and the metabolites helps us better understand how the vaginal microbiota contributes to the host metabolic environment, and identifies potential metabolic biomarkers for vaginal conditions [4]. Here we apply HSICcl and competing methods to test the dependence between the overall vaginal microbiome composition and different metabolic pathways, using data from the Menopause Strategies: Finding Lasting Answers for Symptoms and Health (MsFLASH) Vaginal Health Trial [27].
The MsFLASH trial was a 12-week randomized clinical trial to evaluate the treatment effect of vaginal estradiol vs. placebo on vaginal discomfort in postmenopausal women [27] (see Appendix H.1 for more details). As part of an effort to investigate the mechanism of postmenopausal vaginal discomfort, vaginal microbiota and vaginal fluid metabolites were characterized longitudinally and available in 141 participants at baseline, 4 and 12 weeks [28]. The vaginal microbiome profiles included abundance data of 381 taxa. The metabolome profiles included abundance data of 171 metabolites, which were grouped into 95 metabolic pathways. We apply HSICcl , HSICmean and HSICcat to assess the dependence between the overall vaginal microbiome composition and the metabolites in each pathway, across all 95 pathways. In other words, for each test, we have m = 141, d = 3, X ∈ R381 and Y ∈ Rq , where q is the number of metabolites in a pathway, ranging from 1 to 21 in this data set; 95 tests are performed in total.
Table 2 shows the number of metabolic pathways identified to be associated with the vaginal microbiome composition, at a Bonferroni-corrected significance level α = 0.05/95 = 5.3× 10−4. Due to the close relationship between vaginal microbiota and vaginal metabolites, all methods have identified a considerable number of significant metabolic pathways. Still, HSICcl identifies a larger number of pathways than HSICmean and HSICcat, either with the Gaussian kernel or with the linear kernel. In particular, based on the Gaussian kernel, HSICcl successfully identifies all the significant pathways discovered by HSICmean and HSICcat, and discovers 4 (7) additional pathways compared to HSICcat (HSICmean) (Figure H1). 67 out of 68 pathways discovered by HSICcl using
the linear kernel are also identified by HSICcl using the Gaussian kernel (Figure H2). For this data set, the Gaussian kernel appears to be more powerful in detecting dependence than the linear kernel, indicating a possibly non-linear relationship between certain metabolites and microbial taxa abundances.
We focus on some of the top pathways (with high statistical significance) identified using HSICcl and highlight their biological relevance. The top pathways include multiple metabolic pathways for amino acids. The human vaginal microbiota is dominated by bacteria in the Lactobacillus genus [29], which are known to produce branched-chain amino acids including valine, leucine and isoleucine [30]. All these amino acids are present in our top pathways. In particular, one pathway related to leucine metabolism is only identified by HSICcl but not by HSICmean or HSICcat. Therefore, our finding is consistent with previous studies on bacterial metabolism, confirming the power improvement in using HSICcl for scientific discovery.

heading: 5 Discussion
text: We have introduced a novel kernel-based approach, HSICcl, to evaluate the generalized dependence between two multivariate variables based on cluster-correlated data. Using the previously developed HSIC statistic as our test statistic, we have derived its asymptotic null distribution in the presence of clustered correlation and constructed a statistical test of independence accordingly. We have also established the consistency of the proposed test. Both simulation studies and application to real longitudinal data demonstrate the power gain in using our proposed test, compared to methods based on measurements averaged or concatenated at the cluster level.
One limitation of our framework is that the proposed test only applies to balanced and complete clustered data, which might not be always available in practice. In longitudinal studies, for example, subjects might be followed at different time points from one another (resulting in unbalanced data), or become lost to follow-up (resulting in incomplete data). For incomplete data, one solution is to impute the missing data before applying HSICcl. Further extension of the test for unbalanced or incomplete clustered data will be interesting for future study.
Another limitation is that our proposed test relies on asymptotic results, and the null distribution might not be accurately approximated when the number of clusters is small, which is likely true of many family-based or longitudinal studies. Permutation-based approaches could be a surrogate for HSICcl at small sample sizes (see Appendix F.3), although their computational burden is large compared to HSICcl (see Table G1). Computationally efficient adaptations of HSICcl to small sample sizes, such as those proposed by Lee et al. [31] and Zhan et al. [32], would be another useful extension.
With the continuing emergence of high-dimensional data and the prevalence of cluster-correlated data in different scientific fields, our proposed test is a promising approach to discover novel associations and bring new scientific insights in these settings. Meanwhile, we need to be cautious about potential risks to society that might result from misuse or misinterpretation of our proposed test. For example, confounding is an important factor to consider in genetic and epidemiological studies. A confounder affects both the exposure and the outcome, and could lead to spurious associations between the two variables even if the variables themselves do not have causal relationships. Therefore, as one applies our proposed test to evaluate the association between two variables, it is important to consider the presence of potential confounders and be careful in interpreting the test results. Mistaking certain observed correlation for causation could lead to misinformation in the scientific community and would be especially concerning when the studies being conducted directly influence people’s life.
Acknowledgments and Disclosure of Funding
This work is supported by the National Institutes of Health grant R01-GM129512. The MsFLASH data were provided by the Fred Hutchinson Cancer Research Center (MsFLASH Network), which is supported by the National Institute on Aging grant 5R01-AG048209. We would like to thank the anonymous reviewers for their helpful feedback.
The authors declare no competing interests.

Review summaries:
The paper proposed a modified version of Hilbert-Schmidt Independence Criterion (HSIC) to test the dependence between two multivariate variables. The method is suitable for data that is not iid, i.e. correlation exists between different data points. The paper first demonstrate the asymptotic distributions of HSIC under the null hypothesis (independent) and the condition for the alternative in Theorem 3.2 and 3.3. The approximation of the asymptotic HSIC distribution is provided in proposition 3.4. In the simulation, the proposed method HSICcl shows a better performance than earlier versions of HSICorig both in terms of Type I error and empirical power. In real data, HSICcl identifies more associated pathways than previous versions.  

The paper derives an approximation of the asymptotic null distribution of the Hilbert-Schmidt Independence Criterion under a matched and complete clustering. The null distribution corresponds to independence between the variables.  The samples consist of equally sized clusters each with matched within cluster correlation. The distribution is described as a weighted sum of chi-squared variables, where the weights are approximated by the eigenvalues of a matrix composed from the eigenstructure of the marginal kernel matrices grouped by cluster. Simulations confirm the increased power of the method. Specifically, simulation shows the test is conservative in lower sample size, but increases in power versus other methods (concatenating the samples longitudinally).  Experimental results demonstrate the increased power is able to identify additional metabolite pathways are associated with the vaginal microbiomes composition. 

This submission introduces a calibrated testing procedure for the independence between two variables under non-iid sampling. The sampling is assumed to follow a clustered structure, where samples within clusters are correlated but the clusters are independent from each others. The procedure relies on the Hilbert-Schmidt Independence Criterion (HSIC), and a first contribution is an asymptotic null distribution under the clustered sampling assumption. The test is also shown to be consistent. Experiments on simulated and real data show that the proposed test guarantees a control on the type I error, and has higher power than existing alternatives in a variety of scenarios. 

The authors propose a test of independence for cluster correlated data. The test statistic is still HSIC. However they derive a different asymptotic distribution under cluster correlation.