A Kernel Test for Three-Variable Interactions

Dino Sejdinovic, Arthur Gretton, Wicher Bergsma

Introduction

An important class of high order interactions occurs when the simultaneous effect of two variables on a third may not be additive. In particular, it may be possible that X\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z and Y\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z, whereas \neg\left((X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z\right) (for example, neither adding sugar to coffee nor stirring the coffee individually have an effect on its sweetness but the joint presence of the two does). In addition, study of three-variable interactions can elucidate certain switching mechanisms between positive and negative correlation of two genes expressions, as controlled by a third gene . The presence of such interactions is typically tested using some form of analysis of variance (ANOVA) model which includes additional interaction terms, such as products of individual variables. Since each such additional term requires a new hypothesis test, this increases the risk that some hypothesis test will produce a false positive by chance. Therefore, a test that is able to directly detect the presence of any kind of higher-order interaction would be of a broad interest in statistical modeling. In the present work, we provide to our knowledge the first nonparametric test for three-variable interaction. This work generalizes the HSIC test of pairwise independence, and has as its test statistic the norm of an embedding of an appropriate signed measure to a reproducing kernel Hilbert space (RKHS). When the statistic is non-zero, all third order factorizations can be ruled out. Moreover, this test is applicable to the cases where XX, YY and ZZ are themselves multivariate objects, and may take values in non-Euclidean or structured domains.As the reader might imagine, the situation becomes more complex again when four or more variables interact simultaneously; we provide a brief technical overview in Section 4.3.

One important application of interaction measures is in learning structure for graphical models. If the graphical model is assumed to be Gaussian, then second order interaction statistics may be used to construct an undirected graph . When the interactions are non-Gaussian, however, other approaches are brought to bear. An alternative approach to structure learning is to employ conditional independence tests. In the PC algorithm , a V-structure (two independent variables with directed edges towards a third variable) is detected when an independence test between the parent variables accepts the null hypothesis, while a test of dependence of the parents conditioned on the child rejects the null hypothesis. The PC algorithm gives a correct equivalence class of structures subject to the causal Markov and faithfulness assumptions, in the absence of hidden common causes. The original implementations of the PC algorithm rely on partial correlations for testing, and assume Gaussianity. A number of algorithms have since extended the basic PC algorithm to arbitrary probability distributions over multivariate random variables , by using nonparametric kernel independence tests and conditional dependence tests . We observe that our Lancaster interaction based test provides a strong alternative to the conditional dependence testing approach, and is seen to outperform earlier approaches in detecting cases where independent parent variables weakly influence the child variable when considered individually, but have a strong combined influence.

We begin our presentation in Section 2 with a definition of interaction measures, these being the signed measures we will embed in an RKHS. We cover this embedding procedure in Section 3. We then proceed in Section 4 to define pairwise and three way interactions. We describe a statistic to test mutual independence for more than three variables, and provide a brief overview of the more complex high-order interactions that may be observed when four or more variables are considered. Finally, we provide experimental benchmarks in Section 5.

Matlab code for interaction tests considered in the paper is available at http://www.gatsby.ucl.ac.uk/~gretton/interact/threeWayInteract.htm

Interaction measure

An interaction measure associated to a multidimensional probability distribution PP of a random vector (X1,…,XD)\left(X_{1},\ldots,X_{D}\right) taking values in the product space X1×⋯×XD\mathcal{X}_{1}\times\cdots\times\mathcal{X}_{D} is a signed measure ΔP\Delta P that vanishes whenever PP can be factorised in a non-trivial way as a product of its (possibly multivariate) marginal distributions. For the cases D=2,3D=2,3 the correct interaction measure coincides with the the notion introduced by Lancaster as a formal product

where ∏j=1D′PXij∗\prod_{j=1}^{D^{\prime}}P_{X_{i_{j}}}^{*} is understood as a joint probability distribution of a subvector (Xi1,…,XiD′)\left(X_{i_{1}},\ldots,X_{i_{D^{\prime}}}\right). We will term the signed measure in (1) the Lancaster interaction measure. In the case of a bivariate distribution, the Lancaster interaction measure is simply the difference between the joint probability distribution and the product of the marginal distributions (the only possible non-trivial factorization for D=2D=2), ΔLP=PXY−PXPY\Delta_{L}P=P_{XY}-P_{X}P_{Y}, while in the case D=3D=3, we obtain

For D>3D>3, however, (1) does not capture all possible factorizations of the joint distribution, e.g., for D=4D=4, it need not vanish if (X_{1},X_{2})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}(X_{3},X_{4}), but X1X_{1} and X2X_{2} are dependent and X3X_{3} and X4X_{4} are dependent. Streitberg corrected this definition using a more complicated construction with the Möbius function on the lattice of partitions, which we describe in Section 4.3. In this work, however, we will focus on the case of three variables and formulate interaction tests based on embedding of (2) into an RKHS.

