Sparse Bayesian group factor model for multiple count tables

Fits the sparse Bayesian group factor model (Sp-BGFM) of Zhang, Shen, Chen and Lee by MCMC. The model takes $M$ count tables measured on the same samples (for example bacterial and viral OTU tables from a multi-domain microbiome study) and estimates the interaction structure among features both within and across the tables.

Usage

spbgfm(
  Y,
  J = NULL,
  subject = NULL,
  X = NULL,
  K = 10,
  niter = 1e+05,
  burnin = NULL,
  thin = 5,
  seed = NULL,
  control = spbgfm_control(),
  verbose = TRUE
)

Arguments

Y: Either a list of $M$ count matrices, each n x J_m with samples in rows and features in columns (the domains), or a single n x J matrix of the tables already column-bound, in which case J must be supplied. Counts must be non-negative integers.
J: Integer vector of length $M$ giving the number of features in each table. Required only when Y is a single matrix; inferred from the list otherwise.
subject: Optional length-n vector (factor, integer or character) giving the subject each sample came from, used to share the baseline abundance $\alpha_{smj}$ across repeated samples. Defaults to one subject per sample.
X: Optional n x P covariate design matrix for the mean regression. For a categorical covariate use one indicator column per level and no intercept; the intercept is absorbed into the size factor and baseline abundance. When NULL (default) no covariate effects are estimated.
K: Number of latent factors. Choose a value large enough to capture the covariance; redundant factors are shrunk out through $\tau_k$. See choose_K.
niter: Total number of MCMC iterations. Default 100000.
burnin: Number of burn-in iterations discarded before saving. Default niter/2.
thin: Thinning interval for saved samples. Default 5.
seed: Optional integer seed for reproducible initialisation.
control: A spbgfm_control list of hyperparameters and tuning constants.
verbose: Logical; print progress. Default TRUE.

Value

An object of class "spbgfm": a list with element samples holding posterior draws of Lambda (J x K x nsamp), sig2 (M x nsamp), tau, phi, eta, ri, alpha, xi_alpha, psi_alpha, w_alpha and, when a covariate is supplied, beta (J x P x nsamp); together with runtime, data (dimensions, domain index, subject index and X), prior (the data-driven prior means $\nu^r$ and $\nu^\alpha$), and control.

Summarise the fit with posterior_cor, posterior_Sigma, posterior_beta and posterior_predict.

Details

Counts are modelled directly, without pre-normalisation, as a Dirichlet process mixture of rounded multivariate log-normal kernels. Writing $y^\star_i$ for the latent continuous vector underlying the counts of sample $i$, $$y^\star_i \sim \mathrm{log-N}_J(r_i + \alpha_{s_i} + X_i\beta, \Sigma), \quad y_{i} = \lfloor y^\star_{i} \rfloor,$$ where $r_{im}$ is a sample size factor for table $m$ (model-based normalisation), $\alpha_{s_i m j}$ is a subject-specific baseline abundance, and $\Sigma = \Lambda\Lambda' + V$ is decomposed through a group factor model. The columns of $\Lambda$ are given a Dirichlet-Horseshoe prior, which induces joint sparsity across all $J$ features of all tables and makes the estimate of $\Sigma$ reliable even when $J \gg n$. Both $r$ and $\alpha$ follow mean-constrained Dirichlet process mixtures, which absorb excess zeros and large inter-sample variability.

References

Zhang, S., Shen, Y., Chen, I. A. and Lee, J. (2025). Sparse Bayesian Group Factor Model for Feature Interactions in Multiple Count Tables Data. Journal of the American Statistical Association, 120(550), 723–736. doi:10.1080/01621459.2025.2449721

Examples

# \donttest{
sim <- simulate_spbgfm(n = 20, J = c(30, 10), seed = 1)
fit <- spbgfm(sim$Y, K = 5, niter = 2000, verbose = FALSE)
rho <- posterior_cor(fit)
dim(rho)
#> [1] 40 40
# }