Exact Inference for the Dispersion Matrix

We develop a new and novel exact permutation test for prespecified correlation structures such as compound symmetry or spherical structures under standard assumptions. The key feature of the work contained in this note is the distribution free aspect of our procedures that frees us from the standard and sometimes unrealistic multivariate normality constraint commonly needed for other methods.


Introduction
Let (X 1 ,X 2 , …, X n ) be an iid p-dimensional multivariate sample from an absolutely continuous distribution F with p × p dispersion matrix of X as (1) Inference about the dispersion matrix Σ takes the general form (2) where we assume that Σ 0 is specified in a particular manner, for example, a block diagonal matrix or a spherical type structure or simply an unstructured form.
In general, research and testing methods of this form assume an underlying multivariate normal distribution with associated exact and approximate tests; for example, for a thorough overview and history of this testing problem, see Seber [1] and the references there within. In practice one can safely say that it would be rare that the multivariate normality assumption holds. Hence, we were motivated to develop an exact permutation method approach to this problem. To the best of our knowledge no so-called exact permutation tests have been developed or explored with the exception of the very special case of p = 2 dimensions and testing H 0 : ρ 12 = 0; for example, see Good [2]. Martin [3] provides a bootstrap algorithm for testing H 0 : ρ 12 = ρ 12,0 , which asymptotically can be shown to have the appropriate type I error rate. Unfortunately, the bootstrap methods given by Martin [3] relative to first standardizing the variables and rotating the data so as to transform the problem to the setting of testing H 0 : ρ 12 = 0 do not work in the permutation setting. The permutation test for the case H 0 : ρ 12 = 0 follows by permuting the second column of the n × 2 data matrix (X 1 ,X 2 ) and calculating the test statistic ρ1 2 , where ρ1 2 refers to the standard sample Pearson correlation coefficient, over all permutations. This can be done directly via a computationally expensive algorithm or via the more widely used Monte Carlo techniques. With respect to the Monte Carlo methods we generate B random permutations of the data and denote the permuted value of the test statistic by . Then the one-sided P value for the alternative H 1 : ρ 12 > 0 is given as , where the index i corresponds to a given permutation and I denotes the indicator function. Alternative approaches found in software packages such as SAS PROC FREQ (SAS version 9.3, Cary, NC) utilize hypergeometric probabilities similar to how Fisher's exact test is carried out via treating the fixed data as discrete.
In general permutation testing is most often used for comparing two groups in the context of location differences or other features of distributions such as scale measures. Most of the theoretical work has been done in this setting such as type I error control. For a technical treatment of permutation testing see Romano [4] with respect to a theoretical examination for the behavior of the type I error control for permutation tests under exchangeability versus nonexchangeability conditions. In order to ensure true bounded type I error control in the permutation testing setting either the null hypothesis has to be specified in such a way that exchangeability holds by definition under H 0 or some design feature such as randomization or matching needs to be employed. Commenges [5] studies the more general transformation approach used to preserve exchangeability. Also, see Zhang [6], Huang et al. [7], and Janssen and Pauls [8] with respect to the inflation of type I error rates when comparing means in the two-sample setting along with other types of comparisons. In terms of permutation testing related to correlation structure based hypotheses very little has been accomplished. This paper represents one of the only investigations of this type to date.
In Section 2 we develop the general p-dimensional exact tests given a prespecified covariance structure. Special cases include testing for sphericity, compound symmetry, and block diagonality, to name a few. This presentation is followed by a simulation study in Section 3. We then apply our method in Section 4 to an example involving repeated measures mice weight data.

Exact Tests for Covariance Structures in p Dimensions
The focus of the work in this setting is with respect to two-sided alternatives. In certain instances a subset of these tests with one-sided alternative structure may be constructed. Those tests will not be included as part of this discussion due to the specificity of their applications.

