Gonjo Demonstration

library(TFORGE)

The Gonjo basin anisotropy of magnetic susceptibility (AMS) data from Li et al. (2020) is available with the TFORGE package. This data contains the AMS tensor of 542 specimens along with the depth of each specimen. See ?Gonjo for details. For analysis Li et al. (2020) split the data at depths of 765m and 2154m, leading to three intervals.

Each AMS tensor is a 3 by 3 symmetric matrix with a trace of 1 and non-negative eigenvalues, so we can use test_multiplicity_nonnegative() and test_fixedtrace() to test some hypotheses relevant to AMS data. Because the eigenvalues of the matrices tend to be well away from zero (smallest eigenvalue of all matrices = 0.267), it would likely be okay to use test_multiplicity() too.

Each interval has over 100 specimens. For fast computation in this demonstration we will use only the first 30 specimens of each interval, and only 100 bootstrap resamples.

Preparation

Split the AMS matrices in Gonjo into three groups, keeping the first 30 of each for computation speed in this demonstration.

interval_1 <- Gonjo$matrices[Gonjo$datatable$`really depth` <= 765]
interval_2 <- Gonjo$matrices[Gonjo$datatable$`really depth` <= 2154]
interval_3 <- Gonjo$matrices[Gonjo$datatable$`really depth` > 2154]
samples <- list(`1` = interval_1[1:30], `2` = interval_2[1:30], `3` = interval_3[1:30])

Multiplicity Tests

Following descriptions from Li et al. (2020), here we will test whether population mean of intervals 1, 2, and 3 are oblate, isotropic, or prolate respectively. These are hypotheses about the multiplicity of eigenvalues of the population means.

For interval 1 and 2 the samples are so far from oblate and isotropic, respectively, that weighted bootstrapping (using emplik()) could not be used, and the result is automatically a p-value of 0, corresponding to rejecting the null hypothesis. At a significance level of 0.05, the prolate hypothesis for interval 3 should also be rejected.

Interval 1 oblate:

test_multiplicity_nonnegative(samples[["1"]], mult = c(2, 1), B = 1000)
#> Warning: The proposed null mean appears to be outside the convex hull for at
#> least one sample (emplik() did not converge)
#> $pval
#> [1] 0
#> 
#> $t0
#> [1] 74.30601

Interval 2 isotropic:

test_multiplicity_nonnegative(samples[["2"]], mult = 3, B = 1000)
#> Warning: The proposed null mean appears to be outside the convex hull for at
#> least one sample (emplik() did not converge)
#> $pval
#> [1] 0
#> 
#> $t0
#> [1] 439.724

Interval 3 prolate:

test_multiplicity_nonnegative(samples[["3"]], mult = c(1, 2), B = 1000)
#> $pval
#> [1] 0.018
#> 
#> $t0
#> [1] 20.80563

Because the eigenvalues of the data are so far away from zero (smallest eigenvalue of all matrices = 0.267) the transformation-based bootstrapping used by test_multiplicity() would likely work well too and give the same result. Below we apply test_multiplicity(), and the results lead to the same conclusions as test_multiplicity_nonnegative(). For computational speed in this demonstration we will use only 100 resamples; we recommend using at least 1000 resamples in practise.

test_multiplicity(samples[["1"]], mult = c(2, 1), B = 100)
#> $pval
#> [1] 0
#> 
#> $t0
#> [1] 74.30601
test_multiplicity(samples[["2"]], mult = 3, B = 100)
#> $pval
#> [1] 0
#> 
#> $t0
#> [1] 439.724
test_multiplicity(samples[["3"]], mult = c(1, 2), B = 100)
#> $pval
#> [1] 0
#> 
#> $t0
#> [1] 20.80563

3-Sample Test

Here we can check whether the eigenvalues of the population means are equal. Again the intervals are so different that the convex hulls of these samples do not overlap at the estimated common mean, so equality of eigenvalues is automatically rejected.

test_fixedtrace(samples, B = 1000)
#> Warning in (function (ms, nullmean) : emplik() did not converge, which usually
#> means that the proposed null mean is outside the convex hull of the data
#> Warning: The proposed null mean appears to be outside the convex hull for at
#> least one sample (empirical likelihood weights are zero)
#> $pval
#> [1] 0
#> 
#> $t0
#> [1] 206.5525
#> attr(,"null_evals")
#> [1] 0.3392680 0.3345776 0.3261544
#> attr(,"null_evals_proj")
#> [1] 0.003316612 0.008792391

Confidence Regions

We can also estimate confidence regions for the eigenvalues of the population means. For computational speed in this demonstration we will use only 100 resamples; we recommend using at least 1000 resamples in practise.

s1_CR <- conf_fixedtrace(samples[["1"]],
  alpha = 0.05,
  B = 100, npts = 1000, check = FALSE
)
s2_CR <- conf_fixedtrace(samples[["2"]],
  alpha = 0.05,
  B = 100, npts = 1000, check = FALSE
)
s3_CR <- conf_fixedtrace(samples[["3"]],
  alpha = 0.05,
  B = 100, npts = 1000, check = FALSE
)

These confidence regions are ellipses on a 2D plane. Plotting these confidence regions can be a bit involved - please see the reproducibility documentation of Hingee et al. (2026) for an example using ggplot2.

References

Li, S., van Hinsbergen, D. J. J., Shen, Z., Najman, Y., Deng, C., & Zhu, R. (2020). Anisotropy of Magnetic Susceptibility (AMS) Analysis of the Gonjo Basin as an Independent Constraint to Date Tibetan Shortening Pulses. Geophysical Research Letters, 47(8), e2020GL087531. https://doi.org/10.1029/2020GL087531

Hingee K, Scealy J, Wood A (2026). “Nonparametric bootstrap inference for the eigenvalues of geophysical tensors.” Journal of American Statistical Association. Accepted for publication.