# Permutation test for two multivariate distributions with the same mean, variance, but different covariance/correlation

I am using bivariate normal distributions as an example.

Using independence_test() in R, I can tell if two bivarate distributions have different means. The null hypothesis that their means are the same is rejected in the example below.

library(coin)
library(ggplot2)
library(MASS)

# Generate a vector of two correlated normals.
make_input_data <- function(mu1, mu2, mu3, mu4,
Sigma1, Sigma2,
numSamples) {
out1 <-mvrnorm(n=numSamples,
mu=c(mu1, mu2),
Sigma=Sigma1)
df1 <- data.frame(y1=out1[,1], y2=out1[,2], label='A')
out2 <-mvrnorm(n=numSamples,
mu=c(mu3, mu4),
Sigma=Sigma2)
df2 <- data.frame(y1=out2[,1], y2=out2[,2], label='B')
df <- rbind(df1, df2)
return (df)
}

numSimulations <- 1e5
numSamples <- 1e3
df <- make_input_data(mu1=0.0, mu2=0.0, mu3=1.0, mu4=0.0,
Sigma1=matrix(c(1, 0, 0, 1), ncol = 2),
Sigma2=matrix(c(1, 0, 0, 1), ncol = 2),
numSamples=numSamples)
qplot(y1, y2, data=df, color=label)

# Independence test.
independence_test(y1 + y2 ~ as.factor(label),
distribution=approximate(B=numSimulations-1),
data=df)
---

Approximative General Independence Test

data:  y1, y2 by as.factor(label) (A, B)
maxT = 21.185, p-value < 2.2e-16
alternative hypothesis: two.sided


However, I don't know how to construct a multivariate test to tell if two bivariate distributions have different covariance/correlations.

df <- make_input_data(mu1=0.0, mu2=0.0, mu3=0.0, mu4=0.0,
Sigma1=matrix(c(1, 0.6, 0.6, 1), ncol = 2),
Sigma2=matrix(c(1, -0.6, -0.6, 1), ncol = 2),
numSamples=numSamples)
qplot(y1, y2, data=df, color=label)


The best I can do so far is to propose the following: We run the following 3 univariate tests. Two distributions are deemed the same if all tests pass, because their first and second moments are the same.

• We run a univariate "location" test (e.g., Fisher's test) which tells us if the mean in each dimension $y_i$ is the same.
• Similarly, we run a univariate "scale" test (e.g., Mood test) to see if the variance in each dimension $y_i$ is the same.
• We can compute the correlation of each pair $\text{corr}(y_i, y_j)$ along with a confidence interval. The function cor.test() in R does this.

There are various questions here which are related but which can be answered almost independently:

### Testing for correlation differences

The function cor.test() in base R tests whether there is significant correlation ($H_0: r(y_1, y_2) = 0$) between two variables. What you want to do is to test whether there are correlation differences between two groups ($H_0: r_A(y_1, y_2) = r_b(y_1, y_2)$).

One simple idea for this is to use some transformation of the data so that mean differences on that transformed scale correspond to correlation differences in the original data. A straightforward transformation would be $r_i = (y_{1, i} - \bar y_1)/s_1 \cdot (y_{2, i} - \bar y_2)/s_2$ because up to degrees-of-freedom adjustments the correlation corresponds to the mean of $r_i$.

Using your data-generating process and a standard parametric 2-sample $t$-test we can do:

set.seed(1)
df <- make_input_data(0, 0, 0, 0, matrix(c(1, 0.6, 0.6, 1), ncol = 2),
matrix(c(1, -0.6, -0.6, 1), ncol = 2), 1000)

### Roll your own independence_test interface

Instead of computing these transformations "by hand" and then running a test for location differences on them, you can also define your own transformation function and pass that along to independence_test. This is what all the interfaces like wilcox_test or mood_test in coin do. The underlying ideas are shown in vignette("LegoCondInf", package = "coin") while the computational tools are discussed in more detail in vignette("coin_implementation", package = "coin").