---
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
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.

```{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
```

## 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.

```{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 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.

```{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 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](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: *Types Without Taxa: A
Covariance-Matched-Null Multiverse Test of Categorical versus Continuous
Personality Structure*.