25.3 From Exploration to Hypothesis

"I found seven clusters, named them all, and threw a conference paper at each one. My advisor said that is not how science works. I said that is exactly how science works, just faster."

A Silhouette Score Seeking Validation

1. Internal Cluster Validation

A single retracted paper can cost a lab years of follow-up work and millions in wasted funding, and a frequent trigger is a cluster that looked real on a projection plot but was never tested against independent evidence. The three-layer validation pipeline below exists to prevent exactly that failure.

Your algorithm just split a thousand patients into five clean subgroups and the UMAP projection looks stunning; now prove the clusters are real, without a single ground-truth label to check against. In Section 25.2, we compared clusterings to known labels using the Adjusted Rand Index (ARI), but in genuine exploratory discovery those labels do not exist. You need metrics that evaluate cluster quality from the data alone: internal validation metrics.

Internal validation metrics score how well a clustering partitions the data. They use only the data points and their assigned labels, with no reference to external ground truth. In real discovery settings, you rarely have labels to compare against. These metrics provide your primary quantitative evidence that the structure is geometrically coherent rather than an artifact of the algorithm's assumptions. Each metric compares some measure of within-cluster compactness against between-cluster separation; the silhouette coefficient uses average distances, Density-Based Clustering Validation (DBCV) uses density along minimum spanning trees, and Bayesian Information Criterion (BIC) uses likelihood under a probabilistic model. Use internal validation when you have no labels at all; when partial labels exist, prefer external metrics like ARI or Normalized Mutual Information (NMI), and treat internal scores as complementary evidence.

1.1 The Silhouette Coefficient

What. The silhouette coefficient measures how similar each point is to its own cluster compared to the nearest neighboring cluster. For a point \(\mathbf{x}_i\) assigned to cluster \(C_k\):

where \(a_i\) is the mean distance from \(\mathbf{x}_i\) to all other points in \(C_k\) (intra-cluster distance), and \(b_i\) is the minimum mean distance from \(\mathbf{x}_i\) to points in any other cluster (nearest-cluster distance). The silhouette ranges from \(-1\) (point is in the wrong cluster) to \(+1\) (point is far from neighboring clusters).

Why. The silhouette coefficient is intuitive: \(s_i > 0\) means the point is closer to its own cluster than to any other; \(s_i \approx 0\) means it sits on the boundary between clusters; \(s_i < 0\) means it is probably misassigned. The global silhouette score (mean over all points) summarizes overall clustering quality, while per-cluster silhouette distributions reveal which clusters are well-separated and which are fuzzy.

When. Use silhouette for any distance-based clustering method. It works well for convex clusters (k-means, Gaussian Mixture Models (GMM)) but can be misleading for density-based clusters with irregular shapes, where DBCV (below) is more appropriate.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples

data = load_breast_cancer()
X = StandardScaler().fit_transform(data.data)

# Evaluate k-means for k = 2, 3, 4, 5
results = []
for k in range(2, 6):
    km = KMeans(n_clusters=k, n_init=10, random_state=42)
    labels = km.fit_predict(X)
    sil = silhouette_score(X, labels)
    results.append({"k": k, "silhouette": sil, "inertia": km.inertia_})
    print(f"k={k}: silhouette={sil:.3f}, inertia={km.inertia_:.0f}")

# Per-cluster silhouette analysis for k=2
km2 = KMeans(n_clusters=2, n_init=10, random_state=42)
labels_2 = km2.fit_predict(X)
sample_sil = silhouette_samples(X, labels_2)

for cluster_id in range(2):
    mask = labels_2 == cluster_id
    cluster_sil = sample_sil[mask]
    print(f"\nCluster {cluster_id} (n={mask.sum()}):")
    print(f"  Mean silhouette: {cluster_sil.mean():.3f}")
    print(f"  Min silhouette:  {cluster_sil.min():.3f}")
    print(f"  Fraction < 0:    {(cluster_sil < 0).mean():.3f}")

k=2: silhouette=0.292, inertia=11349
k=3: silhouette=0.276, inertia=9310
k=4: silhouette=0.260, inertia=8074
k=5: silhouette=0.247, inertia=7217

Cluster 0 (n=175):
  Mean silhouette: 0.258
  Min silhouette:  -0.133
  Fraction < 0:    0.051