The implication (3) states that the presence of Lancaster interaction rules out the possibility of any factorization of the joint distribution, but the converse is not generally true; see Appendix C for details. In addition, it is important to note the distinction between the absence of Lancaster interaction and the total (mutual) independence of (X,Y,Z)(X,Y,Z), i.e., PXYZ=PXPYPZP_{XYZ}=P_{X}P_{Y}P_{Z}. While total independence implies the absence of Lancaster interaction, the signed measure ΔtotP=PXYZ−PXPYPZ\Delta_{tot}P=P_{XYZ}-P_{X}P_{Y}P_{Z} associated to the total (mutual) independence of (X,Y,Z)(X,Y,Z) does not vanish if, e.g., (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z, but XX and YY are dependent.

In this contribution, we construct the non-parametric test for the hypothesis ΔLP=0\Delta_{L}P=0 (no Lancaster interaction), as well as the non-parametric test for the hypothesis ΔtotP=0\Delta_{tot}P=0 (total independence), based on the embeddings of the corresponding signed measures ΔLP\Delta_{L}P and ΔtotP\Delta_{tot}P into an RKHS. Both tests are particularly suited to the cases where XX, YY and ZZ take values in a high-dimensional space, and, moreover, they remain valid for a variety of non-Euclidean and structured domains, i.e., for all topological spaces where it is possible to construct a valid positive definite function; see for details. In the case of total independence testing, our approach can be viewed as a generalization of the tests proposed in based on the empirical characteristic functions.

Kernel Embeddings

(Kernel embedding) Let kk be a kernel on Z\mathcal{Z}, and ν∈M(Z)\nu\in\mathcal{M}(\mathcal{Z}). The kernel embedding of ν\nu into the RKHS Hk\mathcal{H}_{k} is μk(ν)∈Hk\mu_{k}(\nu)\in\mathcal{H}_{k} such that ∫f(z)dν(z)=⟨f,μk(ν)⟩Hk\int f(z)d\nu(z)=\left\langle f,\mu_{k}(\nu)\right\rangle_{\mathcal{H}_{k}} for all f∈Hkf\in\mathcal{H}_{k}.

Alternatively, the kernel embedding can be defined by the Bochner integral μk(ν)=∫k(⋅,z) dν(z)\mu_{k}(\nu)=\int k(\cdot,z)\,d\nu(z). If a measurable kernel kk is a bounded function, it is straightforward to show using the Riesz representation theorem that μk(ν)\mu_{k}(\nu) exists for all ν∈M(Z)\nu\in\mathcal{M}(\mathcal{Z}).Unbounded kernels can also be considered, however . In this case, one can still study embeddings of the signed measures Mk1/2(Z)⊂M(Z)\mathcal{M}_{k}^{1/2}(\mathcal{Z})\subset\mathcal{M}(\mathcal{Z}), which satisfy a finite moment condition, i.e., Mk1/2(Z)={ν∈M(Z) : ∫k1/2(z,z) d∣ν∣(z)<∞}\mathcal{M}_{k}^{1/2}(\mathcal{Z})=\left\{\nu\in\mathcal{M}(\mathcal{Z})\,:\,\int k^{1/2}(z,z)\,d|\nu|(z)<\infty\right\} . For many interesting bounded kernels kk, including the Gaussian, Laplacian and inverse multiquadratics, the embedding μk:M(Z)→Hk\mu_{k}:\mathcal{M}(\mathcal{Z})\to\mathcal{H}_{k} is injective. Such kernels are said to be integrally strictly positive definite (ISPD) [27, p. 4]. A related but weaker notion is that of a characteristic kernel , which requires the kernel embedding to be injective only on the set M+1(Z)\mathcal{M}_{+}^{1}(\mathcal{Z}) of probability measures. In the case that kk is ISPD, since Hk\mathcal{H}_{k} is a Hilbert space, we can introduce a notion of an inner product between two signed measures ν,ν′∈M(Z)\nu,\nu^{\prime}\in\mathcal{M}(\mathcal{Z}),

Since μk\mu_{k} is injective, this is a valid inner product and induces a norm on M(Z)\mathcal{M}(\mathcal{Z}), for which ∥ν∥k=⟨⟨ν,ν⟩⟩k1/2=0\left\|\nu\right\|_{k}=\left\langle\left\langle\nu,\nu\right\rangle\right\rangle_{k}^{1/2}=0 if and only if ν=0\nu=0. This fact has been used extensively in the literature to formulate: (a) a nonparametric two-sample test based on estimation of maximum mean discrepancy ∥P−Q∥k\left\|P-Q\right\|_{k}, for samples {Xi}i=1n∼i.i.d.P\left\{X_{i}\right\}_{i=1}^{n}\overset{i.i.d.}{\sim}P, {Yi}i=1m∼i.i.d.Q\left\{Y_{i}\right\}_{i=1}^{m}\overset{i.i.d.}{\sim}Q and (b) a nonparametric independence test based on estimation of ∥PXY−PXPY∥k⊗l\left\|P_{XY}-P_{X}P_{Y}\right\|_{k\otimes l}, for a joint sample {(Xi,Yi)}i=1n∼i.i.d.PXY\left\{\left(X_{i},Y_{i}\right)\right\}_{i=1}^{n}\overset{i.i.d.}{\sim}P_{XY} (the latter is also called a Hilbert-Schmidt independence criterion), with kernel k⊗lk\otimes l on the product space defined as k(x,x′)l(y,y′)k(x,x^{\prime})l(y,y^{\prime}). When a bounded characteristic kernel is used, the above tests are consistent against all alternatives, and their alternative interpretation is as a generalization of energy distance and distance covariance .

In this article, we extend this approach to the three-variable case, and formulate tests for both the Lancaster interaction and for the total independence, using simple consistent estimators of ∥ΔLP∥k⊗l⊗m\left\|\Delta_{L}P\right\|_{k\otimes l\otimes m} and ∥ΔtotP∥k⊗l⊗m\left\|\Delta_{tot}P\right\|_{k\otimes l\otimes m} respectively, which we describe in the next Section. Using the same arguments as in the tests of , these tests are also consistent against all alternatives as long as ISPD kernels are used.

Interaction tests

Notational remarks: Throughout the paper, ∘\circ denotes an Hadamard (entrywise) product. Let AA be an n×nn\times n matrix, and KK a symmetric n×nn\times n matrix. We will fix the following notational conventions: 1\mathbf{1} denotes an n×1n\times 1 column of ones; A+j=∑i=1nAijA_{+j}=\sum_{i=1}^{n}A_{ij} denotes the sum of all elements of the jj-th column of AA; Ai+=∑j=1nAijA_{i+}=\sum_{j=1}^{n}A_{ij} denotes the sum of all elements of the ii-th row of AA; A++=∑i=1n∑j=1nAijA_{++}=\sum_{i=1}^{n}\sum_{j=1}^{n}A_{ij} denotes the sum of all elements of AA; K+=11⊤KK_{+}=\mathbf{1}\mathbf{1}^{\top}K, i.e., [K+]ij=K+j=Kj+\left[K_{+}\right]_{ij}=K_{+j}=K_{j+}, and [K+⊤]ij=Ki+=K+i.\left[K_{+}^{\top}\right]_{ij}=K_{i+}=K_{+i}.

We provide a short overview of the kernel independence test of , which we write as the RKHS norm of the embedding of a signed measure. While this material is not new (it appears in [29, Section 7.4]), it will help define how to proceed when a third variable is introduced, and the signed measures become more involved. We begin by expanding the squared RKHS norm ∥PXY−PXPY∥k⊗l2\left\|P_{XY}-P_{X}P_{Y}\right\|_{k\otimes l}^{2} as inner products, and applying the reproducing property,

where (X,Y)(X,Y) and (X′,Y′)(X^{\prime},Y^{\prime}) are independent copies of random variables on X×Y\mathcal{X}\times\mathcal{Y} with distribution PXYP_{XY}.

Given a joint sample {(Xi,Yi)}i=1n∼i.i.d.PXY\left\{\left(X_{i},Y_{i}\right)\right\}_{i=1}^{n}\overset{i.i.d.}{\sim}P_{XY}, an empirical estimator of ∥PXY−PXPY∥k⊗l2\left\|P_{XY}-P_{X}P_{Y}\right\|_{k\otimes l}^{2} is obtained by substituting corresponding empirical means into (4), which can be represented using Gram matrices KK and LL (Kij=k(Xi,Xj)K_{ij}=k(X_{i},X_{j}), Lij=l(Yi,Yj)L_{ij}=l(Y_{i},Y_{j})),

2 Three-Variable tests

As in the two-variable case, it suffices to derive V-statistics for inner products ⟨⟨ν,ν′⟩⟩k⊗l⊗m\left\langle\left\langle\nu,\nu^{\prime}\right\rangle\right\rangle_{k\otimes l\otimes m}, where ν\nu and ν′\nu^{\prime} take values in all possible combinations of the joint and the products of the marginals, i.e., PXYZP_{XYZ}, PXYPZP_{XY}P_{Z}, etc. Again, it is easy to see that these can be expressed as certain expectations of kernel functions, and thereby can be calculated by an appropriate manipulation of the three Gram matrices. We summarize the resulting expressions in Table 2 - their derivation is a tedious but straightforward linear algebra exercise. For compactness, the appropriate normalizing terms are moved inside the measures considered.

Based on the individual RKHS inner product estimators, we can now easily derive estimators for various signed measures arising as linear combinations of PXYZ,PXYPZ,P_{XYZ},P_{XY}P_{Z}, and so on. The first such measure is an “incomplete” Lancaster interaction measure Δ(Z)P=PXYZ+PXPYPZ−PYZPX−PXZPY\Delta_{(Z)}P=P_{XYZ}+P_{X}P_{Y}P_{Z}-P_{YZ}P_{X}-P_{XZ}P_{Y}, which vanishes if (Y,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X or (X,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y, but not necessarily if (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z. We obtain the following result for the empirical measure P^\hat{P}.

Analogous expressions hold for Δ(X)P^\Delta_{(X)}\hat{P} and Δ(Y)P^\Delta_{(Y)}\hat{P}. Unlike in the two-variable case where either matrix or both can be centered, centering of each matrix in the three-variable case has a different meaning. In particular, one requires centering of all three kernel matrices to perform a “complete” Lancaster interaction test, as given by the following Proposition.

The proofs of these Propositions are given in Appendix A. We summarize various hypotheses and the associated V-statistics in the Appendix B. As we will demonstrate in the experiments in Section 5, while particularly useful for testing the factorization hypothesis, i.e., for (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z\,\vee\,(X,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y\,\vee\,(Y,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X, the statistic ∥ΔLP^∥k⊗l⊗m2\left\|\Delta_{L}\hat{P}\right\|_{k\otimes l\otimes m}^{2} can also be used for powerful tests of either the individual hypotheses (Y,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X, (X,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y, or (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z, or for total independence testing, i.e., PXYZ=PXPYPZP_{XYZ}=P_{X}P_{Y}P_{Z}, as it vanishes in all of these cases. The null distribution under each of these hypotheses can be estimated using a standard permutation-based approach described in Appendix D.

∥ΔLP∥k⊗l⊗m=0\left\|\Delta_{L}P\right\|_{k\otimes l\otimes m}=0 if and only if cov[f(X),g(Y),h(Z)]=0\textrm{cov}\left[f(X),g(Y),h(Z)\right]=0 for all f∈Hkf\in\mathcal{H}_{k}, g∈Hlg\in\mathcal{H}_{l}, h∈Hmh\in\mathcal{H}_{m}.

And finally, we give an estimator of the RKHS norm of the total independence measure ΔtotP\Delta_{tot}P.

Let ΔtotP^=P^XYZ−P^XP^YP^Z\Delta_{tot}\hat{P}=\hat{P}_{XYZ}-\hat{P}_{X}\hat{P}_{Y}\hat{P}_{Z}. Then:

The proof follows simply from reading off the corresponding inner-product V-statistics from the Table 2. While the test statistic for total independence has a somewhat more complicated form than that of Lancaster interaction, it can also be computed in quadratic time.

3 Interaction for D>3𝐷3D>3

Streitberg’s correction of the interaction measure for D>3D>3 has the form

4 Total independence for D>3𝐷3D>3

In general, the test statistic for total independence in the DD-variable case is

A similar statistic for total independence is discussed by where testing of total independence based on empirical characteristic functions is considered. Our test has a direct interpretation in terms of characteristic functions as well, which is straightforward to see in the case of translation invariant kernels on Euclidean spaces, using their Bochner representation, similarly as in [28, Corollary 4].

Experiments

We investigate the performance of various permutation based tests that use the Lancaster statistic ∥ΔLP^∥k⊗l⊗m2\left\|\Delta_{L}\hat{P}\right\|_{k\otimes l\otimes m}^{2} and the total independence statistic ∥ΔtotP^∥k⊗l⊗m2\left\|\Delta_{tot}\hat{P}\right\|_{k\otimes l\otimes m}^{2} on two synthetic datasets where XX, YY and ZZ are random vectors of increasing dimensionality:

where ϵ∼N(0,0.12)\epsilon\sim\mathcal{N}(0,0.1^{2}). Thus, dependence of ZZ on pair (X,Y)(X,Y) is stronger than on XX and YY individually.

In all cases, we use permutation tests as described in Appendix D. The test level is set to α=0.05\alpha=0.05, and we use gaussian kernels with bandwidth set to the interpoint median distance. In Figure 1, we plot the null hypothesis acceptance rates of the standard kernel two-variable tests for X\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y (which is true for both datasets A and B, and accepted at the correct rate across all dimensions) and for X\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z (which is true only for dataset A), as well as of the standard kernel two-variable test for (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z, and the test for (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z using the Lancaster statistic. As expected, in dataset B, we see that dependence of ZZ on pair (X,Y)(X,Y) is somewhat easier to detect than on XX individually with two-variable tests. In both datasets, however, the Lancaster interaction appears significantly more sensitive in detecting this dependence as dimensionality pp increases. Figure 2 plots the Type II error of total independence tests with statistics ∥ΔLP^∥k⊗l⊗m2\left\|\Delta_{L}\hat{P}\right\|_{k\otimes l\otimes m}^{2} and ∥ΔtotP^∥k⊗l⊗m2\left\|\Delta_{tot}\hat{P}\right\|_{k\otimes l\otimes m}^{2}. The Lancaster statistic outperforms the total independence statistic everywhere apart from the Dataset B when the number of dimensions is small (between 1 and 5). Figure 2 plots the Type II error of the factorization test, i.e., test for (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z\,\vee\,(X,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y\,\vee\,(Y,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X with Lancaster statistic with Holm-Bonferroni correction as described in Appendix D, as well as the two-variable based test (which performs three standard two-variable tests and applies the Holm-Bonferroni correction). We also plot the Type II error for the conditional independence test for X\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y|Z from . Under assumption that X\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y (correct on both datasets), negation of each of these three hypotheses is equivalent to the presence of V-structure X→Z←YX\to Z\leftarrow Y, so the rejection of the null can be viewed as a V-structure detection procedure. As dimensionality increases, the Lancaster statistic appears significantly more sensitive to the interactions present than the competing approaches, which is particularly pronounced in Dataset A.

Conclusions

We have constructed permutation-based nonparametric tests for three-variable interactions, including the Lancaster interaction and total independence. The tests can be used in datasets where only higher-order interactions persist, i.e., variables are pairwise independent; as well as in cases where joint dependence may be easier to detect than pairwise dependence, for instance when the effect of two variables on a third is not additive. The flexibility of the framework of RKHS embeddings of signed measures allows us to consider variables that are themselves multidimensional. While the total independence case readily generalizes to more than three dimensions, the combinatorial nature of joint cumulants implies that detecting interactions of higher order requires significantly more costly computation, and is an interesting topic for future work.

References

Appendix A Proofs

Some basic matrix algebra used in this proof is reviewed in Appendix F. The proof of the following simple Lemma directly follows from the results therein.

(K+∘L+∘M)++=(K+⊤∘L+⊤∘M)++=tr(K+∘L+∘M+)=∑a=1nKa+La+Ma+\left(K_{+}\circ L_{+}\circ M\right)_{++}=\left(K_{+}^{\top}\circ L_{+}^{\top}\circ M\right)_{++}=tr(K_{+}\circ L_{+}\circ M_{+})=\sum_{a=1}^{n}K_{a+}L_{a+}M_{a+}

(K+∘L∘M+⊤)++=(KLM)++\left(K_{+}\circ L\circ M_{+}^{\top}\right)_{++}=\left(KLM\right)_{++}

where we used that A++=((K∘M)∘L+)++=((K∘M)L)++,A_{++}=\left(\left(K\circ M\right)\circ L_{+}\right)_{++}=\left(\left(K\circ M\right)L\right)_{++}, and similarly B++=((L∘M)K)++.B_{++}=\left(\left(L\circ M\right)K\right)_{++}. Also, C++=tr(K+∘L+∘M+)C_{++}=tr(K_{+}\circ L_{+}\circ M_{+}) and D++=(LMK)++D_{++}=(LMK)_{++}.

By comparing to the table of V-statistics, we obtain that:

A.2 Proof of Proposition 3

It will be useful to introduce into notation the kernel centered at a probability measure ν\nu, given by:

By expanding the population expression of the kernel norm of the joint under the kernels centered at the marginals, we obtain:

Substituting the definition of the centered kernel in (6), it is readily obtained that

where W=XW=X, YY, or ZZ (individual variable). Without loss of generality, let W=XW=X. Then,

The above is true for any joint distribution PXYZP_{XYZ}, and in particular for the empirical joint, whereby:

A.3 Proof of Proposition 4

By replacing kk, ll, mm with kernels centered at the marginals, we obtain a centered covariance operator Σ(XY)Z\Sigma_{(XY)Z}, for which

Now, consider the supremum of the three-way covariance taken over the unit balls of respective RKHSs:

and thus, ∥ΔLP∥k⊗l⊗m=0\left\|\Delta_{L}P\right\|_{k\otimes l\otimes m}=0 implies sup⁡f,g,hcov[f(X),g(Y),h(Z)]=0\sup_{f,g,h}\textrm{cov}\left[f(X),g(Y),h(Z)\right]=0. Conversely, if cov[f(X),g(Y),h(Z)]=0\textrm{cov}\left[f(X),g(Y),h(Z)\right]=0 ∀f,g,h\forall f,g,h, then Σ(XY)Z[f⊗g]≡0\Sigma_{(XY)Z}\left[f\otimes g\right]\equiv 0 ∀f,g\forall f,g, so the linear operator Σ(XY)Z\Sigma_{(XY)Z} vanishes.

Appendix B The effect of centering

This is no longer true in the three-variable case, where centering of each matrix has a different meaning. Various hypotheses and their corresponding V-statistics are summarized in Table 3. Note that the “composite” hypotheses are obtained simply by an appropriate centering of Gram matrices.

Consider the following simple example with binary variables XX, YY, ZZ with the 2×2×22\times 2\times 2 probability table given in Table 4. It is readily checked that all conditional covariances are equal, so ΔLP=0\Delta_{L}P=0. It is also clear, however, that neither variable is independent of the other two. Therefore, a test for Lancaster interaction per se is not equivalent to testing for the possibility of any factorization of the joint distribution, but our empirical results suggest that it can nonetheless provide a useful surrogate. In other words, while rejection of the null hypothesis ΔLP=0\Delta_{L}P=0 is highly informative and implies that interaction is present and no non-trivial factorization of the joint distribution is available, the acceptance of the null hypothesis should be considered carefully and additional methods to rule out interaction should be sought.

Appendix D Permutation test

A permutation test for total independence is easy to construct: it suffices to compute the value of the statistic (either the Lancaster statistic ∥ΔLP^∥k⊗l⊗m2\left\|\Delta_{L}\hat{P}\right\|_{k\otimes l\otimes m}^{2} or the total independence statistic ∥ΔtotP^∥k⊗l⊗m2\left\|\Delta_{tot}\hat{P}\right\|_{k\otimes l\otimes m}^{2}) on {(X(i),Y(σi),Z(τi))}i=1n\left\{\left(X^{(i)},Y^{(\sigma i)},Z^{(\tau i)}\right)\right\}_{i=1}^{n}, for randomly drawn independent permutations σ,τ∈Sn\sigma,\tau\in S_{n} in order to obtain a sample from the null distribution.

When testing for only one of the hypotheses (Y,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X, (X,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y, or (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z, either with a Lancaster statistic or with a standard two-variable kernel statistic, only one of the samples should be permuted, e.g., if testing for (Y,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X, statistics should be computed on {(X(σi),Y(i),Z(i))}i=1n\left\{\left(X^{(\sigma i)},Y^{(i)},Z^{(i)}\right)\right\}_{i=1}^{n}, for σ∈Sn\sigma\in S_{n}. However, when testing for the disjunction of these hypotheses, i.e., for the existence of a nontrivial factorization of the joint distribution, we are within a multiple hypothesis testing framework (even though one may deal with a single test statistic, as in the Lancaster case). To ensure that the required confidence level α=0.05\alpha=0.05 is reached for the factorization hypothesis, in the experiments reported in Figure 3, the Holm’s sequentially rejective Bonferroni method is used for both the two-variable based and for the Lancaster based factorization tests. Namely, pp-values are computed for each of the hypotheses (Y,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X, (X,Z)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y, or (X,Y)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z using the permutation test, and sorted in the ascending order p(1),p(2),p(3)p_{(1)},p_{(2)},p_{(3)}. Hypotheses are then rejected sequentially if p(l)<α4−lp_{(l)}<\frac{\alpha}{4-l}. The factorization hypothesis is then rejected if and only if all three hypotheses are rejected.

Appendix E Asymptotic behavior

Using terminology from , kernels kk and k′k^{\prime} are said to be equivalent if they induce the same semimetric on the domain, i.e., k(x,x)+k(x′,x′)−2k(x,x′)=k′(x,x)+k′(x′,x′)−2k′(x,x′)k(x,x)+k(x^{\prime},x^{\prime})-2k(x,x^{\prime})=k^{\prime}(x,x)+k^{\prime}(x^{\prime},x^{\prime})-2k^{\prime}(x,x^{\prime}) ∀x,x′\forall x,x^{\prime}. It can be shown that the Lancaster statistic is invariant to changing kernels within the kernel equivalence class, i.e., that

whenever k,k′k,k^{\prime}, l,l′l,l^{\prime} and m,m′m,m^{\prime} are equivalent pairs. From here,

Appendix F Some useful basic matrix algebra

Let AA, BB be n×nn\times n matrices. The following results hold:

[11⊤]ij=1,  ∀i,j[\mathbf{1}\mathbf{1}^{\top}]_{ij}=1,\;\forall i,j, and thus (11⊤)++=n2\left(\mathbf{1}\mathbf{1}^{\top}\right)_{++}=n^{2}

(I−1n11⊤)2=I−1n11⊤.\left(I-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}\right)^{2}=I-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}.

[A1]i=Ai+\left[A\mathbf{1}\right]_{i}=A_{i+}, [1⊤A]j=A+j\left[\mathbf{1}^{\top}A\right]_{j}=A_{+j}

(A11⊤)++=(11⊤A)++=nA++\left(A\mathbf{1}\mathbf{1}^{\top}\right)_{++}=\left(\mathbf{1}\mathbf{1}^{\top}A\right)_{++}=nA_{++}

(αA+βB)++=αA+++βB++\left(\alpha A+\beta B\right)_{++}=\alpha A_{++}+\beta B_{++}

(A11⊤B)++=A++B++\left(A\mathbf{1}\mathbf{1}^{\top}B\right)_{++}=A_{++}B_{++}.

(8): From (4), [A11⊤B]ij=Ai+B+j\left[A\mathbf{1}\mathbf{1}^{\top}B\right]_{ij}=A_{i+}B_{+j}, implying

Now, let KK be a symmetric matrix, and denote H=I−1n11⊤H=I-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top} (the centering matrix). Then:

A∘11⊤=11⊤∘A=AA\circ\mathbf{1}\mathbf{1}^{\top}=\mathbf{1}\mathbf{1}^{\top}\circ A=A

(A∘B)++=tr(AB⊤)\left(A\circ B\right)_{++}=tr(AB^{\top})

For a symmetric matrix KK and any matrix AA, (A∘K+)++=(AK)++\left(A\circ K_{+}\right)_{++}=\left(AK\right)_{++}, (A∘K+⊤)++=(KA)++\left(A\circ K_{+}^{\top}\right)_{++}=\left(KA\right)_{++}

For symmetric matrices KK, LL, (K+∘L+)++=(K+⊤∘L+⊤)++=n(KL)++\left(K_{+}\circ L_{+}\right)_{++}=\left(K_{+}^{\top}\circ L_{+}^{\top}\right)_{++}=n\left(KL\right)_{++}

For symmetric matrices KK, LL, (K+∘L+⊤)++=(K+⊤∘L+)++=K++L++\left(K_{+}\circ L_{+}^{\top}\right)_{++}=\left(K_{+}^{\top}\circ L_{+}\right)_{++}=K_{++}L_{++}.

(4):(A∘K+)++=tr(AK11⊤)=(AK∘11⊤)++=(AK)++.\left(A\circ K_{+}\right)_{++}=tr\left(AK\mathbf{1}\mathbf{1}^{\top}\right)=\left(AK\circ\mathbf{1}\mathbf{1}^{\top}\right)_{++}=\left(AK\right)_{++}. (5): (K+∘L+)++=(K+L)++=(11⊤KL)++=n(KL)++.\left(K_{+}\circ L_{+}\right)_{++}=\left(K_{+}L\right)_{++}=\left(\mathbf{1}\mathbf{1}^{\top}KL\right)_{++}=n\left(KL\right)_{++}. ∎

Denote H=I−1n11⊤H=I-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}. Then:

Let KK and LL be symmetric matrices and consider K∘HLHK\circ HLH. We obtain: