Test of Covariance and Correlation Matrices

TL;DR

This paper introduces a flexible, distribution-free framework for testing covariance and correlation matrices using a generalized cosine measure, applicable to various one-sample and two-sample scenarios with real data and simulations.

Contribution

It develops a novel, distribution-agnostic testing framework for covariance and correlation matrices, including permutation algorithms for common tests, accommodating diverse data distributions and dimensions.

Findings

Performs well in controlling type I error across tests.

Shows high power in detecting covariance and correlation differences.

Applicable to data with different marginal distributions.

Abstract

Based on a generalized cosine measure between two symmetric matrices, we propose a general framework for one-sample and two-sample tests of covariance and correlation matrices. We also develop a set of associated permutation algorithms for some common one-sample tests, such as the tests of sphericity, identity and compound symmetry, and the $K$-sample tests of multivariate equality of covariance or correlation matrices. The proposed method is very flexible in the sense that it does not assume any underlying distributions and data generation models. Moreover, it allows data to have different marginal distributions in both the one-sample identity and $K$-sample tests. Through real datasets and extensive simulations, we demonstrate that the proposed method performs well in terms of empirical type I error and power in a variety of hypothesis testing situations in which data of different sizes and dimensions are generated using different distributions and generation models.

Tables

Table 1. Table 1: Empirical type I error (%) of the sphericity test at the 5% significance level

Sample Size ($n$)	20	40	60	80	20	40	60	80
Dimension ($p$)	$N(0,1)$				Gamma(4, 0.5)
38	4.9	5.2	5.1	5.2	6.2	5.3	4.6	4.7
55	4.8	5.6	5.0	4.5	4.7	5.6	4.2	4.7
89	5.2	5.1	4.5	5.0	5.0	5.2	4.5	4.9
159	4.8	5.2	5.3	4.2	6.0	5.3	4.5	5.1
181	5.2	5.6	5.0	4.4	4.9	4.4	4.7	5.9
331	4.4	4.6	5.0	4.1	4.9	5.5	5.1	4.7
343	5.3	3.4	5.0	5.4	5.8	4.2	5.8	4.2
642	4.8	4.9	5.4	4.8	5.5	4.4	5.6	4.9

Table 2. Table 2: Empirical power (%) of the sphericity test of H 0 : 𝚺 = σ 2 ​ 𝑰 : subscript 𝐻 0 𝚺 superscript 𝜎 2 𝑰 H_{0}:\bm{\Sigma}=\sigma^{2}\bm{I} vs H 1 : 𝚺 = σ 0 2 ​ 𝑰 + σ 1 2 ​ 𝑨 : subscript 𝐻 1 𝚺 superscript subscript 𝜎 0 2 𝑰 superscript subscript 𝜎 1 2 𝑨 H_{1}:\bm{\Sigma}=\sigma_{0}^{2}\bm{I}+\sigma_{1}^{2}\bm{A} at the 5% significance level with 𝑨 = diag ​ ( 𝑰 [ λ ​ p ] , 𝟎 p − [ λ ​ p ] ) 𝑨 diag subscript 𝑰 delimited-[] 𝜆 𝑝 subscript 0 𝑝 delimited-[] 𝜆 𝑝 \bm{A}=\mathrm{diag}(\bm{I}_{\left[\lambda p\right]},\bm{0}_{p-\left[\lambda p\right]}) , σ 0 2 = σ 1 2 = 1 superscript subscript 𝜎 0 2 superscript subscript 𝜎 1 2 1 \sigma_{0}^{2}=\sigma_{1}^{2}=1 , and λ = 0.125 𝜆 0.125 \lambda=0.125 .

Sample Size ($n$)	20	40	60	80	20	40	60	80
Dimension ($p$)	$N(0,1)$				Gamma(4, 0.5)
38	32	66	91	98	17	48	73	90
55	32	73	94	99	16	52	83	95
89	36	83	98	100	17	64	92	98
159	37	85	99	100	21	67	95	100
181	37	86	99	100	20	68	96	100
331	38	90	100	100	20	74	97	100
343	38	88	99	100	21	73	97	100
642	38	89	100	100	22	73	98	100

Table 3. Table 3: Empirical power (%) of the sphericity test H 0 : 𝚺 = σ 2 ​ 𝑰 : subscript 𝐻 0 𝚺 superscript 𝜎 2 𝑰 H_{0}:\bm{\Sigma}=\sigma^{2}\bm{I} vs H 1 : 𝚺 = ( 1 − ρ ) ​ σ 0 2 ​ 𝑰 + ρ ​ σ 1 2 ​ 𝑱 : subscript 𝐻 1 𝚺 1 𝜌 superscript subscript 𝜎 0 2 𝑰 𝜌 superscript subscript 𝜎 1 2 𝑱 H_{1}:\bm{\Sigma}=(1-\rho)\sigma_{0}^{2}\bm{I}+\rho\sigma_{1}^{2}\bm{J} at the 5% significance level with ρ = 0.10 𝜌 0.10 \rho=0.10 , σ 0 2 = 1 superscript subscript 𝜎 0 2 1 \sigma_{0}^{2}=1 and σ 1 2 = 2 superscript subscript 𝜎 1 2 2 \sigma_{1}^{2}=2

Sample Size ($n$)	20	40	60	80	20	40	60	80
Dimension ($p$)	$N(0,1)$				Gamma(4, 0.5)
38	97	100	100	100	93	100	100	100
55	99	100	100	100	97	100	100	100
89	100	100	100	100	99	100	100	100
159	100	100	100	100	100	100	100	100
181	100	100	100	100	100	100	100	100
331	100	100	100	100	100	100	100	100
343	100	100	100	100	100	100	100	100
642	100	100	100	100	100	100	100	100

Table 4. Table 4: Empirical type I error (%) and power (%) of the identity test using block-diagonal model for H 0 : 𝚺 = 𝑰 : subscript 𝐻 0 𝚺 𝑰 H_{0}:\bm{\Sigma}=\bm{I} and H 1 : 𝚺 = 𝑰 4 ⊗ 𝚺 q : subscript 𝐻 1 𝚺 tensor-product subscript 𝑰 4 subscript 𝚺 𝑞 H_{1}:\bm{\Sigma}=\bm{I}_{4}\otimes\bm{\Sigma}_{q} , where 𝚺 q = ( 1 − ρ ) ​ 𝑰 + ρ ​ 𝑱 subscript 𝚺 𝑞 1 𝜌 𝑰 𝜌 𝑱 \bm{\Sigma}_{q}=(1-\rho)\bm{I}+\rho\bm{J} with ρ = 0.15 𝜌 0.15 \rho=0.15

Sample Size ($n$)	$n=100$					$n=4$
Dimension ($p$)	24	32	64	76	92	100	200	300	500	700	1000
	Type I error					Type I error
$N(0,1)$	5.5	5.0	5.5	5.0	4.8	5.3	5.5	4.5	5.0	4.1	4.9
Log-normal(0,1)	5.0	4.2	4.4	4.9	5.4	5.8	4.9	4.4	5.2	5.1	5.0
Student $t_5$	4.2	5.3	4.4	4.8	5.1	4.8	4.9	4.7	5.0	5.3	5.5
Gumbel(10,2)	5.0	4.5	5.2	5.0	4.6	4.9	5.2	4.6	5.2	5.2	5.1
Hybrid	4.6	4.9	4.2	4.4	5.4	4.9	4.5	5.1	5.2	5.1	4.6
	Power					Power
$N(0,1)$	99	100	100	100	100	13	24	31	47	57	70
Log-normal(0,1)	97	100	100	100	100	39	69	82	93	97	99
Student $t_5$	99	100	100	100	100	18	32	43	58	68	80
Gumbel(10,2)	99	100	100	100	100	14	28	40	59	68	79
Hybrid	99	100	100	100	100	22	43	57	74	84	90

Table 5. Table 5: Empirical type 1 error and power of two-sample test H 0 : 𝚺 1 = 𝚺 2 : subscript 𝐻 0 subscript 𝚺 1 subscript 𝚺 2 H_{0}:\bm{\Sigma}_{1}=\bm{\Sigma}_{2} with ρ 1 = ρ 2 = 0.15 subscript 𝜌 1 subscript 𝜌 2 0.15 \rho_{1}=\rho_{2}=0.15 vs. 𝚺 1 ≠ 𝚺 2 subscript 𝚺 1 subscript 𝚺 2 \bm{\Sigma}_{1}\neq\bm{\Sigma}_{2} with ρ 1 = 0.15 subscript 𝜌 1 0.15 \rho_{1}=0.15 and ρ 2 = 0.30 subscript 𝜌 2 0.30 \rho_{2}=0.30

Sample Size ($n$)	$n_1=100$,  $n_2=100$					$n_1=20$, $n_2=20$
Dimension ($p$)	24	32	64	76	92	100	200	300	500	700	1000
	Type I error					Type I error