Cluster 1 (n=394):
  Mean silhouette: 0.307
  Min silhouette:  -0.106
  Fraction < 0:    0.033

1.2 DBCV: Validation for Density-Based Clusters

What. The DBCV index evaluates clusters based on the density connectivity within clusters versus the density separation between clusters. Unlike the silhouette coefficient, DBCV handles clusters of arbitrary shape by measuring density along the minimum spanning tree (the shortest-total-weight tree that connects every point in a cluster without cycles) of each cluster.

Why. Hierarchical DBSCAN (HDBSCAN), a density-based clustering algorithm that extends DBSCAN by building a hierarchy of clusterings at varying density thresholds, finds clusters as dense regions separated by sparse valleys. The silhouette coefficient, which uses mean distances, can penalize elongated clusters even when they are perfectly valid density-based partitions. DBCV respects the density semantics of the clustering by computing a density-based separation that penalizes only clusters whose boundaries cross high-density regions.

How. For each cluster, DBCV computes the minimum density within the cluster's minimum spanning tree (the internal density) and the maximum density between any pair of clusters (the external density). The DBCV score for a cluster is the ratio of separation density to internal density, normalized to [-1, 1]. The global score is the density-weighted mean across clusters.

import hdbscan  # or sklearn.cluster.HDBSCAN in scikit-learn >= 1.3

# HDBSCAN with validity index
clusterer = hdbscan.HDBSCAN(
    min_cluster_size=30,
    min_samples=10,
    gen_min_span_tree=True,  # needed for DBCV
)
labels_hdb = clusterer.fit_predict(X)

# DBCV is available as the relative_validity_ attribute
dbcv = clusterer.relative_validity_
n_clusters = len(set(labels_hdb)) - (1 if -1 in labels_hdb else 0)
n_noise = (labels_hdb == -1).sum()

print(f"HDBSCAN: {n_clusters} clusters, {n_noise} noise points")
print(f"DBCV score: {dbcv:.3f}")
print(f"Cluster persistence scores: {clusterer.cluster_persistence_}")  # persistence: how long each cluster survives across density thresholds

# Compare DBCV across different min_cluster_size values
print("\nDBCV sensitivity to min_cluster_size:")
for mcs in [15, 30, 50, 80]:
    c = hdbscan.HDBSCAN(min_cluster_size=mcs, min_samples=5,
                         gen_min_span_tree=True)
    c.fit(X)
    nc = len(set(c.labels_)) - (1 if -1 in c.labels_ else 0)
    print(f"  min_cluster_size={mcs:3d}: {nc} clusters, "
          f"DBCV={c.relative_validity_:.3f}")

HDBSCAN: 2 clusters, 62 noise points
DBCV score: 0.298
Cluster persistence scores: [0.297 0.118]

DBCV sensitivity to min_cluster_size:
  min_cluster_size= 15: 4 clusters, DBCV=0.221
  min_cluster_size= 30: 2 clusters, DBCV=0.298
  min_cluster_size= 50: 2 clusters, DBCV=0.314
  min_cluster_size= 80: 2 clusters, DBCV=0.281

Internal validation metrics confirm that clusters have good geometric properties, but a geometrically clean partition can still be an accident of which particular data points happened to land in your sample; the next question is whether the same clusters survive when the data are perturbed.

2. Stability Analysis: Bootstrap Resampling

A cluster that disappears when you remove 10% of the data was never really there. Stability analysis tests whether clusters are robust to perturbation by resampling the data and checking whether the same clusters reappear.

Mental Model

Think of stability analysis like checking whether a constellation is real or just a coincidence. If you are stargazing and spot a pattern of five bright stars, you can test whether the pattern is robust by squinting (blurring out some stars) or shifting your viewing angle slightly. If the same constellation is still recognizable each time, the pattern is stable. If squinting makes it vanish or rearrange into something different, you were connecting dots that happened to line up from one vantage point. Bootstrap resampling works the same way: each subsample is a slightly different "viewing angle" on your data, and the Adjusted Rand Index measures how recognizable the original constellation remains. A stable clustering survives dozens of squints; an unstable one was never more than a lucky arrangement.

from sklearn.metrics import adjusted_rand_score
from sklearn.utils import resample

