diff --git a/DESCRIPTION b/DESCRIPTION
index 77bf2fd58..e9e4634ea 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -126,6 +126,7 @@ Suggests:
     loo,
     MASS,
     Matrix,
+    matchednull (>= 0.2.1),
     mclogit,
     mclust,
     metadat,
diff --git a/NAMESPACE b/NAMESPACE
index 3f6718881..6b60377c3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -690,6 +690,7 @@ export(r2_zeroinflated)
 export(rmse)
 export(simulate_residuals)
 export(test_bf)
+export(test_clusters)
 export(test_likelihoodratio)
 export(test_lrt)
 export(test_performance)
diff --git a/NEWS.md b/NEWS.md
index bda1be394..e90648ce6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,10 @@
 
 ## Changes
 
+* New function `test_clusters()` uses the `matchednull` package to test whether
+  a selected cluster count exceeds what matched-null data with the same margins
+  and correlations would produce.
+
 * `check_group_variation()` now returns a numeric effect size of the grouping
   variable's predictive association strength.
 
diff --git a/R/test_clusters.R b/R/test_clusters.R
new file mode 100644
index 000000000..a964b96ce
--- /dev/null
+++ b/R/test_clusters.R
@@ -0,0 +1,155 @@
+#' Test a clustering result against matched nulls
+#'
+#' `test_clusters()` compares the number of clusters selected from observed data
+#' with the numbers selected from matched-null versions of those data. The null
+#' data preserve the observed marginal distributions and correlation structure
+#' but contain no cluster structure by construction. This complements
+#' [`check_clusterstructure()`], which evaluates clustering tendency against a
+#' spatial-uniformity null. The matched-null comparison is implemented with
+#' [`matchednull::matched_null_test()`].
+#'
+#' @param x A numeric matrix or data frame with no missing or infinite values.
+#' @param cluster_function A function that takes a numeric matrix and returns one
+#'   number. By default, `mclust::Mclust()` selects the number of mixture
+#'   components by BIC.
+#' @param iterations Number of matched-null data sets to evaluate.
+#' @param n_max Maximum number of clusters considered by the default
+#'   `mclust::Mclust()` pipeline. Ignored when `cluster_function` is supplied.
+#' @param standardize Logical. If `TRUE`, variables are standardized before the
+#'   observed and matched-null cluster counts are computed.
+#' @param ... Additional arguments passed to
+#'   [`matchednull::matched_null_test()`], such as `copula`, `df`, `probs`,
+#'   `ridge`, or `parallel`.
+#'
+#' @return An object of class `"test_clusters"` and `"matched_null_test"`. It
+#'   contains the observed statistic (`real`), matched-null statistics (`null`),
+#'   null interval (`interval`), one-sided Monte Carlo p-value (`p_exceed`), and
+#'   interval verdict (`within`).
+#'
+#' @details
+#' The default cluster-count pipeline fits Gaussian mixture models for one to
+#' `n_max` components and returns the BIC-selected count. A custom
+#' `cluster_function` receives the same preprocessed data as the default
+#' pipeline and must return one non-missing number.
+#'
+#' The Gaussian matched null tests for structure beyond the observed margins
+#' and correlations; rejection alone does not establish the existence of
+#' discrete types. When tail dependence is plausible, rerun the test with
+#' `copula = "t"` as a sensitivity analysis.
+#'
+#' @examplesIf requireNamespace("matchednull", quietly = TRUE) && requireNamespace("mclust", quietly = TRUE)
+#' \donttest{
+#' set.seed(42)
+#' test_clusters(iris[, 1:4], iterations = 19, n_max = 4)
+#'
+#' # Any scalar-returning clustering pipeline can be tested.
+#' pick_two <- function(data) 2
+#' test_clusters(iris[, 1:4], cluster_function = pick_two, iterations = 19)
+#' }
+#'
+#' @references
+#' Meng, M. (2026). Types Without Taxa: A Covariance-Matched-Null Multiverse
+#' Test of Categorical Versus Continuous Personality Structure. Manuscript
+#' under review. \doi{10.17605/OSF.IO/2EKCG}
+#'
+#' @seealso [`matchednull::matched_null_test()`],
+#'   [`check_clusterstructure()`]
+#' @export
+test_clusters <- function(
+  x,
+  cluster_function = NULL,
+  iterations = 200,
+  n_max = 10,
+  standardize = TRUE,
+  ...
+) {
+  insight::check_if_installed("matchednull", minimum_version = "0.2.1")
+
+  if (!is.logical(standardize) || length(standardize) != 1L || is.na(standardize)) {
+    insight::format_error("`standardize` must be `TRUE` or `FALSE`.")
+  }
+  if (!is.null(cluster_function) && !is.function(cluster_function)) {
+    insight::format_error("`cluster_function` must be a function or `NULL`.")
+  }
+  uses_default <- is.null(cluster_function)
+  if (!.is_positive_integer(iterations)) {
+    insight::format_error("`iterations` must be a positive integer.")
+  }
+  if (uses_default && !.is_positive_integer(n_max)) {
+    insight::format_error("`n_max` must be a positive integer.")
+  }
+
+  x <- .validate_test_clusters_data(x)
+  x <- .standardize_test_clusters_data(x, standardize)
+
+  if (uses_default) {
+    insight::check_if_installed("mclust")
+    mclustBIC <- mclust::mclustBIC
+    max_components <- min(as.integer(n_max), nrow(x) - 1L)
+    cluster_function <- function(data) {
+      mclust::Mclust(
+        data,
+        G = seq_len(max_components),
+        verbose = FALSE
+      )$G
+    }
+  }
+
+  out <- matchednull::matched_null_test(
+    x,
+    cluster_fn = cluster_function,
+    R = as.integer(iterations),
+    ...
+  )
+
+  attr(out, "standardize") <- standardize
+  attr(out, "n_max") <- if (uses_default) max_components else NULL
+  class(out) <- c("test_clusters", class(out))
+  out
+}
+
+
+.validate_test_clusters_data <- function(x) {
+  if (!is.matrix(x) && !is.data.frame(x)) {
+    insight::format_error("`x` must be a numeric matrix or data frame.")
+  }
+  if (is.data.frame(x) && !all(vapply(x, is.numeric, logical(1)))) {
+    insight::format_error("All columns in `x` must be numeric.")
+  }
+
+  x <- as.matrix(x)
+  if (!is.numeric(x)) {
+    insight::format_error("`x` must contain only numeric values.")
+  }
+  if (nrow(x) < 2L || ncol(x) < 1L) {
+    insight::format_error("`x` must contain at least two rows and one column.")
+  }
+  if (anyNA(x) || !all(is.finite(x))) {
+    insight::format_error("`x` must not contain missing or infinite values.")
+  }
+  x
+}
+
+
+.standardize_test_clusters_data <- function(x, standardize) {
+  if (!standardize) {
+    return(x)
+  }
+
+  column_sd <- apply(x, 2L, stats::sd)
+  if (!all(is.finite(column_sd)) || any(column_sd == 0)) {
+    insight::format_error(
+      "`x` cannot contain constant columns when `standardize = TRUE`."
+    )
+  }
+  scale(x)
+}
+
+
+.is_positive_integer <- function(x) {
+  is.numeric(x) &&
+    length(x) == 1L &&
+    !is.na(x) &&
+    x >= 1L &&
+    x %% 1 == 0
+}
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 847ab8bbe..c13d5596f 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -155,6 +155,7 @@ Matuschek
 McElreath
 McKelvey
 Mehrvarz
