Module bhc 

Bayesian Hierarchical Clustering over K clusters with gene-mass sufficient stats.

Post-hoc bottom-up merge of K clusters driven by a pairwise Dirichlet-Multinomial marginal-likelihood Bayes factor under an empirical-Bayes asymmetric Dirichlet prior centered on the pooled empirical marginal bg[g] = (Σ_k T_k,g) / (Σ_k S_k). The prior expects every cluster to look like bg; the BF measures deviation from that baseline. Housekeeping-shaped clusters merge cheaply; clusters with distinct structure stay separate.

Node log marginal under Dir(γ · bg) prior:

f(T, S) = lgamma(γ) − lgamma(γ + S) + Σ_g [lgamma(γ·bg[g] + T_g) − lgamma(γ·bg[g])]

Pairwise log Bayes factor:

log_BF(i, j) = f(T_i + T_j, S_i + S_j) − f(T_i, S_i) − f(T_j, S_j)

Interpreting the scores

log_BF is the log of the ratio of evidences, H₁ (the two clusters came from one underlying distribution) vs H₀ (they came from two distinct ones):

log_BF > 0 — data favor merging (H₁ wins) log_BF < 0 — data favor keeping separate (H₀ wins) log_BF = 0 — model is indifferent (natural cut point)

Magnitude carries BIC-like meaning: at high counts, log_BF ≈ −n · KL(p_i ‖ p_j), where n is the merged sample size and KL is computed on the proportion vectors p_k = T_k / S_k. Rough calibration (symmetric about zero):

|log_BF| < 1 weak — effectively a toss-up 1 ≤ |log_BF| < 3 moderate evidence (≈ e≈2.7× to 20× odds) 3 ≤ |log_BF| < 5 strong (20× to 150× odds) |log_BF| ≥ 5 decisive (>150× odds)

Compare only WITHIN ONE run: absolute magnitudes are not probabilities (the data-dependent normalizer log Γ(y+1) is dropped since it cancels across any partition of the same edges). Cross-run comparisons are not meaningful.

In the bottom-up agglomeration bhc_merge picks the largest log_BF at each step. The scores need not decrease along the merge sequence: a merged node can score higher against a third cluster than the pair that formed it. The crossing point log_BF = 0 is the natural “consensus cut” threshold exposed via bhc_cut(merges, k, 0.0).

Effective sample size is caller-specified. The caller passes effective_size per cluster and the internal (T, S) gets rescaled so that S_eff = effective_size[k], preserving shape T_g / S. Without rescaling, large raw S values (e.g. 10⁴–10⁶ of projection mass on Visium scale, or Σ_n θ_{n,k}·depth_n at 100k-cell scale) make the DM posterior absurdly sharp and blow up log-BF magnitudes. Typical choices:

• link-community: effective_size[k] = edge_count[k] (the natural N).
• topic model: effective_size[k] = (Σ_n θ_{n,k})² / Σ_n θ_{n,k}² (ESS of the soft assignment).

A Heller-Ghahramani DP tree prior can be slotted in later: extend with a BhcConfig { dp_alpha, ... } and replace the pure pairwise log-BF rule with the full posterior r_k recursion.

Runs over K clusters (not N samples). Pairwise log BFs are cached, so only a new node’s K pairs are scored after each merge: K²·M lgamma calls overall, spread over threads. Every node also keeps its best partner among the nodes before it, so picking the best pair is K lookups per merge plus a rescan of the rows whose best partner was just merged away (K² lookups overall when few rows go stale, K³ at worst, no lgamma). The cache holds up to 2K² values: meant for K up to a few thousand. No RNG; the merges match scoring every pair afresh.

Structs

BhcInput

Borrowed view over the per-cluster sufficient statistics needed by bhc_merge.

BhcMerge

One merge event in the dendrogram.

Functions

bhc_cut

Walk the merge tree and collapse children whose log_bf >= cutoff.

bhc_cut_with_sizes

bhc_cut with the effective_size given to bhc_merge: exactly the clusters with effective_size == 0 get label -1, so the only non-empty cluster of a run without merges keeps a label of its own.

bhc_merge

Bottom-up agglomerative merge by pairwise log Bayes factor.