Empirical Likelihood Method for Complete Independence Test on High Dimensional Data

Statistical inference on high dimensional data has gained a wide range of applications in recent years. New techniques generate a vast collection of data sets with high dimensions, for example, trading data from financial market, social network data and biological data like microarray and DNA data. The dimension of these types of data is not small compared with sample size, and typically of the same order as sample size or even larger. Yet classical multivariate statistics usually deal with data from normal distributions with a large sample size $n$ and a fixed dimension $p$, and one can easily find some classic treatments in reference books such as Anderson [1], Morrison [10] and Muirhead [11].

Under multivariate normality settings, the likelihood ratio test statistic converges in distribution to a chi-squared distribution when $p$ is fixed. However, when $p$ changes with $n$ and tends to infinity, this conclusion is no longer true as discovered in Bai et al. [2], Jiang et al. [6], Jiang and Yang [8], Jiang and Qi [7], Qi et al. [16], among others. The results in these papers indicate that the chi-square approximation fails when $p$ diverges as $n$ goes to infinity.

The test of complete independence of a random vector is to test whether all the components of the random vector are mutually independent. In the multivariate normal case, the test of complete independence is equivalent to the test whether covariance matrix is a diagonal matrix, or whether the correlation matrix is the identity matrix.

For more details, we assume $X=(X_{1},\cdots,X_{p})$ is a random vector from a $p$-dimensional multivariate normal distribution $N_{p}(\bm{\mu},\mathbf{\Sigma})$, where $\bm{\mu}$ denotes the mean vector, and $\mathbf{\Sigma}$ is a $p\times p$ covariance matrix. Given a random sample of size $n$ from the normal distribution, $\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n}$, where $\mathbf{x}_{i}=(x_{i1},x_{i2},\cdots,x_{ip})$ for $1\leq i\leq n$, Pearson’s correlation coefficient between the $i$-th and $j$-th components is given by

where $\bar{x}_{i}=\frac{1}{n}\sum\limits_{k=1}^{n}x_{ki}$ and $\bar{x}_{j}=\frac{1}{n}\sum\limits_{k=1}^{n}x_{kj}$ for $1\leq i,j\leq p$. Now we set $\mathbf{R}_{n}=(r_{ij})_{p\times p}$ as the sample correlation coefficient matrix.

The complete independence test for the normal random vector is

where $\bm{\Gamma}$ is the population correlation matrix and $\mathbf{I}_{p}$ is $p\times p$ identity matrix. When $p<n$, the likelihood ratio test statistic for (2) is a function of $|\mathbf{R}_{n}|$, the determinant of $\mathbf{R}_{n}$, from Bartlett [3] or Morrison [10]. In traditional multivariate analysis, when $p$ is a fixed integer, we have under the null hypothesis in (2) that

where $\chi^{2}_{f}$ denotes a chi-square distribution with $f$ degrees of freedom.

When $p=p_{n}$ depends on $n$ with $2\leq p_{n}<n$ and $p_{n}\to\infty$, the likelihood ratio method can still be applied to test (2). The limiting distributions of the likelihood ratio test statistics in this case have been discussed in the aforementioned papers. It is worth mentioning that Qi et al. [16] propose an adjusted likelihood ratio test statistic and show that the distribution of the adjusted likelihood test statistic can be well approximated by a chi-squared distribution whose number of degrees of freedom depends on $p$ regardless of whether $p$ is fixed or divergent.

The limitation of the likelihood ratio test is that the dimension $p$ of the data must be smaller than the sample size $n$. Many other likelihood tests related to the sample covariance matrix or sample correlation matrix have the same problem as the sample covariance matrices are degenerate when $p\geq n$. In order to relax this constraint, a new test statistic using the sum of squared sample correction coefficients is proposed by Schott [17] as follows

Assume that the null hypothesis of (2) holds. Under assumption $\lim\limits_{n\to\infty}p_{n}/n=\gamma\in(0,\infty)$, Schott [17] proves that $t_{np}-\frac{p(p-1)}{2(n-1)}$ converges in distribution to a normal distribution with mean $0$ and variance $\gamma^{2}$, that is,

where $\sigma_{np}^{2}=\frac{p(p-1)(n-2)}{(n-1)^{2}(n+1)}$.

Recently, Mao [9] proposes a different test for complete independence. His test statistic is closely related to Schott’s test and is defined by

It has been proved in Mao [9] that $T_{np}$ is asymptotically normal under the null hypothesis of (2) and the assumption that $\lim\limits_{n\to\infty}p_{n}/n=\gamma\in(0,\infty)$.