Unequal Variance Setting
Now let (X 1 ,X 2 , …, X n ) be an iid p-dimensional multivariate sample from an absolutely continuous distribution F with the first two finite central moments corresponding to each component of X given as E(X i ) = μ x i and , i = 1, 2, …, p. Let Corr(X i , X j ) = ρ ij , i = 1, 2, …, p, j = 1, 2, …, p, (i ≠ j). Furthermore, denote the p × p dispersion matrix of X by (3) where Σ is defined to be a p × p positive definite matrix. We represent the Cholesky decomposition of the p × p matrix Σ as (4) such that A′ ′1 is defined. The Cholesky decomposition is a key component of the permutation test we propose; however it is not unique to the problem; that is, other decomposition methods may yield similar results and alternative solutions. From a practical standpoint the Cholesky decomposition is built in to several statistical softwarepackages, thus making our methodology more feasible for a larger group of practitioners. Now our more general hypothesis of interest takes the form (5) where Σ 0 are the hypothesized value of Σ at (3).
Test Statistic-Let the p × p matrix denote the transpose of A′ −1 with the hypothesized values as given elements from (5). Let the n × p matrix denote the data matrix following transformation. Then the dispersion matrix corresponding to the n × p matrix Z will be a diagonal matrix such that Corr(Z i , X j ) = 0 ∀i < j, j < i if and only if H 0 at (5) holds true. Under these conditions testing H 0 at (5) is equivalent to the test: (6) where the off-diagonal elements of Σ z 0 are equal to 0 under H 0 at (6).
An exact α-level permutation test of H 0 can be defined for (6) by considering the permutation of each column of Z and employing the Pearson correlation coefficient for each combination of columns. Towards this end let us denote the Corr(Z i , Z j ) = ρ z ij with the corresponding Pearson estimator by ρẑ ij . The test statistic of interest in the two-sided case with respect to detecting departures from H 0 at (6)

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript (7) The exactness of the test in terms of the type I error control follows from a straightforward generalization of the form of the dispersion matrix for the 2 × 2 case, where (8) In the 2 × 2 case the off-diagonal elements of the dispersion matrix are given as (9) An examination of the covariance term corresponding to Σ Z|H 1 at (9) clearly indicates that it has the value of 0 if and only if H 0 at (8) is true. When testing H 0 at (8) the Pearson correlation estimate between the transformed variates through , Z 1 , and Z 2 serves to appropriately detect departures from H 0 . Within the permutation testing framework provides an exact α-level test; that is, the covariance of Z 1 and Z 2 is 0 if and only if the correlation of Z 1 and Z 2 is 0.
We resort to a Monte Carlo approximation in order to obtain the P value for testing H 0 defined at (5). The steps for performing the Monte Carlo approximation with respect to estimating the P value are as follows.

2.
Obtain the new random variates Z by applying the transformation to the observed data X. (7).

4.
Permute each column of Z independently such that we have the permuted n × p matrix denoted by Z * .

7.
Calculate the Monte Carlo estimated permutation P value as , where I (·) denotes the indicator function.

Nontransformation Special Cases
In certain special cases we can test specific forms of hypothesis (5) using our permutation approach without specifying a specific subset of the 's or ρ ij 's. One obvious special case  [9] for a description of a likelihood ratio approximation to this test. In this instance we have (10) with unspecified 's under H 0 .
In this case there is no transformation of the data required. An exact α-level permutation test of H 0 can be defined simply by considering the permutation of each column of X and employing the Pearson correlation coefficient for each combination of columns. Towards this end denote the Corr(X i , X j ) = ρ ij with the corresponding Pearson estimator by ρx ij . The test statistic of interest with respect to detecting departures from the diagonal structure is defined as (11) The Monte Carlo estimated permutation P value is calculated similarly as before, where , where I (·) denotes the indicator function.
Another special case where we can have a set of unspecified 's or ρ ij 's is when we may be interested in testing for a block diagonal dispersion matrix structure such that under H 0 at (5) we now have (12) where the partitioned q i × q j matrices are given as (13) The dispersion matrix Σ jj may have different dimensions q j < p(Σq i = p), j = 1, 2, …, b, with unspecified 's and ρ ij 's under H 0 . As in the test for a diagonal dispersion matrix above there is no transformation of the data required. An exact α-level permutation test of H 0 can again be defined simply by considering the permutation of each column of X and employing the Pearson correlation coefficient for each combination of columns. The test statistic of interest with respect to detecting departures from the diagonal structure is a slight modification of the test statistic at (11) defined as (14) where the "off-block" correlation elements at (12), ρ 0,ij = 0, under H 0 and I (·) denote the indicator function. The Monte Carlo estimated permutation p value is calculated similarly as before, where and I (·) denotes the indicator function.

Equal Variance Setting
For the equal variance p-dimensional case we have (15) where we now define the p × p dispersion matrix under H 0 as (16) where Σ 0 is defined to be a p × p positive definite matrix, Var(X i ) = σ 2 , i = 1, 2, …, p, under H 0 and the p × p correlation matrix Γ 0 is given by (17) The Cholesky decomposition of matrix Γ 0 is as (18) In order to test the hypothesis at (15) we utilize the transformation . Note that H 0 at (15) will be sensitive to departures from Γ 0 and unequal marginal variances. Furthermore, explicit values for σ 2 do not need to be specified within this hypothesis testing framework.
The steps for performing the Monte Carlo approximation with respect to estimating the p value are as follows.

2.
Obtain the new random variates Z by applying the transformation to the observed data X.

4.
Permute each column of Z independently such that we have the permuted n × p matrix denoted by Z * .

6.
Repeat steps (4) and (5)  Several other well-known spatial dispersion matrices similar to the spatial power matrix presented above fit within this same framework and will not be presented here.

