Simulation 1: recovering a sparse interaction structure • SpBGFM

Simulation 1 of the paper studies two count tables with $J_1 = 150$ and $J_2 = 50$ features observed on only $N = 20$ samples, one per subject. The true covariance has a sparse factor structure: features 1–25 and 51–75 of domain 1 and features 1–25 of domain 2 load on a single randomly chosen factor, and the remaining features do not interact with anything. About 30% of the counts are zero.

The chunks below are not evaluated when the vignette is built, because the MCMC run takes roughly an hour. Run them interactively to reproduce the results.

Simulating the data

library(SpBGFM)

sim <- simulate_spbgfm(
  n      = 20,
  J      = c(150, 50),
  K_true = 5,
  type   = "block",
  active = list(c(1:25, 51:75), 1:25),   # the configuration used in the paper
  v2     = 0.25,
  seed   = 1
)

mean(do.call(cbind, sim$Y) == 0)   # proportion of zero counts

Choosing K

The paper fixes $K$ empirically: perform PCA on the sample covariance of the log-transformed normalised counts and keep enough components to explain 95% of the total variance.

choose_K(sim$Y)

Simulation 1 uses K = 10. A value that is comfortably larger than the truth is safe: redundant factors are shrunk out through $\tau_k$ .

Fitting the model

fit <- spbgfm(sim$Y, K = 10, niter = 100000, burnin = 50000, thin = 5, seed = 1)
fit

Posterior correlation estimates

posterior_cor() reconstructs $\rho^{mm'}_{jj'}$ from $\Sigma = \Lambda\Lambda' + V$ for every retained draw and returns the posterior median.

rho_hat <- posterior_cor(fit)
rho_true <- sim$cor

ut <- upper.tri(rho_true)
sqrt(mean((rho_hat[ut] - rho_true[ut])^2))   # RMSE reported in Table 1

The heatmap in Fig 4(a) puts the estimate in the upper triangle and the truth in the lower triangle:

panel <- rho_true
panel[ut] <- rho_hat[ut]
image(seq_len(nrow(panel)), seq_len(ncol(panel)), panel,
      zlim = c(-1, 1), xlab = "", ylab = "",
      col = colorRampPalette(c("green", "white", "red"))(64))
abline(h = sim$J[1] + 0.5, v = sim$J[1] + 0.5)

The cross-domain block on its own:

cross <- posterior_cor(fit, domains = c(1, 2))
dim(cross)

Posterior predictive checking

For a new sample with size factor $r^{\mathrm{pred}}_m = 0$ , compare predictive draws with the observed counts rescaled to the same size factor. This reproduces Fig 5.

feat <- c(30, 133, 150 + 31)   # OTUs 30 and 133 of domain 1, OTU 31 of domain 2

pred <- posterior_predict(fit, features = feat, r_pred = 0)
obs  <- normalize_counts(fit, features = feat, r_pred = 0, log = TRUE)

for (q in seq_along(feat)) {
  plot(density(log(pred[, q] + 1)), main = paste("feature", feat[q]),
       xlab = "log(count + 1)")
  points(obs[, q], rep(0, nrow(obs)), pch = 4)
}

Agreement between the predictive density and the empirical distribution of the rescaled counts, including the spike at zero and the multimodality, indicates a good fit.

Reference

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