Add matched-null validation for cluster counts

[haomeng797-ship-it]

Closes easystats/easystats#479.

This adds test_clusters(), which checks whether a selected number of clusters stands out from what the same selection method gives on matched-null "twins" of the data. A twin keeps every variable's exact values, has approximately the same correlations, and has no groups built in. The twins come from the matchednull package.

Usage

• By default the data are standardized and mclust::Mclust() picks the number of components by BIC.
• Any other selection method can be passed as cluster_function, as long as it returns one number. The same method is rerun on every twin, since what's being tested is the selection procedure.
• Options such as the copula family and parallel evaluation go through ... to matchednull::matched_null_test().

Reading the result

The print method reports one of three outcomes:

• No detected excess: the count is inside the twins' range. This isn't evidence of groups, and it isn't evidence against them either.
• Departure from the reference model: the count is above the range. The data have something the twins don't, which may or may not be groups.
• Below the reference range: the count is below the range. This calls for diagnosis and does not, by itself, establish anything about whether groups exist.

The twins are one specific no-groups reference, so this is a diagnostic against that reference. It isn't a general test of whether groups exist. The p-value is a Monte Carlo comparison against twins built from the same data, not an exact p-value. The documentation says both.

Checks

• Every value the clustering function returns, on the real data and on each twin, is validated. Anything other than a single finite number stops with an error.
• tests/testthat/test-test_clusters.R against the installed package: 56 passed, 0 failed, 0 skipped, with matchednull 0.2.1 (CRAN) and with the 0.2.2 patch. This includes a test that runs the default pipeline in a fresh R session without attaching mclust.
• Under devtools::test() with an older performance installed, that fresh-session test skips by design (54 passed, 1 skipped).
• lintr and air format --check are clean on the changed files.

[Copilot]

Pull request overview

Adds a new public helper, test_clusters(), to validate an observed (or user-supplied) cluster count against matched-null reference data generated via matchednull, with a default mclust::Mclust() (BIC) pipeline and support for custom scalar-returning cluster functions.

Changes:

• Introduces test_clusters() (exported) with roxygen documentation and matched-null integration.
• Adds a focused test suite covering custom functions, reproducibility, controls, argument forwarding, and input validation.
• Updates package metadata/documentation (NEWS entry, DESCRIPTION Suggests, WORDLIST, Rd generation, NAMESPACE export).

Reviewed changes

Copilot reviewed 6 out of 7 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
R/test_clusters.R	Adds the new test_clusters() API plus internal validation/standardization helpers.
tests/testthat/test-test_clusters.R	Adds test coverage for the new function’s expected behavior and input checks.
NEWS.md	Documents the new user-facing function in the changelog.
NAMESPACE	Exports test_clusters() as part of the public API.
man/test_clusters.Rd	Generated documentation for test_clusters().
DESCRIPTION	Adds matchednull (>= 0.2.1) to Suggests for the new functionality/tests.
inst/WORDLIST	Adds “Meng” to the spelling whitelist for the new reference entry.

Files not reviewed (1)

• man/test_clusters.Rd: Generated file

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[strengejacke]

Thanks, looks good to me! A version bump (4th digits, to increase dev version) in the DESCRIPTION is missing (or overwritten by merge main into this branch), else I have no other comments. @DominiqueMakowski anything you would like to add?

[haomeng797-ship-it]

Thanks! I bumped the development version to 0.17.1.11.

[DominiqueMakowski]

This is a bit clunky

#' # Any scalar-returning clustering pipeline can be tested.
#' pick_two <- function(data) 2
#' test_clusters(iris[, 1:4], cluster_function = pick_two, iterations = 19)

Rather than having the whole Mclust complication (which is used to get a number), why not streamline the api and just have the n_clusters argument that takes either a number, or a list of numbers to test multiple solutions.

The pipeline could then be like

n <- n_clusters(data)
# this suggests 3
test_clusters(data, n_clusters = c(2, 3, 5))

