Data depth and rank-based tests for HPD matrices

Introduction

Nondegenerate covariance, correlation and spectral density matrices are necessarily sym- metric or Hermitian and positive definite. In (Chau, Ombao, and Sachs 2019), we develop statistical data depths for collections of Hermitian positive definite matrices by exploiting the geometric structure of the space as a Riemannian manifold. Data depth is an important tool in statistical data analysis measuring the depth of a point with respect to a data cloud or probability distribution. In this way, data depth provides a center-to-outward ordering of multivariate data observations, generalizing the notion of a rank for univariate observations.

The proposed data depth measures can be used to characterize most central regions or detect outlying observations in samples of HPD matrices, such as collections of covariance or spectral density matrices. The depth functions also provide a practical framework to perform rank-based hypothesis testing for samples of HPD matrices by replacing the usual ranks by their depth-induced counterparts. Other applications of data depth include the construction of confidence regions, clustering, or classification for samples of HPD matrices.

In this vignette we demonstrate the use of the functions pdDepth() and pdRankTests() to compute data depth values of HPD matrix-valued observations and perform rank-based hypothesis testing for samples of HPD matrices, where the space of HPD matrices can be equipped with several different metrics, such as the affine-invariant Riemannian metric as discussed in (Chau, Ombao, and Sachs 2019) or Chapter 4 of (Chau 2018).

Data depth of HPD matrices with pdDepth()

First, we generate a pointwise random sample of (2,2)-dimensional HPD matrix-valued observations using the exponential map Expm(), with underlying intrinsic (i.e., Karcher or Fréchet) mean equal to the identity matrix diag(2). Second, we generate a random sample of sequences (discretized curves) of (2,2)-dimensional HPD matrix-valued observations, with underlying intrinsic mean curve equal to an array of rescaled identity matrices. For instance, we can think of the first sample as a random collection of HPD covariance matrices, and the second sample as a random collection of HPD spectral matrix curves in the frequency domain.

## Pointwise random sample
library(pdSpecEst); set.seed(100)
X1 <- replicate(50, Expm(diag(2), H.coeff(0.5 * rnorm(4), inverse = T))); str(X1)
#>  cplx [1:2, 1:2, 1:50] 0.7794+0i -0.0314-0.0523i -0.0314+0.0523i ...

## Curve random sample
X2 <- replicate(50, sapply(1:5, function(i) Expm(i * diag(2), H.coeff(0.5 *
                           rnorm(4), inverse = T) / i), simplify = "array")); str(X2)
#>  cplx [1:2, 1:2, 1:5, 1:50] 1.074+0i 0.352+0.147i 0.352-0.147i ...

The function H.coeff() expands an Hermitian matrix with respect to an orthonormal basis of the space of Hermitian matrices equipped with the Frobenius (Euclidean) inner product \(\langle .,. \rangle_F\). By specifying the argument inverse = T, the function computes an inverse basis expansion, transforming a real-valued basis component vector back to an Hermitian matrix.

The function pdDepth() computes the intrinsic data depth of a single HPD matrix (resp. curve of HPD matrices) y with respect to a sample of HPD matrices (resp. sample of curves of HPD matrices) X. The intrinsic data depth is calculated in the space of HPD matrices equipped with one of the following metrics: (i) affine-invariant Riemannian metric (metric = "Riemannian"), (ii) log-Euclidean metric, the Euclidean inner product between matrix logarithms (metric = "logEuclidean"), (iii) Cholesky metric, the Euclidean inner product between Cholesky decompositions (metric = "Cholesky"), (iv) Euclidean metric (metric = "Euclidean") and (v) root-Euclidean metric, the Euclidean inner product between Hermitian square root matrices (metric = "rootEuclidean"). Note that the default choice (affine-invariant Riemannian) has several useful properties not shared by the other metrics. See (Chau, Ombao, and Sachs 2019) and Chapter 4 of (Chau 2018) for more details and additional properties of the available intrinsic depth functions.

Remark: an advantage of substituting, for instance, the Log-Euclidean metric is that the depth computation times may be significantly faster, in particular for large sample sizes. However, the data depths with respect to the Log-Euclidean metric are not general linear congruence invariant according to Chapter 4 of (Chau 2018), and should therefore be used with caution.

## Pointwise geodesic distance, zonoid and spatial depth
point.depth <- function(method) pdDepth(y = diag(2), X = X1, method = method)
point.depth("gdd"); point.depth("zonoid"); point.depth("spatial")
#> [1] 0.4326915
#> [1] 0.7600822
#> [1] 0.7932682