$N(0,1)$	5.5	4.3	4.2	4.8	4.5	4.6	5.5	4.6	5.1	4.8	4.8
Log-normal(0,1)	4.9	5.0	4.8	5.5	5.3	5.8	4.9	4.4	5.2	5.1	5.0
Student $t_5$	5.2	4.4	4.3	4.3	5.3	5.0	5.0	4.4	4.4	5.3	3.9
Gumbel(10,2)	5.0	4.7	5.4	5.2	4.7	5.0	4.8	4.8	5.0	4.8	5.1
Hybrid	4.8	4.4	4.5	5.4	5.6	4.5	5.4	5.0	4.3	5.4	5.0
	Power					Power
$N(0,1)$	46	65	90	92	95	23	30	30	33	34	35
Log-normal(0,1)	9.5	13	30	42	50	24	37	45	47	49	49
Student $t_5$	29	43	81	85	90	18	28	27	30	33	31
Gumbel(10,2)	90	100	100	100	100	100	100	100	100	100	100
Hybrid	37	65	99	99	100	100	100	100	100	100	100

Table 6. Table 6: Empirical type I error (%) and power (%)of the proposed two-sample test

Dimension $p$		50	100	200	300	600	800
Sample 1	Sample 2	Type I error
Gamma(4, 0.5)	Gamma(0.5, $\sqrt{2}$)	100	100	100	100	100	100
Gamma(0.5, $\sqrt{2}$)	Gamma(0.5, $\sqrt{2}$)	5.8	5.7	4.3	5.6	5.0	5.0
Gamma(4, 0.5)	Gamma(4, 0.5)	5.3	5.7	6.0	5.2	5.3	5.5
Sample 1	Sample 2	Power
Gamma(4, 0.5)	Gamma(0.5, $\sqrt{2}$)	100	100	100	100	100	100
Gamma(0.5, $\sqrt{2}$)	Gamma(0.5, $\sqrt{2}$)	13	20	37	51	87	97
Gamma(4, 0.5)	Gamma(4, 0.5)	99	100	100	100	100	100

Table 7. Table 7: Impact of sparse signal and sample signal-to-noise ratio (SSNR) on the proposed two-sample test

Dimension ($p$)		50	100	200	400	800		50	100	200	400	800
	SSNR	Type I error (%)					SSNR	power (%)
$N(0,1)$	0.0 (0.0)*	5.5	5.3	4.4	5.8	5.0	0.0 (0.0)	33	17	10	8	7
$N(2,1)$	2.4 (0.1)	5.5	4.7	4.7	5.0	5.3	2.5 (0.1)	48	27	17	13	9
$N(4,1)$	4.9 (0.1)	5.9	4.7	4.7	5.0	5.3	5.0 (0.1)	99	96	89	80	76
Gamma(5,1)	2.7 (0.1)	4.6	5.3	5.3	5.0	5.3	2.8 (0.1)	55	34	21	15	11
Gamma(10,1)	3.9 (0.1)	4.8	5.3	5.4	5.6	5.2	4.0 (0.1)	89	74	57	42	36
Poisson(5)	2.7 (0.1)	5.1	5.3	4.9	5.0	5.3	2.8 (0.1)	53	35	22	14	10
Poisson(10)	3.8 (0.1)	4.7	5.1	5.2	5.4	4.3	3.9 (0.1)	88	76	57	42	34
Log-normal(0, 0.4)	3.0 (0.1)	4.9	5.0	4.3	5.5	4.8	3.0 (0.1)	58	38	25	17	13
Log-normal(0, 0.3)	4.0 (0.1)	5.1	5.4	4.1	5.8	4.9	4.1 (0.1)	91	77	61	47	39
*Standard deviation of SSNR

Equations

$$cos({\mathbf{M}}_1,{\mathbf{M}}_2)=\frac{⟨f({\mathbf{M}}_1),f({\mathbf{M}}_2)⟩}{\Vert f({\mathbf{M}}_1)\Vert \Vert f({\mathbf{M}}_2)\Vert },$$

$$\left[\begin{array}{ccccc}1.00& 0.50& 0.33& 0.25& 0.20\\ 0.50& 1.00& 0.25& 0.20& 0.17\\ 0.33& 0.25& 1.00& 0.17& 0.14\\ 0.25& 0.20& 0.17& 1.00& 0.12\\ 0.20& 0.17& 0.14& 0.12& 1.00\end{array}\right],$$

$$\mathbf{A}\left(i,j\right)=\{\begin{array}{ll}\frac{1}{i+j-1},& i\ne j,\\ 1,& i=j,\end{array}$$

$$\left[\begin{array}{ccccc}1.00& 0.74& 0.83& 0.54& 0.41\\ 0.74& 1.00& 0.55& 0.60& 0.34\\ 0.83& 0.55& 1.00& 0.28& 0.58\\ 0.54& 0.60& 0.28& 1.00& 0.48\\ 0.41& 0.34& 0.58& 0.48& 1.00\end{array}\right],$$

$${\mathbf{A}}^{\ast }\left(i,j\right)={\mathbf{A}}^{\ast }\left(j,i\right)=\{\begin{array}{ll}\frac{1}{i+j-1}+ϵ,& i\ne j,\\ 1,& i=j,\end{array}$$

$$cos(\mathbf{A},\mathbf{B})=\frac{\mathrm{tr}({\mathbf{A}}^T\mathbf{B})}{\mathrm{tr}({\mathbf{A}}^T\mathbf{A}{)}^{\frac{1}{2}}\mathrm{tr}({\mathbf{B}}^T\mathbf{B}{)}^{\frac{1}{2}}},$$

$$cos(\mathbf{A},\mathbf{B})=\frac{\mathrm{vech}({\mathbf{L}}_1{)}^T\mathrm{vech}({\mathbf{L}}_2)}{\Vert \mathrm{vech}({\mathbf{L}}_1)\Vert \Vert \mathrm{vech}({\mathbf{L}}_2)\Vert },$$

$$cos(\mathbf{A},\mathbf{B})=\frac{v_1^T{v}_2}{\Vert v_1\Vert \Vert v_2\Vert },$$

$$cos(\mathbf{A},\mathbf{B})=\frac{\mathrm{vech}(\mathbf{A}{)}^T\mathrm{vech}(\mathbf{A})}{\Vert \mathrm{vech}(\mathbf{A})\Vert \Vert \mathrm{vech}(\mathbf{B})\Vert },$$

$$cos(\mathbf{A},\mathbf{B})=\frac{{\mathrm{vech}}^{\ast }(\mathbf{A}{)}^T{\mathrm{vech}}^{\ast }(\mathbf{A})}{\Vert {\mathrm{vech}}^{\ast }(\mathbf{A})\Vert \Vert {\mathrm{vech}}^{\ast }(\mathbf{B})\Vert },$$

$$1-cos(\mathbf{S},{\mathbf{\Sigma }}_0),$$

$$1-cos({\mathbf{S}}_1,{\mathbf{S}}_2),$$

$$H_0:\mathbf{\Sigma }={\sigma }^2\mathbf{I}vs.H_1:\mathbf{\Sigma }\ne {\sigma }^2\mathbf{I},$$

$$1-\frac{\mathrm{tr}(\mathbf{S})}{\sqrt{p}\Vert \mathrm{vech}(\mathbf{S})\Vert },$$

$$H_0:\mathbf{\Sigma }=\mathbf{I}vs.H_1:\mathbf{\Sigma }\ne \mathbf{I},$$

$$H_0:\mathbf{R}=\mathbf{I}vs.H_1:\mathbf{R}\ne \mathbf{I},$$

$$1-\frac{\sqrt{p}}{\Vert \mathrm{vech}(\mathbf{S})\Vert }.$$

$$1-\frac{{\mathrm{vech}}^{\ast }(\mathbf{S}{)}^T{\mathrm{vech}}^{\ast }(\mathbf{J})}{\sqrt{0.5p(p-1)}\Vert {\mathrm{vech}}^{\ast }(\mathbf{S})\Vert },$$

$$1-\frac{\mathrm{vech}({\mathbf{S}}_1{)}^T\mathrm{vech}({\mathbf{S}}_2)}{\Vert \mathrm{vech}({\mathbf{S}}_1)\Vert \Vert \mathrm{vech}({\mathbf{S}}_2)\Vert },$$

$$1-\frac{{\mathrm{vech}}^{\ast }({\mathbf{S}}_1{)}^T{\mathrm{vech}}^{\ast }({\mathbf{S}}_2)}{\Vert {\mathrm{vech}}^{\ast }({\mathbf{S}}_1)\Vert \Vert {\mathrm{vech}}^{\ast }({\mathbf{S}}_2)\Vert },$$

$$T=\max \{T_{12},T_{13},\dots ,T_{(K-1)K}\},$$

$$H_0:{\mathbf{\Sigma }}_1={\mathbf{\Sigma }}_2={\mathbf{\Sigma }}_3={\mathbf{\Sigma }}_4={\mathbf{\Sigma }}_5={\mathbf{\Sigma }}_6,$$

$$\mathbf{X}=\mathbf{\Gamma }\mathbf{Z},$$

$$X_{1j}\text{and}{X}_{2j}\text{have the same distribution for all}{j}^{\prime }\text{s}:j=1,\cdots ,p.$$

$$X_{ij}=Z_{ij}+2Z_{i(j+1)},$$

$$X_{2j}=Z_{2j}+2Z_{2(j+1)}+Z_{2(j+2)},$$

$${\mathbf{\Sigma }}_{null}=({\mathbf{\Sigma }}^{\ast }+{\delta }_0\mathbf{I})/(1+{\delta }_0)+{\delta }_1\mathbf{I}.$$

$${\mathbf{\Sigma }}_{alt}=({\mathbf{\Sigma }}^{\ast }+{\delta }_0\mathbf{I})/(1+{\delta }_0)+\mathbf{\Delta }+{\delta }_1\mathbf{I}.$$

$$\Sigma =\left[\begin{array}{cc}{\Sigma }_{{\mathbf{X}}_1{\mathbf{X}}_1}& {\Sigma }_{{\mathbf{X}}_1{\mathbf{X}}_2}\\ {\Sigma }_{{\mathbf{X}}_2{\mathbf{X}}_1}& {\Sigma }_{{\mathbf{X}}_2{\mathbf{X}}_2}\end{array}\right],$$

$$H_0:\Sigma =\left[\begin{array}{cc}{\Sigma }_{{\mathbf{X}}_1{\mathbf{X}}_1}& 0\\ 0& {\Sigma }_{{\mathbf{X}}_2{\mathbf{X}}_2}\end{array}\right].$$

Full text

KEY WORDS: Sphericity test; Identity test; Compound symmetry test; Two-sample test; $K$-sample test

1 Introduction

Testing covariance and correlation matrices is an enduring topic in statistics. In addition to its apparent importance in classical multivariate analysis such as multivariate analysis of variance, Fisher’s linear discriminant analysis, principle component analysis, and repeated measure analysis [Johnson and Wichern (2007], the hypothesis testing of covariance or correlation matrices has a wide range of applications in many scientific fields such as psychometric analysis of behavior difference [Green (1992], co-movement dynamics of credit risk [Annaert, Claes, and De Ceuster (2006], genetic analysis of gene expression [Singh, Febbo, Ross, Jackson, Manola, Ladd, Tamayo, Renshaw, D’Amico, Richie, et al. (2002], dynamic modeling of the correlation structure of the stock market [Djauhari and Gan (2014].

Traditionally, the tests of covariance or correlation matrices are discussed separately in the literature. ?) proposed a method to test if a sample correlation matrix is an identity matrix. ?) proposed an asymptotic $\chi^{2}$ test for the equality of two correlation matrices of multivariate normal distributions. Jennrich’s test has been adopted in several statistical softwares and applied in financial studies. Likelihood ratio based tests for covariance matrices of multivariate normal distributions are well documented in the classical multivariate analysis literature (e.g. ?), Chp. 8). Since the classical likelihood ratio based tests are generally not applicable in high-dimensional data due to the singularity of sample covariance matrices, the recent methodology developments focus on the tests of high-dimensional covariance and correlation matrices, where the dimension $p$ is much bigger than the sample size $n$. Among others, ?) studied the performance of two classical statistics for the one-sample test of sphericity and identity in a high-dimensional setting. ?) adopted the similar test statistics as in ?) but used U-statistics to estimate the population covariance matrix. ?) proposed a two-sample equality test of covariance matrices. ?) provided a summary of some classical and high-dimensional tests of covariance matrices of multivariate normal distributions. ?) discussed the recent developments on estimation and hypothesis tests of covariance matrices in high dimensions from the random matrix theory viewpoint. ?) presented an up-to-date review of methodology developments of the equality test of high-dimensional covariance matrices with sparse signals.

