A reduced rank model for spatial categorical data with many classes

1. Problem Statement

Spatial categorical data (e.g., land cover, soil taxonomy, species distribution) are common in environmental sciences. Modeling these data presents three primary challenges:

• Spatial Dependencies: Locations close together tend to have similar class probabilities.
• Compositional Constraints: Class probabilities must sum to unity (simplex constraint).
• High Dimensionality: As the number of classes ($J$) increases, standard Latent Gaussian Models (LGMs) require inferring a unique spatial effect for every class and their cross-covariances. This leads to a rapid explosion in parameter dimension, causing computational intractability and statistical inefficiency.

Existing Bayesian inference methods for multinomial LGMs, such as the Pólya-Gamma data augmentation (Polson et al., 2013) or logit-beta conjugate priors (Bradley et al., 2019), rely on the separability of the linear predictor. This means the log-odds of each class must depend on a unique, distinct latent effect. This constraint precludes the use of factor models, where multiple classes share a smaller set of latent factors, which is essential for reducing dimensionality in high-class scenarios.

2. Methodology

The authors propose an identifiable reduced-rank spatial multinomial model combined with a novel Gibbs sampler utilizing Laplace approximations.

A. The Reduced-Rank Model

Instead of modeling $J-1$ independent spatial processes, the model assumes the spatial variation is driven by a low-dimensional set of shared latent factors ($u \ll J-1$).

• Latent Structure: The logit vector $\psi(s)$ is modeled as:
  $$\psi(s) = \mu + \Gamma \tilde{\eta}(s)$$
  Where:
  • $\mu$ is a mean vector.
  • $\Gamma$ is a $(J-1) \times u$ factor loading matrix.
  • $\tilde{\eta}(s)$ is a $u$-dimensional multivariate Gaussian process (the shared latent factors).
• Fixed-Rank Process: The latent factors $\tilde{\eta}(s)$ are modeled using a predictive process approach with spatial basis functions and knots to ensure computational feasibility.
• Identifiability: To ensure the factorization is unique, the authors impose constraints: $\Gamma$ is unit-lower triangular, and the precision matrix of the factors is diagonal.

B. Posterior Inference (Gibbs Sampler with Laplace Proposals)

Standard conjugate sampling methods (Pólya-Gamma) fail because the shared factors create dependencies across all class logits, breaking the conditional independence required for those algorithms.

• Approach: The authors employ a Gibbs sampler where the conditional posteriors of the latent spatial effects ($W$) and factor loadings ($\Gamma$) are sampled using Metropolis-Hastings (MH) steps.
• Proposal Distribution: The proposals for the MH steps are generated using Laplace approximations.
  1. Find the mode ($\hat{w}$) of the log-conditional posterior using Newton-Raphson.
  2. Approximate the posterior as a Multivariate Normal distribution centered at the mode with covariance equal to the inverse Hessian ($-H(\hat{w})^{-1}$).
  3. Generate a candidate sample and accept/reject it via the MH ratio.
• Other Parameters:
  • Mean vector $\mu$ and factor matrix $\Gamma$ are sampled similarly using Laplace-MH.
  • Precision parameters ($\omega$) are sampled exactly (Gamma priors lead to Gamma posteriors).
  • Range parameter ($\phi$) is sampled via MH.

C. Prediction

The model supports flexible joint inference:

• Single Class: Probability of a specific class.
• Class Unions: Probability of a subset of classes (e.g., "any Oak species").
• Area Summaries: Probability of occurrence within a specific geographic area (aggregated from grid-level predictions).

3. Key Contributions

1. Reduced-Rank Multinomial Framework: Extends the popular factor model concept to spatial multinomial data, allowing for scalable inference with many classes by reducing the dimension of spatial dependencies.
2. Novel Inference Scheme: Develops a Gibbs sampler using Laplace-approximation proposals to bypass the limitations of Pólya-Gamma augmentation in non-separable factor models. This avoids deterministic quadrature (like INLA) which struggles with high-dimensional factor weights.
3. Identifiability Constraints: Provides a specific, comprehensive parameterization (unit-lower triangular $\Gamma$) to ensure the reduced-rank decomposition is unique.
4. Practical Implementation: Demonstrates the method's scalability on real-world data with 24 classes and 2,000+ observations.

4. Results

Simulation Studies

• Dimension Selection: The study tested model selection using WAIC and PSIS-LOO. While these criteria did not always select the true latent dimension ($u=2$) in every simulation, they consistently selected models with significantly better out-of-sample predictive performance than the full-rank model ($u=4$).
• Laplace Approximation Accuracy: The accuracy of the Laplace proposal (and thus the MH acceptance rate) is heavily dependent on the marginal precision ($\omega$) of the latent factors.
  • High Precision: High acceptance rates (~83%) and accurate posterior approximations.
  • Low Precision: Low acceptance rates (~26%) and severe linear bias in the posterior mean if the MH step is skipped (nested Laplace approximation). The bias arises because low precision leads to large variance in logits, pushing probabilities toward the asymptotes of the softmax function, creating skewed posteriors that a symmetric Gaussian approximation cannot capture well.

Real Data Application (Blue Ridge Mountains)

• Dataset: 2,063 forest inventory plots across North Carolina, classified into 24 dominant tree species (plus "no forest").
• Model Selection: A ternary search was used to find the optimal latent dimension $u$ to minimize WAIC. The optimal dimension was found to be $u=7$.
• Performance:
  • Computation time scaled linearly with $u$ (approx. 0.12s to 1.48s per Gibbs cycle for $u=1$ to $u=15$).
  • The model successfully generated spatially complete predictions for individual species, unions of species (e.g., all Oak classes), and area-level summaries (e.g., probability of Basswood dominance in 10x10km grids).

5. Significance and Discussion

• Scalability: The proposed method makes Bayesian inference feasible for spatial categorical data with many classes, a scenario where standard full-rank LGMs are computationally prohibitive.
• Flexibility: Unlike methods requiring class ordering (e.g., stick-breaking links), this model treats classes symmetrically and allows for joint predictions of arbitrary class subsets.
• Limitations & Future Work:
  • Precision Sensitivity: If the data implies very low marginal precision (high uncertainty) in the latent factors, the Laplace approximation may become inefficient due to skewed posteriors. The authors suggest future work using skew-normal approximations (third-order Taylor expansion) to correct for this bias.
  • Dimension Selection: While effective, selecting the latent dimension $u$ currently requires fitting multiple models. Developing heuristics to avoid repeated fitting for large $J$ is a suggested future direction.
  • Generalizability: The framework is applicable to any multinomial model with known trials, and the authors note it could be adapted for count data (Poisson) with unknown trials.

The paper provides a robust, scalable solution for mapping complex categorical spatial phenomena, bridging the gap between statistical rigor and computational feasibility in high-dimensional settings.