Post-clustering inference after any clustering algorithm

Performs post-clustering inference for the difference between the means of two clusters obtained from a user-specified clustering algorithm, under a general matrix normal model.

The clustering algorithm is given by a user-provided function cl_fun, allowing this method to be applied to a wide range of clustering procedures.

The covariance structure between features Sigma is estimated if not provided, respecting the selective type I error control. The covariance matrix between observations U is assumed known and should have a compound symmetry structure (see details).

Usage

test.clusters.MC(
  X,
  U = NULL,
  Sigma = NULL,
  Y = NULL,
  UY = NULL,
  precUY = NULL,
  clusters,
  cl_fun,
  NC = NULL,
  cl = NULL,
  ndraws = 2000,
  sample_split = FALSE,
  nY = NULL,
  return_Sigma = FALSE,
  return_X_clus = FALSE
)

Arguments

X: A numeric $n \times p$ matrix assumed to arise from a matrix normal distribution $\mathcal{MN}(M, U, \Sigma)$.
U: An $n \times n$ positive-definite matrix describing the dependence structure between the rows of X. If NULL, observations are assumed to be independent and U is set to the identity matrix.
Sigma: A $p \times p$ positive-definite matrix describing the dependence structure between the columns of X. If NULL, Sigma is over-estimated from an auxiliary independent sample Y (in the sense of the Loewner partial order).
Y: If Sigma is NULL, an independent copy of X used to estimate Sigma. It must have the same number of columns as X.
UY: If Sigma is NULL, an $n_Y \times n_Y$ positive-definite matrix describing the dependence structure between the rows of Y. If NULL and precUY is not provided, the identity matrix is used by default.
precUY: The inverse of UY. Providing precUY may improve computational efficiency. If UY is provided but precUY is NULL, it is computed internally.
clusters: Integer vector of length 2 specifying the pair of clusters to compare.
cl_fun: A function implementing the clustering procedure. It must take as input the data matrix X and return an integer vector of length $n$ containing cluster assignments. If required by the clustering algorithm, cl_fun should also take an argument NC specifying the number of clusters to be used.
NC: Integer scalar specifying the number of clusters to be used by cl_fun. If the clustering algorithm determines the number of clusters automatically, this argument should be set to NULL.
cl: Optional integer vector of length $n$ giving a precomputed clustering of X. If provided, cl_fun is not called.
ndraws: Integer. Number of Monte Carlo samples used to approximate the p-value.
sample_split: Logical. Whether to use sample splitting to estimate Sigma when Sigma = NULL. Ignored when Sigma is provided by the user.
nY: Integer. If Y is not provided and sample_split = TRUE, the number of rows of the auxiliary sample Y used to estimate Sigma. If nY is NULL, half of the rows of X are used for estimation. Ignored when Sigma is provided by the user.
return_Sigma: Logical. Whether to include the column covariance matrix used in the test in the returned list. Ignored when Sigma is provided by the user. Default is FALSE.
return_X_clus: Logical. If sample splitting is performed to estimate Sigma, whether to include the data matrix used for clustering in the returned list. Ignored when sample_split = FALSE (as the same data matrix is used for clustering and testing). If further analysis of the retrieved clusters is desired, we recommend setting return_X_clus = TRUE when sample_split = TRUE to avoid confusion. Default is FALSE.

Value

A named list with the following components:

pvalue: The p-value for testing equality of the two selected cluster means.
stat: The observed test statistic.
stderr: Monte Carlo standard error of the estimated p-value.
n_preserved: Number of Monte Carlo samples (out of ndraws) for which preserve.cl returned TRUE, i.e., the perturbed clustering matched the original cluster labels.
clusters: An integer vector of length $n$ giving the cluster membership of each observation.
Sigma: If return_Sigma = TRUE, the column covariance matrix used in the test, either provided by the user or estimated from Y.
X_clus: If return_X_clus = TRUE and sample_split = TRUE, the data matrix used for clustering/testing (i.e., the subsample of X retained for the clustering step).

Details

The method applies to any clustering algorithm that can be expressed through a function cl_fun returning cluster assignments. The selective inference procedure is based on a Monte Carlo approximation of the truncation region induced by the clustering step.

Selective type I error control is guaranteed when the row covariance matrix $\mathbf{U}$ has a compound symmetry (CS) structure, i.e., $$ \mathbf{U} = (a - b)\mathbf{I}_n + b\mathbf{1}_n, $$ for some $a/(n - 1) < b < a$. In practice, the method is robust to moderate deviations from the CS structure. In particular, reliable performance is typically observed when $\mathbf{U}$ remains close to CS, for example in autoregressive (AR(1)), diagonal, banded, or Toeplitz covariance structures, depending on their parameters. We therefore recommend applying the method primarily in settings where $\mathbf{U}$ is known and does not deviate substantially from the compound symmetry structure.

When Sigma = NULL, the column covariance matrix is estimated from an auxiliary independent sample Y. When U = NULL, row-wise independence is assumed.

References

Gao, L. L., Bien, J., and Witten, D. (2022). Selective inference for hierarchical clustering. Journal of the American Statistical Association, 117(540), 2533–2547.

González-Delgado, J., Deronzier, M. Cortés, J., and Neuvial, P. (2023) Post-clustering Inference under Dependence. arXiv.2310.11822.

Examples

n <- 50
p <- 20

# Simulating under the alternative hypothesis
M <- Matrix::Matrix(0, nrow = n, ncol = p)
M[1:floor(n/3),] <- -1
M[(floor(2*n/3)+1):n,] <- 1
Sigma <- stats::toeplitz(seq(1, 0.1, length.out = p))
U <- matrixNormal::I(n)

X <- matrixNormal::rmatnorm(s = 1, M, U, Sigma)
Y <- matrixNormal::rmatnorm(s = 1, M, U, Sigma)

# Example using HDBSCAN clustering (it chooses the number of clusters)
# install.packages("dbscan")

hdbscan_clustering <- function(X, min.occupancy = 5) {
  X.clus <- dbscan::hdbscan(X, minPts = min.occupancy)
  return(X.clus$cluster + 1) 
}

# Precompute clustering
cl_X <- hdbscan_clustering(X)

# Test difference between clusters 3 and 1
test.hdbscan <- test.clusters.MC(
  X = X, U = U, Sigma = Sigma,
  clusters = c(3, 1),
  cl = cl_X,
  cl_fun = hdbscan_clustering,
  NC = NULL,
  ndraws = 500
)
test.hdbscan$pvalue
#> [1] 0.4828108