Technical details and practical advice

What the tests condition on

Write $\hat{\mathcal C}_1$ and $\hat{\mathcal C}_2$ for the two tested clusters and $r = |\hat{\mathcal C}_1| + |\hat{\mathcal C}_2|$ . The tests assume that the rows of the data are independent Gaussian vectors with a common covariance matrix $\boldsymbol\Sigma$ , and they test the null hypothesis that all $r$ observations in the two clusters share one mean.

The deviations of the tested observations from their common mean form an $r \times p$ matrix. Its polar decomposition writes it as $\mathbf{Q}\mathbf{H}$ , where $\mathbf{H}$ is the square root of $r - 1$ times their sample covariance matrix and $\mathbf{Q}$ has orthonormal columns. Under the null hypothesis, $\mathbf{Q}$ is uniformly distributed and independent of $\mathbf{H}$ and of the rest of the data, whatever $\boldsymbol\Sigma$ is. Hotelling’s $T^2$ is a function of $\mathbf{Q}$ alone. So, after fixing $\mathbf{H}$ and the rest of the data, the only thing left random is $\mathbf{Q}$ , whose distribution is known, and the selective $p$ -value is the fraction of the matrices $\mathbf{Q}$ that reproduce the clusters and give a larger $T^2$ than the observed one.

The three p-values

The package computes three selective $p$ -values. All three are exact under the null hypothesis, whatever $\boldsymbol\Sigma$ is; they differ in what they condition on, and the more they condition on, the less power they have.

$p$ -value	Function	Clustering	Conditions on	Power
$p_{\mathrm{polar}}$	`selective_hotelling_test()`, `selective_hotelling_bc()`	any	the clusters returned	highest
$p_{\mathrm{polar,1d}}$	`selective_truncated_F(method = "bc")`	any	also all but one direction of $\mathbf{Q}$	lower
$p_{\mathrm{polar,1d}}^{\mathrm{traj}}$	`selective_truncated_F(method = "exact")`	k-means	also the whole k-means trajectory	lowest

Use $p_{\mathrm{polar}}$ by default. The one-dimensional $p_{\mathrm{polar,1d}}$ samples a single angle instead of the whole matrix $\mathbf{Q}$ , so it is faster, and $p_{\mathrm{polar,1d}}^{\mathrm{traj}}$ has a closed form and needs no sampling at all (the one-dimensional tutorial).

How the main p-value is computed

The set of matrices $\mathbf{Q}$ that reproduce the clusters has no closed form. The package walks through it with a Markov chain: each step rotates two rows of $\mathbf{Q}$ by a random angle, reruns the clustering, and keeps the rotation only if the tested clusters are returned again. The observed $\mathbf{Q}$ is placed at a random position on a path of M + 1 points, and the $p$ -value is its rank among them (Besag and Clifford, 1989). This $p$ -value is exact for any path length M and any number of steps L between points: a short path gives a valid test with less power.

Choosing `M` and `L`

The defaults, M = 5000 and L = 2, are what we use in the paper.
Smaller M keeps the test valid but lowers power. M = 1000 is enough for exploration.
The clustering is rerun L times at each of the M steps, so the computing time is roughly M * L runs of the clustering algorithm.
M and L must be chosen before looking at the data; they must not be tuned on the observed clusters.

Requirements and limitations

Gaussian rows. The guarantee assumes Gaussian data. Heavy tails (for example, $t$ distributions) can inflate the Type I error, especially after average-linkage hierarchical clustering.
Enough observations. The test needs $r \ge p + 2$ , as Hotelling’s $T^2$ does. For high-dimensional data, reduce the dimension first, for example with principal components computed from data that are not tested.
A deterministic clustering function. cl_fun must return the same labels each time it is run on the same data (the tutorial for other clustering algorithms).
The pair to test. The pair should be fixed in advance, or every pair should be tested. Choosing the pair by looking at the data, for example the two clusters that look most different, needs further adjustment.
Large clusters. When the two clusters contain hundreds of observations, each step rotates only two of the $r - 1$ coordinates, and a path of fixed length moves less; the test stays valid but can lose power.

References

Besag, J. and Clifford, P. (1989). Generalized Monte Carlo significance tests. Biometrika, 76(4), 633–642.

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

Gao, L. L., Bien, J., and Witten, D. (2024). Selective inference for hierarchical clustering. Journal of the American Statistical Association, 119(545), 332–342.