There obviously lacks a unified and flexible hypothesis testing framework for practitioners to handle both covariance and correlation matrices without switching between different test statistics or procedures. In this study we propose such a unified and easy-to-use approach for many common tests of covariance and correlation matrices. We exploit a simple geometric property of symmetric matrices to construct test statistics and develop the associated permutation algorithms to compute p-values under a variety of null hypotheses. Our proposed method can handle a variety of one-sample tests including the tests of sphericity, identity and compound symmetry, and two- and $K$-sample equality tests ($K>2$). One advantage of our proposed sphericity and compound symmetry tests is that the proposed method does not involve either the unknown population variance or the intra-class correlation coefficient. Moreover, the proposed method is very flexible in the sense that it does not assume any underlying distributions or data generation models. Through both real data analyses and extensive simulations, we demonstrate that the proposed method can control the type I error at the pre-specified nominal level and provide good power under the alternative in both one-sample and $K$-sample tests, even for high-dimensional data where the dimension $p$ is much bigger than the sample size $n$.

This paper is organized as follows. Section 2 introduces a generalized cosine measure between two symmetric matrices and, based on this notion of generalized cosine, the test statistics for one-sample and two-sample tests are proposed. Section 3 studies the one-sample tests of sphericity, identity and compound symmetry. Section 4 discusses the two-sample equality test of covariance or correlation matrices and extends the method to handle the $K$-sample problem. Section 5 contains extensive simulation studies on the performance of the proposed one-sample and two-sample tests. In addition, Section 5 probes into factors that affect the performance of the proposed methods. The paper concludes with discussions in Section 6.

2 Test statistic

Below is the definition of the generalized cosine measure between two symmetric matrices.

Definition**.**

Let $\mathbb{R}^{p\times p}$ be a vector space of symmetric matrices of size $p\times p$ and $f$ : $\mathbb{R}^{p\times p}$ $\rightarrow$ $\mathbb{R}^{m\times n}$. The angle between two symmetric matrices $\bm{M}_{1}$ and $\bm{M}_{2}$ in $R^{p\times p}$ with respect to $f$ is defined as $arccos(\bm{M}_{1},\bm{M}_{2})$ and

[TABLE]

where $\langle\cdot,\cdot\rangle$ is an inner product and $\|\cdot\|$ is the corresponding norm in $\mathbb{R}^{m\times n}$.

We use two symmetric matrices $\bm{A}$ and $\bm{B}$ to illustrate some possible cosine values with different mappings, their associated inner products and norms. Matrix $\bm{A}$ is

[TABLE]

which is a modified Hilbert matrix defined as follows:

[TABLE]

for $i,j=1,\ldots 5$. Matrix $\bm{B}$ is

[TABLE]

which is a realization of the random matrix $\bm{A^{*}}$, constructed from $\bm{A}$, as follows:

[TABLE]

where $\epsilon$ is generated from a uniform distribution over the interval (0, 0.5).

Case 1: Suppose $f$ is an identity mapping, i.e., $f(\bm{A})=\bm{A}$. We consider the Frobenius inner product $\langle\bm{A},\bm{B}\rangle_{F}=\operatorname{tr}(\bm{A}^{T}\bm{B})$ and its corresponding norm $\|\bm{A}\|_{F}=\operatorname{tr}(\bm{A}^{T}\bm{A})^{\frac{1}{2}}$. By Equation (1), the cosine between $\bm{A}$ and $\bm{B}$ is calculated as

[TABLE]

which is 0.92, corresponding to an angle of 0.41 radian or [math].

Case 2: In this case, we construct a two-step mapping $f$ : $\mathbb{R}^{p\times p}$ $\rightarrow$ $\mathbb{R}^{m}$. In step one we apply the Cholesky decomposition to obtain $\bm{A}=\bm{L}_{1}\bm{L}_{1}^{T}$ and $\bm{B}=\bm{L}_{2}\bm{L}_{2}^{T}$, where $\bm{L}_{1}$ and $\bm{L}_{2}$ are lower triangle matrices with positive diagonal elements. In step two we apply the half-vectorization operator (which vectorizes or stacks only the lower trianguler portion of a symmetric matrix) to $\bm{L}_{1}$ and $\bm{L}_{2}$ separately to obtain two vectors $\operatorname{vech}(\bm{L}_{1})$ and $\operatorname{vech}(\bm{L}_{2})$, each with a length of $m=p(p+1)/2$. To apply Equation (1), we adopt the Euclidean inner product and norm to compute cosine as

[TABLE]

which is 0.87, corresponding to an angle of 0.52 radian or [math]. This mapping is feasible only if $\bm{A}$ and $\bm{B}$ are positive definite matrices.

Case 3: The mapping $f$ can also be constructed by first applying the eigen-decomposition to obtain $\bm{Q}_{1}\bm{\Lambda}_{1}\bm{Q}_{1}^{T}$ for $\bm{A}$ and $\bm{Q}_{2}\bm{\Lambda}_{2}\bm{Q}_{2}^{T}$ for $\bm{B}$, respectively. $\bm{Q}_{1}$ and $\bm{Q}_{2}$ are orthogonal matrices of eigenvectors, while $\bm{\Lambda}_{1}$ and $\bm{\Lambda}_{2}$ are diagonal matrices whose elements are the eigenvalues of $\bm{A}$ and $\bm{B}$, respectively. Provided that not all diagonal elements in $\bm{\Lambda}_{1}$ and $\bm{\Lambda}_{2}$ are zero, we can extract the diagonal elements from $\bm{\Lambda_{1}}$ and $\bm{\Lambda_{2}}$ to create vectors $v_{1}$ and $v_{2}$, respectively. Thus, we construct a mapping $f$ : $\mathbb{R}^{p\times p}$ $\rightarrow$ $\mathbb{R}^{p}$. To apply Equation (1), we adopt the Euclidean inner product and norm to compute the cosine as

