scmix: Bayesian sparse conditional mixture clustering
scmix: Bayesian sparse conditional mixture
clustering
Aqi Dong, Yang-Li Liao, and Danhyang Lee
scmix implements the Bayesian sparse conditional mixture
model of Dong, Liao, and Lee: each mixture component factorizes into a
chain of univariate polynomial regressions, per-component/per-equation
variable selection is sampled (spike-and-slab under a centered
Zellner \(g\)-prior), and an overfitted
sparse mixture selects the number of clusters – all in one run of a
blocked Gibbs sampler. This vignette reproduces, at reduced scale, the
two headline analyses of the accompanying paper.
1. Curved clusters that defeat Gaussian mixtures
Two opposite parabolas: a Gaussian mixture needs many elliptical
components to trace curvature, so BIC over-selects \(K\); the conditional mixture models the
curvature directly.
library(scmix)
n1 <- 400; n2 <- 300
x1 <- rnorm(n1 + n2)
x2 <- c( 0.5 + 1.0 * x1[1:n1]^2,
-2.0 - 0.8 * x1[(n1 + 1):(n1 + n2)]^2) + rnorm(n1 + n2, 0, 0.6)
Y <- cbind(x1, x2)
truth <- rep(1:2, c(n1, n2))
fit <- scmix(Y, K = 7, m = 2, n_iter = 800, burn = 300, seed = 42)
fit
#> Bayesian sparse conditional mixture clustering (scmix)
#> n = 700, p = 2, degree m = 2, order: 1 2
#> clusters (consensus): K = 2 [posterior mode K* = 2]
#> cluster sizes: 398 302
#> runtime: 3.9s
table(consensus = fit$cluster, truth)
#> truth
#> consensus 1 2
#> 1 396 2
#> 2 4 298
A single run with an overfitted cap of \(K
= 7\) recovers the two curved clusters (the chain lengths here
are shortened for vignette build time; the paper uses 1500 sweeps).
What did it learn? Per-cluster regression structure with uncertainty
– the interpretable output that distinguishes sampled selection from
shrinkage:
summary(fit)
#> Bayesian sparse conditional mixture clustering (scmix)
#> n = 700, p = 2, degree m = 2, order: 1 2
#> clusters (consensus): K = 2 [posterior mode K* = 2]
#> cluster sizes: 398 302
#> runtime: 3.9s
#>
#> Per-cluster conditional regressions (terms with inclusion > 0.5 ):
#> Cluster 1 (n = 398):
#> y1 ~ 0.04 (intercept only) [sd 0.97]
#> y2 ~ 0.47 +1.03 y1^2 (P=1.00) [sd 0.64]
#> Cluster 2 (n = 302):
#> y1 ~ -0.09 (intercept only) [sd 1.08]
#> y2 ~ -2.00 -0.81 y1^2 (P=1.00) [sd 0.61]
2. A replication-setting analysis
The paper’s simulation study replicates Melnykov and Wang (2023)
exactly. Here is one draw from their \(p =
3\) supplementary Table S-1 setting (three components, quadratic
conditionals), analyzed the same way:
sim <- scmix_sim_mw3(n = 600, seed = 7) # bundled S-1 generator
fit2 <- scmix(sim$Y, K = 7, m = 2, n_iter = 800, burn = 300, seed = 11)
c(K_consensus = fit2$K, ARI_vs_truth = round(scmix_ari(fit2$cluster, sim$z), 3))
#> K_consensus ARI_vs_truth
#> 3.000 0.969
The consensus partition and \(K^*\)
come from the same single run; no sweep over \(K\), no search over conditioning orders.
The paper’s Section 4 reports the full 100-replicate version of this
analysis (consensus ARI 0.967, \(K^* =
3\) in 100% of replicates), its order-sensitivity study, and the
comparisons with cmbClust, mclust, and the
Bayesian cluster-weighted models of Papastamoulis and Perrakis
(2026).
Options worth knowing
g = "hyper": hyper-\(g\) prior (Liang et al., 2008); removes the
\(\sigma^2\) inflation of the default
unit-information prior at high signal-to-noise, with clustering
conclusions unchanged.gprior_intercept = FALSE (default in
scmix()): centered convention; location-invariant selection
in the linear case.alive_frac in the fit object: per-label occupancy
diagnostics for component-specific summaries. ```