or am i missing somethign

[haomeng797-ship-it]

Thanks for pointing that out. I agree that the pick_two() example is awkward. The function argument is still needed, though, because the same cluster-selection procedure has to run on every matched-null dataset. Supplying fixed values through n_clusters would not specify how a cluster count should be obtained from each null dataset.

I’ll replace the example with a real clustering pipeline and clarify this in the documentation.

[DominiqueMakowski]

You mean the function needs to be aware of the clustering method?

Then why not (I'm just brainstorming I dont have strong opinions) making the function a method that runs on a cluster_analysis() result? So the workflow would go like

n_clusters(data)
Rez <- cluster_analysis(...l

Test_clusters(Rez)
``|

It would actually test against a specific solution 

[haomeng797-ship-it]

Yes that’s what I mean. test_clusters() needs some way to rerun the clustering method on each matched-null dataset, rather than only knowing the number selected for the original data.

I like the idea of making it a method for a cluster_analysis() result. That would probably feel more natural here. We’d just need the result object to retain the method and arguments, so the same analysis can be rerun on the null datasets.

The only thing I want to preserve is that the current test evaluates the selection process, rather than just one fixed solution.
I’ll look at what cluster_analysis() stores and see how the two could fit together!

[DominiqueMakowski]

I’ll look at what cluster_analysis() stores

if need be we can add more attributes to store more info :) we do that quite commonly across the easyverse

[haomeng797-ship-it]

I had a look. cluster_analysis() already stores the method, processed data, model, clusters, and performance, so adding the original arguments should be straightforward.

I think the only catch is a fixed n. If rez was fitted with n = 3, running the same thing on every matched-null dataset would just give 3 each time. For the current test, whatever chooses n needs to run again on each dataset so that there’s actually a null distribution.

So test_clusters(rez) could work if rez also remembers how n was chosen. If you mean testing one fixed solution, we’d probably need to compare something like separation or fit instead. Is that roughly what you were thinking?

[DominiqueMakowski]

do you mean that test_clusters need some alternative scenarios (ie different n solutions)?

If so, would something like:

test_clusters(rez, against = c(1, 2, 5)) work?

---

I reread you message again, and now I am starting to doubt my understanding (sorry i should have done my homework before 😬): do you mean that the matched-null validation is not about testing a given solution, but rather testing a selection procedure? Like if parameters::n_clusters_elbow() returns 3, the idea is to see how clear this selection method's result is on the data vs surrogate clusterless data?

If so... then perhaps test_clusters() should be used within n_clusters()? or as a method for specific cluster number selection methods? I'm confusing myself

[haomeng797-ship-it]

Yes, your second reading is right! 😀 I think my earlier wording mixed up two possible uses of matched_null_test().

For this PR, the statistic is the number of clusters selected. We run the same selector on the real data and on every matched-null dataset. For example, if elbow selects 3 on the real data but usually selects 1 on the null data, then 3 is unusual. If it also often selects 3 or more on the null data, then the result is not much evidence for real clusters.

So there is no need for the user to provide alternative values of n, and against = c(1, 2, 5) is not what I meant.

I think the cleanest API would be something like test_clusters(data, method = "elbow"). It could use the same selection methods as n_clusters(), but remain separate because the simulations can take a while.

Testing a fixed k with something like silhouette width would be another possible test, but that is not what this PR currently does!
Does that sound like the right direction?

[haomeng797-ship-it]

Quick update: this discussion made me go back and check the null generator before changing the API. I found a dependence-matching issue for non-Gaussian margins in matchednull, so I want to finish validating the correction before this is merged.

I'm going to pause this PR for now and come back with the updated calibration. Nothing needed from you in the meantime 🙂

[haomeng797-ship-it]

I opened a separate discussion for the statistical side of this, so we can keep this PR focused on the performance API: easystats/easystats#481

[strengejacke]

I'll update the main branch into this PR, so please pull commits before continue working locally :-)

#945 will also be merged soon.

[haomeng797-ship-it]

Quick update after a couple of weeks of validation: I’ve now settled on correlation-calibrated NORTA as the main direction for the reference generator. The broader checks also uncovered a separate calibration issue for some clustering statistics, so I’m keeping the PR in draft while I work through that part.

The good news is that I’m no longer rethinking the reference-model direction itself. What remains is mainly figuring out what calibration is needed and how cautiously the inferential claims should be stated.

I’ll keep working through that over the next few weeks and update the PR once the results are a bit more settled. Thanks again for the earlier discussion!

[haomeng797-ship-it]

Okay, finally back with this one! It took a lot longer than I expected. Every time I went to finish the API, the stats side turned up something I wanted to sort out first, and I didn't want to hand you something I'd have to walk back later.

The short version: it's ready for you to look at now. I just pushed a round of changes, it works with the matchednull that's already on CRAN, and all the test_clusters tests pass. I'm keeping it as a draft for now, for two reasons I'll explain at the bottom: a small matchednull patch, and a question for you about what kind of tool this should be.

Also, a correction to my last comment: I said I'd probably switch matchednull over to a calibrated generator. I ended up not doing that, and the generator stays exactly as it is on CRAN. What I actually spent the time on was pinning down what a result from this function can and can't tell you, and making the function harder to misuse.

What it does (the way I'd explain it to someone new)

Say you run a cluster analysis and your method picks 2 clusters. The question test_clusters() asks is basically: would the same method also pick 2 if there were no groups in the data?

To check that, it makes "twins" of your data. A twin has exactly the same values on every variable as your real data, so every histogram is identical. The correlations between variables are roughly the same too, and there are no groups built in. Then it runs your cluster-picking method on the real data and on a bunch of twins, and looks at whether your real number stands out.

One thing to keep in mind: the twins are one specific kind of "no groups" data. They can also be missing other ways your variables hang together (anything beyond the kind of dependence a multivariate normal distribution has), so if your data stand out from the twins, it isn't necessarily because of groups.

Here's iris (only 19 twins so the example runs fast; the default is 200):

set.seed(42)
test_clusters(iris[, 1:4], iterations = 19, n_max = 4)

# Matched-null test of the cluster count

  Observed statistic:                 2
  Reference interval:                 [2, 4] (2.5% and 97.5% quantiles of the twins)
  MC exceedance (twins >= observed):  1.000
  Reference model:                    Gaussian copula, 19 matched-null twins

No detected excess: the count lies within the reference interval.

So mclust picks 2 clusters on iris, but it also picks 2 to 4 on the twins, which have no groups in them. The 2 on its own isn't telling you anything the twins wouldn't also tell you.

You get one of three results:

• No detected excess (inside the twins' range, like above). Your method finds just as many clusters in the twins, so this count isn't evidence of groups. It also doesn't mean there are no groups. The test just didn't find anything.
• Departure from the reference model (above the range). Your data have something the twins don't. That could be groups! But it could also be other stuff, like really heavy tails or lopsided relationships between variables. So I'd treat it as "worth a closer look" rather than "the groups are real".
• Below the reference range (below the range). Your method finds fewer clusters in the real data than in the twins. If that happens, it calls for diagnosis (did the models converge, etc.) and does not, by itself, establish anything about whether groups exist.

Where it sits next to the other clustering functions

The way I think about the workflow:

1. check_clusterstructure(): is there any clustering tendency at all? (This one compares against pure random noise, so correlated data can come out as "clusterable" there even without any groups.)
2. n_clusters() or your own method: how many clusters?
3. cluster_analysis(): fit it.
4. test_clusters(): is that number more than you'd get from twins of your data?

So it doesn't replace the Hopkins check. They're asking different questions.

Any method works as long as you wrap it in a function that returns one number. That's why I kept cluster_function as a general argument: the thing being tested is the picking procedure itself, so the exact same procedure has to run again on every twin. Dominique, you were right that pick_two was a weird example, so the docs now use k-means + silhouette instead:

select_k <- function(data) {
  ks <- 2:5
  sil <- vapply(ks, function(k) {
    cl <- stats::kmeans(data, centers = k, nstart = 5)$cluster
    mean(cluster::silhouette(cl, stats::dist(data))[, "sil_width"])
  }, numeric(1))
  ks[which.max(sil)]
}
test_clusters(iris[, 1:4], cluster_function = select_k)

A tidyclust workflow would go in the same way, just wrapped in a function like this.

What's in this push

• It now stops if something goes wrong. If your function returns anything other than a single number (NA, text, two numbers, NULL...), on the real data or on any twin, you get an error instead of a result. Turns out matchednull 0.2.1 quietly drops NA twins from the p-value but still counts them in the total, so a run where most of the twins failed could still print a normal-looking result. Not great! The check lives on the performance side, so it works with the CRAN version too.
• Its own print method, with the three results spelled out.
• Docs that say what a result does and doesn't mean. They also suggest looking at each variable's histogram first with check_multimodal(). If groups are really far apart, it shows up in the histograms. Since the twins copy the histograms exactly, the twins can produce similar selected cluster counts, and then the test may not flag those groups.
• Kept the mclustBIC line Copilot flagged as unused. I know it looks like dead code, but Mclust() actually goes looking for mclustBIC inside whatever function called it, so without that line the default method breaks whenever mclust isn't loaded with library(mclust). I added a test that runs this in a fresh R session, so it'll get caught if someone removes it later.
• Bumped the dev version to 0.18.2.1 and moved the NEWS entry back up under (devel), since the merges from main had left it sitting under 0.18.0.
• Added cluster to Suggests for the example.

With the package installed, all 56 test_clusters tests pass, both with matchednull 0.2.1 from CRAN and with the 0.2.2 patch, and that includes the fresh-session one. I also checked it with two parallel workers.

Two caveats I want to be upfront about (both in the docs now)

• The twins only approximately match the correlations in your data. With skewed variables the difference can be systematic, not just random noise.
• The p-value compares your data against twins that were built from your own data, so it's a Monte Carlo comparison rather than an exact p-value with a guaranteed error rate. I'd read the result as a diagnostic ("does my number of clusters stand out?") and not as a formal significance test.

Neither of these changes how the function runs, just how you should read what it prints.

What's next

I've also made a matchednull 0.2.2 that fixes the NA problem at the source and cleans up the same wording in its docs. I'll send that to CRAN on its own, and once it's through I'll bump the minimum version here.

That patch fixes a software bug, though. The bigger statistical question is still open: how well the p-value is calibrated in general. So I'd really like you to look at this as a diagnostic. It compares a cluster count against one specific no-groups reference, and it can't tell you in general whether groups exist. If that sounds like something you'd want in performance, great! If it needs to be framed differently first, or you're not sure it fits, I'd rather talk about that before this leaves draft.

And if the print output or argument names should look more like the rest of performance, just let me know and I'll change them 🙂

[codecov]

Codecov Report

❌ Patch coverage is 96.26168% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 64.99%. Comparing base (961ba2e) to head (bd8859c).
⚠️ Report is 24 commits behind head on main.

Files with missing lines	Patch %	Lines
R/test_clusters.R	96.26%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #938      +/-   ##
==========================================
+ Coverage   64.07%   64.99%   +0.91%     
==========================================
  Files          94       98       +4     
  Lines        8321     8529     +208     
==========================================
+ Hits         5332     5543     +211     
+ Misses       2989     2986       -3     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[strengejacke]

Merging main into this branch, to resolve conflicts. Next week I'll read your last comment, to get up-to-date on this PR. Thanks a lot for the work so far!