def cluster_stability(X, clusterer_factory, n_bootstrap=50,
                      subsample_frac=0.8, random_state=42):
    """Assess cluster stability via bootstrap resampling.

    Returns mean and std of ARI between full-data and
    subsample clusterings.
    """
    rng = np.random.RandomState(random_state)

    # Reference clustering on full data
    ref_clusterer = clusterer_factory()
    ref_labels = ref_clusterer.fit_predict(X)

    ari_scores = []
    for i in range(n_bootstrap):
        # Subsample without replacement
        n_sub = int(len(X) * subsample_frac)
        idx = rng.choice(len(X), size=n_sub, replace=False)
        X_sub = X[idx]

        # Cluster the subsample
        sub_clusterer = clusterer_factory()
        sub_labels = sub_clusterer.fit_predict(X_sub)

        # Compare with reference labels on the same points
        ari = adjusted_rand_score(ref_labels[idx], sub_labels)
        ari_scores.append(ari)

    return np.mean(ari_scores), np.std(ari_scores)

# Test stability for each algorithm
algorithms = {
    "k-means (k=2)": lambda: KMeans(n_clusters=2, n_init=10,
                                     random_state=None),
    "HDBSCAN": lambda: hdbscan.HDBSCAN(min_cluster_size=30,
                                         min_samples=10),
}

print("Cluster stability (50 bootstrap resamples, 80% subsample):")
for name, factory in algorithms.items():
    mean_ari, std_ari = cluster_stability(X, factory)
    print(f"  {name}: ARI = {mean_ari:.3f} +/- {std_ari:.3f}")

Cluster stability (50 bootstrap resamples, 80% subsample):
  k-means (k=2): ARI = 0.892 +/- 0.034
  HDBSCAN: ARI = 0.831 +/- 0.058

3. Enrichment Analysis: Connecting Clusters to Meaning

Internal validation confirms geometric structure; stability analysis confirms robustness; neither confirms meaning. Enrichment analysis tests whether cluster membership correlates with external variables: metadata, experimental conditions, known labels, or domain-specific annotations.

The core idea is simple: for each cluster, test whether any external variable is overrepresented or underrepresented compared to the dataset as a whole. For categorical variables, use a chi-squared test (which compares observed category frequencies against expected frequencies under the null hypothesis of no association) or Fisher's exact test (a non-approximation alternative to chi-squared that is preferred when sample sizes are small). For continuous variables, use a t-test or Mann-Whitney U test (a non-parametric alternative that compares rank orderings rather than assuming normal distributions). Correct for multiple testing with Bonferroni (which divides the significance threshold by the number of tests) or Benjamini-Hochberg (which controls the expected proportion of false discoveries rather than the probability of any false discovery).

from scipy import stats

def enrichment_analysis(df, cluster_labels, metadata_cols,
                        alpha=0.05):
    """Test whether metadata variables are enriched in clusters.

    Returns a DataFrame of significant associations after
    Bonferroni correction.
    """
    results = []
    n_clusters = len(set(cluster_labels)) - (
        1 if -1 in cluster_labels else 0)
    # Bonferroni correction
    n_tests = n_clusters * len(metadata_cols)

    for cluster_id in range(n_clusters):
        mask = cluster_labels == cluster_id
        for col in metadata_cols:
            values = df[col]
            if values.dtype.name == "category" or values.nunique() < 10:
                # Categorical: chi-squared test
                observed = values[mask].value_counts()
                expected_frac = values.value_counts(normalize=True)
                n_in = mask.sum()
                expected = expected_frac * n_in

                # Align indices
                all_cats = expected_frac.index
                obs = np.array([observed.get(c, 0) for c in all_cats])
                exp = np.array([expected.get(c, 0) for c in all_cats])

                if (exp > 0).all():
                    chi2, p_val = stats.chisquare(obs, exp)
                    effect = dict(zip(all_cats,
                                      (obs / exp).round(2)))
                else:
                    continue
            else:
                # Continuous: t-test
                in_cluster = values[mask].values
                out_cluster = values[~mask].values
                t_stat, p_val = stats.ttest_ind(
                    in_cluster, out_cluster, equal_var=False)
                effect = {
                    "cluster_mean": float(np.mean(in_cluster)),
                    "other_mean": float(np.mean(out_cluster)),
                    "cohens_d": float(
                        (np.mean(in_cluster) - np.mean(out_cluster))
                        / np.sqrt((np.var(in_cluster)
                                   + np.var(out_cluster)) / 2)
                    ),
                }

            p_adj = min(p_val * n_tests, 1.0)  # Bonferroni
            if p_adj < alpha:
                results.append({
                    "cluster": cluster_id,
                    "variable": col,
                    "p_adjusted": p_adj,
                    "effect": effect,
                })

    return results

