📊 statistics

Spatially orthogonal factor models for spatial transcriptomics and remote sensing data

This paper introduces a spatially orthogonal factor model that resolves key limitations of existing methods by enforcing orthogonal loadings, accommodating non-stationary spatial priors, and achieving linear computational complexity, thereby enabling efficient inference of spatial variation in large-scale transcriptomics and remote sensing datasets.

Original authors: Dan Cunha, Lukas M. Weber, Mark A. Friedl, Luis Carvalho

Published 2026-08-27

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Dan Cunha, Lukas M. Weber, Mark A. Friedl, Luis Carvalho

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine trying to understand a complex landscape by looking at a map where every single point holds a secret. In fields ranging from biology to climate science, researchers often face this challenge. They have data collected from thousands of specific locations, and at each spot, there is a rich set of measurements. In biology, these measurements might be the activity levels of thousands of genes inside a slice of brain tissue. In climate science, they could be the greenness of vegetation measured by satellites across an entire continent. The goal is to find the hidden patterns that explain why things look the way they do across space. Scientists have long used a tool called principal component analysis to simplify this chaos. Think of it as a way to find the most important "directions" of change in the data, reducing thousands of measurements into a few key themes. However, standard methods often treat each location as if it were isolated, ignoring the fact that nearby places usually influence one another. When researchers try to force these standard tools to account for geography, they often end up with results that are mathematically messy or computationally impossible to calculate for large datasets.

A team of researchers has developed a new approach to solve these problems, creating a method that respects the geography of the data while keeping the math clean and fast. Their work, which they call a spatially orthogonal factor model, is designed to handle massive amounts of spatial data without losing the crucial information about how nearby points relate to each other. The core idea is to find the main patterns of variation in the data, but to do so in a way that ensures these patterns are distinct from one another and that their shapes change smoothly across the landscape. Unlike older methods that assume the rules of geography are the same everywhere, this new model allows the strength of these connections to vary. In some areas, a pattern might stretch far across the land; in others, it might be very local. The researchers proved that their method can find the best possible answer for these patterns and, crucially, they built a system that can do this calculation quickly, even when dealing with hundreds of thousands of locations.

To test their idea, the team applied it to two very different worlds. First, they looked at a detailed map of a human brain, specifically a region called the dorsolateral prefrontal cortex. This area is known for its layered structure, much like the pages of a book, and scientists have long tried to understand how gene activity changes from one layer to the next. Using data from over 3,500 tiny spots on the brain tissue, the new model successfully identified patterns that matched these known layers. But it also found something extra. It uncovered gene activity patterns that did not follow the neat layers, revealing signals associated with blood and immune cells that were mixed in with the brain tissue. These findings aligned with previous discoveries but were found more clearly and efficiently by the new method. The model also showed that the "reach" of these patterns, or how far a single influence extends, was not the same everywhere; it was strongest in the outermost layer of the brain and varied in other regions.

The second test took the model to a much larger scale: the entire sub-Saharan African continent. Here, the researchers used satellite data measuring the greenness of plants over a single year. This is a difficult task because a single year of data is often too noisy to reveal the true, long-term climate patterns that scientists rely on. Standard methods struggled to find clear patterns in this short, noisy snapshot. However, the new model, by using the spatial relationships between nearby pixels, was able to pull out meaningful seasonal trends. It identified distinct regions where plant life behaves in specific ways, such as the Sahel and the humid tropics. Remarkably, the patterns it found in just one year of data closely matched the major climate zones identified by looking at forty years of historical data. The model even spotted subtle differences in the humid tropics that the standard method missed, showing that it could see the underlying structure of the landscape even when the data was limited.

The success of this work lies in how it handles the mathematics of space. The researchers showed that by treating the spatial connections as a flexible, changing map rather than a fixed rule, they could get better answers. They also developed a way to estimate these changing connections automatically, using a strategy of testing the model on parts of the data it hasn't seen yet to ensure it is learning the right things. This approach allows the model to run on standard computers in a matter of minutes, even for datasets that would take days or weeks to process with older techniques. The result is a tool that helps scientists see the invisible threads connecting places, whether those places are tiny spots on a brain slice or vast stretches of the African savanna, revealing the hidden order in the complexity of our world.

Technical Summary: Spatially Orthogonal Factor Models for Spatial Transcriptomics and Remote Sensing Data

Problem Statement

Spatially indexed data, prevalent in fields such as spatial transcriptomics and remote sensing, often require dimensionality reduction to infer latent modes of spatial variation. While Probabilistic Principal Component Analysis (PCA) is a standard tool, existing methods for incorporating spatial information into this framework face three primary limitations:

Non-Orthogonal Loadings: Standard probabilistic PCA loadings are not orthogonal. Subsequent orthogonalization of these loadings corrupts the original prior spatial information.
Stationarity Assumptions: Current methods typically assume stationarity in the spatial prior, often relying on fixed priors or rule-of-thumb estimators rather than learning the covariance structure from data.
Computational Complexity: Methods incorporating spatial information generally fail to achieve linear-time computational complexity with respect to the number of spatial locations ( $m$ ), often requiring cubic complexity due to eigendecomposition of large covariance matrices.

Furthermore, it remains unclear how the incorporation of spatial information alters maximum likelihood or maximum a posteriori (MAP) inferences in these settings.

Methodology

Model Formulation

