Multivariate Nonparametric Tests of Independence

New test statistics are proposed for testing whether two random vectors are independent. Gieser and Randles, as well as Taskinen, Kankainen, and Oja have introduced and discussed multivariate extensions of the quadrant test of Blomqvist. This article serves as a sequel to this work and presents new multivariate extensions of Kendall's tau and Spearman's rho statistics. Two different approaches are discussed. First, interdirection proportions are used to estimate the cosines of angles between centered observation vectors and between differences of observation vectors. Second, covariances between affine-equivariant multivariate signs and ranks are used. The test statistics arising from these two approaches appear to be asymptotically equivalent if each vector is elliptically symmetric. The spatial sign versions are easy to compute for data in common dimensions, and they provide practical, robust alternatives to normal-theory methods. Asymptotic theory is developed to approximate the finite-sample null distributions as well, as to calculate limiting Pitman efficiencies. Small-sample null permutation distributions are also described. A simple simulation study is used to compare the proposed tests with the classical Wilks test. Finally, the theory is illustrated by an example.


INTRODUCTION
In many problem settings, multiple measurements are made on each experimental unit, resulting in high-dimensional multivariate data. It is often of interest to explore potential relationships among subsets of these measurements. For example, some measurements may represent attributes of psychological characteristics, whereas others may represent attributes of physical characteristics. It may be of interest to determine whether there is a relationship between the psychological and physical characteristics. This requires a test of independence between pairs of vectors, where the vectors potentially have different measurement scales and dimensions.
Let (2)T i is a continuous vector of dimension q. We seek to test H 0 : x (1) i and x (2) i are independent versus H a : x (1) i and x (2) i are dependent. The classical parametric test due to Wilks (1935) is based on i ) T , for s, t = 1, 2. Under the assumption of multivariate normality, Wilks' test is optimal, that is, the most efficient test. Under H 0 with finite fourth moments, −n log W d → χ 2 pq . A nonparametric analog to the Wilks test was given by Puri and Sen (1971). These authors developed a class of tests based on componentwise ranking that uses a test statistic of the form Here the elements of ( p + q) × ( p + q) matrix T are where C si denotes the rank of the sth component of x i among the sth components of all n vectors, J s and J t are arbitrary (standardized) score functions, and T is partitioned in the same manner as in the Wilks test. Under H 0 , −n log S J d → χ 2 pq . Muirhead (1982) examined the effect of the group of transformations {x → Ax + b} on this problem. Here b is a p + q vector and is an arbitrary nonsingular matrix with p × p matrix A 1 and q × q matrix A 2 . The Wilks test is invariant under this group of transformations. Thus its performance does not depend on the variance-covariance structure of either x i . The test given by Puri and Sen is not invariant under this group of transformations. Gieser and Randles (1997) proposed a nonparametric test based on interdirection counts that generalized the univariate ( p = q = 1) quadrant test of Blomqvist (1950) and is invariant under this transformation group. Taskinen, Kankainen, and Oja (2003a) proposed a more practical invariant extension of the quadrant test based on spatial signs that is easy to compute for data in the dimensions (say, <15) commonly encountered in practice. These two invariant quadrant test extensions share certain asymptotic properties, for example, the same asymptotic efficiencies.
In this article we develop tests that generalize the popular univariate tests due to Kendall (1938) and Spearman (1904) to any dimensions (arbitrary p and q). Our tests have advantages over the quadrant test extensions in that they do not require centering (i.e., subtracting a location estimator) and generally have better power properties than the quadrant test extensions. Moreover, the spatial sign versions are easy to compute for data in common dimensions, thus providing intuitive, practical, robust alternatives to multivariate normal-theory methods.
The univariate tests based on Kendall's tau, Spearman's rho, and Blomqvist's quadrant statistic are described in Section 2. Their multivariate analogs based on interdirections are also sketched. In Section 3 the more practical versions of these tests based on spatial signs and ranks are described. Section 4 shows their large-sample properties under the null hypothesis as well as under a sequence of alternatives under mild assumptions (e.g., even the assumption of the existence of first moments can be avoided). Asymptotic efficiencies are given and simulations are used to compare the finite-sample powers of these tests in Section 5. The theory is illustrated by an example in Section 6, and the article concludes with some comments in Section 7.