Very recently, Chang and Qi [4] investigate the limiting distributions for the two test statistics above under less restrictive conditions on $n$ and $p$. Chang and Qi [4] show that (3) is also valid under the general condition that $p_{n}\to\infty$ as $n\to\infty$, regardless of the convergence rate of $p_{n}$. Thus, the normal approximation in (3) based on $t_{np}^{*}$ yields an approximate level $\alpha$ rejection region

where $z_{\alpha}$ is a $\alpha$ level critical value of the standard normal distribution.

Furthermore, Chang and Qi [4] propose adjusted test statistics whose distribution can be fitted by chi-squared distribution regardless of how $p$ changes with $n$ as long as $n$ is large. Chang and Qi’s [4] adjusted test statistics $t_{np}^{c}$ is defined as

Chang and Qi show that

as long as $p_{n}\to\infty$ as $n\to\infty$. Let $\chi^{2}_{f}(\alpha)$ denote the $\alpha$ level critical value of $\chi^{2}_{f}$. Then an approximate level $\alpha$ rejection region based on $t_{np}^{c}$ is given by

One can find more references on test for complete independence in Mao [9] or Chang and Qi [4].

In practice, the assumption of normality for distributions may be violated. Now we assume $X=(X_{1},\cdots,X_{p})$ is a random vector and $X_{1},\cdots,X_{p}$ are identically distributed with distribution function $F$. Given a random sample of size $n$, $\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n}$, where $\mathbf{x}_{i}=(x_{i1},x_{i2},\cdots,x_{ip})$ for $1\leq i\leq n$, are drawn from the distribution of $X=(X_{1},\cdots,X_{p})$, and define Pearson’s correlation coefficients $r_{ij}$’s as in (1). By using the Stein method, Chen and Shao [5] show that (3) holds under some moment conditions of $F$ if $p_{n}/n$ is bounded.

In this paper, we propose to apply empirical likelihood method to the testing problem (2). The empirical likelihood is a nonparametric statistical method proposed by Owen [12, 13], which is originally used to test the mean vector of a population based on a set of independent and identically distributed (i.i.d.) random variables. Empirical likelihood does not require to specify the family of distributions for the data and it possesses some good properties of the likelihood methods.

The rest of the paper is organized as follows. In Section 2, we first introduce a one-sided empirical likelihood method for the mean of a set of random variables with a common mean and then establish the connection between the test of complete independence and the one-sided empirical likelihood method. Our main result concerning the limiting distribution of the one-sided empirical likelihood ratio statistic is also given in Section 2. In Section 3, we carry out a simulation study to compare the performance of the empirical likelihood method and normal approximation based on Schott’s test statistic and chi-square approximation based on Chang and Qi’s adjusted test statistic. In our simulation study, we also apply these methods to some other distributions such as the exponential distributions and mixture of the exponential and normal distributions so as to compare their adaptability to non-normality. The proofs of the main results are given in Section 4.

In this section, we apply the empirical likelihood method to the test of complete independence. First, we assume $X=(X_{1},\cdots,X_{p})$ is a random vector from a $p$-dimensional multivariate normal distribution. Under the null hypothesis of (2), $\{r_{ij}^{2},1\leq i<j\leq p\}$ are random variables from an identical distribution with mean $\frac{1}{n-1}$. As a matter of fact, it follows from Corollary 5.1.2 in Muirhead [11] that $r_{ij}^{2}$ has the same distribution as $T^{2}/(n-2+T^{2})$ under the null hypothesis of (2), where $T$ is a random variable having $t$-distribution with $n-2$ degrees of freedom. $\{r_{ij}^{2},1\leq i<j\leq p\}$ are asymptotically independent if the sample size $n$ is large. We will develop a one-sided empirical likelihood test statistic and apply it to the data set $\{(n-1)r_{ij}^{2},1\leq i<j\leq p\}$, where $p=p_{n}$ is a sequence of positive integers such that $p_{n}\to\infty$ as $n\to\infty$. As an extension, we then consider the case when $X=(X_{1},\cdots,X_{p})$ is a random vector with an identical marginal distribution function $F$ which is not necessarily Gaussian. When the $p$ components of $X$ are independent, we demonstrate that the empirical likelihood method we develop under normality works for general distribution $F$ as well if some additional conditions are satisfied.

In this section, we will consider the following three test statistics for testing complete independence (2), including Schott’s test statistic $t_{np}^{*}$ given in (3), Chang and Qi’s adjusted test statistic $t_{np}^{c}$ defined in (5), and the rescaled empirical likelihood test statistic $\bar{\ell}_{n}$ given in (15). The corresponding rejection regions are given in (4), (6), and (18), respectively. All simulations are implemented by the software R.