## Integrated geodesic distance, zonoid and spatial depth
int.depth <- function(method){ pdDepth(y = sapply(1:5, function(i) i * diag(2),
                                    simplify = "array"), X = X2, method = method) }
int.depth("gdd"); int.depth("zonoid"); int.depth("spatial")
#> [1] 0.751167
#> [1] 0.8631369
#> [1] 0.8825155

By leaving the argument y in the function pdDepth() unspecified, the function computes the data depth of each individual object in the sample array X with respect to the sample X itself.

(dd1 <- pdDepth(X = X1, method = "gdd")) ## pointwise geodesic distance depth
#>  [1] 0.4136872 0.4072360 0.4025387 0.4044413 0.2559378 0.3873209 0.3832640
#>  [8] 0.3002407 0.4034977 0.3358114 0.2597040 0.3120139 0.1967228 0.1876342
#> [15] 0.1922631 0.2483946 0.3696490 0.3386755 0.2068996 0.2058298 0.2020919
#> [22] 0.3757002 0.2523476 0.2142259 0.2864097 0.3353675 0.3192103 0.3354539
#> [29] 0.3210167 0.3004936 0.3990564 0.3107404 0.3665464 0.2987992 0.4306943
#> [36] 0.2448184 0.3170413 0.2921210 0.4100427 0.4292756 0.4257262 0.3645616
#> [43] 0.3668837 0.2471900 0.2340217 0.3391741 0.3632063 0.3623502 0.2293419
#> [50] 0.2719987

(dd2 <- pdDepth(X = X2, method = "gdd")) ## integrated geodesic distance depth
#>  [1] 0.7001805 0.6885243 0.5746611 0.7066090 0.6931021 0.6423633 0.6738678
#>  [8] 0.6197673 0.5654389 0.7123940 0.6469901 0.7158848 0.6487223 0.6243373
#> [15] 0.7075508 0.6356606 0.6654029 0.7162644 0.6380360 0.7091247 0.6595957
#> [22] 0.6520224 0.6627836 0.6932087 0.6783008 0.6798760 0.6490358 0.7024616
#> [29] 0.6324775 0.6889687 0.6728189 0.6927310 0.6464986 0.6492694 0.6591718
#> [36] 0.6874394 0.6743475 0.6365783 0.6644634 0.7112869 0.6935939 0.7361826
#> [43] 0.7209363 0.6161036 0.6952725 0.6470956 0.6250286 0.7059682 0.7125312
#> [50] 0.7113942

A center-to-outward ordering of the individual objects in the sample of (curves of) HPD matrices is then obtained by computing the data depth induced ranks, with the most central observation (i.e., highest depth value) having smallest rank and the most outlying observation (i.e., lowest depth value) having largest rank.

(dd1.ranks <- rank(1 - dd1)) ## pointwise depth ranks
#>  [1]  4  6  9  7 37 11 12 31  8 22 36 28 48 50 49 39 14 21 45 46 47 13 38
#> [24] 44 34 24 26 23 25 30 10 29 16 32  1 41 27 33  5  2  3 17 15 40 42 20
#> [47] 18 19 43 35

(dd2.ranks <- rank(1 - dd2)) ## integrated depth ranks
#>  [1] 14 21 49 11 18 40 26 47 50  6 38  4 36 46 10 43 28  3 41  9 31 33 30
#> [24] 17 24 23 35 13 44 20 27 19 39 34 32 22 25 42 29  8 16  1  2 48 15 37
#> [47] 45 12  5  7

## Explore sample X1
head(order(dd1.ranks)) ## most central observations
#> [1] 35 40 41  1 39  2
rev(tail(order(dd1.ranks))) ## most outlying observations
#> [1] 14 15 13 21 20 19
X1[ , , which(dd1.ranks == 1)] ## most central HPD matrix
#>                       [,1]                  [,2]
#> [1,]  0.9407902+0.0000000i -0.0154483+0.1752587i
#> [2,] -0.0154483-0.1752587i  1.2710700+0.0000000i
X1[ , , which(dd1.ranks == 50)] ## most outlying HPD matrix
#>                     [,1]                [,2]
#> [1,]  1.847918+0.000000i -1.288255+1.075924i
#> [2,] -1.288255-1.075924i  2.490724+0.000000i

We can compare the most central HPD matrix above with the intrinsic median of the observations under the affine-invariant metric obtained with pdMedian(). Note that the intrinsic sample median maximizes the geodesic distance depth in a given sample of HPD matrices. The point of maximum depth (i.e., deepest point) with respect to the intrinsic zonoid depth is the intrinsic sample mean, which can be obtained with pdMean(). For more details, see (Chau, Ombao, and Sachs 2019) or Chapter 4 of (Chau 2018).