ANALOGS TO KENDALL'S TAU AND SPEARMAN'S RHO
Consider the problem of measuring dependence within the components of bivariate vectors. Let (x (1) i , x (2) i ) T , i = 1, . . . , n, be a random sample from a bivariate, continuous population. Using the univariate sign function S(x) = sign(x) = −1, 0, 1 as x <, =, > 0, the sign of the univariate x i , µ (2) , and S (2) i analogously, the popular nonparametric measures of dependence are now conveniently expressed. They are Blomqvist's quadrant statistic, ij , and Spearman's rho, The test statistics are thus covariances (or correlations) between centered signs, signs of the pairwise differences, and centered ranks.
When correlating univariate pairs, the most interpretable feature is the magnitude of the correlation, that is, the square of the correlation coefficient. Multivariate measures of correlation provide analogs to where the averages are computed over all possible indices. We can now interpret S (1) i j , for example, as the cosines of the angles between x (1) j . In the univariate cases, the cosine values are +1 or −1. Randles (1989) used so-called "interdirection proportions" to estimate the cosines of angular distance between two vectors relative to the positions of the other vectors. In the current context, for example, the cosine of the angle between x (1) are on opposite sides of the hyperplane. (Here i * and j * would not use any of {i, j, i or j }.) Defining similar proportions among x (2) vector differences yields τ 2 1 = ave cos πp (1) (i, j; i , j ) cos πp (2) (i, j; i , j ) and ρ 2 1 = ave cos πp (1) (i, j; i , j ) cos πp (2) (i, k; i , k ) as natural multivariate analogs to Kendall's τ 2 and Spearman's ρ 2 . Also here the averages are computed over all possible indices. Kendall's tau and Spearman's rho analogs provide natural counterparts to the quadrant test analog.
Q 2 1 = ave cos π p (1) (i; i ) cos π p (2) (i; i ) , which was studied by Gieser and Randles (1997), where, for example, p (1) (i; i ) is the fraction of hyperplanes formed by p − 1 vectors x are on opposite sides of the hyperplane, where i * is not equal to either i or i . Here the observations are centered on some affineequivariant location estimator µ (1) based on x n , for example, the Oja (1983) median or the transformationretransformation spatial median (Hettmansperger and Randles 2002). The Kendall and Spearman analogs do not require centering on a location estimator.
Although Q 2 1 , τ 2 1 , and ρ 2 1 are meaningful extensions of the quadrant test statistic, Kendall's tau and Spearman's rho, they are difficult to compute if either p or q exceeds 2. More practical extensions are described in the next section.

TESTS BASED ON SPATIAL SIGNS AND RANKS
Intuitive, computationally convenient extensions of Kendall's and Spearman's statistics can be created using notions of spatial signs and spatial centered ranks as developed by Möttönen and Oja (1995). To make the resulting correlation measures affine invariant, the data points are transformed before spatial signs and ranks are formed. Consider a sample x 1 , . . . , x n , where each vector is of dimension k. The data points are transformed into where µ is an affine-equivariant location estimator and V is a shape matrix estimate.
The spatial sign function is defined as where x = (x T x) 1/2 is the (Euclidean) length of the vector x. Applying the spatial sign function to the transformed data points produces standardized sign vectors S i = S(z i ) and S ij = S(z i − z j ). The standardized rank vector is then R i = ave j S ij . Because the S ij 's and R i 's are based on differences, they do not depend on µ at all, unless µ plays a role in the computation of V. The following lemma is proved in Appendix A.
Lemma 1. The standardized sign vectors S i and S ij and standardized rank vector R i are affine equivariant (and location invariant) in the sense that if S * i , S * ij , and R * i are calculated from

The multivariate extension of Kendall's tau is created by forming sign vectors S
(1) ij based on the first components where I p is the p × p identity matrix. This is the transformation studied by Tyler (1987) computed only on differences x (also see Dümbgen 1998). With this choice, the The multivariate Kendall's tau squared is then where · 2 = Tr(· T ·) is the squared matrix norm. τ 2 2 is not the same as τ 2 1 , but like τ 2 1 it can be described as a correlation of cosines. Unlike τ 2 1 , the statistic τ 2 2 is easily computed for data in common dimensions.
The multivariate extension of Spearman's rho uses R (1) i based on the first components x (1) n transformed by a shape matrix V (1) chosen so that (See Visuri, Oja, and Koivunen 2000 for the development of these rank covariance matrices.) With this choice, the rank vectors R (1) i do not depend on µ (1) in any way. With analogous descriptions of the q-dimensional rank vectors R (2) i , the multivariate Spearman's rho squared is then It is not the same as ρ 2 1 , but it has a similar interpretation as a correlation of cosines. The statistic ρ 2 2 is easy to compute for data in common dimensions. This property makes it much more practical than ρ 2 1 . The quadrant statistic of Taskinen et al. (2003a) The estimates µ (1) and V (1) are the transformation-retransformation spatial median and Tyler's shape estimate as proposed by Hettmansperger and Randles (2002). With analogous descriptions of the q-dimensional standardized signs S i , the test of Taskinen et al. test is based on the statistic Again, Q 2 2 is much easier to compute than Q 2 1 and thus is more practical.