For sample size $n=20,50,100$ and dimension $p=10,20,50,100$, we apply the three test statistics to each of five distributions for 10000 iterations to obtain the empirical sizes and the empirical powers of the tests. We set the nominal type I error $\alpha=0.05$. The five distributions include the normal, the uniform over $[-1,1]$, the exponential, the mixture of the normal and exponential distributions, and the sum of normal and exponential distributions.

To control the dependence structure, we introduce a covariance matrix $\bm{\Gamma}_{\rho}$ defined by

which is also a correlation matrix. In our simulation study, we generate random samples from the distribution of a random vector $X=(X_{1},\cdots,X_{p})$ with covariance matrix $\bm{\Gamma}_{\rho}$ or correlation matrix $\bm{\Gamma}_{\rho}$. For details, see the five distributions described below. For all distributions we consider, the observations have independent components when $\rho=0$ and positively dependent components when $\rho>0$. We choose very small values for $\rho$ such as $\rho=0.02$ and $0.05$. When the value of $\rho$ is large, the resulting powers for all three methods will be too close to $1$, and the comparison is meaningless. Therefore, based on $10,000$ replicates, the sizes for three test statistics are estimated when $\rho=0$, and their powers are estimated when $\rho=0.02$ and $0.05$. All results are reported in Tables 1 to 5.

The observations are drawn from a multivariate normal random vector $X=(X_{1},\cdots,X_{p})$ with mean $\bm{\mu}=(0,\cdots,0)$ and variance matrix $\bm{\Gamma}_{\rho}$ specified in (19). The results on the empirical sizes and powers are given in Table 1.

We first generate $p+1$ i.i.d. random variables $Y_{0},Y_{1},\cdots,Y_{p}$ from Uniform $(-1,1)$ distribution, then set $X_{i}=\frac{\sqrt{\rho}}{\sqrt{1-\rho}}Y_{0}+Y_{i}$, $i=1,2,\cdots,p$. It is easy to verify that random vector $X=(X_{1},\cdots,X_{p})$ has mean $\bm{\mu}=(0,\cdots,0)$ and correlation matrix $\bm{\Gamma}_{\rho}$ as defined in (19). The results on the empirical sizes and powers are given in Table 2.

We generate $p+1$ i.i.d. random variables $Y_{0},Y_{1},\cdots,Y_{p}$ from the unit exponential distribution, then define $X_{i}=\frac{\sqrt{\rho}}{\sqrt{1-\rho}}Y_{0}+Y_{i}$, $i=1,2,\cdots,p$. The random vector $X=(X_{1},\cdots,X_{p})$ has a correlation matrix $\bm{\Gamma}_{\rho}$ as defined in (19) for $\rho\in[0,1)$. The results on the empirical sizes and powers are given in Table 3.

The random vector $X=(X_{1},\cdots,X_{p})$ is sampled from a mixture of the normal and exponential distributions which is with 90% probability from the multivariate normal with mean $\bm{\mu}=(1,\cdots,1)$ and covariance matrix $\bm{\Gamma}_{\rho}$ given in (19) and with 10% probability from a random vector $(Y_{1},\cdots,Y_{p})$ where $Y_{1},\cdots,Y_{p}$ are i.i.d. unit exponential random variables. The results on the empirical sizes and powers are given in Table 4.

The random vector $X=(X_{1},\cdots,X_{p})$ is a weighted sum of two independent random vectors, $U$ and $V$, $X=U+0.01V$, where $U$ is from a multivariate normal distribution with mean $\bm{\mu}=(0,\cdots,0)$ and covariance matrix $\bm{\Gamma}_{\rho}$ defined in (19), and $V=(Y_{1},\cdots,Y_{p})$ with $Y_{i}$’s being i.i.d. unit exponential random variables. The results on the empirical sizes and powers are given in Table 5.

From the simulation results, the empirical sizes for all three tests are close to $0.05$ which is the nominal type I error we set in the simulation, especially when both $n$ and $p$ are large. Test statistic $t_{np}^{c}$ has the smallest size in most cases, and it is a little bit conservative sometimes. The size of $\bar{\ell}_{n}$ is between that of $t_{np}^{c}$ and $t_{np}^{*}$ and both $\bar{\ell}_{n}$ and $t_{np}^{*}$ are comparable for most combinations of $n$ and $p$.

