Expand description
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, so scores are monotonically non-increasing along the merge sequence.
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. Picking the best pair still scans the cached scores, K² lookups per merge (K³ overall, no lgamma), and 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.