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).
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+5.55e-17i 0.352+1.47e-01i 0.352-1.47e-01i ...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 ⟨., .⟩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.8825155By 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.7113942A 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 44 34
#> [26] 24 26 23 25 30 10 29 16 32 1 41 27 33 5 2 3 17 15 40 42 20 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 17 24
#> [26] 23 35 13 44 20 27 19 39 34 32 22 25 42 29 8 16 1 2 48 15 37 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+1.387779e-17i -0.0154483+0.1752587i
#> [2,] -0.0154483-1.752587e-01i 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.000000iWe 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+5.079487e-17i -0.03426934+3.226425e-02i
#> [2,] -0.03426934-3.226425e-02i 1.15446038+5.421011e-19i
pdDepth(y = med.X1, X = X1, method = "gdd") ## maximum out-of-sample depth
#> [1] 0.4409941
max(dd1) ## maximum in-sample depth
#> [1] 0.4306943pdRankTests()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-14To 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.001444632The 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.003768685To 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