Bayesian sparse principal coordinates analysis with delta-tolerant linear approximation for microbiome data

1. Problem Statement

Principal Coordinates Analysis (PCoA) is the standard exploratory tool for visualizing beta-diversity in microbiome studies, particularly when using ecologically meaningful dissimilarity measures like Bray–Curtis or Hellinger distances. However, classical PCoA has a critical limitation:

• Lack of Interpretability: The axes in PCoA are derived from the sample-to-sample dissimilarity matrix, not directly from the original taxa (features). Consequently, it is difficult to identify which specific taxa drive the observed ordination patterns.
• The Gap: Researchers need a method that preserves the geometric structure of the distance-based ordination while providing a sparse, interpretable mapping back to the original taxa.
• Challenge: Unlike Principal Component Analysis (PCA), where coordinates are linear combinations of centered variables, PCoA coordinates arise from a double-centered distance matrix. This relationship is not always linear, especially with non-Euclidean distances, making a direct sparse linear surrogate difficult to justify without quantifying the approximation error.

2. Methodology: Bayesian Sparse PCoA (BSPCoA)

The authors propose BSPCoA, a post-hoc framework that approximates the leading PCoA coordinates using a sparse linear surrogate in the original taxa space.

A. The Linear Surrogate Framework

Let $X \in \mathbb{R}^{n \times p}$ be the preprocessed feature matrix (taxa) and $Z_k \in \mathbb{R}^{n \times k}$ be the matrix of the first $k$ leading PCoA coordinates. The goal is to find a coefficient matrix $B \in \mathbb{R}^{p \times k}$ such that:
$$Z_k \approx XB$$
If $Z_k$ can be well-approximated by $XB$, the columns of $B$ represent the "loadings" of the taxa on the ordination axes.

B. Delta-Tolerance Diagnostic

To address the fact that non-Euclidean distances may not admit an exact linear representation, the authors introduce a $\delta$-tolerance diagnostic:

• Relative Reconstruction Error ($\delta$): $\delta(B) = \frac{\|Z_k - XB\|_F}{\|Z_k\|_F}$.
• Explanation Index (ExI): The cosine of the angle between the surrogate and the true coordinates.
• Significance: $\delta^*$ (the minimum possible error) quantifies the irreducible loss of forcing a linear model on the ordination. If $\delta^*$ is small, a sparse linear interpretation is well-supported.

C. Bayesian Formulation

The method treats $Z_k$ as a "pseudo-response" in a multivariate regression model $Y = XB + E$. To induce row sparsity (selecting entire taxa across all axes), they employ Three-Parameter Beta Normal (TPBN) global-local shrinkage priors on the rows of $B$.

• Prior Structure: A hierarchical model where each row $b_j$ of $B$ has a local scale parameter $\psi_j$, which is shrunk by a global parameter $\tau$.
• Specifics: The authors utilize the Horseshoe prior (a special case of TPBN) to aggressively shrink irrelevant taxa coefficients to zero while preserving large signals.
• Algorithm: An iterative Gibbs sampling procedure (Algorithm 1) is used to estimate $B$ and update the alignment matrix $A$ (via Procrustes analysis) to handle the rotational ambiguity of PCoA.

D. Theoretical Guarantees

The paper establishes posterior concentration results for $\delta$-tolerant surrogates.

• Theorem 2: The posterior distribution concentrates around the set of linear surrogates that approximate the PCoA geometry within the tolerance level $\delta$.
• Corollary 1: In the special case of Euclidean distance (where $Z_k$ is exactly linear in $X$), the method reduces to standard high-dimensional Sparse PCA (SPCA) with established convergence rates.

3. Key Contributions

1. BSPCoA Framework: A novel post-hoc method that bridges the gap between distance-based ordination (PCoA) and feature-level interpretation (sparse loadings).
2. $\delta$-Tolerance Diagnostic: A quantitative metric to assess when a linear taxon-level explanation is valid for a given distance metric, preventing over-interpretation of non-linear structures.
3. Theoretical Extension: Extends Bayesian posterior concentration theory to "approximate" surrogates, acknowledging that non-Euclidean ordinations may not have exact linear representations.
4. Scalability: Introduces a subsample projection strategy for large datasets. Instead of decomposing the full $n \times n$ distance matrix, the model is fitted on a small subsample ($m \ll n$) and projected onto the full dataset, drastically reducing computational cost.

4. Results

Simulation Studies

• Euclidean Benchmark: BSPCoA successfully recovers the same sparse support and explained variance as frequentist SPCA when the distance is Euclidean, validating the Bayesian implementation.
• Microbiome Simulations (Dirichlet-Multinomial):
  • Data was generated with varying sparsity levels ($\alpha = 0.5$ and $\alpha = 0.1$) and Bray-Curtis distances.
  • Performance: BSPCoA achieved clustering accuracy (BM-ACC) comparable to classical PCoA (often 100% in sparse settings) while providing a sparse set of explanatory taxa.
  • Robustness: Even when the surrogate explained a smaller proportion of the total eigenvalue mass, the clustering performance remained high, indicating that the selected taxa captured the dominant group structure effectively.

Real Data Application: Hadza Gut Microbiome

• Dataset: Seasonal variation in the Hadza hunter-gatherer gut microbiome (Smits et al., 2017).
• Findings:
  • BSPCoA produced an ordination highly aligned with classical PCoA (using Hellinger distance).
  • Clustering: The clustering accuracy for wet vs. dry seasons improved from 0.559 (PCoA) to 0.609 (BSPCoA), suggesting the shrinkage removed noise.
  • Interpretability: The method identified a parsimonious set of taxa driving the seasonal variation, whereas classical PCoA offered no direct taxon-level coefficients.

5. Significance

• Interpretability in Ecology: Microbiome studies often rely on complex distance metrics (Bray-Curtis, UniFrac) that obscure the underlying biological drivers. BSPCoA provides a rigorous way to link these complex distances back to specific taxa.
• Handling Non-Linearity: By introducing the $\delta$-tolerance, the method avoids the pitfall of forcing a linear interpretation where none exists, offering a diagnostic tool for researchers to gauge the validity of their feature selection.
• Computational Efficiency: The subsample projection strategy makes Bayesian ordination feasible for large-scale microbiome datasets (thousands of samples), overcoming the $O(n^3)$ bottleneck of classical PCoA eigendecomposition.
• Statistical Rigor: The use of Bayesian shrinkage provides natural uncertainty quantification (posterior intervals) for the selected taxa, a feature often missing in frequentist sparse PCA implementations.

In summary, BSPCoA offers a statistically sound, computationally efficient, and interpretable solution for modern microbiome beta-diversity analysis, allowing researchers to visualize complex community structures while pinpointing the specific taxa responsible.