(med.X1 <- pdMedian(X1)) ## intrinsic sample median
#>                         [,1]                    [,2]
#> [1,]  0.90753401+0.00000000i -0.03426934+0.03226425i
#> [2,] -0.03426934-0.03226425i  1.15446038-0.00000000i


pdDepth(y = med.X1, X = X1, method = "gdd") ## maximum out-of-sample depth
#> [1] 0.4409941

max(dd1) ## maximum in-sample depth
#> [1] 0.4306943

Rank-based tests for HPD matrices with pdRankTests()

The null hypotheses of the available rank-based hypothesis tests in the function pdRankTests() are specified by the argument test and can be one of the following:

• "rank.sum": homogeneity of distributions of two independent samples of HPD matrices (resp. sequences of HPD matrices).
• "krusk.wall": homogeneity of distributions of more than two independent samples of HPD matrices (resp. sequences of HPD matrices).
• "signed-rank": homogeneity of distributions of independent paired or matched samples of HPD matrices.
• "bartels": exchangeability (i.e. randomness) within a single independent sample of HPD matrices (resp. sequences of HPD matrices).

Below, we construct several simulated examples for which: (i) the null hypotheses listed above are satisfied; and (ii) the null hypotheses listed above are not satisfied. Analogous to the previous section, we generate pointwise random samples (resp. random samples of sequences) of (2,2)-dimensional HPD matrix-valued observations, with underlying geometric mean equal to the identity matrix (resp. sequence of scaled identity matrices).

Let us first consider simulated examples of the intrinsic Wilcoxon rank-sum test (test = "rank.sum") and intrinsic Kruskal-Wallis test (test = "krusk.wall").

## Generate data (null true)
data1 <- array(c(X1, replicate(50, Expm(diag(2), H.coeff(0.5 * rnorm(4), inverse = T)))), dim = c(2, 2, 100)) ## pointwise HPD sample
data2 <- array(c(X2, replicate(50, sapply(1:5, function(i) Expm(i * diag(2), H.coeff(0.5 * rnorm(4), inverse = T) / i), simplify = "array"))), dim = c(2, 2, 5, 100)) ## HPD curve sample

## Generate data (null false)
data1a <- array(c(X1, replicate(50, Expm(diag(2), H.coeff(rnorm(4), inverse = T)))), dim = c(2, 2, 100)) ## pointwise HPD scale change
data2a <- array(c(X2, replicate(50, sapply(1:5, function(i) Expm(i * diag(2), H.coeff(rnorm(4), inverse = T) / i), simplify = "arra"))), dim = c(2, 2, 5, 100)) ## HPD curve scale change

## Rank-sum test
pdRankTests(data1, sample_sizes = c(50, 50), "rank.sum")[1:4] ## null true (pointwise)
#> $test
#> [1] "Intrinsic Wilcoxon rank-sum"
#> 
#> $p.value
#> [1] 0.1037495
#> 
#> $statistic
#> [1] 1.626941
#> 
#> $null.distr
#> [1] "Standard normal distribution"

pdRankTests(data2, sample_sizes = c(50, 50), "rank.sum")[2] ## null true (curve)
#> $p.value
#> [1] 0.9834998

pdRankTests(data1a, sample_sizes = c(50, 50), "rank.sum")[2] ## null false (pointwise)
#> $p.value
#> [1] 6.958285e-11

pdRankTests(data2a, sample_sizes = c(50, 50), "rank.sum")[2] ## null false (curve)
#> $p.value
#> [1] 1.020181e-15

## Kruskal-Wallis test
pdRankTests(data1, sample_sizes = c(50, 25, 25), "krusk.wall")[1:4] ## null true (pointwise)
#> $test
#> [1] "Intrinsic Kruskal-Wallis"
#> 
#> $p.value
#> [1] 0.1443239
#> 
#> $statistic
#> [1] 3.87139
#> 
#> $null.distr
#> [1] "Chi-squared distribution (df = 2)"

pdRankTests(data2, sample_sizes = c(50, 25, 25), "krusk.wall")[2] ## null true (curve)
#> $p.value
#> [1] 0.1526552

pdRankTests(data1a, sample_sizes = c(50, 25, 25), "krusk.wall")[2] ## null false (pointwise)
#> $p.value
#> [1] 5.495634e-10

pdRankTests(data2a, sample_sizes = c(50, 25, 25), "krusk.wall")[2] ## null false (curve)
#> $p.value
#> [1] 1.033775e-14

