1
$\begingroup$

Background. I need to compare two correlation matrices: $R_1$ is observed matrix and $R_2$ is expected matrix. The original question is here. In the paper (Modarres and Jernigan, 1993) two different tests were proposed. The quadratic form test is more simple one.

A QUADRATIC FORM TEST

Let $Vech(R_i)$ denote the sample correlation vector obtained from the $i$-th sample correlation matrix, $R_i$, by placing its subdiagonal elements underneath each other, by columns. Let $P_i$ denote the $i$-th population correlation matrix. It is well known (Steiger and Hakstian, 1982) that $$\sqrt{n_i} Vech(R_i - P_i)$$ is asymptotically multivariate normal with covariance matrix $\Gamma$ whose elements depend on the fourth order moments.

For testing the null hypothesis $H_0: P_1 = P_2 = . . . = P_k = P$, where $P$ is unspecified, against the general alternative $H_1: P_i \neq P_j$, for at least one $i \neq j$, the $k$-sample test statistic is given by

$$Q_m=\sum_{\alpha=1}^k n_\alpha Vech(R_\alpha-\bar{R})'\hat{\Gamma}^{-1}Vech(R_\alpha-\bar{R})$$ where $\hat{\Gamma}$ is any consistent estimator of $\Gamma$ and $\bar{R}$ is the weighted average of the correlation matrices. This test has an asymptotic $\chi^2$ distribution with $(k - 1)p(p- 1)/2$ degrees of freedom under the null hypothesis.

In the two sample case, the test statistics $Q$ is given by $$Q=\frac{n_1n_2}{n_1+n_2}Vech(R_1-R_2)'\hat{\Gamma}^{-1}Vech(R_1-R_2).$$ To estimate $\Gamma$ we use a pooled estimate of the covariance of correlations from each population. That is, $\hat\Gamma = [(n_1 - 1)\hat{\Gamma}_1 + (n_2 - 1)\hat{\Gamma}_2]/(n_1 + n_2 - 2)$ where $\hat{\Gamma}_1$, and $\hat{\Gamma}_2$ are obtained by evaluating the Steiger-Hakstian expression at $R_1$ and $R_2$, respectively.

I don't find the $R$ package with realization of the quadratic form test. My code is below.

Edit 1. I have increased the value of n to 10000 in the corSample function.

library(fungible)
set.seed(40)
n <- 4
k <- 2 # k-sample test statistic
p <- n # number of variables
rand_vec <- rnorm(n);rand_vec
n1 <- 100; n2 <- 100

# correlation matrix on empirical data
R1 <- matrix(c(
1.00, 0.51, 0.44, 0.22,
0.51, 1.00, 0.36, 0.21,
0.44, 0.36, 1.00, 0.26,
0.22, 0.21, 0.26, 1.00), n, n)

# moments of empirical data
skew_vec = c(-0.254, -0.083, 0.443, -0.017) 
kurt_vec = c(6.133,   4.709, 6.619,  4.276)
# random terms
rand_vec <- rnorm(n)
X1 <- monte1(seed = 123, nvar = n, nsub = n1, cormat = R1,
skewvec = skew_vec, 
kurtvec = kurt_vec)$data #; X1

# generate a sample correlation from population 'R1' with n = 10000
R2 <- corSample(R1, n = 10000)$cor.sample

X2 <- monte1(seed = 1234, nvar = n, nsub = n2, cormat = R2, 
            skewvec = skew_vec + rand_vec,
            kurtvec = kurt_vec + rand_vec)$data#; X2

# vectorization operator 
delta <- row(R1) - col(R2)
vR1 <- as.vector(t(R1[delta > 0])) 
vR2 <- as.vector(t(R2[delta > 0]))

# Create ADF Covariance Matrix of Correlations by Steiger-Hakstian expression
G1 <- adfCor(X1)#; G1
G2 <- adfCor(X2)#; G2

G    <- ((n1 - 1)*G1 + (n2 - 1)*G2)/(n1 + n2 - 2)
Ginv <- MASS::ginv(G)#; round(Ginv%*%G, 3)

Q    <- n1*n2/(n1 + n2) * ((vR1 - vR2) %*% Ginv) %*% (vR1 - vR2); Q  
# 3.797144

The hypothesis $H_1: P_1 \neq P_2$ is the complex hypothesis.

# This test has an asymptotic x2 distribution with (k - 1)p(p- 1)/2 degrees # of freedom under the null hypothesis.

alpha = 0.05

left  <- qchisq(alpha/2, df=(k - 1)*p*(p-1)/2); left
# 1.237344
right <- qchisq(1- alpha/2, df=(k - 1)*p*(p-1)/2); right
# 14.44938

You can see that $\chi^2_{\alpha=0.025, df=6}=1.237344<Q=3.797144 <\chi^2_{\alpha=0.975, df=6}=14.44938$. And I should not to reject the null hypothesis.

Question. Is there an alternative way to verify correctness of the calculations?

References.

Reza Modarres & Robert W. Jernigan (1993) A robust test for comparing correlation matrices, Journal of Statistical Computation and Simulation, 46:3-4, 169-181.

Steiger. J. H. and Hakstian, A. R. (1982). 'The asymptotic distribution of elements of a correlation matrix: Theory and application.' British Journal of Mathematical and Statistical Psychology 35, 208-215.

$\endgroup$
  • $\begingroup$ One (natural) answer to your question is "yes: apply the calculations to simulated data." $\endgroup$ – whuber Sep 2 '18 at 14:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.