Post-clustering inference after hierarchical clustering

Performs post-clustering inference for the difference between the means of two clusters obtained by hierarchical agglomerative clustering, under a general matrix normal model.

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).

Supported linkage criteria (as in Gao et al. 2022) are "single", "average", "centroid", "ward.D", "median", "mcquitty", and "complete".

Usage

test.clusters.hc(
  X,
  U = NULL,
  Sigma = NULL,
  Y = NULL,
  UY = NULL,
  precUY = NULL,
  NC,
  clusters,
  linkage = "average",
  hcl = NULL,
  ndraws = 2000,
  sample_split = FALSE,
  nY = NULL,
  return_Sigma = FALSE,
  return_X_clus = FALSE,
  dismat = NULL
)

Arguments

X: An $n \times p$ data 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 $n \times n$ 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. Supplying precUY may improve computational efficiency. If UY is provided but precUY is NULL, precUY is computed internally by matrix inversion.
NC: Integer. Number of clusters in the partition obtained by cutting the hierarchical clustering tree.
clusters: Integer vector of length 2 containing the labels of the two clusters to compare. Its entries must belong to 1:NC.
linkage: Character string specifying the linkage criterion used in hierarchical clustering. Must be one of "single", "average", "centroid", "ward.D", "median", "mcquitty", or "complete". When hcl is provided, the linkage is inferred from hcl$method and this argument is overridden. If linkage is supplied explicitly alongside hcl but does not match hcl$method, a warning is emitted and the method stored in hcl takes precedence.
hcl: An optional precomputed hierarchical clustering object of class "hclust", as returned by fastcluster::hclust() or stats::hclust() on the squared Euclidean distance matrix of X. When supplied, the clustering step is skipped and the linkage criterion is inferred from hcl$method, overriding the linkage argument. Passing a precomputed hcl is useful when testing several cluster pairs from the same clustering, as it avoids recomputing the dendrogram each time. The number of leaves in hcl must equal nrow(X); an error is raised otherwise. Ignored (with a warning) when sample_split = TRUE, because sample splitting changes X after the clustering would have been computed.
ndraws: Integer. Number of Monte Carlo samples used to approximate the p-value when linkage = "complete". Ignored otherwise.
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.
dismat: An optional precomputed squared Euclidean distance object of class "dist", as returned by stats::dist(X, method = "euclidean")^2. When supplied alongside hcl, both the distance object and the dendrogram computations are skipped, which is useful when testing several cluster pairs from the same clustering. When supplied without hcl, a warning is emitted and the dendrogram is computed from dismat. When hcl is supplied without dismat and the linkage is not "complete", a warning is emitted and dismat is recomputed internally. Ignored (with a warning) when sample_split = TRUE, because sample splitting changes X after the distance object would have been computed.

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 when linkage = "complete". This component is omitted otherwise.
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. Only present when linkage = "complete".
S: The truncation set used to compute the p-value when linkage != "complete". This component is omitted otherwise.
hcl: An integer vector of length $n$ giving the cluster membership of each observation in the partition with NC clusters.
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 and testing (i.e., the subsample of X retained after sample splitting).

Details

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.

For linkage criteria other than "complete", the p-value is computed using an exact characterization of the truncation region, adapted from clusterpval. For "complete" linkage, the truncation region is not computed explicitly and the p-value is approximated by Monte Carlo. In this case, the Monte Carlo procedure can be automatically parallelized using the future package (see examples).

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 null hypothesis
M <- Matrix::Matrix(0, nrow = n, ncol = p) 
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)

# Hierarchical clustering with average linkage and known Sigma
test.hc <- test.clusters.hc(
  X = X, U = U, Sigma = Sigma,
  NC = 3, clusters = sample(1:3, 2),
  linkage = "average"
)
test.hc$pvalue
#> [1] 0.6363383

# Testing two more cluster pairs reusing the same dendrogram
hcl_obj <- fastcluster::hclust(dist(X)^2, method = "average")
test.hc12 <- test.clusters.hc(
  X = X, U = U, Sigma = Sigma,
  NC = 3, clusters = c(1, 2), hcl = hcl_obj
)
#> Warning: 'hcl' is provided but 'dismat' is not. The distance matrix will be recomputed internally. Pass 'dismat' as well to avoid recomputing it when testing multiple cluster pairs.
test.hc13 <- test.clusters.hc(
  X = X, U = U, Sigma = Sigma,
  NC = 3, clusters = c(1, 3), hcl = hcl_obj
)
#> Warning: 'hcl' is provided but 'dismat' is not. The distance matrix will be recomputed internally. Pass 'dismat' as well to avoid recomputing it when testing multiple cluster pairs.

# Hierarchical clustering with complete linkage and estimated Sigma
test.hc <- test.clusters.hc(
  X = X, U = U, Sigma = NULL, Y = Y,
  NC = 3, clusters = sample(1:3, 2),
  linkage = "complete"
)
#> Sigma not provided: plugging an over-estimate.
#> Clustering with complete linkage. Monte-Carlo approximation of the p-value.
test.hc$pvalue
#> [1] 0.07380339

# Hierarchical clustering with complete linkage, estimated Sigma, and parallelization
if (FALSE) { # \dontrun{
library(future)
# Set parallelization plan (adjust workers according to your machine):
plan(multisession, workers = 4)
test.hc <- test.clusters.hc(
  X = X, U = U, Sigma = NULL, Y = Y,
  NC = 3, clusters = sample(1:3, 2),
  linkage = "complete", ndraws = 500
)
test.hc$pvalue
plan(sequential) # Reset to avoid leaving open connections
} # }