[TABLE]

which is 0.93, corresponding to a angle of 0.37 radian or approximately [math].

Case 4: Since a $p\times p$ symmetric matrix is completely determined by its lower triangular elements together with the symmetry, we construct the mapping $f$ by applying the half-vectorization operator on $\bm{A}$ and $\bm{B}$ directly. The cosine between $\bm{A}$ and $\bm{B}$ can be obtained from Equation (1) as follows:

[TABLE]

which is 0.94, corresponding to a angle of 0.35 radian or [math].

For correlation matrices, the computation can be simplified by introducing a modified half-vectorization operator $\operatorname{vech}^{*}(\cdot)$ from $\operatorname{vech}(\cdot)$ by excluding the diagonal elements of the matrix. Suppose $\bm{A}$ and $\bm{B}$ are two correlation matrices, the cosine can be computed by Equation (1) as

[TABLE]

which is 0.95, corresponding to 0.31 radian or [math].

In summary, we have considered four different mappings to compute the generalized cosine between the two modified Hilbert matrices $\bm{A}$ and $\bm{B}$. The mapping based on the Cholesky decomposition gives the lowest cosine value of 0.87. The mappings based on the Frobenius inner product and the eigen-decomposition give the similar cosine values of 0.92 and 0.93, respectively. When the redundant information is removed, the half-vectorization and the modified half-vectorization operators give the largest cosine values of 0.94 and 0.95, respectively. In this study we choose the half-vectorization operator, $\operatorname{vech}(\cdot)$, for covariance matrices and the modified half-vectorization operator, $\operatorname{vech}^{*}(\cdot)$, for correlation matrices, because these two operators completely remove the redundancy in symmetric matrices and are very easy to compute. Moreover, these two operators demonstrate better statistical power in our pilot simulation studies. The cosine value computed from Equation (1) measures the similarity between two symmetric matrices. When this value is one, the two matrices are identical. We exploit this property to propose new test statistics for one-sample and two-sample tests of covariance and correlation matrices.

One-sample test statistic**.**

Let $\bm{S}$ be the sample covariance (correlation) matrix with the corresponding population covariance (correlation) matrix $\bm{\Sigma}$. The test statistic for the equality of $\bm{\Sigma}$ and $\bm{\Sigma}_{0}$, which is either known or specified up to some unknown parameters, is

[TABLE]

where the cosine is computed according to Equation (1) under a pre-determined mapping $f$ with a properly chosen inner product and norm.

Two-sample test statistic**.**

Let $\bm{S}_{1}$ and $\bm{S}_{2}$ be the sample covariance (correlation) matrices with the corresponding population covariance (correlation) matrices $\bm{\Sigma}_{1}$ and $\bm{\Sigma}_{2}$, respectively. The test statistic for the equality of $\bm{\Sigma}_{1}$ and $\bm{\Sigma}_{2}$ is

[TABLE]

where the cosine is computed according to Equation (1) under a pre-determined mapping $f$ with a properly chosen inner product and norm.

It is worth pointing out that the two newly proposed test statistics exploit a simple geometric property of symmetric matrices and require no distributional assumptions of data. We further adopt a permutation approach to obtain the null distributions of the two test statistics and compute p-values under the null hypotheses of different tests.

3 One-sample test

Suppose a matrix, $\bm{D}_{n\times p}=\{x_{ij}\}$, $i=1,\ldots,n$ and $j=1,\ldots,p$, contains data collected on $n$ subjects and each with $p$ covariates. Its $i$th row is a realization of random vector $\bm{X}=(X_{1},X_{2},\ldots,X_{p})^{T}$, which has an unknown distribution with a covariance matrix $\bm{\Sigma}$ and a correlation matrix $\bm{R}$. In the one-sample problem one may use the data matrix $\bm{D}_{n\times p}$ to test if the covariance matrix $\bm{\Sigma}$ is proportional to an identity matrix, known as the sphericity test, or to test if $\bm{\Sigma}$ or $\bm{R}$ equals to an identity matrix, known as the identity test.

3.1 Sphericity test

The hypotheses in the sphericity test for a covariance matrix are as follows:

[TABLE]

where $\sigma^{2}$ is the unknown population variance under the null. In this case, $\bm{\Sigma}_{0}$ in the one-sample test statistic (9) is $\sigma^{2}\bm{I}$ and $\bm{S}$ is the sample covariance matrix of the data matrix $\bm{D}_{n\times p}$. To test the sphericity of a covariance matrix, using $\operatorname{vech}(\cdot)$, the Euclidean inner product and norm, the test statistic (9) can be shown to be

[TABLE]

where $\operatorname{tr}(\cdot)$ denotes the trace. The derivation of (12) is provided in Appendix I. Note that the test statistic (12) does not depend on the unknown parameter $\sigma^{2}$.

The sphericity test (11) is a combination of two tests [Fujikoshi, Ulyanov, and Shimizu (2011]: 1) $H_{0}:\bm{\Sigma}=\sigma^{2}\bm{I}$ vs. $H_{1}:\bm{\Sigma}=\mathrm{diag}(\sigma_{1}^{2},\sigma_{2}^{2},\ldots,\sigma_{p}^{2})$ and 2) $H_{0}:\bm{\Sigma}=\mathrm{diag}(\sigma_{1}^{2},\sigma_{2}^{2},\ldots,\sigma_{p}^{2})$ vs. $H_{1}:\bm{\Sigma}\neq\mathrm{diag}(\sigma_{1}^{2},\sigma_{2}^{2},\ldots,\sigma_{p}^{2})$. To handle these two tests simultaneously, we propose Algorithm 1 to compute p-values under the null of sphericity. This algorithm first permutes the elements in each row of $\bm{D}_{n\times p}$ to generate a permuted data matrix $\bm{D}^{\dagger}_{n\times p}$, then the elements in each column of $\bm{D}^{\dagger}$ to generate $\bm{D}^{\ddagger}_{n\times p}$. The permutation on each row of the data matrix is to test if the diagonal elements in $\bm{\Sigma}$ are identical; the permutation on each column is to test if the off-diagonal elements of $\bm{\Sigma}$ are zero. Algorithm 1 assumes that the elements in the random vector $\bm{X}$ are exchangeable in the sense that $\bm{\Sigma}$ is invariant to the permutations of the elements in $\bm{X}$ under the null hypothesis of (11). This algorithm requires that the elements in $\bm{X}$ have the same distribution so that the elements in each row of $\bm{D}_{n\times p}$ can be permuted.

