A Kernel-based Test of Independence for Cluster-correlated Data
Authors: Hongjiao Liu
Anna M Plantinga
Yunhua Xiang
Michael Wu
Venue: NeurIPS 2021 Poster (2021)
Keywords: multivariate analysis
kernel methods
independence test
correlated data

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 cluster-correlated 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.

A Kernel-based Test of Independence for Cluster-correlated Data
Hongjiao Liu
Anna M. Plantinga
Department of Biostatisties University Washington liuhj@us edu
Department of Mathematics and Statistics
Williams College
amp9@villiams. edu
Yunhua Xiang Department of Biostatistics University of Washington xiangyh@uw. edu
Michael C. Wu
Public Health Sciences Division
Fred Hutchinson Cancer Research Center mcru@fredhutch.com org
Abstract
The Hilbert-Schmidt Independence Criterion (HSIC) is a powerful kernel-based statistic for assessing the generalized dependence between two multivariate vari- ables. However, independence testing based on the HSIC not directly possible for correer-oorrlltte data. Such correlation pattern among the observations arises in many practical situations, c.g., family- and longitudinal data, and requires proper accommodation. Therefore, propose novel HSIC-based independence test to evaluate the dependence between two multivariate variables based on cluster- correlated data. Using the previously proposed empirical HSIC we derive its asymptotic distribution under the null hypothesis independence between the variables but the presence sample correlation. Based on both simulation studies and real data analysis, we show that, with clustered data, our approach effectively controls type error and has higher statistical power than competing methods.
test statistic,
Introduction
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 gene and group of traits that likely share common genetic mechanism 21 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 particular metabolic pathway [3] []. Such multivariate analyses aggregate information across variables and often more powerful than univariate analyses. Meanwhile, correlated observations arise in many practical situations. Family-basee designs are common in genetic studies [5], where multiple family members are recruited together into 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.
variety of parametric and parimearamttric methods has been proposed to study the association between one or multiple exposure variables and 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 single response vector and then apply the
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

usual methods that account for clustered data, such as generalized estimating equations (GEE) 7 and random effects models 99 10]. Such approaches generally apply now-dimensional setting [1]] and are subject to limitations typical of parametric and semi-parametric methods. Random- effect-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 limited number of dependence patterns, and are not sufficient as independence tests.
Here we base our approach on the Hilbert-Schmidt Independence Criterion (HSIC), non parametric kernel-based measure for assessing the generalized dependence between two multivariate, and potentially -igh-dimensional, variables [1]. By mapping the two variables into reproducing kernel Hilbert spaces (RKHS's), the population HSIC can be viewed as measure of maximized covariance between functions in the two RKHS's. When the RKHS's s being used are characteristic [I3], the population HSIC is zero if and only the two variables are independent. This measure makes no assumption the distributions of the variables the nature of the dependence.
The original HSIC-based independence test 1E] applies 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 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. [99]: They analyzed the association between multiple genetic variants and 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 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.
graphical model and deriving the test
In this work, we present the first HSIC-based independence test for correl--orrellted 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 clustered correlation among observations. We also examine the behavior of the test statistic under the alternative hypothesis and establish the consistency 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 range of scenarios. In an application to longitudinal microbiome- metabolite 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.
our test. Furthermore, we provide a way to approximate the asymptotic
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. 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 summarize our work, discuss the limitations of our proposed test and provide a conclusion.
2 Background
In this section, we introduce our assumption on cluster-correlated data and give an overview of the HSIC statistic.

---