+Meng
 Merkle
 Methoden
 Michalos
diff --git a/man/test_clusters.Rd b/man/test_clusters.Rd
new file mode 100644
index 000000000..abf5d2e56
--- /dev/null
+++ b/man/test_clusters.Rd
@@ -0,0 +1,81 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/test_clusters.R
+\name{test_clusters}
+\alias{test_clusters}
+\title{Test a clustering result against matched nulls}
+\usage{
+test_clusters(
+  x,
+  cluster_function = NULL,
+  iterations = 200,
+  n_max = 10,
+  standardize = TRUE,
+  ...
+)
+}
+\arguments{
+\item{x}{A numeric matrix or data frame with no missing or infinite values.}
+
+\item{cluster_function}{A function that takes a numeric matrix and returns one
+number. By default, \code{mclust::Mclust()} selects the number of mixture
+components by BIC.}
+
+\item{iterations}{Number of matched-null data sets to evaluate.}
+
+\item{n_max}{Maximum number of clusters considered by the default
+\code{mclust::Mclust()} pipeline. Ignored when \code{cluster_function} is supplied.}
+
+\item{standardize}{Logical. If \code{TRUE}, variables are standardized before the
+observed and matched-null cluster counts are computed.}
+
+\item{...}{Additional arguments passed to
+\code{\link[matchednull:matched_null_test]{matchednull::matched_null_test()}}, such as \code{copula}, \code{df}, \code{probs},
+\code{ridge}, or \code{parallel}.}
+}
+\value{
+An object of class \code{"test_clusters"} and \code{"matched_null_test"}. It
+contains the observed statistic (\code{real}), matched-null statistics (\code{null}),
+null interval (\code{interval}), one-sided Monte Carlo p-value (\code{p_exceed}), and
+interval verdict (\code{within}).
+}
+\description{
+\code{test_clusters()} compares the number of clusters selected from observed data
+with the numbers selected from matched-null versions of those data. The null
+data preserve the observed marginal distributions and correlation structure
+but contain no cluster structure by construction. This complements
+\code{\link[=check_clusterstructure]{check_clusterstructure()}}, which evaluates clustering tendency against a
+spatial-uniformity null. The matched-null comparison is implemented with
+\code{\link[matchednull:matched_null_test]{matchednull::matched_null_test()}}.
+}
+\details{
+The default cluster-count pipeline fits Gaussian mixture models for one to
+\code{n_max} components and returns the BIC-selected count. A custom
+\code{cluster_function} receives the same preprocessed data as the default
+pipeline and must return one non-missing number.
+
+The Gaussian matched null tests for structure beyond the observed margins
+and correlations; rejection alone does not establish the existence of
+discrete types. When tail dependence is plausible, rerun the test with
+\code{copula = "t"} as a sensitivity analysis.
+}
+\examples{
+\dontshow{if (requireNamespace("matchednull", quietly = TRUE) && requireNamespace("mclust", quietly = TRUE)) withAutoprint(\{ # examplesIf}
+\donttest{
+set.seed(42)
+test_clusters(iris[, 1:4], iterations = 19, n_max = 4)
+
+# Any scalar-returning clustering pipeline can be tested.
+pick_two <- function(data) 2
+test_clusters(iris[, 1:4], cluster_function = pick_two, iterations = 19)
+}
+\dontshow{\}) # examplesIf}
+}
+\references{
+Meng, M. (2026). Types Without Taxa: A Covariance-Matched-Null Multiverse
+Test of Categorical Versus Continuous Personality Structure. Manuscript
+under review. \doi{10.17605/OSF.IO/2EKCG}
+}
+\seealso{
+\code{\link[matchednull:matched_null_test]{matchednull::matched_null_test()}},
+\code{\link[=check_clusterstructure]{check_clusterstructure()}}
+}
diff --git a/tests/testthat/test-test_clusters.R b/tests/testthat/test-test_clusters.R
new file mode 100644
index 000000000..926ea8cbe
--- /dev/null
+++ b/tests/testthat/test-test_clusters.R
@@ -0,0 +1,89 @@
+test_that("test_clusters supports a custom cluster function", {
+  skip_if_not_installed("matchednull", minimum_version = "0.2.1")
+
+  set.seed(12)
+  x <- matrix(rnorm(120), ncol = 3)
+  out <- test_clusters(
+    x,
+    cluster_function = function(data) 2,
+    iterations = 9
+  )
+
+  expect_s3_class(out, "test_clusters")
+  expect_s3_class(out, "matched_null_test")
+  expect_identical(out$real, 2)
+  expect_length(out$null, 9)
+  expect_true(out$within)
+  expect_true(attr(out, "standardize"))
+  expect_null(attr(out, "n_max"))
+  expect_output(print(out), "Matched-null test")
+
+  out_t <- test_clusters(
+    x,
+    cluster_function = function(data) 2,
+    iterations = 3,
+    copula = "t",
+    df = 4
+  )
+  expect_identical(out_t$copula, "t")
+  expect_identical(out_t$df, 4)
+})
+
+
+test_that("test_clusters is reproducible from a seed", {
+  skip_if_not_installed("matchednull", minimum_version = "0.2.1")
+
+  set.seed(3)
+  x <- matrix(rnorm(180), ncol = 3)
+  cluster_summary <- function(data) sum(data[1:5, 1])
+
+  set.seed(11)
+  out1 <- test_clusters(x, cluster_summary, iterations = 9)
+  set.seed(11)
+  out2 <- test_clusters(x, cluster_summary, iterations = 9)
+
+  expect_identical(out1$null, out2$null)
+  expect_gt(stats::var(out1$null), 0)
+})
+
+
+test_that("default mclust pipeline distinguishes null and positive controls", {
+  skip_if_not_installed("matchednull", minimum_version = "0.2.1")
+  skip_if_not_installed("mclust")
+
+  set.seed(42)
+  x_null <- matrix(rnorm(150 * 4), 150, 4) %*%
+    chol(diag(4) * 0.5 + 0.5)
+  set.seed(7)
+  null_result <- test_clusters(x_null, iterations = 9, n_max = 4)
+  expect_true(null_result$within)
+  expect_identical(attr(null_result, "n_max"), 4L)
+
+  set.seed(42)
+  group <- sample.int(2, 150, replace = TRUE)
+  x_positive <- matrix(rnorm(150 * 4), 150, 4)
+  positive_correlation <- chol(matrix(c(1, 0.9, 0.9, 1), 2, 2))
+  negative_correlation <- chol(matrix(c(1, -0.9, -0.9, 1), 2, 2))
+  x_positive[group == 1, 1:2] <- x_positive[group == 1, 1:2] %*% positive_correlation
+  x_positive[group == 1, 3:4] <- x_positive[group == 1, 3:4] %*% positive_correlation
+  x_positive[group == 2, 1:2] <- x_positive[group == 2, 1:2] %*% negative_correlation
+  x_positive[group == 2, 3:4] <- x_positive[group == 2, 3:4] %*% negative_correlation
+
+  set.seed(7)
+  positive_result <- test_clusters(x_positive, iterations = 9, n_max = 4)
+  expect_gt(positive_result$real, positive_result$interval[2])
+})
+
+
+test_that("test_clusters validates its inputs", {
+  skip_if_not_installed("matchednull", minimum_version = "0.2.1")
+
+  x <- matrix(rnorm(40), ncol = 2)
+  expect_error(test_clusters(1:10), "matrix or data frame")
+  expect_error(test_clusters(data.frame(x = 1:3, y = letters[1:3])), "numeric")
+  expect_error(test_clusters(cbind(x, NA_real_)), "missing or infinite")
+  expect_error(test_clusters(cbind(x, 1)), "constant columns")
+  expect_error(test_clusters(x, cluster_function = 1), "must be a function")
+  expect_error(test_clusters(x, iterations = 1.5), "positive integer")
+  expect_error(test_clusters(x, n_max = 0), "positive integer")
+})