To apply the manifold Wilcoxon signed-rank test (test = "signed-rank"), we generate paired observations for independent trials (or subjects) by introducing trial-specific random effects, such that the paired observations in each trial share a trial-specific geometric mean. Note that for such data the intrinsic Wilcoxon rank-sum test is no longer valid due to the introduced sample dependence.

## Trial-specific means
mu <- replicate(50, Expm(diag(2), H.coeff(0.1 * rnorm(4), inverse = T)))

## Generate paired samples X,Y
make_sample <- function(null) sapply(1:50, function(i) Expm(mu[, , i], pdSpecEst:::T_coeff_inv(ifelse(null, 1, 0.5) * rexp(4) - 1, mu[, , i])), simplify = "array")

X3 <- make_sample(null = T) ## refernce sample 
Y3 <- make_sample(null = T) ## null true
Y3a <- make_sample(null = F) ## null false (scale change)

## Signed-rank test
pdRankTests(array(c(X3, Y3), dim = c(2, 2, 100)), test = "signed.rank")[1:4] ## null true
#> $test
#> [1] "Intrinsic Wilcoxon signed-rank"
#> 
#> $p.value
#> [1] 0.2025809
#> 
#> $statistic
#>   V 
#> 505 
#> 
#> $null.distr
#> [1] "Wilcoxon signed rank test with continuity correction"

pdRankTests(array(c(X3, Y3a), dim = c(2, 2, 100)), test = "signed.rank")[2] ## null false
#> $p.value
#> [1] 0.001444632

The intrinsic signed-rank test also provides a valid procedure to test for equivalence of spectral matrices of two (independent) multivariate stationary time series based on the HPD periodogram matrices obtained via pdPgram(). In contrast to other available tests in the literature, this asymptotic test does not require consistent spectral estimators or resampling of test statistics, and therefore remains computationally efficient for higher-dimensional spectral matrices or a large number of sampled Fourier frequencies.

## Signed-rank test for equivalence of spectra
## vARMA(1,1) process: Example 11.4.1 in (Brockwell and Davis, 1991)
Phi <- array(c(0.7, 0, 0, 0.6, rep(0, 4)), dim = c(2, 2, 2))
Theta <- array(c(0.5, -0.7, 0.6, 0.8, rep(0, 4)), dim = c(2, 2, 2))
Sigma <- matrix(c(1, 0.71, 0.71, 2), nrow = 2)
pgram <- function(Sigma) pdPgram(rARMA(2^10, 2, Phi, Theta, Sigma)$X)$P ## HPD periodogram

## Null is true
pdRankTests(array(c(pgram(Sigma), pgram(Sigma)), dim = c(2, 2, 2^10)), test = "signed.rank")[2]
#> $p.value
#> [1] 0.4047506

## Null is false
pdRankTests(array(c(pgram(Sigma), pgram(0.9 * Sigma)), dim = c(2, 2, 2^10)), test = "signed.rank")[2]
#> $p.value
#> [1] 0.003768685

To conclude, we demonstrate the intrinsic Bartels-von Neumman test (test = "bartels") by generating an independent non-identically distributed sample of HPD matrices with a trend in the scale of the distribution across observations. In this case, the null hypothesis of randomness breaks down.

## Null is true
data3 <- replicate(200, Expm(diag(2), H.coeff(rnorm(4), inverse = T))) ## iid HPd sample
data4 <- replicate(100, sapply(1:5, function(i) Expm(i * diag(2), H.coeff(rnorm(4), inverse = T) / i), simplify = "array")) ## iid HPD curve

## Null is false
data3a <- sapply(1:200, function(j) Expm(diag(2), H.coeff(((200 - j) / 200 + j * 2 / 200) * rnorm(4), inverse = T)), simplify = "array") ## pointwise trend in scale
data4a <- sapply(1:100, function(j) sapply(1:5, function(i) Expm(i * diag(2), H.coeff(((100 - j) / 100 + j * 2 / 100) * rnorm(4), inverse = T) / i), simplify = "array"), simplify = "array") ## curve trend in scale

## Bartels-von Neumann test
pdRankTests(data3, test = "bartels")[1:4] ## null true (pointwise)
#> $test
#> [1] "Intrinsic Bartels-von Neumann"
#> 
#> $p.value
#> [1] 0.2641266
#> 
#> $statistic
#> [1] 1.116691
#> 
#> $null.distr
#> [1] "Standard normal distribution"

pdRankTests(data4, test = "bartels")[2] ## null true (curve)
#> $p.value
#> [1] 0.7612759

pdRankTests(data3a, test = "bartels")[2] ## null false (pointwise)
#> $p.value
#> [1] 1.612058e-05

pdRankTests(data4a, test = "bartels")[2] ## null false (curve)
#> $p.value
#> [1] 0.000732231