2.1 General setting
Let Pxy be probability measure defined on a sample space X' V where both X and can be multi dinessoiimesi..... Let Px and Py be the marginal distributions on X' and y respectively. The variables X and are statistically independent Pxy = Px (equivalently, we can write X -. Y).
We consider a sample of clustered data {(Xj. K))) drawn from Pxy, where the pattern of clustered correlation is balanced and complete:
Assumption 2.1. The observations (X1, Yi),... , (Xn, Yn) are identically distributed according to Pxy, and can be divided into m clusters of fixed size (i.e., n = md). In particular, the m. clusters
{[xxa=a-1 Yai-at1), (Xai, Yas) )
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 {(X,, X,)
2.2 Hilbert-Schmidt Independence Criterion
We briefly review the Hilbert-Schmidt Independence Criterion (HSIC) proposed by Gretton et al. (2005a) [1]. The HSIC measures the generalized dependence between two variables X and Y, by embedding X and Y into reproducing kernel Hilbert spaces (RKHS's) s) and maximizing the covariance between functions of X and 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 X' -> R. and let Hy be an RKHS on y with associated kernel funetion ky : : Following Gretton et al. (2007) [14), the cross covariovv operator Cxy Hy - Hx can be defined such that, for any Hx and g E Hy,
RR
(f. CxYg)x = Exy (f(x) - -xXf((x)) ((((())))(())>>)) - -xyl)))))))) = Cov (f())
As shown by Gretton et al. (2005b) 200, the operator norm (i.c., the largest singular value) of CXY. Cov (f(X).((Y)) is zero if and only if X IL Y. given that the kernels kx and ky are universal (see Definition 5 of (20). In this sense, ||xxy|l is a measure of independence between X and Y.
defined by ||Cxy]l := sup x9E Hy lo Stlllll ll
The largest singular value of Cxy becomes zero when the sum of all squared singular values, denoted as the squared Hilbert-Schmidt norm [1]. is zero. Therefore, the squared Fillber-Schhidd norm of CXY, ||Cxy s. is also an independence criterion. This measure is defined as the population HSIC, which can be expressed conveniently in terms of kernel functions:
HSIC(Pxx) = |Cxx = E (X, ky YY)
+ E Exxx [xx x') Eyy Jay (Y,Y") - 2xxy Ex [lex (X, By [ky [yy
where X' is an independent copy of X. It is obvious that, if X is independent from Y, then we have USIC(Pxyy)
Furthermore, for certain characteristic 1x and by 13. HSIC( only if X TL Y. Example characteristic kernels include Gaussian kernels and Laplacian kerness 21
if and
To estimate the population HSIC from sample {(X,, Y,)) the empirical HSIC can be used:
HSIC(P.) ==
1
n kx (Xi, K,)y
Σ k.j.
n
nd
2
kx XX), (Y,, Yg)-
n
Define the kernel matrices Kx and Ky such that the (i, element of Kx x is kx (X, Xy) and the (ih element of Ky is by (Y, Y,). Then the empirical HSIC can also be written in terms of Kx and Ky:
HSIC(P.) n² tr HKy),
3

---

where H = 1117 is a centering matrix.
Both Gretton et al. (2007) [1] and Zhang et al. (2012) 22 have derived the asymptotic distribution
of HSTC(P;) under the null hypothesis of independence between X and Y as a weighted sum i.i.d. In Section S examine the asymptotic behavior of HS))(P,) based on cluster-correlated observations. It turns out that the null distribution in this case still weighted sum of chi-square variables, where the weights are now modified.
of chi-square variables, when the observations
3 HSIC for cluster-correlated data
Based on the clustered data setting in Section 2.1 we aim to test the null hypothesis Ho XX T Y using the empirical HSIC statistic. We first define some useful parameters and statistics.
Assume that the kernel matrices Kx and Ky defined in Section 2.2 are positive semi-definite. We focus our attention on the centered kernel matrices: Kx := HK x H and Kr := HKy H. Let kx and ky be the centered kernel functions derived from 4x and ky, with associated RKHS's H. and My, respectively. Note that the empirical HSIC can be written as HSIC(P,) = , tr ir(Kxe Ky).
Let 7x, be the r-th largest eigenvalue and ux.r = (ux,,(X)),....... (X,)) be the r-th eigen- vector of K x- Similarly, we define the eigenvalues Ty,""s 's and eigenvectors uy, 's for Ky. On the other hand, let Ax,- be the r-th largest eigenvalue of the kernel kx with respect to Px. with associated eigenfunction (.). such that f hx (zr, x (x))xx,( ((d')d((('))) (d')= = (x,)),,((x). Similarly, we define the eigenvalues Xy,,'s 's and eigenfunctions oy,,'s 's for the kernel ky with respect to Py.
For any fixed R E N, let let,, = E let (x) == v Ej=ı (x, (.e,xx))un (X,) Define ky,n and gy,r similarly. The upcoming theorems will rely the following assumption:
(x) For each where 7x, >
"x.,r
Assumption 3.1. Suppose that (X, < oo and (Y, Y'I < cc. Assume that, for each R E N. the classes Cx == fxr + (ex - kx.n)
x' E X3 and Cy == {y + y' EV} are Px-Donsker and Py- Py-DDonser (23). respectively. Further assume that, for each r, the functions x + x and + gy)((y) converge uniformly in probability as m DG, with their limit functions in L22 (x)) and I (2(P)), respectively.
(by ky, (),
In general, Assumption 3.1 ensures that the data-dependent eigenvalues and (elements of) cigenvectors of the kernel matrices K x and Ky converge in probability eigenvalues and eigenfunctions of the kernels xx and ky. We show in Appendix A. that the Donsker class condition in Assumption 3.] holds for Gaussian kernels. Now we can establish the asymptotic distribution of HSIC Pn)) under Hg : X IL Y based on clustered data.
Theorem 3.2. Suppose that, for two multivariate random variables X and Y, we have centered kernels xox and by with discrete eigenvalues. Suppose that Assumption 2.1 and Assumption 3.1 hold. Under the null hypothesis Ho X 11 Y, we have
HSIC(P.) -r d
Σ li + co,
where zi's ss are i..dd. standard normal variables, and l, are the solutions the eigenvalue problem
d
1
Σ E
(Σ Vxxp)rapxx,(XXxxrr.. dr.pg
d
p.g=1
for some double sequence ( 100 1r,am
E IR.
For the kernel function kx, the corresponding centered kernel function lx is: (x(x,x)) = Ax (x(xxx)) Ex [ex (X, æ²)] } Ex (x, X') ] Exx (xxx//xxxi (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 Kx and Ky in the presence of clustered data. We then adopt a strategy similar to that of Zhang et al, (2012) [22]: The test statistic n HSIC can be expressed as sum of squared terms, Dr. Q? where the terms Qr 's depend on eigenvalues and cigenvectors of the kernel matrices. We could show that Qr 's are asymptotically jointly normal with mean zero under Hg. and the asymptotic variances and covariances of these terms depend on eigenvalues and eigenfunctions of the kernels k-x and ky.
As result, the asymptotic distribution of HSIC(P,) under Ho is weighted sum of chi-square variables. In particular, we require the number of clusters, m, to sufficiently large. Knowing the null distribution of the test statistic enables to construct a statistical test at a given significance level. To examine the power of the proposed test, we further explore behavior of the test statistic when the null hypothesis is violated. The next theorem states the asymptotic behavior of HSIC(P,) under the alternative hypothesis H1 X Y.
Theorem 3.3. Suppose that, for two multivariate random variables and Y. we have centered
kernels kx and by with discrete eigenvalues. Suppose that Assumption 2.1 and Assumption 3.1 hold. If
there exists some r, N such that (X)x,
then
n. -t 3 00 as m + D0.
00.
nI
When kx and by are characteristic kernels, Condition equivalent H1 X Y.
Here Condition (2 is a sufficient condition for X V If 11 Y, then E( (oyy,((Y) = El (Y))
for all r, E N; a contrapositive, (2) implies X K Y. When characteristic kernels are used, based on the definition and property of the population HSIC, we can show that X 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 statistic smaller than d d'HSIC(P,) that converges in probability positive constant, which results n HSIC(P,) = md HSIC(Pp) going to infinity, as the number of clusters (m)) goes infinity. When the test statistic goes infinity, the rejection rate of the would we have established the consistency of the proposed
approach one. Hence, based on Theorem 3.3 test.
In practice, the weights 4 in [I of Theorem 3.2 are unknown and we need to estimate them with empirical counterparts. In similar spirit Theorem 4 of Zhang et al. (2012) (22). the following proposition provides an approximation for the asymptotic null distribution of HSIC(P,) and allows for independence testing in clustered data.
Proposition 3.4. Assume that the conditions in Theorem 3.2 hold. To test the mull hypothesis
Ho X IL Y at significance level a, compare the statistic n HSIC(P.) = Lrr((( Ky) against the (1 a)- -eqquannee of the distribution of
T
where 2, are i.i.d. standard normal variables and i 's are eigenvalues of VVT, with V [0,,-.. "ūm]- Each vector i; is obtained by vectorizing (i.e., stacking the columns of) the matrix M,, whose (r, 8)tth entry
Mirs
(XX,,..xxx, j=di-d
The proof of Proposition 3.4 provided in Appendix D, where we show that T has the same asymptotic distribution as n HSIC(P,) under Hg. Note that the cigenvalues of V an n? x n² the
matrix, are the same as the eigenvalues of İTV V. an m X matrix. In practice, we can calculate
cigenvalues of V instead to avoid excessive computational burden.
5

---

The distribution of T, which is mixture of chi-square variables, can be efficiently approximated by Davies' exact method (24] This method shown work well in previous studies 251 that have statistical tests based on a mixture of chi-square distributions.
4
Experiments
In this section, we conduct simulation studies and real data analysis to demonstrate the performance of our proposed test.
4.1
Simulation studies
4.1.1
Methods
We consider longitudinal data setting, where a set of exposure variables X E IP and a set of outcome variables Y E R2 are measured on m subjects observations {(X,. Y,) are grouped into m clusters with cluster size d = 3. To introduce correlation across different time points as well across different variables, we use a Kronecker product-based covariance structure [7]. which has often been used to model multivariate longitudinal data
time points. In other words, the
The general simulation setting is as following. For each cluster, let Tij denote the i-th variable X measured the j-th time point, for i I,... P and j 2, Let gig be defined similarly. Within each cluster, let (I11, 3 12, 113, Fpp, Typ; x23)² N5 x lap; )x), where Σx Rx Ra. with
1
PX
Px
PX
px
Pe
Pe
Rx-
R.
Pe
1
Pe
2 Pc
Pe
1
PX
PX
pxp
Here the Kronecker product. Marginally, we have imposed an exchangeable correlation structure Rxx across the p variables in X and an AR(1) correlation structure R&t across the three time points. The correlations between distinct variables at different time points products of the marginal correlations: e.g., Corr([11, 422) PXPe-
We simulate situation where single exposure (say, the r-th variable in X) affects multiple outcomes, with different effect sizes on different outcomes. Within each cluster, we use the model:
(111. V12, V13. ... M41. Mo2. Vy3) ... (z,(( = ( f (2))), ( (f Pyf(2,2), (, +e,
where 8., with s 1,... is the effect size of the chosen exposure on the s-th outcome, and (0, ),), with Σy Ry & R.s Ry is the correlation matrix for the q variables in Y. In simulations, we set q 20, = 0.5 and consider various levels of within-cluster correlation: Pe 0.3, 0.5 0.7. We also let Ry Rx.
Type I error simulation To evaluate the type I error rate (rejection rate under Ho), we let B1 Bg = 0, so that the null hypothesis H, X is true. We perform both the proposed HSIC test with proper accommodation for clustered correlation (HSICa). and the original HSIC test without any adjustment (HSICarig) as in [22] and [.. These two methods are applied to the data
{(X,, Yy) at the observation level. We consider = 500, 1000 or 1500 clusters and calculate empirical type error rates each setting based on 1000 simulated data sets.
From Model B3. both X and Y have multivariate normal distributions under Ho. While we focus on normal data here, additional Type 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 n proportion (1) = 10%, 20%, 30%, 40%) of outcomes Y depend on that exposure (with nonzero B3 (). We let the function f(x) take two forms: f(z) (Power Scenario 1) and f(x) = log((x 4) (Power

---

Scenario 2). For s = 1, ... ,7. the effect sizes B 's are generated from a Uniform (0, /25/m) distribution under Power Scenario 1, and from Uniform(O, /10/m) under Power Scenario 2.
In the power simulation, we perform HSICa 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 == -3 Ej=1 Tij.".. 3 Σj=1 "pp)² and Y = (3 Eー=1 Uij," .. .3 Ej= Daj)" and then perform the original HSIC test based on {(x; Y i=1- In the second method (HSICeat). follow the strategy of Rudra al.
[19] and concatenate the observations different time
points for each cluster: We consider the new variables X := (#11, 12, 213." 22p, T2p,, Ipss) Fpp) and Y : (111, 912, M/3, Myl, M42, 843)² and then perform the original HSIC test based on {(X;", p )
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 ky (xx,212)) 22) = ky (21, cxp( - 22 113//T and the linear kernel kx (21, 22) = ky (21, 22) = 21 22- For the Gaussian kernel, the shape parameter is chosen 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 255 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.
4.1.2 Results
Table T shows the empirical type I error rates of HSIC and HSICa for normal data. The type I error rate of HSICorig is inflated in each setting, where the inflation becomes greater as the within- cluster correlation (pe) increases. In contrast, HSICa has well-controlled type I error rate across all levels of within-cluster correlation. Using the linear kernel, HSICa has type I error rates close to the nominal a in all situations. Using the Gaussian kernel, HSICa conservative when the number of clusters is moderate (m 500), but its type I error rate gets close to the nominal a at a larger sample size (m 1500). For non-normal data (Figure F1 in Appendix), the pattern similar: HSICa is able to control the type error rate, either with the Gaussian kernel or with the linear kernel.
HSICa 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 maaussiannkerel-baasee in the weighted sum of chi-square variables (as the weights depend on cigenvalues), which might aggregate more finite-sample errors and make the test statistic converge slower to the asymptotic distribution.
HSICa thus involves more terms
Figure I shows the empirical power of HSICa and the two competing methods under Power Scenario 1. In all situations, HSICa has a higher power than both HSICmenn and HSICcatt regardless of the level of within-cluster correlation the kernel being used. In addition, the power of HSIC, improves quickly as 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 comparable performance as the linear kernel.
Figure 2 shows the empirical power under Power Scenario 2, where X and Y have non-linear relationship. Similar to Power Scenario 1, for both the Gaussian kernel and the linear kernel, HSICa achieves a higher power than the competing methods under all levels of within-cluster correlation. When compared between kernels, HSICa based on the Gaussian kernel is more powerful than HSICa based the linear kernel in each setting, showing the advantage of the Gaussian kernel as characteristic kernel to detect general dependence patterns.
Overall, both Power Scenario 1 and 2 show the considerable power gain of HSICa over analyzing data at the cluster level. We also note that, the power gain of HSICa decreases as the within-cluster correlation increases. This is expected since there will be less pronounced information loss in
7

---

Table 1: Empirical type error rate of HSI urig and HSICa at nominal level a for normal data under simulation.
Gaussian kernel
Linear kernel
m
Pe
HSIC
HSICa
orig
HSICa
0.05
500
0.3
0.119
0,024
0.068
0,047
0.5
0.603
0.031
0.141
0.044
0.7
1.000
0.030
0.330
0.044
1000
0.3
0.115
0.029
0.070
0.043
0.5
0.591
0.034
0.166
0.054
0.7
1.000
0.034
0.348
0.043
1500
0.3
0.113
0.047
0.082
0.053
0.5
0.608
0.043
0.145
0.053
0.7
1,000
0.044
0.352
0.052
0.01
500
0.3
0.018
0.005
0.021
0.013
0.5
0.190
0.005
0.035
0.011
0.7
0.998
0.009
0.114
0,008
1000
0.3
0.022
0.006
0.015
0.010
0.5
0.180
0.010
0.050
0.008
0.7
0.999
0.010
0.111
0.010
1500
0.3
0.027
0.007
0.019
0.010
0.5
0.209
0.008
0.047
0.010
0.7
0.999
0.008
0.117
0.009
averaging concatenating the data the cluster level observations within a cluster are highly correlated.
Pc 0.3 Gaussian
Ps 0.5, Gaussian
Pc 0.7 Gaussian
1.0
1.0
1.0
0.8
0.8
Power
Power
Power
0.4
0.4
0.4
0.2
0.2
10%
20%
30%
40%
10%
20%
30%
40%
10%
20%
30%
40%
Pc 0.3 Linear
Pc
Linear
Pc
Linear
1.0
1.0
1.0
0.8
0.8
Power
0.6
Power
0.4
0.4
0.4
0.2
0.2
0.2
10%
20%
30%
40%
10%
20%
30%
40%
10%
20%
30%
40%
Method
HSICmean
HSICCat
HSICJ
Figure 1: Empirical power of HSICa and competing methods nominal level 0.05 under Power Scenario 1. The x-axis represents the proportion of variables Y that are associated with X. The top shows results based the Gaussian kernel, and the bottom row shows results based the linear kernel.
While the above results are based on fixed cluster size, we have also investigated the effect of cluster size on performance HSICa (Appendix F.2). When the number of clusters (m) and the level of within-cluster correlation (pe) are fixed, type error control is similar for different cluster sizes (Figure F2-F3), suggesting that the convergence speed the test statistic under Ho is likely not affected by cluster size. However, larger cluster size tends result in higher statistical

---

Pc 0.3 Gaussian
Ps 0.5. Gaussian
Pc 0.7 Gaussian
1.0
1.0
1.0
0.8
0.8
Power
0.6
Power
0.5
Power
0.4
0.4
0.4
0.2
0.2
10%
20%
30%
40%
10%
20%
30%
40%
10%
20%
30%
40%
Pc 0.3. Linear
Pc
Linear
Pc
Linear
1.0
1.0
1.0
0.8
0.8
Power
0.6
Power
0.5
0.4
0.4
0.4
0.2
0.2
0.2
10%
20%
30%
40%
10%
20%
30%
40%
10%
20%
30%
40%
Method
HSICmean
HSICaat
HSICa
Figure 2: Empirical power of HSICa and competing methods nominal level 0.05 under Power Scenario 2. The x-axis represents the proportion of variables Y that are associated with X. The top row shows results based the Gaussian kernel, and the bottom row shows results based on the linear kernel.
power (Figure F4), possibly due to an increase in the overall sample size, which allows for additional information gain.
4.2 Application to real data
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 the host metabolic environment, and identifies potential metabolic biomarkers for vaginal conditions [4]. Here we apply
HSICa 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 (MsFLASSH Vaginal Health Trial 27
The MsFLASH trial was 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.I for more details), As part of an effort to investigate the mechanism of postmenopausal vaginal discomfort, vaginal mierobiota 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 HSICa HSICmean and HSIC to assess the dependence between the overall vaginal microbiome composition and metabolites in each pathway, across all 95 pathways. In other words, for each test, we have = 141, 3, X 381 and Y IR9, where g is the number of metabolites in pathway, ranging from to 21 in this data set; 95 tests performed in total.
Table shows the number of metabolic pathways identified be associated with the vaginal microbiome composition, at mnnferroni-correcter significance level a 0.05/95 /95 = 5.3 10-4. methods
Due the close relationship between vaginal microbiota and vaginal metabolites, have identified considerable number of significant metabolic pathways. Still, HSICa identifies a larger number pathways than HSICmean and HSICaat. either with the Gaussian kernel or with the linear kernel. In particular, based on the Gaussian kernel, HSICa successfully identifies all the significant pathways discovered by HSICmenn and HSIC and discovers (7) additional pathways compared to HSICcat (HSICmean) (Figure H1). 67 out of 68 pathways discovered by HSICa using

---

the linear kernel are also identified by HSICa 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 possibly non-linear relationship between certain metabolites and microbial taxa abundances.
Table 2: Number metabolic pathways identified to be associated with the vaginal microbiome composition based on the MsFLASH data set (a 5.3 x 10 10=4).
Kernel
HSICmean
HSIC
HSICa
Gaussian
68
71 64
75
Linear
64
68
We focus on some of the top pathways (with high statistical significance) identified using HSICa 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 Lactobacilltu 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 HSICa but not by HSICmean or HSICcat. Therefore, our finding is consistent with previous studies bacterial metabolism, confirming the power improvement in using HSICa for scientific discovery.
Discussion
We have introduced novel kernel-based approach, HSICa, evaluate the generalized dependence between two multivariate variables based on correla-e data. Using the previously developed HSIC statistic 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 HSIC Further extension of the test for unbalanced or incomplete clustered data will be interesting for future study.
Another limitation 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-basee approaches could be a surrogate for HSICa at small sample sizes (see Appendix F.3), although their computational burden is large compared to HSICa (see Table G1). Computationally efficient adaptations of HSICa to small sample sizes, such as those proposed by Lee et al. (11 and Zhan et al. 32). would be another useful extension.
With the continuing emergence of sigh-dimensional data and the prevalence of musstercorrrelaed data in different scientific fields, our proposed test is 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 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 's life.
10

---

Acknowledgments and Disclosure of Funding
This work supported by the National Institutes of Health grant ROI-GMI29512. 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.
References
[1] K Alaine Broadaway, David J Cutler, Richard Duncan, Jacob L Moore, Erin Ware, Min A Jhun, Lawrence F Bielak, Wei Zhao, Jennifer A Smith, Patricia A Peyser, et al. A statistical approach for testing pross-phenotype effects of rare variants. The American Journal of Human Genetics, 98(3):525-540, 2016.
[2] Xiang Zhan, Ni Zhao, Anna Plantinga, Timothy A Thornton, Karen N Conneely, Michael P Epstein, and Michael Wu. Powerful genetic association analysis for common rare variants with high-dimensional structured traits. Genetics, 206 1790, 2017.
[3] Alessia Visconti, Caroline I Le Roy, Fabio Rosa, Niccolò Rossi, Tiphaine Martin, Robert Mohney, Weizhong Li, Emanuele de Rinaldis, Jordana T Bell, J Craig Venter, et al. Interplay between the human gut microbiome and host metabolism. Nature Communications, 10(1):1-10, 2019.
[4] Amy McMillan, Stephen Rulisa, Mark Sumarah, Jean M Macklaim, Justin Renaud, Jordan Bisanz, Gregory B Gloor, and Gregor Reid. A multi-platform metabolomics approach identifies highly specific biomarkers of bacterial diversity in the vagina of pregnant and non-pregnant women. Scientific Reports, 5(1):1-14, 2015.
[5] Jurg Ott, Yoichiro Kamatani, and Mark Lathrop. Family-based designs for genome-wide association studies. Nature Reviews Genetics, 12(7):465 474, 2011.
[6] Jos WR Twisk. Applied longitudinal data analysis for epidemiology: practical guide. Cambridge University Press, 2013.
[7] Andrzej T Galecki. General class of covariance structures for two repeated factors in longitudinal data analysis. Communications in Statistics-Theory and Methods, 23(11):3105 3119, 1994.
[8] Liam M O' OrBrien and Garrett M Fitzmaurice. Analysis of longitudinal multiple-source binary data using generalized estimating equations. Journal of the Royal Statistical Society: Series (Applied Statistics), 53(1):177-193, 2004.
[9] Gregory Reinsel. Estimation and prediction a multivariate random effects generalized linear model. Journal of the American Statistical Association, 79(386):406 414, 1984.
[10] Laurel A Beckett, Daniel Tancredi, and RS Wilson. Multivariate longitudinal models for complex change processes. Statistics in Medicine, 23(2):231-239, 2004.
[11] Geert Verbeke, Steffen Fieuws, Geert Molenberghs, and Marie Davidian. The analysis of multivariate longitudinal data: review. Statistical Methods in Medical Research, 23(1):42-59, 59, 2014.
[12] Arthur Gretton, Olivier Bousquet, Alex Smola, and Bernhard Schölkopf. Measuring statistical dependence with hilbert-schmidt norms. In International Conference on Algorithunic Learning Theory, pages 63-77. Springer, 2005a.
[13] Bharath K Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert RG Lanckriet. Hilbert space embeddings and metrics on probability measures. The Journal of Machine Learning Research, 11:1517-1561, 2010.

---

[14] Arthur Gretton, Kenji Fukumizu, Choon Hui Teo, Le Song. Bernhard Schölkopf, Alexander J Smola, et al. kernel statistical test of independence. In NIPS, volume 20, pages 585-592. Citeseer, 2007.
[15] Xinhua Zhang, Le Song, Arthur Gretton, Alexander J Smola, et al. Kernel measures of independence for non-iid data. In NIPS, volume 22, pages 1937-1944, 2008.
[16] Kacper Chwialkowski and Arthur Gretton. A kernel independence test for random processes. In International Conference on Machine Learning, pages 1422-1430. PMLR, 2014.
[17] Guochang Wang, Wai Keung Li, and Zhu. New HSIC-based tests for independence between two stationary multivariate time series. arXiv preprint urXiv:/804.090 2018.
[18] Seth R Flaxman, Daniel B Neill, and Alexander Smola. Gaussian processes for independence tests with non-iid data in causal inference. ACM Transactions on Intelligent Systems and Technology (TIST), 7(2):123,, 2015.
[19] Pratyaydipta Rudra, K Alaine Broadaway, Erin B Ware, Min A Jhun, Lawrence Bielak, Wei Zhao, Jennifer A Smith, Patricia A Peyser, Sharon LR Kardia, Michael P Epstein, et al. Testing cross-phenotype effects of rare variants in longitudinal studies of complex traits. Genetic Epidemiology, 42(4):320-332, 2018.
[20] Arthur Gretton, Ralf Herbrich, Alexander Smola, Olivier Bousquet, Bernhard Schölkopf, et al. Kernel methods for measuring independence. The Journal of Machine Learning Research, 6:2075-2133 2005b.
[21] Kenji Fukumizu, Arthur Gretton, Xiaohai Sun, and Bernhard Schölkopf. Kernel measures of conditional dependence. In NIPS, volume 20, pages 489-496, 2007.
[22] Kun Zhang, Jonas Peters, Dominík Janzing, and Bernhard Schölkopf. Kernel-based conditional independence test and application in causal discovery. In Proceedings of the Twenty-Sevnnth Conference on Uncertainty in Artificial Intelligence, 1AI' page 804-813, Arlington, Virginia,
USA, 2011. AUAI Press.
[23] Aad W Van der Vaart. Asymptotic statistics. Cambridge University Press, 1998.
[24] Robert B Davies. The distribution of linear combination of x² random variables. Journal of the Royal Statistical Society: Series C (Applied Statisties), 29(3):323-33) 3333-333, 1980.
[25] Michael C Wu, Seunggeun Lee, Tianxi Cai, Yun Li, Michael Boehnke, and Xihong Lin. Rare- variant association testing for sequencing data with the sequence kernel association test. The American Journal of Human Genetics, 89(1):22-93, 2011.
[26] Sujatha Srinivasan, Martin T Morgan, Tina L Fiedler, Danijel Djukovic, Noah G Hoffman, Daniel Raftery, Jeanne M Marrazzo, and David N Fredricks. Metabolic signatures of bacterial vaginosis. MBio, 6(2), 2015.
[27] Caroline M Mitchell, Susan D Reed, Susan Diem, Joseph C Larson, Katherine M Newton, Kristine E Ensrud, Andrea Z LaCroix, Bette Caan, and Katherine A Guthrie. Efficacy of vaginal estradiol or vaginal moisturizer vs placebo for treating postmenopausal vulvovaginal symptoms: a randomized clinical trial. JAMA Internal Medicine, 178(5):681-699) 2018.
(28) Caroline M Mitchell, Nanxun Ma, Alissa J Mitchell, Michael C Wu, DJ Valint, Sean Proll, Susan D Reed, Katherine Guthrie, Andrea Z Lacroix, Joseph C Larson, et al. Association be- tween postmenopausal vulvovaginal discomfort, vaginal microbiota, and mucosal inflammation. American Journal of Obstetrics and Gynecology, 2021.
[29] Steven S Witkin and Iara M Linhares. Why do lactobacilli dominate the human vaginal microbiota? BJOG: An International Journal of Obstetrics & Cynaecology, 126,4);60--611, 2017.
[30] Camilla Ceccarani, Claudio Foschi, Carola Parolin, Antonietta Antuono. Valeria Gaspari, Clarissa Consolandi, Luca Laghi, Tania Camboni, Beatrice Vitali, Marco Severgnini, et al. Diversity of vaginal microbiome and metabolome during genital infections. Scientific Reports, 91):1-12, 2019.
12

---

[31] Seunggeun Lee, Mary J Emond, Michael J Bamshad, Kathleen C Barnes, Mark Rieder, Deborah Nickerson, ESP Lung Project Team, David C Christiani, Mark M Wurfel, Xihong Lin, et al. Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies. The American Journal of Human Genetics, (1 237, 2012.
[32] Xiang Zhan, Anna Plantinga, Ni Zhao, and Michael C Wu. A fast small-sample kernel indepen- dence test for microbiome community- association analysis. Biometrics, 7354)31453-1463, 2017.
13


Reviews:
[{"review_id": "tzhLKHkB7O", "reviewer": "Reviewer_8fxs", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper seems well wrtitten. I don't have many comments about the writing  or the structure.\n\nUnfortunately I believe the problem/contribution to be rather narrow. I am not sure if cluster correlated independence tests are impactful/ how an important this problem is.\n\nIn the real data experiment, could you write explictly what is X, what is Y, and what are their dimensions.\n\nI DID NOT CHECK PROOFS IN THE APPENDIX. I think 14 pages of technical appendix is challenging to completely verify. I believe it would be stronger if most of the proof could be given in the main text, or some insightful high-level discussion/sketch could be given that suggests the proof is correct.\n\nExperimental evidence that the results are correct are currently not very strong I believe. For the Gaussian kernel the p-values are not that close to alpha (although they are smaller so in terms of type 1 error it seems ok). Could you experimentally demonstrate the p-values are uniformly distributed under the null? Or take larger samples?\n\nWhile it is nice there is a check fo the type 1 error experimentally in table 1, I'm not sure how convincing one gaussian distribution is to conclude from this that the type 1 error rate is correct. \n\nOverall view: \n\nThe paper is generally well written and structured. The experiments are well organized and described. I'm a bit in doubt between weak accept/weak reject. The proofs are all in the appendix and I was not able to check them due to length. More insightful proof sketch that suggests the analysis is correct/ more extensive type 1 experiments could help in that respect. Main question for me is if cluster correlated independence tests are relevant/impactful enough for neurips.\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "-", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "pW6dt5SN0A", "reviewer": "Reviewer_WwvK", "summary": "This submission introduces a calibrated testing procedure for the\nindependence between two variables under non-iid sampling. The\nsampling is assumed to follow a clustered structure, where samples\nwithin clusters are correlated but the clusters are independent from\neach others. The procedure relies on the Hilbert-Schmidt Independence\nCriterion (HSIC), and a first contribution is an asymptotic null\ndistribution under the clustered sampling assumption. The test is also\nshown to be consistent. Experiments on simulated and real data show\nthat the proposed test guarantees a control on the type I error, and\nhas higher power than existing alternatives in a variety of scenarios.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The submission is clear and well written. The problem of non-iid data\nis an important and (in my opinion) under-addressed problem in\nstatistics and machine learning, making this work a useful and\ninteresting contribution. This contribution is significant as the\nproposed results are non-trivial. The experiments are thorough and\nillustrate the proposed results well.\n\nI would appreciate some clarification on the effect of d (the size of\nthe clusters) on the proposed results. The sampling assumption\nconstrains all clusters to be of equal size d, so that the sample size\nn is m*d where m is the number of clusters. All asymptotics are given\nin m->inf (and therefore n->inf) and these asymptotics don't seem to\ndepend on d. Is it possible to say something about non-asymptotic\nproperties of the statistic that would vary with the cluster size?\nDoes this size affect the convergence speed? Experiments with varying\ncluster sizes would also be welcome. With d=1 one recovers the iid\nsetting, and all experiments are done with d=3. Larger d will always\nincrease the dependence. Would the performance (calibration, power) of\nthe proposed procedure be affected in this case?\n\nI also have a few minor comments:\n\n- I think that the statement of Proposition 3.4 could be made more\nprecise. It currently states that under some assumptions, \"we can\ncompare the statistics [...] against the (1-alpha) quatile of the\ndistribution [...]\". I would find it more clear if the convergence\nresult was stated in the Proposition (and the consequence on practical\ntesting could be discussed after the Proposition).\n\n- line 129, the \"centered kernel functions\" could be defined more\n  clearly. The only centering defined so far has been on empirical nxn\n  Gram matrices, not on functions.\n\n- lines 234-236, I couldn't understand how X** differed from X. Both\n  are defined as (x_11, x_12, x_13, ..., x_p1, x_p2n x_p3). If they\n  are identical, this should be stated clearly. And in this case,\n  HSIC_cat is equivalent to HSIC_orig?\n\n- line 254, it could be useful to provide intuition as to why HSIC_cl\n  with Gaussian kernel tends to be conservative.\n\n- line 254, it is not clear from the table that the type I error rate\n  gets closer to the nominal alpha when m=1000, the difference with\n  m=500 is tiny. Did you observe a stronger improvement with larger m?\n\n- Figure 1: power is provided at nominal level for all methods. Are\n  the baseline methods correctly calibrated?\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors adequately addressed the limitations and potential negative societal impact of their work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "VKbIKih30wi", "reviewer": "Reviewer_wjAG", "summary": "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. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Originality** This is a useful tool for dependence testing in cases of non-i.i.d. instances. The paper's contribution is useful. The paper builds from Zhang et al. taking the assumption of equivalent correlation within each cluster structure to derive the novel asymptotic distribution.\n\n**Quality** The technical description and experiment are statistically sound. The discussion regarding small sample methods (non-asymptotic) could be expanded beyond the reference to Zhan et al. (2017). \n\nIn terms of comparisons, it seems permutation based tests could derive surrogate null distributions to handle the cases of missing data and unequal cluster size.  Especially for the simulated data, I would like to see a comparison versus permutation based methods. On the real data, it seems that the  stringent significance threshold (based on Bonferroni correction) would require a large number of permutations, but it could be useful to compare.\n\n**Clarity** The paper is well-organized and well-written. \n\n**Significance** A solid statistical test for dependence has the potential to impact scientific discovery in exploratory work such as biomarker discovery. Longitudinally, clustered data is increasing in importance in biomedical health applications.  A properly documented and released code may positively impact biostatistical practice. \n \n### After authors' response.\nI have read the authors' response and the other reviews. I maintain my overall score. I hope the authors can report computational complexity (cubic in the sample size to determine asymptotic distribution for the Gaussian kernel) could be mentioned in addition to the wall clock time. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "A statement of the computational complexity and its comparison to the computation required to run permutation tests with HSIC would be useful. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3.5"}, "extra_text": {}}, {"review_id": "JQBiVTbU97F", "reviewer": "Reviewer_myr6", "summary": "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. \n", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I have not be able to go through the theoretical part and just comment on the real experiment part. It will be nice to justify the identified pathways are biologically meaningful, especially the newly identified pathways. In line 284, it says “The original study found no additional benefit of estradiol or moisturiser over placebo.” How is this conclusion linked with results provided in this paper? Also, it says “In a second study …..”, what is the main conclusion of that paper? Is the main conclusion supported by this paper’s results?\n\nI think this test will be interesting if its predictions could be justified in the biological setting.\n\n---------------------------\nI am convinced by the authors' reply that the extra detected pathways are biologically sensible and I will raise my rating for this paper. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "yes", "review_text": "", "extra_scores": {"ethical_concerns": "no", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "7"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a test of independence for cluster correlated data based on HSIC and derive an asymptotic distribution for HSIC under cluster correlation.
In evaluations, the test is slightly conservative.
The test is applied to a microbiome study.
The paper is well organized and well written.
After discussion with the authors, only minor issues remain.

However, the reviewers agree that this is a good contribution and should be accepted.
Decision: Accept (Poster)