coverage_correlation

The covercorr package implements the coverage correlation coefficient introduced in the paper Coverage correlation: detecting singular dependencies between random variables. The coverage correlation coefficient is a nonparametric measure of statistical association designed to detect dependencies concentrated on low-dimensional structures within the joint distribution of two (or more) random variables or vectors. Based on Monge-Kantorovich ranks and geometric coverage processes, this statistic quantifies the extent to which the joint distribution concentrates on a singular subset with respect to the product of the marginals. The coverage correlation coefficient is distribution-free, admits an analytically tractable asymptotic null distribution, and can be computed efficiently, making it well-suited for uncovering complex, potentially nonlinear associations in large-scale testing.

This vignette covers the original pairwise function coverage_correlation(), and then the functions added in version 1.1.0: the unified front-end covercorr(), deterministic reference grids via coverage_correlation_grid(), and joint dependence among many variables via coverage_correlation_K().

library(covercorr)

The pairwise coefficient

Example 1

In this example, we demonstrate how to use coverage_correlation() with a simple simulation. We compute the coverage correlation coefficient between two one-dimensional Normal random variables, X and Y, and then vary the strength of their relationship to observe how the statistic changes.

set.seed(1)
n <- 100
X <- rnorm(n)
Y <- rnorm(n)
result <- coverage_correlation(Y, X, visualise = TRUE)

result
#> 
#> Coverage correlation
#> 
#> statistic = 0.03622657, p-value = 0.1151168
#> n = 100, dimension = 2, method = exact

In the example above, X and Y are independent.

The parameter visualise defaults to FALSE, but setting it to TRUE produces a plot that illustrates the intuition behind the coverage correlation coefficient. The coverage correlation coefficient first transforms X and Y into their Monge-Kantorovich ranks, denoted by X_rank and Y_rank, which are uniformly distributed on \([0, 1]\). The plot displays the pairs \((X_{\text{rank}_i}, Y_{\text{rank}_i})\) along with cubes of volume \(n^{-1}\).

Inside the function coverage_correlation(), we compute \(V_n\), the total uncovered area after taking the union of these cubes. The coverage correlation coefficient is then defined as \[ \kappa_n^{X, Y} := \frac{V_n - e^{-1}}{1 - e^{-1}}. \]

The function returns an object of class "covercorr" (a list) with the following elements:

Printing the result (as above) uses the print method for the "covercorr" class; use str(result) to see the full list structure.

By default, method = "auto". In this mode, if the total dimension of X and Y (i.e., ncol(X) + ncol(Y), treating vectors as one-dimensional) is at most 6, the method is set to "exact"; otherwise, it uses "approx".

Next we can see how the result changes as we introduce dependence between X and Y.

set.seed(2)
n <- 100
X <- rnorm(n)
Z <- rnorm(n)
rho <- 0.9
Y <- rho * X + sqrt(1 - rho^2) * Z
result <- coverage_correlation(Y, X, visualise = TRUE)

result
#> 
#> Coverage correlation
#> 
#> statistic = 0.3591986, p-value = < 2.2204e-16
#> n = 100, dimension = 2, method = exact

You may notice parts of some cubes appearing at the corners of the plot. This happens because we treat \([0, 1]^2\) as a torus. If a cube centered at one of the rank points lies partially outside \([0, 1]^2\), we wrap it around so that the plot reflects this topology.

Example 2

The coverage correlation coefficient can handle multidimensional random vectors as well.

set.seed(3)
n <- 100
p <- 2
X <- matrix(rnorm(p * n), ncol = p)
Y <- matrix(0, nrow = n, ncol = p)
Y[, 1] <- X[, 1]^2
Y[, 2] <- X[, 1] * X[, 2]
result <- coverage_correlation(Y, X)
result
#> 
#> Coverage correlation
#> 
#> statistic = 0.273216, p-value = < 2.2204e-16
#> n = 100, dimension = 4, method = exact

In this case we cannot visualise the whole plot as X and Y are not one-dimensional.

Example 3

In the example below, X and Y are independent and 2-dimensional. We set the method parameter equal to "approx".

set.seed(4)
n <- 50
p <- 2
X <- matrix(rnorm(p * n), ncol = p)
Y <- matrix(rnorm(p * n), ncol = p)
result <- coverage_correlation(Y, X, method = "approx")
result
#> 
#> Coverage correlation
#> 
#> statistic = 0.03919494, p-value = 0.0682784
#> n = 50, dimension = 4, method = approx
#> Monte Carlo standard error = 0.02595504

Deterministic reference grids

The Monge-Kantorovich rank transformation requires a set of reference points. By default coverage_correlation() generates these randomly with runif(), which means repeated calls on the same data can give slightly different results. When both inputs are one-dimensional, the natural reference set is the deterministic uniform grid \(\{1/n, 2/n, \ldots, 1\}\), which makes the computation reproducible. coverage_correlation_grid() uses this grid by default.

set.seed(5)
n <- 100
X <- rnorm(n)
Y <- sin(3 * X) + rnorm(n) * 0.05

# Deterministic: repeated calls give identical results
r1 <- coverage_correlation_grid(Y, X)
r2 <- coverage_correlation_grid(Y, X)
c(r1$stat, r2$stat)
#> [1] 0.3624634 0.3624634

You can also supply your own reference points through u and v, for example a low-discrepancy sequence, as long as they match the shape of the corresponding input.

