An exact test for a change in covariance structure

The problem

Did the dependence structure change between two periods? The likelihood-ratio statistic is Box’s M, and its null distribution has traditionally been handled by asymptotic approximations that are poor exactly where the question is interesting: small samples and moderate dimension.

set.seed(1)
p <- 3
XA <- matrix(rnorm(60 * p), 60, p)
XB <- matrix(rnorm(50 * p), 50, p) %*% diag(c(1, 1.6, 0.7))
cov_M(XA, XB)
#> [1] 35.19953

The null is covariance-free

The statistic vanishes when nothing differs, and it is invariant to a common change of basis:

cov_M(XA, XA)
#> [1] 0
A <- matrix(c(2, 0.3, -1, 0.1, 1.4, 0.2, 0, 0.5, 1), 3, 3)
c(original = cov_M(XA, XB), transformed = cov_M(XA %*% A, XB %*% A))
#>    original transformed 
#>    35.19953    35.19953

That invariance is the point. Because the null law does not depend on the unknown common covariance, it can be calibrated once at the identity and used at every covariance, with no nuisance parameter to estimate.

An exact p-value

cov_pexact evaluates the null tail by inverting the exact characteristic function, so the test is exact rather than approximate. It agrees with simulation from the covariance-free null, which is an independent route:

set.seed(3)
null <- cov_null(p = 3, nA = 40, nB = 40, B = 20000, seed = 5L)
q <- unname(quantile(null, c(0.5, 0.9)))
rbind(simulated_tail = c(0.5, 0.1),
      exact_tail = vapply(q, function(m) cov_pexact(m, 39, 39, 3), 0))
#>                     [,1]       [,2]
#> simulated_tail 0.5000000 0.10000000
#> exact_tail     0.5013875 0.09727948

Applied to the data above, the scale change in one coordinate is detected:

cov_test(XA, XB, method = "exact")$p.value
#> [1] 6.38278e-06

and two samples from the same law are not:

set.seed(11)
cov_test(matrix(rnorm(180), 60, 3), matrix(rnorm(180), 60, 3), method = "exact")$p.value
#> [1] 0.8665873

Differential networks

The same machinery gives a log-domain differential network, \(D = \log \mathrm{cov}(X_B) - \log \mathrm{cov}(X_A)\). It is symmetric, antisymmetric under swapping the samples, and inversion-invariant, so it says the same thing whether you think in covariances or precisions:

D <- cov_logdiff(XA, XB)
round(D, 4)
#>         [,1]   [,2]    [,3]
#> [1,]  0.3851 0.0774 -0.0867
#> [2,]  0.0774 1.2499  0.1094
#> [3,] -0.0867 0.1094 -0.9849
c(antisymmetry = max(abs(D + cov_logdiff(XB, XA))), symmetry = max(abs(D - t(D))))
#> antisymmetry     symmetry 
#> 0.000000e+00 1.387779e-17

A max-type statistic in the style of Cai, Liu and Xia is also provided for the sparse alternative:

clx_stat(XA, XB)
#> [1] 11.64884

Implementation

Everything above runs in a shared C back-end with its own random number generator, so results do not depend on R’s RNG state, and the same sources are bound from Python.