Post-clustering inference after k-means clustering

Performs post-clustering inference for the difference between the means of two clusters obtained by k-means 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).

Usage

test.clusters.km(
  X,
  U = NULL,
  Sigma = NULL,
  Y = NULL,
  UY = NULL,
  precUY = NULL,
  NC,
  clusters,
  itermax = 10,
  tol = 1e-06,
  km_at_cl = NULL,
  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 $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. Providing precUY may improve computational efficiency. If UY is provided but precUY is NULL, it is computed internally.
NC: Integer scalar giving the number of clusters.
clusters: Integer vector of length 2 specifying the pair of clusters to compare. Entries must belong to 1:NC.
itermax: Integer. Maximum number of iterations for the k-means algorithm, passed to KmeansInference::kmeans_inference.
tol: Numeric tolerance parameter controlling convergence of the k-means algorithm, passed as tol_eps.
km_at_cl: An optional precomputed output of KmeansInference::kmeans_inference(). When supplied, the clustering and truncation-set computation steps are skipped and the precomputed partition and interval are used directly. This is useful when testing the same cluster pair under several covariance assumptions, as it avoids redundant calls to the inference algorithm. The object must be compatible with X and NC: the number of cluster labels must equal nrow(X), the number of distinct clusters must equal NC, and all labels in clusters must be present in the partition. Ignored (with a warning) when sample_split = TRUE, because sample splitting changes X after the clustering would have been computed.
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 in the returned list the column covariance matrix Sigma used in the test. When return_Sigma = TRUE, the returned matrix is included whether Sigma was provided by the user or estimated internally. 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.
km: An integer vector of length $n$ giving the cluster membership of each observation returned by the k-means algorithm.
S: The truncation set used to compute the p-value.
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 (subsample of X) used for clustering and testing.

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.

The clustering step is performed using the function KmeansInference::kmeans_inference, which provides both the clustering partition and the truncation region associated with the selection event. The p-value is computed using an exact characterization of the truncation region for the Euclidean norm, following Chen and Witten (2022), and then mapped to the corresponding statistic under the general matrix normal model.

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.

Chen, Y. T., and Witten, D. M. (2023). Selective inference for k-means clustering. Journal of Machine Learning Research, 24(152), 1-41.

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)

# k-means with known Sigma
test.km <- test.clusters.km(
  X = X, U = U, Sigma = Sigma,
  NC = 3, clusters = sample(1:3, 2)
)
test.km$pvalue
#> [1] 0.941819

# k-means with over-estimation of Sigma
test.km <- test.clusters.km(
  X = X, U = U, Sigma = NULL, Y = Y,
  NC = 3, clusters = sample(1:3, 2)
)
#> Sigma not provided: plugging an over-estimate.
test.km$pvalue
#> [1] 0.7318592