We illustrate the proposed sphericity test and Algorithm 1 using the bfi data in the R package pysch [Revelle (2018]. The bfi data contain twenty-five personality self reported items from 2800 subjects who participated in the Synthetic Aperture Personality Assessment (SAPA) web based personality assessment project. These twenty-five items are organized by five putative factors: Agreeableness, Conscientiousness, Extraversion, Neuroticism, and Openness. Each subject answered these personality questions using a six point response scale ranging from 1 (Very Inaccurate) to 6 (Very Accurate). After removing the missing data, we apply the proposed sphericity test and Algorithm 1 to test the sphericity of the $25\times 25$ covariance matrix of personality measurements of 2436 subjects. We obtain a p-value of 0.01 based on 100 permutations using Algorithm 1. Therefore, the null hypothesis of sphericity is rejected at the significance level of 0.05. We will further investigate the empirical type I error and power of our proposed sphericity test in Section 5.

3.2 Identity test

The identity test of a covariance matrix can be formulated as [Ledoit and Wolf (2002]

[TABLE]

which includes the more general null hypothesis $H_{0}:\bm{\Sigma}=\bm{M}_{0}$, where $\bm{M}_{0}$ is a known covariance matrix, since we can replace $\bm{M}_{0}$ by $\bm{I}$ and multiply data by $\bm{M}_{0}^{-1/2}$.

Because $\bm{I}$ in the hypothesis (13) and $\sigma^{2}\bm{I}$ in the sphericity hypothesis (11) only differ by a constant $\sigma^{2}$ and the test statistic (12) does not depend on $\sigma^{2}$, we can adopt the test statistic (12) to test the hypothesis (13).

For the identity test of a correlation matrix,

[TABLE]

for which, the test statistic (12) can be simplified into

[TABLE]

The derivation of (15) is provided in Appendix I.

Since the identity test mainly concerns the off-diagonal elements of a covariance or correlation matrix, we propose Algorithm 2, which permutes the elements in each column of the data matrix $\bm{D}_{n\times p}$, to compute p-values under the null. Permuting the elements in columns can disrupt the covariance structure, but it does not change the sample variance of each column of the data matrix, i.e. the diagonal elements of the sample covariance matrix remain unchanged in different permutations. This implies that, for a given $\bm{D}_{n\times p}$, $\operatorname{tr}(\bm{S})$ remains a constant in (12) in different permutations. For a correlation matrix $\bm{S}$, we have $\operatorname{tr}(\bm{S})$=$p$. Since the column-wise permutations can sufficiently disassociate every pair of random variables, and hence, disrupt the covariance structure of the $p$ random variables, the off-diagonal elements of the sample covariance matrix of the permuted data are much closer to zero than the original sample covariance (correlation) matrix under the alternative. Therefore, the sample covariance and correlation matrices with a good portion of non-zero off-diagonal elements have larger values of the statistics of (12) and (15), respectively, to reject the null hypothesis of identity. Moreover, Algorithm 2 does not require each column of a data matrix to have the same distribution.

In addition to the Pearson’s correlation matrix, the test statistic defined in (15) and its associated Algorithm 2 can handle the Spearman’s and the Kendall’s rank correlation matrices, because this test statistic is based on the cosine measure between a sample correlation matrix $\bm{S}$ and $\bm{I}$, and does not require any particular types of correlation matrices.

We use the bfi data as described in section 3.1 to illustrate the proposed identity test and its associated Algorithm 2 for correlation matrices. For a Pearson’s correlation matrix of twenty-five personality self reported items, Algorithm 2 produces a p-value of 0.01 based on 100 permutations. The Bartlett’s identity test implemented in the same R package gives a p-value of zero. We also apply the identity test to the Spearman’s rank and Kendell’s rank correlation matrices. Algorithm 2 produces a p-value of 0.01 for both types of correlation matrices based on 100 permutations . Therefore, we reject the identity of these three types of correlation matrices. We will further study the performance of our proposed identity test in Section 5.

3.3 Compound symmetry test

The hypothesis of compound symmetry is $H_{0}:\bm{\Sigma}=\sigma^{2}[(1-\rho)\bm{I}+\rho\bm{J}]$ for a covariance matrix or $H_{0}:\bm{R}=[(1-\rho)\bm{I}+\rho\bm{J}]$ for a correlation matrix, where $\bm{J}$ is a square matrix of ones and $\rho>0$ is the intra-class correlation coefficient. In this case, $\Sigma_{0}$ in the one-sample test statistic (9) is $\sigma^{2}[(1-\rho)\bm{I}+\rho\bm{J}]$ for a covariance and $[(1-\rho)\bm{I}+\rho\bm{J}]$ for a correlation matrix, respectively. Using the modified half-vectorization operator, $\operatorname{vech}^{*}(\cdot)$, a unified test statistic for testing the compound symmetry of both covariance and correlation matrices can be obtained from (9) as

[TABLE]

where $\bm{S}$ is either a sample covariance or correlation matrix. The derivation of (16) is provided in Appendix II. Note that this test statistic does not depend on the unknown parameters $\sigma^{2}$ and $\rho$.

Suppose a data matrix $\bm{D}_{n\times p}$ contains rows that are the realizations of a $p\times 1$ random vector $\bm{X}=(X_{1},\ldots,X_{p})^{T}$, in which all $X_{i}$’s have the same distribution. Under the null hypothesis of compound symmetry, it is true that $Var(X_{i})=\sigma^{2}$, for all $i$’s and $Cov(X_{i},X_{j})=\sigma^{2}(1-\rho)$ for all pairs of $i,j$, provided $i\neq j$. This implies that the random permutations of the elements in $\bm{X}$ do not alter its covariance matrix under the assumption of compound symmetry. Using this fact, we propose Algorithm 3 to compute p-values for the test statistic (16) by randomly permuting the elements in each row of the data matrix. For this algorithm to be applicable, each column of the data matrix must have the same distribution.

For an illustrative purpose, we consider the Cork data (?), Table 6.21). Four cork borings, one in each of the four cardinal directions, were taken from each of twenty-eight trees and their weights were measured in integers. The goal is to test the compound symmetry of the covariance matrix of these weight measurements taken in the four cardinal directions. Using the modified half-vectorization operator $\operatorname{vech}^{*}(\cdot)$ and the Euclidean inner product and norm, we obtain a cosine of 0.99 (0.998) between the sample covariance (correlation) matrix and the compound symmetry, indicating a strong similarity. We apply our proposed compound symmetry test (16) to this dataset and, using Algorithm 3, obtain a p-value of 0.099 for the compound symmetry test of the sample covariance matrix and 0.069 for the sample correlation matrix based on 100 permutations. Hence, the null hypothesis of compound symmetry can not be rejected at the significance level of 0.05. ?), however, rejected the compound symmetry based on a $\chi^{2}$-approximation to a likelihood ratio test using a multivariate normal distribution (?), Example 7.2.3). This discrepancy is due to their inappropriate adoption of the multivariate normal distribution, since three of the four marginal distributions failed to pass the normality test. The test statistic (16) and Algorithm 3, on the other hand, require no distributional assumptions and are robust to the underlying distribution of Cork borings weights.

4 $K$-sample test ($K>2$)

In this section we will first discuss the test statistics and their associated algorithms for the two-sample problem ($K$=2), then extend the method to handle the $K$-sample problem. Let $\bm{D}_{n_{1}\times p}$ and $\bm{D}_{n_{2}\times p}$ be two data matrices with sample covariance (correlation) matrices $\bm{S}_{1}$ and $\bm{S}_{2}$, respectively. Suppose the rows of $\bm{D}_{n_{1}\times p}$ are the realizations of a $p\times 1$ random vector $\bm{X}_{1}=(X_{11},\ldots,X_{1p})^{T}$ and the rows of $\bm{D}_{n_{2}\times p}$ are of random vector $\bm{X}_{2}=(X_{21},\ldots,X_{2p})^{T}$. Let $\bm{\Sigma}_{1}$ ($\bm{R}_{1}$) be the covariance (correlation) matrix of $\bm{X}_{1}$ and $\bm{\Sigma}_{2}$ ($\bm{R}_{2}$) be the covariance (correlation) matrix of $\bm{X}_{2}$.

In the two-sample problem, the null hypotheses are $H_{0}:\bm{\Sigma}_{1}=\bm{\Sigma}_{2}$ for covariance matrix and $H_{0}:\bm{R}_{1}=\bm{R}_{2}$ for correlation matrix. We adopt $\operatorname{vech}(\cdot)$, the Euclidean inner product and norm in the two-sample test statistic (10) to test the equality of covariance matrices, leading to a test statistic

[TABLE]

and $\operatorname{vech}^{*}(\cdot)$, the Euclidean inner product and norm to test the equality of correlation matrices, yielding a test statistic

[TABLE]

where $\bm{S}_{1}$ and $\bm{S}_{2}$ are the sample covariance or correlation matrices.

We propose Algorithm 4 to compute p-values for these test statistics under the null hypothesis of equality. In this algorithm two data matrices $\bm{D}_{n_{1}\times p}$ and $\bm{D}_{n_{2}\times p}$ are stacked to form a new data matrix $\bm{D}_{n\times p}$ ($n=n_{1}+n_{2}$). In each permutation the rows of $\bm{D}_{n\times p}$ are randomly permuted to generate a permuted data matrix $\bm{D}^{*}_{n\times p}$, which is then split into two data matrices of $\bm{D}_{n_{1}\times p}^{*}$ and $\bm{D}_{n_{2}\times p}^{*}$ to compute the test statistic (17) or (18). Algorithm 4 assumes that $X_{1i}$ in $\bm{X}_{1}$ and $X_{2i}$ in $\bm{X}_{2}$ have the same distribution for all $i$’s. One advantage of our proposed two-sample test is that $X_{1i}$ and $X_{1j}$, $i\neq j$, need not to have the same distribution.

Under the null hypothesis of equality $\bm{\Sigma}_{1}=\bm{\Sigma}_{2}=\bm{\Sigma}$, the sample covariance matrices $\bm{S}_{1}^{*}$ of $\bm{D}_{n_{1}\times p}^{*}$ and $\bm{S}_{2}^{*}$ of $\bm{D}_{n_{2}\times p}^{*}$ can be considered as the realizations of matrix $\bm{\Sigma}$. The rationale behind Algorithm 4 is that the cosine value between $\bm{S}_{1}^{*}$ and $\bm{S}_{2}^{*}$ is similar to that of $\bm{S}_{1}$ and $\bm{S}_{2}$ under the null hypothesis and the permutations provide a good control of the type I error. Under the alternative, the repeated random-mixing rows of $\bm{D}_{n_{1}\times p}$ and $\bm{D}_{n_{2}\times p}$ produce $\bm{S}_{1}^{*}$ and $\bm{S}_{2}^{*}$ such that the cosine value between the two is bigger than that of $\bm{S}_{1}$ and $\bm{S}_{2}$, therefore the test statistics (17) and (18) have good power to reject the null at a pre-determined significance level.

