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.
