Fit the Zero-Inflated Multivariate Rounded Log-Normal (ZI-MLN) model

Runs the Markov chain Monte Carlo (MCMC) sampler for the ZI-MLN model of Zhang et al. (2023). The model infers the interaction (covariance) structure between microbial features (OTUs) directly from an observed count table, accommodating overdispersion, compositionality, subject-level random effects, covariate effects on abundance, and excess (structural) zeros.

Usage

ZI_MLN(
  Y,
  X = NULL,
  m = NULL,
  M = NULL,
  K = 10,
  niter = 20000,
  burnin = floor(niter/2),
  seed = 3,
  a.phi = 1/2,
  a.tau = 1,
  b.tau = 1/50,
  a.sig = 3,
  b.sig = 3,
  a.vs2 = 3,
  b.vs2 = 3,
  acc.target = 0.44,
  Lr = 8,
  a.phi.r = 3,
  a.w = 5,
  b.w = 5,
  ur2 = 1,
  L.theta = 15,
  a.phi.theta = 1,
  a.w.theta = 5,
  b.w.theta = 5,
  u2.theta = 1,
  m_kappa = 0,
  sig_kappa = 3,
  m_beta = 0,
  sig_beta = 25
)

Arguments

Y: Integer count matrix, n samples (rows) by J OTUs (columns). Raw counts; no normalization is required. Zeros are allowed and modelled.
X: Optional numeric covariate matrix, n by P, one row per sample. Used both in the mean model for abundance and (with an added intercept) in the probit model for the zero-inflation probability. NULL (default) fits the no-covariate model.
m: Integer vector of length n giving the subject/group index of each sample. For example c(1, 1, 2, 3) means samples 1-2 come from subject 1, sample 3 from subject 2 and sample 4 from subject 3. Defaults to seq_len(nrow(Y)) (each sample is its own subject).
M: Number of distinct subjects/groups. Defaults to length(unique(m)).
K: Factor dimension of the low-rank covariance decomposition $\Sigma = \Lambda\Lambda' + \sigma^2 I$. Fixed at a moderately large value (default 10).
niter: Total number of MCMC iterations (default 20000).
burnin: Number of initial iterations discarded as burn-in. Draws are stored for iterations > burnin. Defaults to floor(niter / 2).
seed: Random seed for reproducibility (default 3).
a.phi: Dirichlet-Laplace concentration $a_\phi$ controlling sparsity of the covariance (default 1/2). Smaller values shrink more.
a.tau, b.tau: Gamma shape/rate hyperparameters for $\tau_k$ (defaults 1 and 1/50).
a.sig, b.sig: Inverse-gamma hyperparameters for the idiosyncratic variance $\sigma^2$ (defaults 3, 3).
a.vs2, b.vs2: Inverse-gamma hyperparameters for the subject random-effect variance $u_s^2$ (defaults 3, 3).
acc.target: Target acceptance rate for the adaptive Metropolis update of $\phi_j$ (default 0.44).
Lr, a.phi.r, a.w, b.w, ur2: Hyperparameters of the mean-constrained mixture-of-mixtures prior on the sample-size factors $r_i$: number of mixture components, Dirichlet concentration, two Beta parameters and the component variance.
L.theta, a.phi.theta, a.w.theta, b.w.theta, u2.theta: Corresponding hyperparameters of the mixture-of-mixtures prior on the OTU-size factors $\alpha_j$.
m_kappa, sig_kappa: Prior mean and variance for the probit zero-inflation coefficients $\kappa$ (defaults 0, 3).
m_beta, sig_beta: Prior mean and variance for the abundance regression coefficients $\beta$ (defaults 0, 25). Ignored when X = NULL.

Value

A list of length niter - burnin. Each element is one posterior draw, itself a list with components:

Lambda: J by K factor loading matrix.
sig2: idiosyncratic variance $\sigma^2$.
vs2: subject random-effect variance $u_s^2$.
phi,tau.k: Dirichlet-Laplace parameters.
eta: n by K latent factor scores.
ri,thetaj: sample-size factors $r_i$ and OTU-size factors $\alpha_j$ (only $r_i + \alpha_j$ is identifiable).
sij,smj: sample- and subject-level random effects.
beta: J by P abundance regression coefficients (present only when X is supplied).
delta: n by J presence indicators; 1 = OTU present (drawn from the count distribution), 0 = structural/absent zero. (Note this is the complement of the $\delta_{ij}$ in the paper.)
eps.ij: n by J matrix of estimated absence probabilities $\epsilon_{ij}$ (the paper's zero-inflation probability).
kappa: J by P + 1 probit zero-inflation coefficients (intercept in the first column).

Common posterior summaries: the marginal OTU correlation matrix is cov2cor(tcrossprod(draw$Lambda) + diag(draw$vs2 + draw$sig2, J)), averaged over draws.

Details

A single, unified interface handles both the no-covariate and the with-covariate cases: when X = NULL the covariate regression terms are dropped and the zero-inflation component reduces to an intercept-only probit model, reproducing the "without covariate" sampler.

References

Zhang, S., Shen, Y., Chen, I. A. and Lee, J. (2023). Bayesian Modeling of Interaction between Features in Sparse Multivariate Count Data with Application to Microbiome Study. The Annals of Applied Statistics, 17(3). doi:10.1214/22-AOAS1690

Examples

# \donttest{
sim <- simulate_zimln(n = 30, J = 30, K = 3, seed = 1)
fit <- ZI_MLN(sim$Y, m = sim$m, M = sim$M, niter = 400, burnin = 200)
# posterior mean marginal correlation
J <- ncol(sim$Y)
rho <- Reduce(`+`, lapply(fit, function(d)
  cov2cor(tcrossprod(d$Lambda) + diag(d$vs2 + d$sig2, J)))) / length(fit)
# }