We use the flea beetles data (?), Table 5.5) to illustrate our proposed two-sample equality test. This dataset contains four measurements of two species of flea beetles, Haltica oleraces and Haltica cardorum. To test the equality of two covariance matrices of the four measurements of the two flea beetle species, we apply the proposed two-sample test (17) and obtain a p-value of 0.37 using Algorithm 4 with 100 permutations. Assuming normality and using $\chi^{2}$- and $F$- approximation, the $Box^{\prime}s$ $M$-$test$ gives a similar scale p-value of 0.56 [Rencher and Christensen (2012]. In this case our proposed equality test (17) agrees with the likelihood ratio based $Box^{\prime}s$ $M$-$test$, and the null hypothesis of the equality of two covariance matrices is not rejected at the significance level of 0.05.

To test the multivariate equality of several covariance or correlation matrices, we consider $H_{0}:\bm{\Sigma}_{1}=\bm{\Sigma}_{2}=\cdots=\bm{\Sigma}_{K}$ for the equality of multiple covariance matrices and $H_{0}:\bm{R}_{1}=\bm{R}_{2}=\cdots=\bm{R}_{K}$ for the equality of multiple correlation matrices. We propose the test statistic

[TABLE]

where $T_{ij}$ is the test statistic (17) or (18) for the pairwise two-sample comparison of $\bm{\Sigma}_{i}$ ($\bm{R}_{i}$) and $\bm{\Sigma}_{j}$ ($\bm{R}_{j}$) for all possible unique pairs of $K$ populations. Let $\bm{D}_{n_{1}\times p},\bm{D}_{n_{2}\times p},\ldots,\bm{D}_{n_{K}\times p}$ be the data matrices from $K$ populations. Algorithm 5 provides a permutation approach to compute p-values under the null hypothesis of the multivariate equality of covariance or correlation matrices.

We use the Rootstock data (?), Table 6.2) to illustrate our proposed $K$-sample test. Eight apple trees from each of the six rootstocks were measured on four variables: 1) trunk girth at 4 years; 2) extension growth at 4 years; 3) trunk girth at 15 years; 4) weight of the tree above ground at 15 years. The null hypothesis is

[TABLE]

for the equality of six covariance matrices of the six rootstocks. The likelihood ratio based $Box^{\prime}s$ $M$-$test$ yields a p-value of 0.71 and 0.74 based on $\chi^{2}$- and $F$- approximation [Rencher and Christensen (2012], respectively. We apply the multivariate equality test (19) to this dataset and obtain a p-value of 0.99 using Algorithm 5 with 100 permutations. For this dataset, our method agrees with the $Box^{\prime}s$ $M$-$test$, and the null hypothesis of the multivariate equality of these six covariance matrices is not rejected at the significance level of 0.05.

5 Simulation studies

In this section we carry out extensive simulations to investigate the empirical type I error and power of our proposed one-sample and two-sample tests. In the one-sample setting we investigate the proposed sphericity test and identity test. In the two-sample setting, we investigate the proposed equality test of two covariance matrices. We design and conduct these simulation studies with the backdrop of high dimensions. We also probe into factors that may affect the performance of the proposed methods. All simulations are conducted in the R computing environment [R Core Team (2018].

5.1 One-sample tests

5.1.1 Sphericity test

In this section we evaluate the empirical type I error and power of our proposed sphericity test. We adapt the simulation design and model of ?) for a $p$-dimensional random vector $\bm{X}=\left(X_{1},\ldots,X_{p}\right)^{T}$ as

[TABLE]

where $\bm{\Gamma}$ is a $p\times m$ constant matrix with $p\leq m$ and $\bm{Z}=(Z_{1},\cdots,Z_{m})^{T}$ is a $m$-dimensional random vector with $Var(\bm{Z})=\bm{I}$. We further let $\bm{\Gamma}\bm{\Gamma}^{T}=\bm{\Sigma}$= $Var(\bm{X})$. The elements $Z_{i}$’s in $\bm{Z}$ are IID random variables with a pre-specified distribution. It is worth noting that Model (21) is more general than Model (2.4) in ?), where they additionally require that $E(Z_{i})=0$ and a uniform bound for the 8th moment for all $Z_{i}$’s.

We consider two distributions for $Z_{i}$: 1) $N(0,1)$; 2) Gamma(4,0.5). ?) also considered these two distributions but forcing the mean of Gamma(4,0.5) to be zero in order to meet their requirement of data generation model for their sphericity test.

To evaluate the empirical type I error, we set $\bm{\Gamma}=\bm{I}$ to simulate data under the null with $\bm{\Sigma}=\bm{I}$. To evaluate the power, following ?), we consider two forms of alternatives: 1) $\bm{\Sigma}=\sigma_{0}^{2}\bm{I}+\sigma_{1}^{2}\bm{A}$, where $\bm{A}=\mathrm{diag}(\bm{I}_{\lfloor\lambda p\rfloor},\bm{0}_{p-\lfloor\lambda p\rfloor})$, $\sigma_{0}^{2}=\sigma_{1}^{2}=1$ and $\lfloor\cdot\rfloor$ is the floor function. In this case $\bm{\Gamma}=\text{diag}\left(\sqrt{2}\bm{I}_{\lfloor\lambda p\rfloor},\bm{I}_{p-\lfloor\lambda p\rfloor})\right)$ and $\bm{Z}$ is a $p\times 1$ random vector with elements that are IID $N(0,1)$ or Gamma(4,0.5) random variables. 2) $\bm{\Sigma}=(1-\rho)\sigma_{0}^{2}\bm{I}+\rho\sigma_{1}^{2}\bm{J}$ with $\sigma_{0}^{2}=1$ and $\sigma_{1}^{2}=2$. In this case $\bm{\Gamma}=\text{diag}\left(\sqrt{1-\rho}\bm{I},\sqrt{2}\bm{1}\right)$ and $\bm{Z}$ is a $(p+1)\times 1$ random vector with elements that are IID $N(0,1)$ or Gamma(4,0.5) random variables. We choose the challenging cases in ?) by setting $\lambda=0.125$ in the first alternative and $\rho=0.1$ in the second, where the proposed tests of ?) had lower power in their simulation study.

We apply the sphericity test (12) to a collection of simulated datasets with the varying samples sizes and dimensions as considered in ?), and compute p-values with 100 permutations using Algorithm 1. We replicate each simulation setting 2000 times and reject the null hypothesis at the significance level of 0.05. Table LABEL:tab:sph_null shows that the proposed sphericity test produces p-values compatible with the nominal level of 0.05 under the null. Table LABEL:tab:sph_pw1 presents the empirical power of our proposed sphericity test under the first alternative, where for a fixed dimension, the empirical power increases quickly with the sample size for both distributions. For a fixed sample size, except for the case with the sample size of twenty, the empirical power increases with the dimension. Table LABEL:tab:sph_p2 shows that our proposed sphericity test demonstrates superb empirical power, almost 100$\%$, under the second alternative.

5.1.2 Identity test

In this section we evaluate the performance of the proposed identity test (12) for covariance or correlation matrices. We design a block-diagonal model for a $p\times 1$ random vector $\bm{X}=(\bm{X}_{1}^{T},\bm{X}_{2}^{T},\bm{X}_{3}^{T},\bm{X}_{4}^{T})^{T}$, where $\bm{X}_{i}$ is a $q\times 1$ sub-vector, $i=1,2,3,4$, such that $Var(\bm{X})=\bm{\Sigma}$, $Var(\bm{X}_{i})=\bm{\Sigma}_{i}$ and $Cov(\bm{X}_{i},\bm{X}_{j})=\bm{0}$ for $i\neq j$. The covariance matrix in this block-diagonal model takes the form of $\bm{\Sigma}=\text{diag}(\bm{\Sigma}_{1},\cdots,\bm{\Sigma}_{4})$, where $\bm{\Sigma}_{i}$ is a $q\times q$ matrix for all $i$’s and $\bm{\Sigma}$ has dimensions of $p\times p$ with $p=4q$. We first simulate each $\bm{X}_{i}=\bm{L}\bm{U}_{i}$, where $\bm{U}_{i}$ is a $q\times 1$ random vector and $\bm{L}$ is a lower triangle matrix from the Cholesky decomposition $\bm{\tilde{\Sigma}_{q}}=\bm{L}\bm{L}^{T}$, where $\bm{\tilde{\Sigma}_{q}}$ has a structure of $(1-\rho)\bm{I}+\rho\bm{J}$. Thus $\bm{\Sigma}_{i}=\bm{L}\bm{\Lambda}_{i}\bm{L}^{T}$, where $\bm{\Lambda}_{i}$ is a diagonal matrix whose diagonal elements are the variances of the elements of $\bm{U}_{i}$. Then these four $\bm{X}_{i}$’s are stacked to obtain the random vector $\bm{X}$.

We consider two configurations for $\bm{X}$. In the first configuration, the four $\bm{U}_{i}$’s are IID and the elements in all $\bm{U}_{i}$’s have one of the following four distributions: 1) $N(0,1)$; 2) Log-normal(0,1); 3) Student $t_{5}$; 4) Gumbel(10, 2) to generate the random vector $\bm{X}$. In the second configuration, the four $\bm{U}_{i}$’s are independently, but not identically, distributed and each $\bm{U}_{i}$ has a different distribution. The second configuration is referred as a hybrid configuration. The reason to choose these four distributions is to investigate the impact of different types of distributions, such as asymmetric distributions, heavy-tail distributions, and extreme value distributions, on the performance of the proposed identity test (12) and its associated Algorithm 2.

