The technological development in genomic sequencing has enabled researchers to unveil the wide variability of bacteria presented within different locations of the body. It is necessary to better understand both environmental and host genetic factors impact the composition of the microbiome to improve disease management. Several analytic approaches are introduced to summarize and assess the single or multiple OTUs using different computational algorithms. Motivated by the multivariate nature of microbiome data with hierarchical taxonomic clusters, counts that are often skewed and zero inflated, we propose a Bayesian latent variable methodology to jointly model multiple operational taxonomic units within a single taxonomic cluster. This novel method can incorporate both negative binomial and zero-inflated negative binomial responses, and can account for serial and familial correlations. This method can help discover both genetic and environmental factors that influence the microbiome.