1
$\begingroup$

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)

enter image description here

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.
$\endgroup$
1
$\begingroup$

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)
df$r <- (df$y1 - mean(df$y1))/sd(df$y1) * (df$y2 - mean(df$y2))/sd(df$y2)
t.test(r ~ label, data = df)
##  Welch Two Sample t-test
## 
## data:  r by label
## t = 23.31, df = 1996.8, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.091220 1.291708
## sample estimates:
## mean in group A mean in group B 
##       0.5983141      -0.5931498 

Note that the means in the two groups almost exactly recover the true simulated correlations. Not surprisingly, the result is highly signiciant.

Rank-based tests

If you want to use a robust version of this test you could compute r based on rank(y1) and rank(y2). This would then correspond to a Spearman correlation (up to degrees-of-freedom adjustments).

Permutation tests

If you don't trust the central limit theorem in this test you can, of course, approximate the corresponding permutation test using coin.

set.seed(1)
independence_test(r ~ label, data = df, distribution = approximate(9999))
##  Approximative General Independence Test
## 
## data:  r by label (A, B)
## Z = 20.673, p-value < 2.2e-16
## alternative hypothesis: two.sided

However, in this situation the normal approximation appears to be working perfectly well.

Testing for scale differences

To test for scale differences you can, in principle, also compute analogous transformations, e.g., $s_{1, i} = (y_{1, i} - \bar y_1)^2$, so that the corresponding mean yields the variance (up to degrees-of-freedom adjustments). Alternatively, one could again use a rank-based version of this transformation which is what the Mood test does.

df$s1 <- (df$y1 - mean(df$y1))^2
df$s2 <- (df$y2 - mean(df$y2))^2
set.seed(1)
independence_test(s1 ~ label, data = df, distribution = approximate(9999))
##  Approximative General Independence Test
## 
## data:  s1 by label (A, B)
## Z = -0.11806, p-value = 0.9121
## alternative hypothesis: two.sided

And analogously for s2.

Combining tests for location, scale, and correlation differences

If you want to test for differences in all five directions (two locations, two scales, one correlation), one simple approach would be to conduct each test and reject if the maximum $t$ or $z$ statistic is large. The $p$-values or critical values should be Bonferroni-adjusted to account for five separate tests being conducted.

You can do this "by hand" using any tests of your choice. In coin you can also conduct these in one go using the default teststat = "maximum":

set.seed(1)
it <- independence_test(y1 + y2 + s1 + s2 + r ~ label, data = df,
  distribution = approximate(9999))
it
##  Approximative General Independence Test
## 
## data:  y1, y2, s1, s2, r by label (A, B)
## maxT = 20.673, p-value < 2.2e-16
## alternative hypothesis: two.sided

Note that this yields the correlation statistic from above as the maximum test statistic. You have also access to the five separate $z$ statistics:

statistic(it, "standardized")
##          y1         y2         s1        s2        r
## A 0.3894535 -0.5173593 -0.1180622 0.3051142 20.67332

which clearly shows that the correlation differences are the reason for rejection here.

Instead of the teststat = "maximum" you can also use a teststat = "quadratic" which corresponds to a $\chi^2$ type test statistic. The former will typically have higher power if only a single of the five parameters differs between the groups. The latter will typically perform better if several parameters change.

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").

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.