We also consider two scenarios for the sample size: 1) a low-dimensional case with $n>p$; 2) a high-dimensional case with $n\ll p$. Using the block-diagonal model, we generate data under the null with $\rho=0$ to evaluate the empirical type I error and under the alternative with $\rho=0.15$ to evaluate the power. Using Algorithm 2, we compute p-values based on 100 permutations and reject the null hypothesis at the significance level of 0.05. We replicate each simulation setting 2000 times. Table LABEL:tab:idt summarizes the results of the empirical type I error and power of our proposed identity test. It is obvious that the p-values are compatible with the nominal level of 0.05 under the null. The empirical power is close to 100% when the sample size is bigger than the dimension of the data. When the sample size is much smaller than the dimension, the empirical power increases with the increase of the dimension for this block-diagonal model, though there are some variations across different distributions.

5.2 Two-sample test

In this section we study the empirical type I error and power of our proposed two-sample test (17) for the equality of two covariance matrices. We consider two dimensionality configurations: 1) $\min(n_{1},n_{2})>p$ and 2) $\max(n_{1},n_{2})\ll p$, where $n_{1}$ and $n_{2}$ are the two sample sizes and $p$ is the dimension. We adopt the block-diagonal model as in Section 5.1.2 to simulate data in these two dimensional configurations. Under the null hypothesis we set $\rho=0.15$ in both samples and under the alternative we set $\rho_{1}=0.15$ for sample one and $\rho_{2}=0.30$ for sample two. We apply Algorithm 4 with 100 permutations to compute p-values and reject the null hypothesis at the significance level of 0.05. We replicate each simulation setting 2000 times.

Table LABEL:tab:two shows that the proposed equality test maintains a good control of the empirical type I error at the nominal level of 0.05 under the null hypothesis regardless of the underlying distributions and data dimensions.

Note that 75% of the elements in the covariance matrix are zero under the block-diagonal model. This implies that at least 75% of the elements of $\bm{\Sigma}_{1}$ and $\bm{\Sigma}_{2}$ are the same regardless of the value of $\rho$, which makes it rather difficult to test the equality of two covariance matrices. In fact, for the same data, the two-sample $Box^{\prime}s$ $M$-$test$, which approximates the two-sample likelihood ratio test based on a multivariate normal distribution, produces empirical power: 30%, 30%, 66%, 96% and 100% for the sample sizes of $n_{1}$=$n_{2}$ =100 and the dimensions $p$= 24, 32, 64, 76, 92, respectively. Clearly, the empirical power is not high even for the $Box^{\prime}s$ $M$-$test$ based on a correct model. The results of our proposed two-sample tests are reported in Table LABEL:tab:two. For the multivariate normal distribution, our proposed method, on average, outperforms the two-sample $Box^{\prime}s$ $M$-$test$. The model with log-normal(0,1) gives the lowest power in the first dimensionality configuration where $n_{1}$=$n_{2}$=100 ($>p$). The model with Gumbel(10, 2) gives the highest power among the four distributions in both configurations. Furthermore, it is interesting to notice that the model with hybrid distributions gives the second highest power in both dimensionality configurations. We discover that the empirical power of our proposed two-sample test is sensitive to the signal-to-noise ratio (SNR) of the underlying distribution, defined as the reciprocal of the coefficient of variation. The model with Gumbel(10,2) has the highest empirical power accompanied with the highest theoretical SNR of 1.70, followed by the model with log-normal(0,1) which has the second highest power with a theoretical SNR of 0.35. The models with $N(0,1)$ and Student $t_{5}$ have similar low empirical power and both have the theoretical SNR of zero. We will further investigate the impact of SNR in Section 5.3.

5.3 Factors affecting the performance of the proposed two-sample test

In this section, we probe into some factors that may affect the performance of our proposed two-sample test (17) and its associated Algorithm 4. Define a random vector $\bm{X}_{1}=(X_{11},\ldots,X_{1p})^{T}$ with $Var(\bm{X}_{1})=\bm{\Sigma}_{1}$ for the population one and $\bm{X}_{2}=(X_{21},\ldots,X_{2p})^{T}$ with $Var(\bm{X}_{2})=\bm{\Sigma}_{2}$ for the population two. We recall the assumption required for the proposed two-sample test as

[TABLE]

We first investigate the impact of violating assumption (22) on the empirical type I error. Secondly, we use distributions with different SNR’s to examine the sensitivity of the proposed two-sample test to the change of SNR. Finally, we study the power of the proposed two-sample test under a sparse alternative as in ?), where the two covariance matrices only differ in a very few and fixed number of off-diagonal elements. In this section we focus on the scenarios where the sample size $n$ is smaller than the dimension $p$.

5.3.1 Impact of assumption (22)

To investigate the impact of violating assumption (22) on our proposed two-sample test, we adapt the moving average models as in ?). To generate $\bm{X}_{1}$ and $\bm{X}_{2}$ under the null hypothesis, we use the model

[TABLE]

where $j=1,\cdots,p$ and $i=1,2$ for samples one and two, respectively. The two sequences $\{Z_{1k}\}$ and $\{Z_{2k}\}$, $k=1,\cdots,p+1$, consist of IID random variables. Under the alternative, we generate the sample one using Model (23) and sample two using the following model

[TABLE]

where the sequence $\{Z_{2j}\}$ consists of $(p+2)$ IID random variables.

To examine the impact of violating assumption (22) under the null hypothesis, we set $\{Z_{1k}\}$ to be IID Gamma(4, 0.5) random variables and $\{Z_{2k}\}$ to be IID Gamma(0.5, $\sqrt{2}$) random variables in Model (23). We choose the smallest sample size of twenty as used in ?). We replicate the simulation 2000 times. In each replicate we compute p-values using 100 permutations in Algorithm 4 and reject the null hypothesis at the significance level of 0.05. Table LABEL:tab:tab3 shows that the violation of assumption (22) can result in incorrect empirical type I errors in our proposed two-sample test. On the other hand, when the two samples are generated using the same distribution, the proposed two-sample test produces correct empirical type I errors compatible with the nominal level of 0.05.

5.3.2 Impact of SNR

To investigate the impact of SNR on the power of our proposed two-sample test, we generate $\bm{X}_{1}$ using Model (23) in the sample one and $\bm{X}_{2}$ using Model (24) in sample two. The sequences $\{Z_{1k}\}$ and $\{Z_{2k}\}$ have the same distribution of either Gamma(4, 0.5) or Gamma(0.5, $\sqrt{2}$). The two Gamma distributions have the same variance but different means, which imply that they have different SNR’s. Table LABEL:tab:tab3 shows that the model with Gamma(4, 0.5) has much higher empirical power than that of the model with Gamma(0.5, $\sqrt{2}$). This is in line with our observation in Section 5.2 that higher SNR’s result in higher power in our proposed two-sample test, since the model with Gamma(4, 0.5) has a SNR that is twice bigger than that of the model with Gamma(0.5, $\sqrt{2}$). To investigate the impact of SNR further, we consider an additional simulation study with two samples generated using the distributions of $N(5,1)$ and $N(5,5^{2})$, respectively, and the former has a SNR that is five times bigger than the latter. We use $n_{1}=n_{2}=20$ for the sample sizes and $p=50$ for the dimension. Simulation results in Table LABEL:tab:tab3 show that the empirical power of the model with $N(5,1)$ is about 100%, whereas the power reduces significantly to 27% using the sample generated by the model with $N(5,5^{2})$.

5.3.3 Impact of sparseness

To investigate the impact of sparse signals, we adopt a model from the supplemental simulation study of [Cai, Liu, and Xia (2013]. Under the null hypothesis, we set $\bm{\Sigma}_{1}=\bm{\Sigma}_{2}=\bm{\Sigma}_{null}$, where

[TABLE]

Under the alternative, we set $\bm{\Sigma}_{1}=\bm{\Sigma}_{null}$ and $\bm{\Sigma}_{2}=\bm{\Sigma}_{alt}$, where

[TABLE]

In models (25) and (26), $\bm{\Sigma^{*}}=(\sigma_{ij}^{2})$ is a symmetric matrix, in which $\sigma_{ii}=1$ and $\sigma_{ij}=0.5*\text{Bernoulli}(1,0.05)$ for $i<j$. $\delta_{0}$ = $|\lambda_{\text{min}}(\bm{\Sigma}^{*})|$ + 0.05 and $\delta_{1}=|\text{min}\{\lambda_{\text{min}}((\bm{\Sigma}^{*}+\delta_{0}\bm{I})/(1+\delta_{0})),\lambda_{\text{min}}((\bm{\Sigma}^{*}+\delta_{0}\bm{I})/(1+\delta_{0})+\bm{\Delta})\}|+0.05$, where $\lambda_{\text{min}}(\cdot)$ returns the smallest eigenvalue of the matrix which it applies to. $\bm{\Delta}$ is a square matrix containing thirty-two entries of 0.9 and all the other entries are zero. The locations of sixteen nonzero entries are randomly selected in the lower triangle portion and the other sixteen nonzero entries are located in the upper triangle portion accordingly to make $\bm{\Delta}$ symmetric. Regardless of the dimension of the covariance matrix, the number of nonzero entries in $\bm{\Delta}$ remains thirty-two to keep the signals sparse. Note that the matrix $\bm{\Delta}$ changes from one replicate to another in this simulation study.

Through the Cholesky decomposition we obtain $\bm{\Sigma}_{null}=\bm{\Gamma}_{1}\bm{\Gamma}_{1}^{T}$ and $\bm{\Sigma}_{alt}=\bm{\Gamma}_{2}\bm{\Gamma}_{2}^{T}$. Under the null we generate $\bm{X}_{1}=\bm{\Gamma}_{1}\bm{U}$ and $\bm{X}_{2}=\bm{\Gamma}_{1}\bm{U}$ to evaluate the empirical type I error, where $\bm{U}$ is a random vector whose elements are IID random variables. Under the alternative we generate $\bm{X}_{1}=\bm{\Gamma}_{1}\bm{U}$ and $\bm{X}_{2}=\bm{\Gamma}_{2}\bm{U}$ to evaluate the power of our proposed two-sample test.

Following ?), we set the sample size to thirty for each sample. We use Algorithm 4 to compute p-values with 100 permutations. We replicate each simulation setting 2000 times. We also study several distributions for the elements in $\bm{U}$. Due to the sampling variations in small samples, the sample SNR’s may differ from their corresponding theoretic SNR’s. We therefore consider the sample signal-to-noise (SSNR) computed using the sample mean and the sample standard deviation. Table LABEL:tab:title presents the influence of sparse signals and SSNR’s on the performance of our proposed two-sample test and its associated Algorithm 4. The results show that, under the null hypothesis, our approach produces p-values compatible with the nominal level of 0.05, regardless of the dimension of data and the level of SSNR’s. In terms of the empirical power, our method is not sensitive to different distributions but to the level of SSNR’s: the higher SSNR is, the higher the empirical power is. Moreover, Table LABEL:tab:title shows that distributions with similar SSNR’s tend to have similar power.