As we expect, the powers of all three test statistics become higher as $p$ grows larger. The increase in $n$ also brings about an increase in power, but not as much as the increase in $p$ does, because $\frac{p(p-1)}{2}$ is the number of $r_{ij}^{2}$’s involved in the test. All test statistics achieve high power when $\rho=0.05$. Three test statistics result in comparable powers in general, although the power of Chang and Qi’s test statistic is occasionally a little bit less than the other two test statistics. These differences may be due to the fact that Chang and Qi’s test statistic maintain a lower type I error.

In summary, in this paper, we have developed the one-sided empirical likelihood method and proposed the rescaled empirical likelihood test statistic for testing the complete independence for high dimensional random vectors. The rescaled empirical likelihood test statistic performs very well in terms of the size and power and can serve as a good alternative to the existent test statistics in the literature.

For every $t\in(0,1)$, set $\omega_{ti}=t\omega_{1i}+(1-t)\omega_{2i}$, $i=1,\cdots,N$. Then $\omega_{ti}>0$, $\sum\limits^{N}_{i=1}\omega_{ti}=1$, $\sum\limits^{N}_{i=1}\omega_{ti}y_{i}=t\mu_{1}+(1-t)\mu_{2}\in H$. Since $\log x$ is strictly concave in $(0,\infty)$, we have

and at least one of the inequalities is strict, i.e, $\log(\omega_{ti})>t\log\omega_{1i}+(1-t)\log\omega_{2i}$ for some $i$, since $\mu_{1}\neq\mu_{2}$ implies $(\omega_{11},\omega_{12},\cdots,\omega_{1N})\neq(\omega_{21},\omega_{22},\cdots,\omega_{2N})$. Therefore, we get

which implies

proving (20).

When $\mu=\bar{y}$, an obvious solution to (9) is $\lambda=0$. Since the solution to (9) is unique, we see that $\omega_{i}=N^{-1}$, and thus, $L(\bar{y})=N^{-N}$. We also notice that

The last step is obtained by using the Lagrange multipliers. We omit the details here. Therefore, we conclude that $L(\bar{y})=\sup\limits_{\mu}L(\mu)=N^{-N}$. $\Box$

Define $\sigma_{n}^{2}=\frac{2(n-2)}{n+1}$ and $S_{n}^{2}=\frac{1}{N}\sum\limits_{1\leq i<j\leq p}\big{(}(n-1)r_{ij}^{2}-1\big{)}^{2}$. Review that $N=p(p-1)/2$. We have $\sigma_{np}^{2}=N\sigma_{n}^{2}/(n-1)^{2}$. Since the distribution of $y_{j}$’s depends on $n$, $\{y_{j},\,1\leq j\leq N\}$ forms an array of random variables.

If the following three conditions are satisfied: (i). $\frac{1}{\sigma_{n}}\max\limits_{1\leq j\leq N}|y_{j}-1|=o_{p}(N^{1/2})$ as $n\to\infty$; (ii). $\frac{1}{N\sigma_{n}^{2}}\sum\limits^{N}_{j=1}(y_{j}-1)^{2}\overset{p}{\to}1$ as $n\to\infty$; (iii). $\displaystyle{\frac{\sum\limits^{N}_{j=1}y_{j}-N}{\sqrt{N\sigma_{n}^{2}}}}\overset{d}{\to}N(0,1)$ as $n\to\infty$, equivalently, in term of $r_{ij}^{2}$’s,

1. (C1).

   $\displaystyle\frac{1}{\sigma_{n}}\max_{1\leq i<j\leq p}|(n-1)r_{ij}^{2}-1|=o_{p}(N^{1/2})$ as $n\to\infty$;

2. (C2).

   $\displaystyle\frac{1}{\sigma_{n}^{2}}S_{n}^{2}\overset{p}{\to}1$ as $n\to\infty$;

3. (C3).

   $\displaystyle z_{n}:=\frac{\sum\limits_{1\leq i<j\leq p}(n-1)r_{ij}^{2}-N}{\sqrt{N\sigma_{n}^{2}}}\overset{d}{\to}N(0,1)$ as $n\to\infty$,

we can follow the same procedure as in Owen [14] or use Theorem 6.1 in Peng and Schick [15] to conclude that

where $\ell(1)$ is defined in (10) with $\mu_{0}=1$. Again, by using condition (C3), we have

as $n\to\infty$, where $\ell_{n}$ is defined in (13), proving Theorem 1.

Now we will verify conditions (C1), (C2) and (C3). (C3) has been proved by Chang and Qi [4] as we indicate below equation (3).

Assume $(i,j)$ is a pair of integers with for $1\leq i<j\leq p$. It is proved in Schott [17] that