LIMITING DISTRIBUTIONS
To establish a large-sample distribution for the statistics described in Section 3, we first examine the statistics behavior under the null hypothesis condition that x (1) is independent from x (2) . Assume now that x (1) has a continuous elliptically symmetric distribution with density function of the form where µ 1 denotes the location of the distribution and 1 defines the shape of its contours. Similarly, assume that x (2) has an elliptically symmetric distribution with density In the following we also need to assume that the location and shape estimates used in the standardization are √ n-consistent. Note that if, for example, the standardization for τ 2 2 and ρ 2 2 is based on Dümbgen's (1998) shape estimate, then no assumptions on the existence of first moments of f and g are needed.
Theorem 1. Under the null hypothesis with elliptically sym- Also, small-sample simulations showed only minor differences in the performance of the statistics within each pair. Thus, in what follows we exclusively examine the more practical versions τ 2 2 , ρ 2 2 , and Q 2 2 . Theorem 2. Under the null hypothesis assumptions described earlier, the limiting distributions of npq and npqQ 2 2 are chi-squared distributions with pq degrees of freedom. The constants c 2 F and c 2 G depend on the marginal distributions F and G.
The constants c 2 F and c 2 G depend on the underlying distributions 1 and 2 , that is, the distributions of the respective standardized component radius. When 1 = I, the centered rank and c 2 G is defined analogously. Appendix B provides expressions of c 2 F at selected distributions. To use τ 2 2 and ρ 2 2 as test statistics for detecting correlation, obviously the constants c 2 F and c 2 G must be consistently estimated, because they depend on the respective marginal distributions. Fortunately, doing this is not difficult, because convergent estimates are readily available.
To consider and compare the efficiencies of different test statistics, we next derive the limiting distributions of τ 2 2 , ρ 2 2 , and Q 2 2 under certain contiguous alternative sequences. We use similar sequences of alternatives as those used by Gieser and Randles (1997) and Taskinen et al. (2003a). Let x (1) i and x (2) i be independent with spherical marginal densities exp{− 1 ( x (1) )} and exp{− 2 ( x (2) )} and write, for some choices of M 1 and M 2 , The optimal likelihood ratio test statistic for testing H 0 against H is then We need the general assumption that under the null hypothesis, where This should be checked separately in each case. Appendix A provides the formula for k i . Under this assumption, the sequence of alternatives is then contiguous to the null hypothesis (LeCam's first lemma).
Theorem 3. For max( p, q) > 1, the limiting distributions of [npq/4c 2 F c 2 G ]τ 2 * 2 and [npq/c 2 F c 2 G ]ρ 2 * 2 are noncentral chi-squared distributions with pq degrees of freedom and noncentrality parameter The limiting distribution of npqQ 2 * 2 is a noncentral chi-squared distribution with pq degrees of freedom and noncentrality parameter (For the case where p = q = 1, see Taskinen et al. 2003b.) Using these noncentrality parameters in asymptotic efficiencies requires some clarifying comments. We are considering efficiencies for the dependent alternatives in the "directions" determined by fixed M 1 and M 2 . The alternative distributions are not elliptical (except in the normal case), but the sequence is contiguous to an elliptical null distribution. Next, note that the efficiencies are of the same type h 1 M 1 + h 2 M T 2 2 , where h 1 and h 2 depend on the marginal spherical distributions and on the test used. If the marginal distributions are the same (which implies that h 1 = h 2 ), then the efficiencies do not depend on M 1 and M 2 at all. We say that the test is "asymptotically unbiased" if the noncentrality parameter is positive. Our tests are clearly not asymptotically unbiased for choices M 1 and M 2 such that h 1 M 1 = −h 2 M T 2 . The only exception is again the case where both marginal distributions are multivariate normal; in that case M 1 = −M T 2 implies uncorrelateness as well as independence. Finally, note that if the marginal distributions are elliptical, then, due to affine invariance, the efficiencies are given by h 1 In the next section our efficiency comparisons use M 1 = M T 2 . As test statistics for testing the null hypothesis of independence, none of the τ 2 2 , ρ 2 2 , or Q 2 2 has an unconditional distribution-free property for small n. However, their permutation null distributions, and hence permutation p values, are easy to generate. The permutation principle dictates that all possible n! pairings of the given set of n first-component x (1) vectors with the set of n second-component x (2) vectors are equally likely. Because computation of marginal standardized signs and ranks does not depend on the pairing, the null permutation distributions can be generated by finding where α = (α 1 , . . . , α n ) T is drawn at random from the n! permutations of the integers (1, . . . , n).

LIMITING AND FINITE-SAMPLE EFFICIENCIES
We next compare τ 2 2 and ρ 2 2 with Wilks' test W through limiting efficiencies and simple simulation studies. (For comparisons between quadrant statistics and W, see Gieser andRandles 1997 andTaskinen et al. 2003a.) Because the limiting distributions are of the same type, χ 2 pq , the efficiency comparisons may be based on noncentrality parameters only. The efficiency comparisons are now made in the multivariate normal distribution, t distribution, and contaminated normal distribution cases. For simplicity, we assume that M 1 = M T 2 . Because W has a limiting noncentral chi-squared distribution with pq degrees of freedom and noncentrality parameter δ 2 M 1 + M T 2 2 , the asymptotic efficiencies are given simply by with d 1 and d 2 as given in Theorem 3. For a k-dimensional multivariate normal distribution with cdf , . For a k-dimensional standardized t-distribution with ν degrees of freedom, .
The resulting limiting efficiencies for selected degrees of freedom and dimensions are listed in Table 1. Note that because limiting multinormality of the regular covariance matrix holds if the fourth moments of the underlying distribution are finite, W has a limiting distribution only when ν > 4. When the underlying distribution is multivariate normal (ν = ∞), the Wilks test is the best, but τ 2 2 and ρ 2 2 are very competitive with it. As the underlying population becomes heavy-tailed (i.e., ν gets smaller), τ 2 2 and ρ 2 2 are better than the Wilks test. Next consider the contaminated normal distribution with cdf , and limiting efficiencies for = .1 and for selected values of c are listed in Table 2. Here the superior performance of the nonparametric tests is obvious. A simulation study was performed to compare the finitesample powers of τ 2 2 , ρ 2 2 , and W. (For comparisons between quadrant statistics, see Taskinen et al. 2003a.) Independent x (1) and x (2) samples of sizes n = 30 and 60 were generated from a multivariate standard normal distribution, from a standardized t distribution with ν 1 = ν 2 = 5, and from a contaminated normal distribution with = .1 and c = 6. The transformation in (1) with M 1 = M T 2 = I was used for chosen values of δ, to introduce dependence into the model. The tests were applied, and the process was replicated 1,500 times. The critical values used in test construction were based on the chi-squared approximations to the null distributions as given in Theorem 1 (with estimated c 2 F and c 2 G ). The empirical powers for p = q = 3 are given in Figure 1, and those for p = q = 5 are given in Figure 2. First, consider the simulation results in the case where p = q = 3. In the multivariate normal case, W is slightly better than the other tests. In the t distribution case considered, the tests are almost equally powerful, and in the contaminated normal case, τ 2 2 and ρ 2 2 perform much better than W. As n = 30, ρ 2 2 seems to be slightly more powerful than τ 2 2 , but as n increases, no significant differences between tests can be seen. As p = q = 5 and n = 30, ρ 2 2 is clearly more powerful than τ 2 2 , but again, for large n, the tests are almost equally powerful. The size of ρ 2 2 is very close to the designated size of .05 in all cases. For small sample sizes, the size of τ 2 2 is often slightly below .05, and for heavy-tailed distributions, W varies widely above .05.
Because interdirection-based τ 2 1 and ρ 2 1 are computationally very intensive, we compared the empirical powers of τ 2 1 and ρ 2 1 with those of τ 2 2 and ρ 2 2 for n = 15 and p = q = 2 only. As shown in Figure 3, both Spearman's rho statistics are clearly better than the two Kendall's tau statistics, as was expected. However, no significant differences can be seen between τ 2 1 and τ 2 2 or between ρ 2 1 and ρ 2 2 .

A REAL DATA EXAMPLE
McNaughton and Davies (1987) studied the effect of an aerobic conditioning program on cholesterol levels. Several physical and blood-related variables of 12 subjects were measured before and after an aerobic conditioning program. Table 3 lists changes in the respiratory variables vital capacity (VC), forced expiratory volume (FEV), and maximum oxygen uptake (VO 2 ) and in the cholesterol variables total cholesterol (TC), triglycerides (TG), and high-density lipoprotein cholesterol (HDL). We wish to test whether the respiratory changes are independent of cholesterol changes.
The p values given by different tests were calculated using the chi-squared approximations to the null distributions. For the original dataset, the p values for W, τ 2 2 , ρ 2 2 , and Q 2 2 are .022, .028, .018, and .025. To illustrate the robustness of τ 2 2 , ρ 2 2 , and Q 2 2 , the first VO 2 measurement was moved from −25 to 50 and the p values were recalculated. The results are shown in Figure 4. When the first VO 2 measurement decreases, the p value given by W increases substantially, whereas only small changes can be seen in the p values of other tests.
Finally, the permutation p values were obtained for the original dataset using 1,500 permutation samples. The p values for W, τ 2 2 , ρ 2 2 , and Q 2 2 are .069, .011, .015, and .013. Note that even for this very small sample size (n = 12), the permutation p values of τ 2 2 , ρ 2 2 , and Q 2 2 are quite similar to those given by the chi-squared approximation, and the approximation for ρ 2 2 is extremely good. More variation is seen in the p value of W.

FINAL REMARKS
In this article we have developed new nonparametric, intuitively appealing tests for testing whether two random vectors are independent. The tests generalize the popular univariate tests due to Kendall (1938) and Spearman (1904). In the first approach, the interdirections of Randles (1989) are used in test construction. It is remarkable that in this approach, no explicit location vector and scatter matrix estimates are needed. In the second approach, the extensions based on the spatial signs and ranks of the standardized observations are introduced. Natural shape matrix estimates for the standardization are proposed. The test statistics arising from these two approaches appear to be asymptotically equivalent if each vector is elliptically symmetric.
The new tests have an appealing invariance property and a convenient limiting chi-squared null distribution. The tests also have excellent asymptotic and finite-sample efficiencies. The spatial sign and rank versions are easy to compute for data in any dimension. As demonstrated in the example, the tests are resistant to outliers, and the p values for the conditionally distribution-free permutation test are also easily generated. To conclude, the spatial sign and rank versions provide intuitively appealing, practical, robust alternatives to multivariate normaltheory methods.

Proof of Lemma 1
Using V * = k[Tr(A VA T )] −1 A VA T and µ * = A µ + b, it is straightforward to show that The transformation matrix P is clearly orthogonal. Similarly, S * ij = P S ij and R * i = P R i .

Proof of Theorem 1
To shorten the notation in what follows, we write x j . Due to the affine invariance of the tests, we consider only the case with spherical marginal densities. First, Gieser and Randles (1997) showed that Q 2 1 and 1 n i are asymptotically equivalent. The same asymptotic representation is found for Q 2 2 in the proof of theorem 1 of Taskinen et al. (2003a). Second, under the null model, the interdirection proportions p (1) (i, j; i , j ) are U-statistics with bounded kernels and expected values and similarly for p (2) (i, j; i , j ). Thus As was done by Gieser and Randles (1997), we can then show that the multivariate Kendall's τ 2 1 is asymptotically equivalent to and similarly for Spearman's ρ 2 1 . The proof is then completed by finding the same asymptotic representations for τ 2 2 and ρ 2 2 ; see the proof of Theorem 2.

Proof of Theorem 2
The Limiting Distribution of τ 2 2 . As in the proof of theorem 1 of Taskinen et al. (2003a), it can be shown that under H 0 ,  (2)T ij which completes the proof.
The Limiting Distribution of ρ 2 * 2 . Proceeding as before, we get Finally, note that 2E G (q G (r)r) = E( x ij ) and 2E F (q F (r)r) = E( x (1) ij ). The limiting distribution of Q 2 * 2 was derived by Taskinen et al. (2003a).

APPENDIX B: EXPRESSIONS FOR c 2 F AT CERTAIN MULTIVARIATE DISTRIBUTIONS
The constants c 2 F for standard multivariate normal distribution F = , multivariate t distribution F = t ν and multivariate contaminated normal distribution F = ,c are and c 2 ,c = 1 2π where Expressions for the multivariate normal and t distribution cases were derived by Möttönen, Oja, and Tienari (1997). The expression for the contaminated normal case is found in a similar way. Constants for selected distributions and dimensions are listed in Table B.1. [Received February 2003. Revised December 2004