--- title: "Testing a cluster count with matched nulls" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Testing a cluster count with matched nulls} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") has_mclust <- requireNamespace("mclust", quietly = TRUE) ``` ## The question A clustering method applied to continuous data will report some number of clusters whether or not any exist. So when an analysis announces that a dataset contains, say, four types, the number may describe the population, or it may be what the method would report for any data with the same margins and correlations. `matchednull` compares the reported number with what a specified reference model of exactly such data produces. ## The idea: a null twin `copula_null()` builds a synthetic twin of your data that keeps everything innocent and removes the one thing in question: * every variable's distribution is preserved **exactly** (the twin reuses the real values, reshuffled); * the correlation matrix is reproduced approximately: the observed correlations serve as the latent correlation of the copula, and restoring the non-normal margins by rank reordering moves the output Pearson correlations, so they come out close but not exact; * all remaining dependence is Gaussian, so the twin contains **no cluster structure by construction**. If your pipeline finds as many clusters in the twins as in the real data, there is no detected excess under this reference, which is not evidence that the clusters are absent. ```{r} library(matchednull) set.seed(1) x <- matrix(rnorm(400 * 3), 400, 3) %*% chol(matrix(c(1, .5, .3, .5, 1, .4, .3, .4, 1), 3, 3)) twin <- copula_null(x) all(sort(twin[, 1]) == sort(x[, 1])) # margins: identical round(cor(x) - cor(twin), 2) # correlations: close, not exact ``` ## The test `matched_null_test()` takes your data and **your own pipeline**, wrapped as a function that returns one number (typically the selected number of clusters). It runs the identical pipeline on the real data and on `R` twins, and asks whether the real answer stands out. The `p_exceed` it reports is a plug-in Monte Carlo exceedance quantity under a reference fitted to the same data, not an exact finite-sample p-value. ```{r, eval = has_mclust} suppressPackageStartupMessages(library(mclust)) pick_k <- function(d) Mclust(d, G = 1:4, modelNames = "VVV", verbose = FALSE)$G # 1. Typeless data: the test should stay quiet. set.seed(7) matched_null_test(x, pick_k, R = 30) ``` ```{r, eval = has_mclust} # 2. Two genuine types, hidden in the dependence structure # (identical margins, opposite correlation orientation). set.seed(42) z <- sample(2, 400, replace = TRUE) X <- matrix(rnorm(400 * 4), 400, 4) L1 <- chol(matrix(c(1, .85, .85, 1), 2, 2)) L2 <- chol(matrix(c(1, -.85, -.85, 1), 2, 2)) X[z == 1, 1:2] <- X[z == 1, 1:2] %*% L1 X[z == 1, 3:4] <- X[z == 1, 3:4] %*% L1 X[z == 2, 1:2] <- X[z == 2, 1:2] %*% L2 X[z == 2, 3:4] <- X[z == 2, 3:4] %*% L2 set.seed(7) matched_null_test(X, pick_k, R = 30) ``` The first call reports *no detected excess*: whatever clustering the pipeline reports is reproduced by twins with no types, which is not the same as evidence that types are absent. The second reports a *departure from the reference model*: the grouping lives in structure the twins cannot carry. That is a departure from a specified reference, not yet a finding of types; what kind of structure it is remains to be settled. ## Plugging in your own pipeline `cluster_fn` can wrap anything: a k-means heuristic, a published typology's exact workflow, or any scalar measure of clustering strength. The null is defined at the level of the data, not the pipeline. Two practical notes: * **Set a seed before calling** — reproducibility is deliberately left to the caller. * **Check the margins first.** The one regime the null cannot flag is types so separated that they are visible in the margins and covariance themselves (strongly bimodal variables); inspect the margins for pronounced multimodality before relying on the count test, and tie-break granular (Likert-type) scales before formal unimodality tests. ## A sensitivity variant: heavier tails A departure from the Gaussian twins means only *structure beyond margins and correlations*, not types. Heavy-tailed dependence, where extreme values across variables arrive together, also exceeds a Gaussian null, and real questionnaire and clinical data are heavy-tailed more often than not. To separate the two readings, rerun the test with t-copula twins: same margins, same correlations, but tails that co-move. ```{r, eval = FALSE} matched_null_test(x, pick_k, R = 200, copula = "t", df = 8) # moderate tails matched_null_test(x, pick_k, R = 200, copula = "t", df = 3) # heavy tails ``` The Gaussian twins remain the test. A t copula could not take that role, since it would absorb into the reference the tail structure the test is meant to detect, and that is exactly what makes it useful as a sensitivity variant: a departure that survives twins which have already matched tail co-movement away cannot be attributed to tail co-movement. Run this variant on a departure the Gaussian run has already shown; if the Gaussian twins were not exceeded there is no excess for the variant to probe, and the question does not arise. `df` is set by you and never estimated from the data, and both values should be reported rather than the one that suits the conclusion. A result quoted without its Gaussian outcome is not a matched-null result. The ladder reads: a result that exceeds the Gaussian twins *and* the t twins is hard to attribute to tails; a result the t twins reproduce is consistent with tail dependence rather than types. The t variant does not by itself establish or rule out clusters. Whether a surviving departure is *categorical* is a separate question, answered by the taxometric checks rather than by any reference model. ## Going further Why the construction works, what the positive controls show, how the method relates to earlier reference-distribution tests, and the one regime where the test declines to answer are all worked through in the [How it works](https://haomeng797-ship-it.github.io/matchednull/articles/how-it-works.html) article on the package website, which is kept current between releases. ## Reference The method, its positive controls, and its false-positive calibration are described in the accompanying paper: Meng (2026), *A Matched-Null Test of Whether the Clusters a Pipeline Reports Are Real, With an Application to Personality Types*, manuscript under review. Its personality application was preregistered as *Types Without Taxa* ().