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

plot(fit, Y)

Scatter plot of the two curved clusters colored by the fitted consensus partition

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