Simulation Study
In this section we examine the test at (19), where we specify ρ 0 and the form of the correlation structure Γ at (17), for example, compound symmetry. Our simulation study for the p × p case will utilize a p-variate standardized multivariate normal distribution with p = 5 and a special case mixing the marginal distributions across normal, exponential, and uniform forms. Again, differing location and scale doe not vary the general conclusions. In terms of our simulation study we set the null value of ρ 0 = 0, 0.5, 0.9 under a compound symmetry assumption and a first-order autoregressive assumption, where Γ will take the forms: 1. compound symmetry: 2. first-order autoregressive: Note that the special case ρ 0 = 0 is the same for both covariance structures and is only presented once. It should also be noted that under the assumption of multivariate normality testing H 0 : Σ = σΓ 0 (0) under the compound symmetry or first-order autoregressive structure is a special case (equal variance assumption) of the well-known test for "complete independence"; for example, see Mudholkar et al. [9]. Under nonnormality we are essentially testing the "complete uncorrelated" case. In this special case the methods presented here are the first exact methods developed for tackling this particular hypothesis.
In terms of large sample theory around similar results see Jiang [10] and Xiao and Wu [11].
For our simulation study we used 1000 replicates for our study at n = 10, 20, 30, 40 and set α = 0.05. The covariance structure was the same under H 0 and H 1 for this set of simulations.
The results are contained in Figures 1, 2, 3, 4, and 5. As anticipated we see the expected results of appropriate type I error control and monotone power functions increasing in either direction about the null value for ρ. The range of ρ under the alternative was dictated by the constraint that Γ 0 (ρ 0 ) is defined to be positive definite.
For the sake of example we modified our simulation and took ρ 0 = 0.5 with marginals given by X 1 , X 2 ~ N(0, 1), , and X 5 exp(1, −1) with Γ 0 (ρ 0 ) assumed to be compound symmetric under H 0 and under H 1 . The results are shown in Figure 6 and as we can see they do not differ dramatically from Figure 2 assuming multivariate normality, thus illustrating the flexibility and nonparametric nature of our methodology.
As an additional result we studied the power under the correctly specified ρ under H 0 with Γ 0 (ρ 0 ) differing in structure. For this study we set ρ 0 = 0.5, 0.9 with Γ 0 (ρ 0 ) set to compound symmetry under H 0 and Γ 0 (ρ 0 ) set to the first-order autoregressive structure under H 1 . In other words what is the power to detect a correlation structure different from the null structure given that ρ 0 is the true correlation. At α = 0.05 the power to detect a different correlation structure under the alternative for ρ 0 = 0.5 and 0.9 at n = 10, 20, 30, 40 and α = 0.05 was 0.245, 0.539, 0.761, and 0.894 and 0.520, 0.858, 0.951, and 0.998, respectively.

Example
As an illustration of our method we will use phenotypic weight data from n = 16 mice as contained in Table 1 from a recent unpublished study conducted within Roswell Park Cancer Institute. The estimated correlation matrix is provided in Table 2. The respective sample variance estimates were , and . For example, suppose we were interested in testing (26) for both Γ(·) having the compound symmetry structure or the first-order autoregressive structure as defined in (23). In this instance the test corresponding to the above hypothesis under a compound symmetry correlation structure yielded a Monte Carlo estimated P value <0.0001 (B = 10, 000). While the test corresponding to the above hypothesis under the firstorder autoregressive correlation structure yielded a Monte Carlo estimated P value = 0.002 (B = 10, 000). For this example this provides some measure of evidence that the correlation structure does not fit the compound symmetry structure and that the first-order autoregressive structure assuming ρ = 0.95 may be more appropriate. Similarly, the test for diagonality under the equal variances assumption (sphericity), which does not assume a value for ρ 0 , yielded a MonteCarlo estimated P value <0.0001 (B = 10, 000). Note that we may be rejecting H 0 at some specified level α under at least one of 3 scenarios: unequal marginal variances, ρ ≠ ρ 0 or Γ ≠ Γ 0 .
Given our overall P value from above which was 0.002, we may wish to examine in further detail what is driving us to reject H 0 . In this case we can examine specific submatrices of the dispersion matrix of . For this example we could test H 0 : Σ = σΓ 0 (.95) using days 0, 1, and 2, only or any other combinations of days such that the appropriate correlation substructure is extracted from the original hypothesized values for Γ. For our example subtest we get P = 0.32 indicating no strong evidence against a first order autoregressive "substructure" with equal variances and ρ 0 = 0.95. If we add day 3, our P value=0.04, indicating either the correlation structure may be misspecified at this point or the variance is different. Note that further work relative to the multiple comparison problem of subtests and their relative correlation is needed. This is simply an exploratory approach to this issue relative to the example at hand.

Concluding Remarks
In this note we provided a method for exact testing around specific covariance structures. We employed the Cholesky decomposition for this purpose. It was noted by a reviewer that other decomposition methodologies may lead to extensions of this methodology, which we will consider in terms of future work.

Author Manuscript
Hutson et al. Page 18 Table 1 Mice weights (grams).  Table 2 Estimated correlation matrix for mouse weight example data.