Inside-out cross-covariance for spatial multivariate data

As the spatial features of multivariate data are increasingly central in researchers’ applied problems, there is a growing demand for novel spatially-aware methods that are flexible, easily interpretable, and scalable to large data. We develop inside-out cross-covariance (IOX) models for multivariate spatial likelihood-based inference. IOX leads to valid cross-covariance matrix functions which we interpret as inducing spatial dependence on independent replicates of a correlated random vector. The resulting sample cross-covariance matrices are “inside-out” relative to the ubiquitous linear model of coregionalization (LMC). However, unlike LMCs, our methods offer direct marginal inference, easy prior elicitation of covariance parameters, the ability to model outcomes with unequal smoothness, and flexible dimension reduction. As a covariance model for a q-variate Gaussian process, IOX leads to scalable models for noisy vector data as well as flexible latent models. For large n cases, IOX complements Vecchia approximations and related process-based methods based on sparse graphical models. We demonstrate superior performance of IOX on synthetic datasets as well as on colorectal cancer proteomics data. An R package implementing the proposed methods is available at github.com/mkln/spiox.

💡 Research Summary

The paper introduces a novel cross‑covariance construction for multivariate spatial data called Inside‑Out Cross‑Covariance (IOX). Traditional approaches, especially the Linear Model of Coregionalization (LMC), model multivariate spatial dependence by linearly combining a set of univariate correlation functions. While computationally convenient, LMC suffers from several drawbacks: indirect marginal interpretation, difficulty in prior elicitation, inability to handle variables with differing smoothness, challenges incorporating nugget effects, and limited scalability when the number of variables (q) or locations (n) is large.

IOX addresses these issues by defining the cross‑covariance matrix function directly from q univariate correlation functions ρ₁,…,ρ_q and a positive semidefinite matrix Σ that captures inter‑variable dependence. A set of reference locations S (typically the observed sites) is chosen. For each outcome i, the matrix ρ_i(S) is factorized as L_i L_iᵀ (lower‑triangular Cholesky). The mapping h_i(ℓ)=ρ_i(ℓ,S) ρ_i(S)⁻¹ projects any location ℓ onto the space spanned by S. The IOX cross‑covariance is then

C_{ij}(ℓ,ℓ′)=σ_{ij} h_i(ℓ) L_i L_jᵀ h_j(ℓ′)ᵀ + ξ_{ij}(ℓ,ℓ′),

where σ_{ij} are entries of Σ and ξ_{ij} adds a nugget term when ℓ=ℓ′.

Two key theoretical results are proved. First, the marginal covariance C_{ii} reduces to σ_{ii} ρ_i(ℓ,ℓ′) whenever at least one of the two locations belongs to S (or they coincide). Consequently, the marginal parameters are directly those of the chosen univariate correlation functions, making prior specification and interpretation straightforward. Second, the cross‑covariance is bounded above by σ_{ij}; when ρ_i=ρ_j the bound is attained at zero distance, i.e., C_{ij}(ℓ,ℓ)=σ_{ij}. This property mirrors the role of Σ in multivariate Matérn models but without imposing the complex validity constraints that those models require.

IOX’s construction permits each outcome to have its own range, smoothness, and nugget, thereby handling unequal smoothness and non‑stationarity naturally. Non‑stationarity arises because the covariance depends on the reference set S; as S grows, the model approximates the original stationary correlation, but for prediction locations outside S the covariance is effectively a predictive‑process covariance, which is inherently non‑stationary. Dimension reduction is achieved by imposing low‑rank structure on Σ or clustering outcomes, allowing the method to scale to high‑dimensional q.

Scalability to large n is obtained by embedding IOX within Vecchia‑type approximations. For each ρ_i, only m nearest neighbours are used to construct a sparse Cholesky factor L_i, reducing storage and computation to O(n m²). The resulting sparse directed acyclic graph (DAG) can be incorporated into Bayesian hierarchical models, enabling efficient Gibbs or Metropolis‑in‑Gibbs samplers for posterior inference.

Empirical evaluation includes synthetic experiments with up to 40 000 locations and 3–4 outcomes, where IOX outperforms LMC and full multivariate Matérn models in predictive accuracy, log‑likelihood, and computational time. A real‑world application to colorectal cancer proteomics data (hundreds of proteins measured on thousands of tissue samples) demonstrates that IOX can capture biologically meaningful cross‑protein spatial patterns while providing interpretable estimates of range, smoothness, and nugget for each protein.

The authors release an R package “spiox” (github.com/mkln/spiox) implementing likelihood‑based inference, Vecchia approximations, and posterior sampling for IOX models. In summary, IOX offers a flexible, interpretable, and scalable alternative to existing multivariate spatial covariance models, bridging the gap between full‑rank Gaussian processes and computationally tractable approximations, and opening new avenues for high‑dimensional spatial analysis in ecology, epidemiology, and omics studies.