6 Discussion

In this study we introduce a generalized cosine measure between two symmetric matrices and, based on this geometry, we propose new test statistics for one-sample and two-sample tests of covariance and correlation matrices. These test statistics and their companioned algorithms for computing p-values can be applied to a variety of hypothesis tests of covariance and correlation matrices. In the one-sample setting, we implement the test statistics and algorithms for the hypothesis tests of sphericity, identity and compound symmetry. We also propose new test statistics and algorithms to test the equality of two or multiple covariance or correlation matrices. We demonstrate the effectiveness of these tests and algorithms through several real datasets and extensive simulation studies, where data are generated using a variety of models and distributions adopted in the previous literature.

In order for the proposed algorithms to work properly, the number of permutations is important and usually can be determined by the pre-specified significance level. For instance, we recommend at least one hundred permutations for the significance level of 0.05, for this number can provide the required precision to two decimal places. Moreover, the extensive simulations show that one hundred permutations appear to be sufficient to control the empirical type I error and provide good power for a variety of models and distributions at the significance level of 0.05. If a greater significance level is desirable, for example 0.001 in a multiple testing problem, then at least one thousand permutations is recommended. Furthermore, the dimension and the sample size of data also impact on the number of recommended permutations. For algorithms associated with the proposed one-sample tests, the dimension plays a more prominent role than the sample size, since the number of permutations is often determined by the number of columns in the data matrix, i.e., the dimension of data. We show in Table LABEL:tab:idt that the proposed algorithm can provide a correct empirical type I error for data with sample size as small as four in identity test, provided that the dimension of data is sufficiently large. In the $K$-sample test the sample size may put a restriction on the number of possible permutations, since the algorithms for the $K$-sample tests randomly shuffle rows of the stacked data matrix. For example, in the two-sample test, the minimum recommended sample size for the combined sample is ten with the smaller sample containing at least three subjects, for $10\choose 3$ is 120 and exceeds the number of permutations recommended for the significance level of 0.05. Generally for a $K$-sample test, we recommend that ${n\choose n_{(1)}}\geq k_{0}$, where $n_{(1)}$ is the smallest sample size of the $K$ samples and $k_{0}$ is the required number of permutations for a pre-specified significance level.

Our proposed test statistics are general and flexible so that they can be easily extended to handle other complex situations. As an example, we outline one possible extension as follows. It is often of interest to test the independence of several normal random vectors as in ?) or more generally uncorrelation of multiple random vectors. Suppose a random vector $\bm{X}$ can be partitioned into $k$ subvectors of $\bm{X}_{1},\ldots,\bm{X}_{k}$ with the length of $p_{1},\ldots,p_{k}$, respectively. Below we let $k=2$ for an illustrative purpose. Let the covariance matrix of $\bm{X}$ be

[TABLE]

where $\Sigma_{\bm{X}_{i}\bm{X}_{j}}$ is the covariance matrix of $\bm{X}_{i}$ and $\bm{X}_{j}$. The null hypothesis of uncorrelation between $\bm{X}_{1}$ and $\bm{X}_{2}$ is

[TABLE]

To test this null hypothesis, we specify $\bm{\Sigma}_{0}$ in the one-sample test statistic (9) as

[TABLE]

where $\bm{J}_{i}$ is a $p_{i}\times p_{i}$ matrix of ones for $i=1,2$, respectively. Using the mapping $\operatorname{vech}(\cdot)$, the Euclidean inner product and norm, it can be shown that test statistic (9) in this case can be expressed as

[TABLE]

where $\bm{S}$ is the sample covariance matrix. Below we outline an associated algorithm to compute the p-value of (29) under the null hypothesis of uncorrelation. Analogous to the identity test, we consider permutations of the elements in each column of the data matrix. In this case the data matrix has two “group-column”’s: the first “group-column” consists of the first $p_{1}$ columns in the data matrix and the second “group-column” includes the remaining $p_{2}$ columns. We then permute the elements in each of the “group-column”’s to compute the p-value under the null, for this permutation procedure retains the sample covariance matrices for each of the two sub-vectors, while disrupting the covariance between the two. In future study we plan to extend the proposed tests and permutation algorithms to handle other covariance structures, for example, those frequently adopted in repeated measure data analysis such as AR(1) and Toeplitz.

7 Acknowledgement

We would like to thank the computing facility SHARCNET (Ontario, Canada) for a partial computation support.

Appendix

Appendix I Derivation of sphericity test statistic (12) and identity test statistic (15)

Let $\bm{S}$ be a $p\times p$ covariance matrix and $\bm{I}$ a $p\times p$ identity matrix. For the sphericity test, $\bm{\Sigma}_{0}$ in the one-sample test statistic (9) is $\sigma^{2}\bm{I}$. We adopt the mapping $\operatorname{vech}(\cdot)$, the Euclidean inner product and norm and apply Equation (1) to compute $cos(\bm{S},\sigma^{2}\bm{I})$ as follows:

[TABLE]

where $s_{ii}^{2}$, $i=1,\ldots,p$, are the diagonal elements of $\bm{S}$. Then the sphericity test statistic (12) follows from $1-cos(\bm{S},\sigma^{2}\bm{I})$.

For the identity test of a correlation matrix $\bm{S}$, the $\bm{\Sigma}_{0}$ in the one-sample test statistic (9) is $\bm{I}$. We adopt the mapping $\operatorname{vech}(\cdot)$, the Euclidean inner product and norm and apply Equation (1) to compute $cos(\bm{S},\bm{I})$, which leads to equation (30). Since $\bm{S}$ is a sample correlation matrix and $\operatorname{tr}(\bm{S})$=$p$, the $cos(\bm{S},\bm{I})$ can be simplified into

[TABLE]

then the identity test statistic (15) follows from $1-cos(\bm{S},\bm{I})$.

Appendix II Derivation of the compound symmetry test statistic (16)

Let $\bm{S}$ be a $p\times p$ sample covariance or correlation matrix and $\bm{J}$ a $p\times p$ matrix of one’s. For the compound symmetry test, $\bm{\Sigma}_{0}$ in the one-sample test statistic (9) is $\sigma^{2}[(1-\rho)\bm{I}+\rho\bm{J}]$ for the covariance matrix and $[(1-\rho)\bm{I}+\rho\bm{J}]$ for the correlation matrix, where $\rho>0$ is the intra-class coefficient coefficient. We first derive the test statistic for the covariance matrix and show that it does not depend on the unknown parameters of $\sigma^{2}$ and $\rho$.

We adopt the modified half-vectorization $\operatorname{vech}^{*}(\cdot)$, the Euclidean inner product and norm to compute $cos(\bm{S},\sigma^{2}[(1-\rho)\bm{I}+\rho\bm{J}])$ according to Equation (1) as follows:

[TABLE]

which does not involve the unknown parameters of $\sigma_{2}$ and $\rho$. The test statistic (16) for the covariance matrix $\bm{S}$ follows from $1-cos(\bm{S},\sigma^{2}[(1-\rho)\bm{I}+\rho\bm{J}])$.

Suppose $\bm{S}$ is a sample correlation matrix, we adopt the modified half-vectorization $\operatorname{vech}^{*}(\cdot)$, the Euclidean inner product and norm to compute $cos(\bm{S},[(1-\rho)\bm{I}+\rho\bm{J}])$ according to Equation (1). Following the above calculation for the sample covariance matrix, we can show that $cos(\bm{S},[(1-\rho)\bm{I}+\rho\bm{J})])$ yields the same equation (32) and the test statistic (16) follows from $1-cos(\bm{S},[(1-\rho)\bm{I}+\rho\bm{J}])$.