The unified interface: covercorr()

covercorr() is a single front-end that dispatches to the appropriate routine based on its input. It accepts a pair of variables, a matrix or data frame whose columns are treated as variables, or a list of variables. The reference argument chooses between random reference points ("random", the default) and the deterministic grid ("deterministic").

set.seed(6)
n <- 100
X <- rnorm(n)
Y <- sin(3 * X) + rnorm(n) * 0.05

# Pairwise, random reference points
covercorr(X, Y)
#> 
#> Coverage correlation
#> 
#> statistic = 0.4009571, p-value = < 2.2204e-16
#> n = 100, dimension = 2, method = exact

# Pairwise, deterministic grid
covercorr(X, Y, reference = "deterministic")
#> 
#> Coverage correlation
#> 
#> statistic = 0.3661019, p-value = < 2.2204e-16
#> n = 100, dimension = 2, method = exact

When given a matrix or data frame, each column is treated as a separate variable. A two-column input is handled as the pairwise case; more columns are treated jointly as several variables.

set.seed(7)
n <- 100
X <- rnorm(n)
Y <- sin(3 * X) + rnorm(n) * 0.05
covercorr(cbind(X, Y))
#> 
#> Coverage correlation
#> 
#> statistic = 0.37483, p-value = < 2.2204e-16
#> n = 100, dimension = 2, method = exact

Joint dependence among many variables

coverage_correlation_K() generalises the coefficient from a pair to K variables, supplied as a list. It tests the null hypothesis of mutual independence among all the inputs.

set.seed(8)
n <- 100
X1 <- rnorm(n)
X2 <- sin(3 * X1) + rnorm(n) * 0.05
X3 <- X1 * X2 + rnorm(n) * 0.1
coverage_correlation_K(list(X1, X2, X3))
#> 
#> Coverage correlation
#> 
#> statistic = 0.4638644, p-value = < 2.2204e-16
#> n = 100, dimension = 3, method = exact

The same dependence can be tested through covercorr() by passing a list, or a matrix with more than two columns.

set.seed(9)
n <- 100
X1 <- rnorm(n)
X2 <- sin(3 * X1) + rnorm(n) * 0.05
X3 <- X1 * X2 + rnorm(n) * 0.1
covercorr(cbind(X1, X2, X3))
#> 
#> Coverage correlation
#> 
#> statistic = 0.439989, p-value = < 2.2204e-16
#> n = 100, dimension = 3, method = exact

A note on the K-variable fixed-grid case

A deterministic-grid p-value is available for the pairwise case and for the two-variable fixed-grid case. For three or more variables with a fixed grid, the null distribution has not been derived, so coverage_correlation_K_grid() (and covercorr() with reference = "deterministic" on three or more variables) returns pval = NA and reports pval_available = FALSE, while still computing the statistic itself.

set.seed(10)
n <- 100
X1 <- rnorm(n)
X2 <- sin(3 * X1) + rnorm(n) * 0.05
X3 <- X1 * X2 + rnorm(n) * 0.1
res <- coverage_correlation_K_grid(list(X1, X2, X3))
#> No p-value is available for the K-variable (K > 2) fixed-grid case: the null distribution has not been derived. Returning pval = NA.
res$stat
#> [1] 0.4222658
res$pval
#> [1] NA
res$pval_available
#> [1] FALSE

Visualising the joint density on real data

The examples above use simulated data so that the dependence structure is known. The package also ships a real dataset, CD8T, containing single-cell RNA-seq gene expression measurements for 9369 fetal CD8+ T cells across 1000 genes (Suo et al., 2022). We can use it to illustrate visualise_density(), which gives a picture of where in the unit square a pair of variables departs from independence.

visualise_density() takes the two Monge-Kantorovich rank vectors and estimates their joint density with a toroidally-wrapped two-dimensional kernel density estimate (the same torus topology used in the cube plots above, so that mass near an edge is not lost). It then plots the density on a log scale relative to the uniform density: regions redder than the mid-tone are more concentrated than they would be under independence, bluer regions are less concentrated, and a flat mid-toned image corresponds to independence.

For illustration we take two genes; you can substitute any pair of column indices. We first compute the coverage correlation coefficient for the pair, then place the raw scatter plot and the rank-density plot side by side.

data(CD8T)
n <- nrow(CD8T)

i <- 1
j <- 2
x <- CD8T[, i]
y <- CD8T[, j]

# The coefficient and its p-value for this pair
covercorr(x, y)
#> 
#> Coverage correlation
#> 
#> statistic = 0.0132117, p-value = 1.24082e-05
#> n = 9369, dimension = 2, method = exact

# Monge-Kantorovich ranks, using random reference points
set.seed(1)
x_rank <- MK_rank(as.matrix(x), runif(n))
y_rank <- MK_rank(as.matrix(y), runif(n))

op <- par(mfrow = c(1, 2))
plot(x, y, pch = ".", xlab = colnames(CD8T)[i], ylab = colnames(CD8T)[j])
visualise_density(x_rank, y_rank)

par(op)

Single-cell expression data are sparse, so many cells share the value zero for a given gene; MK_rank() breaks these ties at random, which is why the scatter plot shows banding along the axes. Pairs that the coverage correlation test flags as dependent typically show clear ridges of concentrated (red) mass in the density plot, whereas independent pairs give an essentially flat image.