# Add diagnosis as metadata for enrichment analysis
import pandas as pd
df = pd.DataFrame(data.data, columns=data.feature_names)
df["diagnosis"] = pd.Categorical(data.target_names[data.target])

km = KMeans(n_clusters=2, n_init=10, random_state=42)
labels = km.fit_predict(X)

enrichments = enrichment_analysis(
    df, labels, ["diagnosis", "mean radius", "mean concavity"]
)

print("Significant enrichments (Bonferroni-corrected p < 0.05):")
for e in enrichments:
    print(f"\n  Cluster {e['cluster']}, {e['variable']}:")
    print(f"    p_adjusted = {e['p_adjusted']:.2e}")
    for k, v in e['effect'].items():
        print(f"    {k}: {v}")

Significant enrichments (Bonferroni-corrected p < 0.05):
  Cluster 0, diagnosis:
    p_adjusted = 3.41e-42
    benign: 0.20
    malignant: 2.35

  Cluster 0, mean radius:
    p_adjusted = 1.18e-67
    cluster_mean: 17.12
    other_mean: 12.55
    cohens_d: 1.64

  Cluster 0, mean concavity:
    p_adjusted = 4.92e-51
    cluster_mean: 0.14
    other_mean: 0.05
    cohens_d: 1.47

Notice that the k=2 silhouette score was only 0.292, barely above the "no structure" threshold, yet enrichment analysis revealed a malignant-diagnosis association at p = 3.41 x 10^-42 with Cohen's d of 1.64, where Cohen's d measures the difference between two group means in units of their pooled standard deviation (values above 0.8 are conventionally considered large effects). A metric that measures geometric separation nearly missed a partition with overwhelming biological significance.

4. LLM-Assisted Hypothesis Generation

Enrichment analysis tells you what distinguishes a cluster. The next step, generating a why hypothesis, has traditionally required domain expertise and creative thinking. Large language models offer a new tool for this step: given a structured summary of cluster characteristics, an LLM can propose candidate explanations, suggest follow-up experiments, and connect patterns to existing scientific literature.

This is not automated discovery. The LLM does not know your data; it knows language patterns that correlate with scientific reasoning. Treat the hypotheses it generates as suggestions for human evaluation, not as conclusions. But as hypothesis generators, LLMs are remarkably productive, especially when given well-structured input.

def build_cluster_summary(df, labels, features, metadata_cols):
    """Build a structured text summary of cluster properties
    suitable for LLM hypothesis generation."""
    summaries = []
    for cluster_id in sorted(set(labels)):
        if cluster_id == -1:
            continue
        mask = labels == cluster_id
        n = mask.sum()

        lines = [f"## Cluster {cluster_id} (n={n}, "
                 f"{n/len(labels)*100:.1f}% of data)"]

        # Feature statistics
        lines.append("\n### Key Feature Statistics:")
        for feat in features:
            vals = df.loc[mask, feat]
            all_vals = df[feat]
            z = (vals.mean() - all_vals.mean()) / all_vals.std()
            if abs(z) > 0.5:  # only report distinctive features
                direction = "higher" if z > 0 else "lower"
                lines.append(
                    f"- {feat}: mean={vals.mean():.3f} "
                    f"(dataset mean={all_vals.mean():.3f}, "
                    f"{abs(z):.1f} SD {direction})"
                )

        # Metadata composition
        lines.append("\n### Metadata Composition:")
        for col in metadata_cols:
            if df[col].dtype.name == "category":
                composition = df.loc[mask, col].value_counts(
                    normalize=True)
                lines.append(f"- {col}: " + ", ".join(
                    f"{k}={v:.1%}" for k, v in composition.items()
                ))

        summaries.append("\n".join(lines))

    return "\n\n".join(summaries)

# Build the summary
mean_features = [c for c in df.columns if c.startswith("mean")]
summary = build_cluster_summary(
    df, labels, mean_features, ["diagnosis"]
)
print(summary[:1200])  # preview first portion

