Simulation 2: an arbitrary covariance and a categorical covariate • SpBGFM

Simulation 2 keeps $J_1 = 150$ and $J_2 = 50$ but makes the problem harder in two ways. The true covariance is no longer generated from a factor model: it is an arbitrary sparse correlation matrix produced by the vine method of Lewandowski, Kurowicka and Joe (2009), rescaled by feature-specific variances drawn from $\mathrm{Unif}(1, 1.5)$ . And a binary categorical covariate is added, with two samples per subject (one under each level), giving $N = 40$ samples from $S = 20$ subjects. About 45% of the counts are zero.

Simulating the data

The baseline abundance mixture weights are taken from the empirical zero rates of the skin microbiome data, so that the simulated tables resemble the real ones.

library(SpBGFM)
data(skin)

zero_rate <- c(colMeans(skin$Y$bacteria == 0), colMeans(skin$Y$virus == 0))
J <- c(150, 50)
set.seed(2)
zp <- c(sample(zero_rate[seq_len(skin$J[1])], J[1], replace = TRUE),
        sample(zero_rate[skin$J[1] + seq_len(skin$J[2])], J[2], replace = TRUE))

sim <- simulate_spbgfm(
  n         = 40,
  J         = J,
  type      = "vine",
  covariate = TRUE,
  beta_zero = 0.8,
  zero_prop = zp,
  seed      = 2
)

mean(do.call(cbind, sim$Y) == 0)
table(sim$subject)          # two samples per subject
head(sim$X)                 # (1, 0) and (0, 1) alternating

Fitting with a covariate

Pass the subject index so that the baseline abundance $\alpha_{smj}$ is shared across the two samples of a subject, and the design matrix X. Use one indicator column per level and no intercept: the intercept is absorbed into the size factor and baseline abundance.

fit <- spbgfm(sim$Y, subject = sim$subject, X = sim$X,
              K = 15, niter = 100000, thin = 5, seed = 2)
fit

Interaction structure

rho_hat <- posterior_cor(fit)
ut <- upper.tri(sim$cor)
sqrt(mean((rho_hat[ut] - sim$cor[ut])^2))

Covariate effects

posterior_beta() with a contrast gives the difference between two levels, which is what Fig 7 plots. The interval estimates that exclude zero are flagged in the significant column.

b <- posterior_beta(fit, contrast = c(1, 2))
head(b)
sum(b$significant)

truth <- sim$beta[, 1] - sim$beta[, 2]
plot(truth, b$median, xlab = "true beta_1 - beta_2",
     ylab = "posterior median", pch = 19, cex = 0.6)
segments(truth, b$lower, truth, b$upper, col = "grey")
abline(0, 1, lty = 2)

Predictive checking under each condition

feat <- c(12, 32, 150 + 11)
pred_a <- posterior_predict(fit, features = feat, r_pred = 0, x = c(1, 0))
pred_b <- posterior_predict(fit, features = feat, r_pred = 0, x = c(0, 1))

for (q in seq_along(feat)) {
  plot(density(log(pred_a[, q] + 1)), main = paste("feature", feat[q]),
       xlab = "log(count + 1)")
  lines(density(log(pred_b[, q] + 1)), lty = 2)
}

The shift between the solid and dashed densities is the covariate effect on the predicted counts for that feature.

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