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 only the shape of the data under that
method. matchednull separates the two.
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 preserved to within sampling error;
- 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, those clusters were the data’s shape, not its people.
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
#> [1] TRUE
round(cor(x) - cor(twin), 2) # correlations: close
#> [,1] [,2] [,3]
#> [1,] 0.00 0.01 0.05
#> [2,] 0.01 0.00 -0.03
#> [3,] 0.05 -0.03 0.00The 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.
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)
#> Matched-null test (30 Gaussian null twins)
#> real statistic: 1
#> null interval: [1, 1]
#> p (real >= nulls): 1
#> verdict: null-like (within the twins' interval)
# 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)
#> Matched-null test (30 Gaussian null twins)
#> real statistic: 2
#> null interval: [1, 1]
#> p (real >= nulls): 0.032
#> verdict: exceeds the null (beyond margins + covariance)The first call returns null-like: whatever clustering the pipeline reports is reproduced by twins with no types. The second returns exceeds the null: the grouping lives in structure the twins cannot carry.
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.
Stress-testing an exceedance: heavier tails
“Exceeds the null” under Gaussian twins licenses 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.
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 tailsThe Gaussian twins remain the test. A t copula could not take that
role, since it would absorb into the null the tail structure the test is
meant to detect, and that is exactly what makes it useful as a second,
deliberately blinded null: an exceedance that survives twins which have
already matched tail co-movement away cannot be attributed to tail
co-movement. Run this rung only on an exceedance the Gaussian run has
already licensed; if the Gaussian twins were not exceeded, the count is
what margins and correlations alone produce 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 verdict 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 was tail dependence, not types. Either way the verdict is sharper than what a single null could give. Whether a surviving exceedance is categorical is a separate question, answered by the taxometric checks rather than by any null.
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 article on the package website, which is kept current between releases.