def generate_hypotheses_prompt(cluster_summary, domain_context):
    """Create a prompt for LLM hypothesis generation from
    cluster analysis results."""
    prompt = f"""You are a scientific data analyst examining
clustering results from a biomedical dataset. Based on the
cluster characteristics below, generate 3 testable hypotheses
that could explain the observed patterns.

For each hypothesis:
1. State the hypothesis clearly and specifically
2. Explain what evidence in the clusters supports it
3. Propose one experiment or analysis to test it
4. Rate your confidence (low/medium/high) and explain why

Domain context: {domain_context}

Cluster Analysis Results:
{cluster_summary}

Generate your hypotheses:"""
    return prompt

# Example: create the prompt (would be sent to an LLM API)
domain = ("Breast cancer cell nuclei measurements from "
          "digitized fine needle aspirate images. Features "
          "describe nuclear size, shape, and texture.")
prompt = generate_hypotheses_prompt(summary, domain)

# In production, you would call:
# from anthropic import Anthropic
# client = Anthropic()
# response = client.messages.create(
#     model="claude-sonnet-4-20250514",
#     max_tokens=2000,
#     messages=[{"role": "user", "content": prompt}]
# )
# hypotheses = response.content[0].text

print("Prompt length:", len(prompt), "characters")
print("Ready for LLM API call.")

With validated clusters in hand and LLM-generated candidate explanations on the table, the remaining challenge is organizational: how do you move from a promising pattern to a result you can trust enough to publish? Figure 25.5 illustrates the five-stage pipeline that enforces this separation.

5. The Exploration-to-Confirmation Pipeline

The complete pipeline from exploration to hypothesis follows a strict separation between discovery and confirmation. This separation, rooted in the scientific method covered in Chapter 2, prevents the most dangerous mistake in data science: presenting an exploratory finding as a confirmed result. Figure 25.3.1 illustrates this exploration-to-confirmation pipeline.

1. Explore (this chapter): cluster, project, visualize, identify patterns.
2. Validate (this section): internal metrics, stability, enrichment.
3. Hypothesize (LLM-assisted): generate candidate explanations.
4. Design (see Chapter 46): plan a confirmatory experiment on fresh data.
5. Confirm: test the hypothesis with pre-registered analysis (where the hypothesis, statistical tests, and success criteria are publicly declared before seeing the new data) on the new data.

Steps 1 through 3 use the same data. Step 5 must use different data. If you test your hypothesis on the data that generated it, you are very likely to find support for it, because the hypothesis was tailored to that data's idiosyncrasies. This is not a technicality; it is the difference between science and storytelling. In short: A pattern is not a finding until it survives validation, stability, enrichment, and replication on data it has never seen.

6. Pitfalls in Exploratory Discovery

The most common ways exploratory analysis goes wrong fall into five categories.

Pitfall 1: Clustering in the embedding. As noted in Section 25.2, Uniform Manifold Approximation and Projection (UMAP) distorts distances. Clusters found in a UMAP projection may not exist in the original space. Always cluster in the original (or Principal Component Analysis (PCA)-reduced) space and use UMAP only for visualization.

Pitfall 2: Treating cluster labels as ground truth. Cluster labels are hypotheses, not facts. A cluster labeled "cell type A" is a group of observations that the algorithm found similar. Whether they are genuinely the same cell type requires external validation (enrichment analysis, marker gene expression, functional assays).

Pitfall 3: Ignoring the noise label. HDBSCAN assigns ambiguous points to a noise cluster (label -1). Silently dropping these points biases all downstream analysis toward the cleanest, most prototypical examples. Report how many points are noise and analyze their properties.

Interpretation Traps

Pitfall 4: Over-splitting. More clusters always reduce intra-cluster variance (in the limit, \(N\) clusters with one point each have zero variance). Metrics like BIC and DBCV penalize complexity, but no metric can substitute for the question: "Does this additional cluster correspond to a scientifically meaningful distinction?"

Pitfall 5: Feature leakage into cluster interpretation. If you cluster on 30 features and then report that the clusters differ on those same 30 features, you have said nothing. Enrichment analysis must use variables that were not part of the clustering input. Otherwise, the differences are tautological.