The authors propose Spatially Orthogonal Factor Models (SOFM), which parameterize the model directly with orthogonal loadings ( $U$ ). Unlike traditional approaches that place spatial dependence on latent factors, SOFM places spatial dependence on the loadings matrix $U$ . The model is specified as:
$y_i | U, L, \sigma^2 \sim N(0, UL^2U' + \sigma^2 I_m)$
where $y_i$ is the observation vector at location $i$ , $L$ contains singular values, and $U$ is an orthogonal matrix representing spatial patterns.

Prior Derivation

To address the lack of orthogonality in standard priors, the authors derive the sampling distribution for the Singular Value Decomposition (SVD) of a matrix where the first $k$ singular values are unique and the remaining $m-k$ are repeated. This distribution, derived from James (1964) and adapted for PCA constraints, serves as the prior for $U, L, \sigma^2$ .

Proposition 1: Establishes that the conditional prior of $U$ given $L, \sigma^2, \Sigma$ is a matrix-Fisher-Bingham distribution, which is conjugate to the likelihood.
Theorem 1: Demonstrates that the unconditional MAP estimator for the orthogonal loadings $U$ corresponds to the eigenvectors of $(S + \frac{1}{n}\Sigma)$ , where $S$ is the empirical covariance matrix and $\Sigma$ is the prior spatial covariance matrix.

Estimation Algorithm (MM-EM)

Direct eigendecomposition of $(S + \frac{1}{n}\Sigma)$ is computationally infeasible for large $m$ (cubic complexity). To resolve this, the authors develop a Minorization-Maximization within Expectation-Maximization (MM-EM) algorithm:

E-step: Computes posterior expectations of latent states.
M-step: Uses a minorization-maximization approach to update $U$ . By exploiting the convexity of the objective function and the Procrustes solution (Schönemann, 1966), the algorithm maximizes a surrogate function of the form $\text{tr}(M^T U)$ .
Complexity: This approach achieves linear-time computational complexity with respect to the number of spatial locations.

Non-Stationary Covariance Estimation

To address the stationarity limitation, the authors propose a validation framework to estimate a non-stationary spatial covariance function (based on Paciorek and Schervish, 2006).

Held-Out Locations: A set of spatial locations is held out during training. The model treats outcomes at these locations as latent.
Optimization: The algorithm estimates local length scales and sill parameters by maximizing the held-out predictive likelihood.
Smoothing: Local estimates are generalized across the domain using two-dimensional P-splines to ensure spatial smoothness of the covariance parameters.

Key Contributions

Theoretical Derivation: Derivation of the sampling distribution for SVD under repeated eigenvalue constraints, enabling a conjugate prior for orthogonal loadings.
MAP Characterization: Proof that the MAP estimate for orthogonal loadings is the eigendecomposition of the sum of empirical and prior spatial covariances.
Scalable Algorithm: Development of an MM-EM algorithm that reduces computational complexity from cubic to linear with respect to spatial locations.
Non-Stationary Estimation: A novel validation strategy for estimating non-stationary spatial covariance parameters using held-out data and P-spline smoothing.
Applications: Demonstration of the method on two distinct datasets: human brain spatial transcriptomics and continental-scale remote sensing.

Results

Case Study 1: Spatial Transcriptomics (DLPFC)

Applied to postmortem human brain tissue from the dorsolateral prefrontal cortex (DLPFC):

Laminar and Non-Laminar Patterns: The model successfully identified spatial factors corresponding to cortical layers (laminar structure) as well as non-laminar patterns.
Biological Relevance: Factors captured genes associated with blood (e.g., HBB) and immune cell populations (e.g., IGKC), consistent with previous findings by Maynard et al. (2021).
Performance: The SOFM model outperformed standard PCA and a supervised "Laminae-PCA" in capturing non-laminar spatial variability and explained a higher proportion of variance in the empirical covariance for the first few factors.
Efficiency: The algorithm ran in under 90 seconds for the full gene set (33,538 genes) on 32 cores.

Case Study 2: Remote Sensing (Sub-Saharan Africa)

Applied to NDVI data from the GIMMS dataset over sub-Saharan Africa:

Non-Stationarity: Estimated length scales varied spatially, showing stronger dependence in the Sahel and sub-humid tropics compared to transition zones.
Short-Term vs. Long-Term: When applied to a single year of data (2022), SOFM identified spatial patterns (specifically in the eastern humid tropics and Mozambique) that were present in a 40-year PCA model but missed by a single-year PCA model.
Significance: This demonstrates that spatial regularization allows the model to recover underlying spatial factors from limited temporal data that standard PCA cannot detect.

Significance and Claims

The paper claims that SOFM resolves the trade-off between incorporating spatial information and maintaining orthogonal, interpretable factors. By deriving a specific prior for orthogonal loadings, the method avoids the corruption of spatial information seen in post-hoc orthogonalization. The linear-time complexity makes the method scalable to large spatial datasets (e.g., continental scales or high-resolution tissue sections).

The authors emphasize that their approach to estimating non-stationary covariance via held-out likelihood maximization is a novel contribution to the application areas of spatial transcriptomics and environmental science. The results suggest that spatially orthogonal factors can reveal biologically and environmentally relevant modes of variation that are obscured by standard dimensionality reduction techniques, particularly when data is sparse in the temporal or observational domain.

The paper concludes modestly, noting that while the current framework focuses on spatial dependence in loadings, future work could explore temporal dependence in latent variables or evolving loadings, as well as formal inferential procedures for gene association in transcriptomics.