SpaceBF: spatial coexpression analysis using Bayesian fused approaches in spatial omics datasets

TL;DR

SpaceBF is a new Bayesian method for analyzing spatial co-expression patterns in spatial omics data, revealing cancer-related molecular interactions with high accuracy.

Contribution

Introduces SpaceBF, a Bayesian fused modeling framework for robust spatial co-expression analysis in spatial omics.

Findings

SpaceBF outperforms existing methods in specificity and power in simulations.

The method identifies cancer-relevant molecular interactions in real spatial omics datasets.

SpaceBF preserves spatial structure while allowing for large edge-specific differences.

Abstract

Advances in spatial omics enable measurement of genes (spatial transcriptomics) and peptides, lipids, or N-glycans (mass spectrometry imaging) across thousands of locations within a tissue. While detecting spatially variable molecules is a well-studied problem, robust methods for identifying spatially varying co-expression between molecule pairs remain limited. We introduce SpaceBF, a Bayesian fused modeling framework that estimates co-expression at both local (location-specific) and global (tissue-wide) levels. SpaceBF enforces spatial smoothness via a fused horseshoe prior on the edges of a predefined spatial adjacency graph, allowing large, edge-specific differences to escape shrinkage while preserving overall structure. In extensive simulations, SpaceBF achieves higher specificity and power than commonly used methods that leverage geospatial metrics, including bivariate Moran’s I and Lee’s L. We also benchmark the proposed prior against standard alternatives, such as intrinsic conditional autoregressive and Matérn priors. Applied to spatial transcriptomics and proteomics datasets, SpaceBF reveals cancer-relevant molecular interactions and patterns of cell–cell communication (e.g., ligand–receptor signaling), demonstrating its utility for principled, uncertainty-aware co-expression analysis of spatial omics data.

Figures

Cutaneous melanoma data analysis. (A) Annotated H&E-stained image. (B) Cell types based on gene expression. (C) Minimum spanning tree (MST) capturing the spatial structure. (d) Clustering of spots based on centered and scaled estimates of slope surfaces of 53 statistically significant LR pairs. (E) The 3 main spatial patterns of the estimated surfaces. (F) Enrichment of LR interactions in 3 major cell types, with LR names arranged and color-coded according to their respective patterns. (G) The first 2 columns show the expression of 3 LR pairs. The third column displays the centered and scaled slope surfaces. In the fourth column, insignificant spot-level slope estimates are grayed. Figure 1A, https://zenodo.org/records/8215682 (reference 230), it has CC4 license.

cSCC data analysis. (A) H&E-stained image. (B) MST capturing the spatial structure. (C) Spatial clusters obtained using the BayesSpace package. (d) The 3 main spatial patterns of the estimated surfaces. (E) Bipartite graphs between Type 1 and Type 2 keratins based on their spatial pattern. (F) Study of the binding between the Type 1 keratin KRT17 and 3 different Type 2 keratins, with each slope surface exhibiting a unique spatial pattern. The insignificant spot-level slope estimates are grayed in the last column. Figure 2A, https://zenodo.org/records/8215682 (reference 231), it has CC4 license.

DCIS data analysis. (A) H&E-stained image. (B) Patterns of standardized co-expression (slope) of 66 peptide pairs (from 12 peptides). (C) Peptide description and dendrogram corresponding to the patterns. Mean slope estimates are presented as a heatmap on the right. (d) Clustering of spots based on the slope surfaces. (E) Scaled expression of the peptide 1125 forming the yellow-bordered module in the dendrogram. (F) Scaled expression of 3 peptides belonging to the same module (top row) and their spatial co-expression with peptide 1125 (bottom row).

Comparison of global tests under the simulation design 1. (A) Heatmap of the spatial covariance matrix H with an exponential kernel and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $L^2$\end{document} distance, for varying values of the lengthscale parameter l. (B) Simulated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\bf X} ^m$\end{document} based on Equation (2), for a fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\bf X} ^{m^{\prime }}$\end{document} but different values of the constant slope \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\beta _1^{mm^{\prime }}(s) = \nu$\end{document}. (C) Performance of the methods in terms of Type 1 error (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu = 0$\end{document}) and power (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu \ne 0$\end{document}). The dotted line represents the significance level 0.05.

Comparison of global tests under simulation design 2 for lengthscale l between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\lbrace 0.6, 1.8, 3.6, 7.2\rbrace$\end{document}.

(A) Power comparison of spatial priors under simulation design 2 for lengthscale l between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\lbrace 1.8, 3.6, 7.2\rbrace$\end{document}. (B) MSE comparison of spatial priors under simulation design 3, linear partition boundary. (C) MSE comparison of spatial priors under simulation design 3, circular boundary. In panel (A), sdmTMB models are omitted due to recurrent convergence issues.

SVC simulation based on linear partition boundary from the “Simulation design 3” section with effect size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu = 4$\end{document}. The boundary pattern changes based on the parameter r: (A) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r = 1$\end{document}, more zeros on the lower left, and (B) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r = 2$\end{document}, more negative slope values on the lower left.

SVC simulation based on circular boundary from the “Simulation design 3” section with effect size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu = 4$\end{document}. (A) A smaller circle with positive \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu$\end{document} values, radius \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r = 4$\end{document}, and (B) a larger circle with positive \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu$\end{document} values, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r = 6$\end{document}.

Graphical summary of the proposed approach.

Run-time comparison of SpaceBF, with different priors: the horseshoe GMRF and ICAR, with varying spatial adjacency graphs.

Tables

Table 1. Comparison of the methods in terms of the underlying assumptions.

Method	Central metric or concept	Global and local	Test type	Potential sensitivity
		co-expression tests
SpaceBF	Spatially varying	Global and local	Exact	Spatial adjacency graph
	coefficients model
	with NB distribution
MERINGUE	Bivariate Moran’s I	Global	Permutation test	Spatial adjacency graph
SpatialDM	Bivariate Moran’s I	Global and local	Permutation or	Lengthscale in the
			exact test	kernel weight matrix
LIANA+	Bivariate Moran’s I	Global and local	Permutation test	As above
	and cosine similarity
Voyager	Lee’s L	Global and local	Permutation test	Spatial adjacency graph
SpaGene	kNN network and	Global but with local	Permutation test	Choice of k
	earth mover’s distance	interaction scores
PearsonCorr	Pearson correlation	Global	Permutation or	None
			exact test
SpatialCorra	Spatial kernel-weighted	Global and local	Permutation test	Lengthscale in the
	sample correlation			kernel weight matrix
Copulaccia	Bivariate Poisson	Global but with local	Permutation test	Spatial adjacency graph
	distribution along an	interaction scores		and fixed correlation
	adjacency graph			term across locations

Table 2. Methods, graphs, priors, and limitations. HS can be placed on any spatial adjacency graph; denser graphs (e.g., Delaunay or large-k kNN) increase the risk of over-smoothing. HS-Del and HS-kNN (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k{=}3$\end{document}) probe denser and intermediate backbones than HS-MST. ICAR variants are included for comparison. The sdmTMB Matérn models differ by mesh density; a coarser mesh leads to a smoother spatial field.

Method	Graph	Prior	Limitations
HS-MST	Minimum spanning tree (MST)	Horseshoe (HS)	May under-smooth within dense micro-domains
HS-Del	Delaunay triangulation (Del)	HS	Risk of over-smoothing and spurious association; approximate local scale update
HS-kNN	k-nearest neighbors (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k{=}3$\end{document})	HS	Sensitive to k: larger k increases density and risk of over-smoothing
ICAR-Del	Delaunay triangulation	ICAR (special case of HS)	Risk of diffusing boundaries on dense graphs; tends to attenuate sharp slope changes.
ICAR-kNN	k-nearest neighbors	ICAR	Similar boundary diffusion; performance depends on k
sdmTMB-Matérn1	Mesh (cutoff \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $=1$\end{document}; denser)	Matérn (SPDE)	Risk of diffusing boundaries; mesh tuning adds burden, Laplace approximation may not converge
sdmTMB-Matérn2	Mesh (cutoff \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $=1.5$\end{document}; coarser)	Matérn (SPDE)	Above problems and risk of oversmoothing

Equations

$$\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} X^{m}(s_k) &=& \beta ^{mm^{\prime }}_{0}(s_k) + X^{m^{\prime }}(s_k) \beta _{1}^{mm^{\prime }} (s_k) + C(s_k)^\top \alpha _m + \epsilon (s_k), \\ &&\quad k = 1, \ldots , n,\\ \end{eqnarray*}\end{document}$$

$$\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \eta _m(s_k) &=& \beta ^{mm^{\prime }}_{0}(s_k) + X^{m^{\prime }}(s_k) \beta _{1}^{mm^{\prime }} (s_k) + C(s_k)^\top \alpha _m \\ p(X^{m}(s_k)|\psi _m(s_k), r_m) &\propto & (1 - \psi _m(s_k))^{r_{m}} \psi _m(s_k)^{X^{m}(s_k)}, \\ &&\psi _m(s_k) = \frac{\exp (\eta _m(s_k))}{1 + \exp (\eta _m(s_k))}, \end{eqnarray*}\end{document}$$

$$\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \pi (\pmb {\beta }^{mm^{\prime }}_0|\ldots ) &\propto & \prod _{{i \in E}} \exp \left( -\frac{\lambda _0}{\sigma }\Delta \beta ^{(0)}_{i}\right), \\ &&\pmb {\beta }^{mm^{\prime }}_0 = \big(\beta ^{mm^{\prime }}_{0}(s_1), \ldots , \beta ^{mm^{\prime }}_{0}(s_n)\big)^\top , \\ \pi (\pmb {\beta }^{mm^{\prime }}_1|\ldots ) &\propto & \prod _{{i \in E}} \exp \left( -\frac{\lambda _1}{\sigma }\Delta \beta ^{(1)}_{i} \right), \\ &&\pmb {\beta }^{mm^{\prime }}_1 = \big(\beta ^{mm^{\prime }}_{1}(s_1), \ldots , \beta ^{mm^{\prime }}_{1}(s_n)\big)^\top , \end{eqnarray*}\end{document}$$

$$\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \Delta \beta ^{(0)}_{i}\mid \Lambda ^{2}_{0i},\tau _{0}^{2},\sigma ^{2} &\sim & N \big (0,\ \Lambda ^{2}_{0i}\, \tau _{0}^{2}\, \sigma ^{2}\big ),\Lambda _{0i}\sim C^{+}(0,1),\\ &&\quad\tau _{0}\sim C^{+}(0,1),\\ \Delta \beta ^{(1)}_{i}\mid \Lambda ^{2}_{1i},\tau _{1}^{2},\sigma ^{2} &\sim & N \big (0,\ \Lambda ^{2}_{1i}\, \tau _{1}^{2}\, \sigma ^{2}\big ),\Lambda _{1i}\sim C^{+}(0,1),\\ &&\quad\tau _{1}\sim C^{+}(0,1), \end{eqnarray*}\end{document}$$

$$\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \pi\big (\boldsymbol{\beta }^{mm^{\prime }}_j \mid \Lambda _j,\tau _j^{2},\sigma ^{2}\big ) &\propto & \big(\tau _j^{2}\sigma ^{2}\big)^{-\ell _{\mathrm{rank}}/2}\, \exp \left\lbrace -\frac{1}{2\, \tau _j^{2}\sigma ^{2}} \sum _{i=1}^{p}\frac{\big (\Delta \beta ^{(j)}_{i}\big )^{2}}{\Lambda ^{2}_{ji}}\right\rbrace \\ &=& \big(\tau _j^{2}\sigma ^{2}\big)^{-\ell _{\mathrm{rank}}/2}\, \exp \left\lbrace -\frac{1}{2\, \tau _j^{2}\sigma ^{2}}\, {\boldsymbol{\beta }^{mm^{\prime }}_j}^{\top } L(\Lambda _j)\boldsymbol{\beta }^{mm^{\prime }}_j\right\rbrace ,\\ \end{eqnarray*}\end{document}$$

Full text

Introduction

Technological advances in spatial omics [1–3] have enabled in situ profiling of varying molecules, including genes (via spatial transcriptomics [ST]) [4–7], lipids, or peptides (using mass spectrometry imaging [MSI]) [8–11], and immune proteins (through multiplex immunofluorescence [mIF]) [12–15], within tissues. The technologies offer distinct yet complementary biological insights, differing in spatial resolution and the number of detectable molecules (throughput). For example, the next-generation sequencing-based ST platform Visium (from 10X Genomics) [16] offers transcriptome-wide gene-expression profiling (throughput \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sim 20,000\end{document} ) at a 55 µm spot-level resolution. MALDI MSI-based platforms (from Bruker Daltonics [17] and others) offer profiling of different types of molecules, such as peptides, lipids, nucleotides, proteins, metabolites, and N-glycans (throughput ∼50–1,000) at 10 µm spot-level resolution. The mIF platform PhenoCycler (from Akoya Biosciences) [18] enables protein profiling (throughput \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sim 40\end{document} ) at a 0.6 µm cellular resolution. Despite these differences, the underlying data structure remains largely consistent across technologies and platforms, comprising a collection of spatial locations (from single or multiple samples) with observed expression or intensity of various molecules. Consequently, common biostatistical questions arise, centering on the spatial dynamics of molecules within the complex tissue or tumor microenvironment (TME) [19–23].

In the context of ST datasets, identifying spatially variable genes (SVGs), i.e., the genes exhibiting spatially structured expression patterns across the tissue, has gained significant attention [24–38]. It enables critical downstream analyses such as discovering potential biomarkers and defining tissue regions that influence cellular differentiation and function [39–42]. Analogously, for mIF or imaging mass cytometry (IMC) datasets, innovative methods [43–51] have been proposed to understand the spatial distribution of immune cell types (defined by binarizing the expression profile of immune proteins) across the TME. Building upon this univariate framework, which typically analyzes one molecule at a time, another widely investigated problem has been spatial domain detection, i.e., deconvolving the tissue into distinct, spatially contiguous neighborhoods based on multivariate gene expression (ST) [52–63] or immune cell type composition (mIF) [64–69]. It aids in mapping the molecular and functional landscape of tissues, elucidating disease progression, and guiding targeted therapies [70–72]. While some of the referenced methods can be adapted for use with MSI datasets, it is important to underscore the lack of sophisticated spatial functionalities of the existing bioinformatics toolboxes [73–76].

While univariate and multivariate spatial analyses have garnered significant attention, a critical intermediate task remains underexplored: bivariate spatial co-expression analysis of molecular pairs at both “local” (spot/cell-specific) and “global” (tissue-wide) levels, aimed at precisely characterizing the spatial interaction or binding pattern of any 2 molecules throughout the tissue plane. To emphasize the importance of such an analysis, we review the concepts of cell–cell communication (CCC) [77–80]. CCC is a fundamental biological process through which cells exchange information via direct contact or signaling molecules (ligands) binding to receptor molecules present on the same or different cells. It regulates essential biological functions, including tissue development [81] and immune responses [82], and its disruption has been implicated in the onset and progression of cancer [83]. Autocrine, juxtacrine, and paracrine signaling are 3 major pathways of CCC [84]. In autocrine signaling, ligands released by a cell bind to receptors on the same cell, while in juxtacrine and paracrine signaling, the ligands target adjacent and nearby cells. The study of ligand–receptor interactions (LRI), which involves identifying gene pairs (ligands and receptors) that show coordinated upregulation or downregulation across groups of cells, has become a fundamental approach for inferring CCC from single-cell RNA sequencing (scRNA-seq) datasets [85–93]. However, these approaches are prone to false positive interactions due to the lack of spatial context in scRNA-seq datasets, treating distant cell pairs similarly to nearby ones [94–96], which potentially leads to an overestimation of juxtacrine and paracrine signaling. ST datasets offer a natural avenue for improvement by enabling spatially constrained LRI analysis.

A limited number of tools exist for spatial LRI analysis or, more broadly, for assessing bivariate spatial co-expression of molecules in ST or MSI datasets. It should be emphasized that bivariate co-expression can manifest in 2 ways: (a) joint over- or under-expression within the same cells (correlation) and (b) joint over- or under-expression in neighboring cells (cross-correlation [97]). Some relevant methods include MERINGUE [98], Giotto [99], SpaGene [100], SpaTalk [101], SpatialDM [102], CellChat V2 [103], LIANA+ [104], and Copulacci [105]. We skip the approaches that jointly analyze multiple LR pairs [106, 107]. Methods such as MERINGUE, Giotto, and SpaTalk provide only a global summary of spatial co-expression across a tissue, whereas others also offer local (spot/cell-specific) estimates. Let the standardized expression of 2 genes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (m, m^{\prime })\end{document} be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^m(s)\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^{m^{\prime }}(s)\end{document} at location s for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} s \in \lbrace s_1, \ldots , s_n\rbrace\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^m = (X^m(s_1), \ldots , X^m(s_n))^\top\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^{m^{\prime }} = (X^{m^{\prime }}(s_1), \ldots , X^{m^{\prime }}(s_n))^\top\end{document} . For a global summary of spatial co-expression, MERINGUE and SpatialDM leverage a popular geospatial metric termed the bivariate Moran’s I ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} I_{BV}\end{document} ) [108, 109], interpreted as the Pearson correlation between one variable and the spatial lagged version of the other [110–112]. Mathematically, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} I_{BV} \propto (X^m)^\top WX^{m^{\prime }}\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} W = [[w_{k_1k_2}]]\end{document} is the spatial weight matrix that controls the spatial lagging. As W, MERINGUE uses a binary adjacency matrix based on the Delaunay triangulation [113] of the spatial locations ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} w_{k_1k_2} = 1\end{document} if locations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (s_{k_1}, s_{k_2})\end{document} are connected, or 0 otherwise). SpatialDM uses a kernel covariance matrix or Gram matrix [114] based on the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} L^2\end{document} distance between locations ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} w_{k_1k_2} = k_l(|s_{k_1}-s_{k_2}|^2)\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} k_l\end{document} is a kernel function with lengthscale parameter l [115]). For local estimates of spatial correlation, SpatialDM considers the bivariate local Moran’s I ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} I_{BV}^{local}(s)\end{document} ) based on the local indicators of spatial association approach [116]. The LIANA+ toolbox implements SpatialDM and introduces a similar spatially weighted cosine similarity index. Of note, a newer package named Voyager [117] considers Lee’s L statistic [110], which has a slightly different formulation than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} I_{BV}\end{document} . A critical yet often overlooked aspect of ST data analysis is that gene expression, measured in terms of unique molecular identifier (UMI) count, is inherently a discrete random variable (RV). However, the above methods assume normality upon a variance-stabilizing transformation [118–120], which may obscure true signals and have been widely criticized both within the ST literature [26, 55, 121, 122] and in broader contexts [123–126]. Addressing this issue, Copulacci models a pair of genes as bivariate Poisson-distributed RVs, with their correlation in spatially adjacent cells captured using a Gaussian copula [127]. For inference, these methods typically rely on a permutation test [128].

Bivariate Moran’s I ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} I_{BV}\end{document} ) and Lee’s L statistic, as implemented in MERINGUE, SpatialDM, LIANA+, and Voyager, are primarily recommended as exploratory metrics for assessing cross-correlation rather than as rigorous hypothesis testing tools [97, 129], in traditional spatial statistical literature. In simulation studies (see the “Simulation studies” section), we have shown that even when 2 variables are independently simulated with certain spatial covariance structures, the unmodeled spatial autocorrelation introduces a confounding effect on the bivariate association, leading to significantly inflated Type 1 error rates. A similar issue is well documented, as extensive literature highlights the limitations of using simple Pearson correlation to assess dependencies between 2 variables in the presence of spatial autocorrelation [130–134]. By extension, since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} I_{BV}\end{document} and Lee’s L are both fundamentally based on Pearson correlation between spatially lagged variables, they may be susceptible to similar pitfalls. In addition, we show that although the asymptotic mean of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} I_{BV}\end{document} tends to 0 under the null hypothesis of independence, its asymptotic variance can be large when the marginal spatial autocorrelation patterns of the two molecules are aligned. Consequently, in real datasets, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} I_{BV}\end{document} may take spuriously large values even under independence. Further, these spatially weighted association indices, being model-free, are unable to seamlessly adjust for cell-level covariates such as cell type, a limitation also present in Copulacci. As a side note, mapping to the aforementioned CCC pathways, Pearson correlation between ligand and receptor can be interpreted as a proxy for autocrine signaling, while cross-correlation may reflect a combination of juxtacrine and paracrine signaling.

We approach the bivariate spatial co-expression detection as a generalized linear regression problem, modeling a molecule m as the outcome and the other molecule \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} m^{\prime }\end{document} as the predictor (see the “Methods” section). For ST datasets, gene expression or UMI count is modeled as an overdispersed negative binomial (NB)-distributed RV [135], while an alternative Gaussian model is considered for continuous cases. The regression coefficients, both intercept ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{0}(s)\end{document} ) and slope ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s)\end{document} ), are assumed to vary across locations (s) exhibiting spatial dependency. Known as the spatially varying coefficients (SVCs) model [136], this framework provides exceptional flexibility and precision in capturing locally changing co-expression patterns through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s)\end{document} . A large positive \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s)\end{document} suggests strong positive co-expression at location s, i.e., joint up- or down-regulation, whereas a large negative value indicates avoidance or repulsion. The average of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s)\end{document} ’s, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overline{\beta 1^{mm^{\prime }}} = \sum k \beta 1^{mm^{\prime }}(s_k)/n\end{document} , provides a summary of the global co-expression pattern. Similar models have been widely used in fields such as disease mapping [137, 138], econometrics [139, 140], ecological studies [141, 142], and neuroimaging research [143, 144]. In the Bayesian paradigm, the spatial dependency between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta 0^{mm^{\prime }}(s)\end{document} ’s and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta 1^{mm^{\prime }}(s)\end{document} ’s is typically modeled using a conditional autoregressive (CAR) [145–147] or Gaussian process (GP) priors [97, 148, 149]. In contrast, we introduce a locally adaptive spatial Gaussian Markov random field (GMRF) prior [150] based on the concepts of fusion penalties [151–153] and horseshoe prior [154–156], extending a related work in the frequentist setup [157]. Briefly, the prior incorporates the spatial similarity between 2 adjacent locations, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (s{k_1}, s{k_2})\end{document} , by encouraging \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} |\beta ^{mm^{\prime }}{0}(s{k_1}) - \beta ^{mm^{\prime }}{0}(s{k_2})| \approx 0\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} |\beta ^{mm^{\prime }}{1}(s{k_1}) - \beta ^{mm^{\prime }}{1}(s{k_2})| \approx 0\end{document} . We parameterize spatial adjacency primarily with the minimum spanning tree (MST) [158–160], following Li et al. [157], which provides cycle-free, globally economical connectivity by minimizing total edge weight. In practice, we find that modestly denser graphs, e.g., k-nearest neighbor with small k, can yield improved performance. We evaluate the proposed method SpaceBF against established approaches under realistic simulation scenarios, demonstrating high specificity and power. For broader applicability, we also benchmark our prior against the standard intrinsic CAR (ICAR) and a stochastic partial differential equation (SPDE)-based Matérn prior [161, 162], where SpaceBF consistently outperforms both alternatives. SpaceBF is applied to 3 real datasets: (a) an ST dataset on cutaneous melanoma [39] for spatial LRI analysis, (b) an ST dataset on cutaneous squamous cell carcinoma (cSCC) [163] for keratin-interaction analysis, and (c) a spatial proteomics dataset on ductal carcinoma in situ (DCIS) from the Medical University of South Carolina (MUSC) for peptide co-localization analysis.

Result

Real data analysis

We use the MST as the spatial adjacency graph for the real datasets in the main text and provide complementary kNN-based results in the Supplementary Material.

Melanoma ST dataset

We analyzed a cutaneous melanoma dataset [39] from a long-term survivor (10+ years), collected using the ST technology [4], comprising 293 spots, each 100 µm in size and at a 200 µm center-to-center distance. There are 16,148 genes, forming 1,180 known ligand–receptor (LR) pairs as available from CellChatDB [88]. There are 3 major pathologist-annotated regions as seen in the histology image (Fig. 1A), collected from Thrane et al. [39], and 6 major cell types (Fig. 1B) predicted using the RCTD [164] package based on overall gene expression [102]. After filtering out genes with extremely low expression ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} < 0.2 \times 293 \approx 59\end{document} reads), 161 LR pairs remain, which were examined using our method SpaceBF, without adjusting for any covariates. To briefly summarize the SpaceBF workflow (see Fig. 9), it first constructs an MST based on the spatial coordinates of the spots (Fig. 1C). Then, for every LR pair: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (m^{\prime }, m)\end{document} , it considers Equation (2) with the receptor expression as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^m(s_k)\end{document} and the ligand expression as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^{m^{\prime }}(s_k)\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} s_k\end{document} representing a spot. Following parameter estimation via a Markov Chain Monte Carlo (MCMC) procedure, the framework performs 2 hypothesis tests to assess the significance of spatial co-expression at both global and local levels (see the “Hypothesis testing” section). Using the global test in this dataset, SpaceBF identified 53 LR pairs at a significance level of 0.05 (33 at an false discovery rate (FDR) of 0.1). The estimated slope surface \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _1^{mm^{\prime }}(s_k)\end{document} of different LR pairs exhibits distinct patterns. To highlight these differences, we classify the detected LR pairs into 3 major patterns (Fig. 1E) based on hierarchical clustering [165] of the standardized vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{\boldsymbol \beta } _1^{mm^{\prime }} = (\beta _1^{mm^{\prime }}(s_1) - \overline{\beta _1^{mm^{\prime }}}, \ldots , \beta _1^{mm^{\prime }}(s_n) - \overline{\beta _1^{mm^{\prime }}})^\top /\sigma _{\beta }^{mm^{\prime }}\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overline{\beta _1^{mm^{\prime }}} = \sum _k \beta _1^{mm^{\prime }}(s_k)/n\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sigma _{\beta }^{mm^{\prime }}\end{document} are the tissue-wide average and the SD of estimated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _1^{mm^{\prime }}(s_k)\end{document} ’s, respectively. Twenty LR pairs follow pattern 1, while 22 and 11 LR pairs correspond to patterns 2 and 3, respectively. Similarly, the spots are grouped into 4 clusters based on the spot-level vectors of slopes corresponding to the 53 detected LR pairs (Fig. 1D). It is evident that clusters 1 and 3 correspond to the melanoma region, while clusters 2 and 4 loosely correspond to the stroma and lymphoid regions, respectively. Returning to the LR patterns, in Fig. 1F, the LR pairs are arranged sequentially from pattern 1 to 3, highlighting the enrichment of their interaction in 3 major cell types. For example, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sum _{k \in \text{B/T cells}} \beta _1^{mm^{\prime }}(s_k)\end{document} represents the enrichment within B/T cells relative to the average enrichment \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overline{\beta _1^{mm^{\prime }}}\end{document} and scaled by the SD. The levels “highest,” “medium,” and “lowest” indicate the degree of enrichment, with “highest” corresponding to the greatest or most positive enrichment and so on. The majority of LR pairs following pattern 1 exhibit higher or more positive interaction in B/T cells within the lymphoid region (some in CAF cells) and more negative interaction (avoidance or repulsion) in the melanoma region or cells. Pattern 2 mostly corresponds to LR pairs with the highest enrichment in CAF cells, while pattern 3 clearly corresponds to the pairs with the highest enrichment in melanoma cells. Next, we investigate the biological relevance of the estimated slope surfaces for a selected set of LR pairs. The LR pair (IGF2, IGF1R) [166] corresponds to pattern 1 and demonstrates a negative association overall, with an estimated average slope of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overline{\beta _1^{mm^{\prime }}} = -0.212\end{document} , and the P-value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} = 0.024\end{document} , which is consistent with a visual inspection (Fig. 1G). It could indicate a lack of binding between these genes, which would be a generally favorable factor for the survivor [167]. Setting the insignificant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _1^{mm^{\prime }*}(s_k)\end{document} values to 0 based on the local test, the negative interaction found in the melanoma region has the highest credibility. The pair (PTPRC, CD22) [168] follows pattern 2, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overline{\beta _1^{mm^{\prime }}} = 0.422\end{document} and P-value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} 6.28 \times 10^{-6}\end{document} . PTPRC, also known as CD45, is a facilitator of T-cell receptor (TCR) and B-cell receptor (BCR) signaling [169], while CD22 is primarily an inhibitor of BCR signaling [170]. Their overall positive co-expression, particularly in the lymphoid region, is likely associated with a balanced B-cell regulation, helping to prevent autoimmunity and promoting lymphoid growth in other regions as part of the immune response. The final LR pair we discuss is (SPP1, CD44) [171], which follows pattern 3, exhibiting a highly positive overall co-expression with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overline{\beta _1^{mm^{\prime }}} = 0.79\end{document} and P-value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} 1.03 \times 10^{-6}\end{document} . This strong interaction displays a decreasing gradient from the melanoma region to the lymphoid region, which aligns with its known role in dysregulated cytoskeletal remodeling [172], facilitating melanoma cell invasion into surrounding tissues.

Cutaneous melanoma data analysis. (A) Annotated H&E-stained image. (B) Cell types based on gene expression. (C) Minimum spanning tree (MST) capturing the spatial structure. (d) Clustering of spots based on centered and scaled estimates of slope surfaces of 53 statistically significant LR pairs. (E) The 3 main spatial patterns of the estimated surfaces. (F) Enrichment of LR interactions in 3 major cell types, with LR names arranged and color-coded according to their respective patterns. (G) The first 2 columns show the expression of 3 LR pairs. The third column displays the centered and scaled slope surfaces. In the fourth column, insignificant spot-level slope estimates are grayed. Figure 1A, https://zenodo.org/records/8215682 (reference 230), it has CC4 license.

cSCC ST dataset

We analyzed a cSCC dataset [163] on a patient sample with a histopathologic subtype of “moderately differentiated” cSCC [173]. The dataset was collected using the ST technology with 621 spots, each of size 110 µm and a center-to-center distance of 150 µm. There are 16,643 genes, of which 45 are keratins (14 after filtering low-count genes, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} < 0.2 \times 621 \approx 124\end{document} reads). These keratins can be classified into 2 types: (1) Type 1, which includes KRT10, KRT14-KRT17, and KRT23, and (2) Type 2, which includes KRT1, KRT2, KRT5, KRT6A, KRT6B, KRT6C, KRT78, and KRT80. The keratins pair together to form intermediate filaments, providing structural support to epithelial cells [174]. In the context of cSCC and other carcinomas, keratins are emerging as highly significant targets for therapeutic intervention [175–177]. Of note, some of the keratins belong to the GO term “keratinocyte differentiation” (GO:0030216) and were reported to exhibit strong spatial correlation in an earlier work [178] involving the same dataset. We utilized SpaceBF to investigate the binding between Type 1 and Type 2 keratins, resulting in a set of 48 keratin pairs. In the histology image (Fig. 2A), the deep blue areas at the top and left sides correspond to tumor regions, while the whitish region at the bottom represents a non-tumor region possibly composed of keratinized layers and stroma [179]. However, the tumor and non-tumor regions are not clearly delineated, a feature characteristic of moderately differentiated cSCC, though the spatial clusters obtained using the BayesSpace package [52] on the transcriptome-wide gene expression profile (Fig. 2C) partially elucidate this distinction. We emphasize that these clusters are shown for visualization only and are not used in our analysis. The constructed MST is shown in Fig. 2B. Using the global test, SpaceBF identified 39 keratin pairs at a significance level of 0.05 (41 at an FDR \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} < 0.1\end{document} ), suggesting that most pairs bind to each other, albeit to varying degrees. Similar to the earlier analysis, we classify the detected slope surfaces into 3 major patterns (Fig. 2D) based on hierarchical clustering of the standardized vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{\boldsymbol \beta } _1^{mm^{\prime }*}\end{document} . We represent the keratin pairs as bipartite graphs between Type 1 and Type 2 keratins under each pattern (Fig. 2E). One important observation is that the Type 2 keratins KRT6A, KRT6B, and KRT6C are closely related isoforms of keratin 6 [180] and therefore tend to be strongly co-expressed. Consequently, their spatial association patterns with a given Type 1 keratin are expected to be similar, consistent with the patterns recovered by SpaceBF. For instance, the slope surfaces of KRT10 with KRT6A, 6B, and 6C all align with pattern 1, while the slope surfaces of KRT16 with KRT6A, 6B, and 6C all correspond to pattern 2. This consistency underscores the reliability of SpaceBF in identifying true local patterns. In Fig. 2F, we present the estimated slopes for KRT17, which is a well-established therapeutic target in various cancers [181–183], binding with 3 Type 2 keratins: KRT80 (pattern 1, P-value = 0.006), KRT78 (pattern 2, P-value = 0.03), and KRT6B (pattern 3, P-value = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} 5.59 \times 10^{-6}\end{document} ). Notably, the average slope estimates for KRT17-KRT80 and KRT17-KRT78 interactions are small ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \approx 0.1\end{document} ), whereas for KRT17-KRT6B, the average slope is substantially higher at 0.82, with the highest local estimates observed mostly in tumor regions. These trends are also evident from the individual expression profiles provided in Fig. 2F. Although the expression patterns of KRT80 and KRT78 appear similar, a closer examination reveals that KRT80 exhibits a thicker band of expression on the left, specifically within the tumor regions. This distinction contributes to the difference in co-expression patterns of KRT17-KRT80 and KRT17-KRT78. As previously noted, both association levels are low, also indicated by the small number of significant spots identified by the local test, 72 and 49, respectively.

cSCC data analysis. (A) H&E-stained image. (B) MST capturing the spatial structure. (C) Spatial clusters obtained using the BayesSpace package. (d) The 3 main spatial patterns of the estimated surfaces. (E) Bipartite graphs between Type 1 and Type 2 keratins based on their spatial pattern. (F) Study of the binding between the Type 1 keratin KRT17 and 3 different Type 2 keratins, with each slope surface exhibiting a unique spatial pattern. The insignificant spot-level slope estimates are grayed in the last column. Figure 2A, https://zenodo.org/records/8215682 (reference 231), it has CC4 license.

DCIS proteomics dataset

We analyzed a single-sample DCIS dataset collected using the MALDI-MSI spatial proteomics platform, as part of an ongoing study at the MUSC aimed at defining the proteomic landscape of DCIS and invasive breast cancer (IBC), in terms of collagen peptides and immune cell types. DCIS is marked by the abnormal growth of malignant epithelial cells confined to the breast’s milk ducts, without invading the surrounding stromal tissue. While prognosis is excellent, around 20–40% of diagnosed DCIS progress to IBC [184, 185]. Understanding proteomic co-localization within the extracellular matrix (ECM) of a DCIS tissue is crucial for assessing progression risk and predicting therapeutic response, as the ECM plays a key role in regulating tumor cell proliferation, migration, and survival [186]. In this dataset, there are 5,548 tissue spots and 12 ECM peptides, whose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \big({{12}\atop{2}}\big)\end{document} pairwise interactions are of our interest. As the peptide expression is continuous-valued, we used the Gaussian model of SpaceBF for this analysis (Equation 1). Seven of the peptides are derived from the COL1A1 gene, while the remaining peptides originate from COL1A2, COL3A1, and FN1 (Fig. 3C). From the histology image (Fig. 3A), the stromal ECM can be identified by the light pink staining of fibrous connective tissue, while epithelium regions are highlighted in deep blue. Using hierarchical clustering, we group the standardized slope surfaces ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{\boldsymbol \beta } _1^{mm^{\prime }*}\end{document} ) into 3 patterns (Fig. 3B), and the spots into 3 clusters based on the spot-level vectors of standardized slopes (Fig. 3D). While the differences among the 3 patterns are subtle, the spot clusters are well defined and spatially distinct: the red cluster (cluster 3) aligns with stromal regions, the green cluster (cluster 2) with epithelial regions, and the light blue cluster (cluster 1) is a mixture of both. It is important to note that the MSI image has substantially lower resolution compared to the histology image, making one-to-one correspondence between the 2 inherently challenging. Patterns 1 and 2 (Fig. 3B), which visually resemble each other, both suggest strong co-localization of the associated peptide pairs in the stroma. This is expected, as all of these peptides are known to constitute the stromal ECM. The tree diagram in Fig. 3C shows the hierarchical relationships between the peptide pairs, with their patterns indicated on the right. The module highlighted by the yellow box includes pairs involving peptide 1125 (from COL1A2) and 7 other peptides. From Fig. 3E and F (top row), peptides 1125, 1212, 1386, and 1681 (from the module) show pronounced co-expression in the stromal region. Correspondingly, the estimated slope surfaces (Fig. 3F, bottom row) for the pairs (1125, 1212), (1125, 1386), and (1125, 1681) all fall under pattern 1, but the association strength is notably higher for (1125, 1212): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overline{ \mathrm{\boldsymbol \beta } _1^{mm^{\prime }}} = 0.92\end{document} , compared to 0.72 and 0.68 for the other 2. Although the existing literature on these interactions is limited, the findings will inform future comparative analyses of ECM compositions across DCIS subtypes and stages of progression [187, 188].

DCIS data analysis. (A) H&E-stained image. (B) Patterns of standardized co-expression (slope) of 66 peptide pairs (from 12 peptides). (C) Peptide description and dendrogram corresponding to the patterns. Mean slope estimates are presented as a heatmap on the right. (d) Clustering of spots based on the slope surfaces. (E) Scaled expression of the peptide 1125 forming the yellow-bordered module in the dendrogram. (F) Scaled expression of 3 peptides belonging to the same module (top row) and their spatial co-expression with peptide 1125 (bottom row).

Simulation studies

Simulation design 1: comparison between global methods under linear association

We consider the spatial coordinates ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} n = 293\end{document} ) from the cutaneous melanoma dataset. In simulation design 1, 1 NB-distributed RV, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} is generated using a Gaussian copula with a spatial covariance matrix H based on an exponential kernel (for varying lengthscale l) and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} L^2\end{document} distance. Another NB-distributed RV, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m}\end{document} is then generated using the NB model from Equation (2) with a constant slope \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta 1^{mm^{\prime }}(s) = \nu\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{\boldsymbol \beta } 0^{mm^{\prime }}\end{document} simulated using a GP model [97] with the spatial covariance matrix H. More details on the design are provided in the “Simulation design 1” section. From Fig. 4A, we notice how the structure of H changes as the lengthscale l varies. The off-diagonal elements of H ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} H{k_1k_2}\end{document} ) can range between 0 and 1. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} l = 0.6\end{document} , only the nearest locations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (k_1, k_2)\end{document} exhibit high \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} H{k_1k_2}\end{document} with most of the other values being close to 0. In contrast, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} l = 18\end{document} , the majority of location pairs have high \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} H_{k_1k_2}\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \approx 1\end{document} ), inducing an exceptionally strong spatial autocorrelation in both variables. To visibly understand how \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu\end{document} might affect the relationship between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m}\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} , in Fig. 4B, we show the spatial expression of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m}\end{document} for the same \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} but 3 different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lbrace -0.75, 0, 0.75\rbrace\end{document} . It is somewhat evident that nonzero \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu\end{document} ’s result in a visibly positive or negative association, while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu = 0\end{document} produces a random pattern of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^m\end{document} . In Fig. 4C, we show the Type 1 error ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu = 0\end{document} ) and power ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu \ne 0\end{document} ) comparison of the different methods, including SpaceBF, for 3 values of the lengthscale l. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} l = 3.6\end{document} , both variables exhibit considerable spatial autocorrelation, yet SpaceBF maintains the correct Type 1 error. In contrast, all other methods suffer from inflated Type 1 errors. Notably, simple Pearson correlation, while still inflated, performs better in controlling Type 1 error compared to methods based on bivariate Moran’s I or Lee’s L. Although SpaGene does not rely on these traditional metrics, it still fails to control Type 1 error. The issue becomes more pronounced as l increases. Notably, Lee’s \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} L\end{document} exhibits the highest inflation in the majority of cases. SpaceBF also retains a high detection power throughout all 3 cases. Although the power declines slightly for the largest l, as expected, due to a decrease in effective sample size from increased spatial autocorrelation. In summary, the simulation effectively demonstrates the specificity and power of our method.

Comparison of global tests under the simulation design 1. (A) Heatmap of the spatial covariance matrix H with an exponential kernel and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $L^2$\end{document} distance, for varying values of the lengthscale parameter l. (B) Simulated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\bf X} ^m$\end{document} based on Equation (2), for a fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\bf X} ^{m^{\prime }}$\end{document} but different values of the constant slope \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\beta _1^{mm^{\prime }}(s) = \nu$\end{document}. (C) Performance of the methods in terms of Type 1 error (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu = 0$\end{document}) and power (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu \ne 0$\end{document}). The dotted line represents the significance level 0.05.

Comparison of global tests under simulation design 2 for lengthscale l between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\lbrace 0.6, 1.8, 3.6, 7.2\rbrace$\end{document}.

Simulation design 2: comparison between global methods under non-linear association

In simulation design 2, we generate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ( {\bf X} ^{m}, {\bf X} ^{m^{\prime }})\end{document} jointly as bivariate spatially correlated NB-distributed RVs. This setup is more complex than the previous one, as the association is non-linear and driven by the Kronecker product-based spatial covariance structure (see the “Simulation design 2” section and Fig. 5). The methods, except SpaceBF, perform poorly in terms of the Type 1 error for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} l \ge 1.8.\end{document} When the spatial autocorrelation is the weakest ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} l = 0.6\end{document} ), Pearson correlation performs well as expected, while spatially weighted indices still show a slight inflation. SpaceBF achieves controlled Type 1 error and steady detection power across varying l’s. As earlier, the power decreases as the effective sample size decreases. Together, these 2 simulation designs demonstrate SpaceBF’s robustness under complex data generation processes. Finally, we argue that incorporating \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}_{0}(s)\end{document} in SpaceBF accounts for the spatial autocorrelation of variable m, thereby mitigating bias in the association analysis of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (m, m^{\prime })\end{document} . As previously noted, Pearson correlation is already recognized to be suboptimal in such scenarios, and spatially weighted indices, essentially Pearson correlation between spatially lagged variables, thus also remain susceptible to spurious detections. Recall that MERINGUE uses a binary spatial weight matrix W based on the Delaunay triangulation, while SpatialDM uses a continuous spatial weight matrix having a similar form as H for a particular choice of the lengthscale. Lee’s L is based on a binary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} k-\end{document} NN network in our study. Hence, the performance of these methods could be sensitive to choices of the spatial weight matrices, i.e., different networks or lengthscale l values. SpatialDM and SpaGene focus solely on the joint over-expression of molecules, neglecting joint under-expression. As a result, for most values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu < 0\end{document} , these methods show almost no detection power.

Simulation design 3: comparison between spatial priors under SVC framework

While the spatial horseshoe (HS) prior is introduced on an MST, it can be placed on any spatial backbone (e.g., Delaunay or kNN graphs), albeit with a potential risk of oversmoothing. This simulation study evaluates how graph choice affects HS performance. A Delaunay network is substantially denser than an MST, whereas a kNN network can serve as a middle ground for small k. In Fig. 6, HS-MST denotes HS on the MST (the original SpaceBF setting used in previous simulations and applications), HS-Del denotes HS on the Delaunay graph, and HS-kNN denotes HS on a kNN graph with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} k=3\end{document} . As noted in the “Methods” section, the ICAR prior is a special case of the HS prior; we therefore include ICAR-Del and ICAR-kNN for comparison. For completeness, we also consider a SPDE [161]-based NB SVC model implemented in the efficient R package sdmTMB [162], which uses a Matérn prior: sdmTMB-Matérn1 uses a denser mesh (cutoff \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} =1\end{document} ), and sdmTMB-Matérn2 uses a coarser mesh (cutoff \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} =1.5\end{document} ), see the Supplementary Material for a visual comparison.

(A) Power comparison of spatial priors under simulation design 2 for lengthscale l between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\lbrace 1.8, 3.6, 7.2\rbrace$\end{document}. (B) MSE comparison of spatial priors under simulation design 3, linear partition boundary. (C) MSE comparison of spatial priors under simulation design 3, circular boundary. In panel (A), sdmTMB models are omitted due to recurrent convergence issues.

We first revisit the simulation design 2 to assess global performance under a nonlinear model. As shown in Fig. 6A, all priors perform similarly across lengthscales, with a modest inflation in Type I error for HS-Del and ICAR-Del. This aligns with the oversmoothing tendency of denser backbones, which can spuriously propagate a small local positive association (e.g., confined to one corner of the tissue) across the entire domain. The overall similarity among priors is unsurprising in this particular simulation setting: with co-expression effectively constant over space, the additional flexibility of HS (an adaptive precision matrix) is not meaningfully exercised. Notably, sdmTMB-based methods failed to converge in several scenarios (and are therefore omitted here), especially when the outcome contained a high proportion of zeros. We attribute these failures to the near-vanishing curvature of the NB log-likelihood under heavy zero counts, which potentially destabilizes the Laplace-based frequentist optimization; in such cases, a Bayesian approach via R-INLA [189] may be preferable.

Next, we assess each prior’s ability to recover spatially varying slopes under the linear and circular boundary designs in the “Simulation design 3” section. In Fig. 6B, HS-kNN achieves the lowest mean squared error (MSE) across both linear-boundary settings (both values of r) and for both effect sizes ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu =2\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu =4\end{document} ). The gain is most pronounced at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu =4\end{document} , where its MSE is roughly half of the competing methods. By contrast, HS-MST, sdmTMB-Matérn1, and sdmTMB-Matérn2 generally perform the worst. For the circular boundary (Fig. 7C), HS-Del and HS-kNN perform similarly and outperform the remaining approaches. Overall, the MSE results favor the HS prior over standard ICAR models on generic spatial graphs and over the Matérn (SPDE-based) alternative. Because HS-MST is consistently outperformed by HS-kNN, we recommend using a denser graph such as kNN in typical applications. Moreover, Figs 6C and 7C show that HS-kNN best preserves sharp slope changes at boundaries, whereas ICAR methods (especially ICAR-Del) and sdmTMB methods (notably sdmTMB-Matérn2 with a coarser mesh) tend to blur these transitions, highlighting the HS prior’s ability to capture subtle co-expression shifts in the TME. Finally, while HS-Del can occasionally perform well, it may oversmooth and produce spurious detections, cautioning against the use of overly dense graphs such as Delaunay.

SVC simulation based on linear partition boundary from the “Simulation design 3” section with effect size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu = 4$\end{document}. The boundary pattern changes based on the parameter r: (A) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r = 1$\end{document}, more zeros on the lower left, and (B) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r = 2$\end{document}, more negative slope values on the lower left.

Discussion

We have developed a rigorous framework for studying spatial co-expression of a pair of molecules in the context of ST and MSI datasets, at both global (tissue-wide) and local (cell/spot-specific) levels. Existing tools mostly rely on 2 exploratory geospatial metrics, namely, bivariate Moran’s I [108] and Lee’s L [110], which lead to highly spurious association inference as demonstrated by our simulation studies. Our proposed approach, SpaceBF, builds on the widely used SVCs model [136], effectively capturing spatial autocorrelation and locally varying co-expression patterns. We introduce a spatial GMRF prior inspired by the fused horseshoe [156], with a global scale that regulates smoothness along graph edges and edge-specific local scales that allow large discontinuities to escape shrinkage. Setting the local scales to unity recovers the standard ICAR prior. We elucidate the prior’s theoretical properties and evaluate its performance against alternative spatial priors. Integrating the prior within both a Gaussian linear regression model and a more complex NB regression framework [135], SpaceBF is broadly applicable across various analytical contexts and data types.

We conduct a comprehensive evaluation of the proposed method under challenging simulation scenarios, demonstrating its ability to maintain well-controlled Type 1 error rates alongside strong detection power. In 3 real-world applications, 2 ST datasets and 1 MSI dataset, the method exhibits robust performance in identifying biologically meaningful molecular interactions, including LR signaling, keratin binding, and peptide co-localization. Notably, the analysis of the cutaneous melanoma sample reveals spatially variable patterns of CCC, as assessed through LR interactions. The LR pairs exhibit coordinated over- or under-expression within spatially distinct tissue regions (identified from histology) or specific cell types (inferred from transcriptome-wide gene expression). This level of granular understanding may provide critical insights for developing novel, targeted tissue-specific therapies in broader clinical settings [190–192].

Using the MST as the spatial graph offers several benefits: (i) uniqueness, removing the need to tune additional graph hyperparameters (e.g., GP lengthscales [97]); (ii) reduced computational burden via an exceptionally sparse precision matrix; and (iii) exact Gibbs updates for local horseshoe scales. In our simulations with spatial autocorrelation generated from a GP with an exponential kernel and varying lengthscales (but a domain-constant slope), the MST performs well, underscoring its robustness. Nonetheless, restricting the spatial structure to a single spanning tree can exclude salient edges [193], yielding noisier local slope estimates and overly sharp transition boundaries when coefficients vary spatially. In practice, a moderately denser graph, such as a kNN network with a small k, often achieves a better trade-off between computational efficiency and appropriate smoothness, as observed in our simulations. A more principled avenue could be to treat the spanning tree as unknown and update it iteratively within the model [194]. While we leverage the spam package [195] for fast sparse Cholesky factorization, overall complexity is graph-structure dependent (e.g., near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} O(n)\end{document} on trees/MSTs and typically around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} O(n^{3/2})\end{document} time for 2D planar/kNN graphs) [196]. As future work, we will pursue MCMC-free, variational-inference-based estimation to improve scalability [197, 198]. We have focused on pairwise analyses thus far; extending to joint modeling will follow prior works [199, 200]. Although we have primarily used SpaceBF in ST and MSI datasets, it could also be useful in mIF or IMC datasets where the molecular outcome of interest is generally immune cell types. To study cell type co-localization in such cases, one could split an mIF image into regular grids and count how many cells of 2 types \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (m, m^{\prime })\end{document} fall into each grid. Assuming that the grid centers are the locations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} s_k\end{document} ’s, a spatial graph can be constructed, and the spatial cell counts \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^m(s_k)\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^{m^{\prime }}(s_k)\end{document} can be analyzed using our framework.

Methods

Gaussian and NB regression models

We assume a single sample or image with n spots/cells. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^m(s_k)\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^{m^{\prime }}(s_k)\end{document} denote the expression of a pair of molecules m and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} m^{\prime }\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} C(s_k)\end{document} be a vector of p covariates observed at spot/cell location \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} s_k\end{document} , for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} k \in \lbrace 1, \ldots , n\rbrace\end{document} . For example, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} C(s_k)\end{document} can be the cell Type 1 indicator or a vector of cell type proportions [164]. We consider the following Gaussian SVCs model [136]:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} X^{m}(s_k) &=& \beta ^{mm^{\prime }}_{0}(s_k) + X^{m^{\prime }}(s_k) \beta _{1}^{mm^{\prime }} (s_k) + C(s_k)^\top \alpha _m + \epsilon (s_k), \\ &&\quad k = 1, \ldots , n,\\ \end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{0}(s_k)\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s_k)\end{document} denote spatially varying intercept and slope, respectively, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \alpha m\end{document} is a fixed effect vector, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \epsilon (s_k)\end{document} is an independent error term. To interpret the model, a significantly positive \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s_k)\end{document} implies that the molecules \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (m, m^{\prime })\end{document} co-express at the location \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} s_k\end{document} , while a significantly negative value suggests avoidance. Intuitively, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}_{0}(s_k)\end{document} accounts for spatial autocorrelation of molecule m. More discussion on the underlying bivariate spatial process is provided in the Supplementary Material. For a count-valued \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^{m}(s_k)\end{document} (e.g., genes in the ST datasets), we consider a spatially varying NB distribution [135] as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{ NB}(\psi m(s_k), r{m})\end{document} with the failure probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \psi _m(s_k)\end{document} modeled as [201]

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \eta _m(s_k) &=& \beta ^{mm^{\prime }}_{0}(s_k) + X^{m^{\prime }}(s_k) \beta _{1}^{mm^{\prime }} (s_k) + C(s_k)^\top \alpha _m \\ p(X^{m}(s_k)|\psi _m(s_k), r_m) &\propto & (1 - \psi _m(s_k))^{r_{m}} \psi _m(s_k)^{X^{m}(s_k)}, \\ &&\psi _m(s_k) = \frac{\exp (\eta _m(s_k))}{1 + \exp (\eta _m(s_k))}, \end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} p(.|.)\end{document} denotes the conditional probability mass function and the dispersion parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_m (> 0)\end{document} is assumed to be constant across locations. To explain how this framework effectively models overdispersion: as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_m \rightarrow \infty\end{document} , it reduces to a Poisson model; in contrast, as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_m \rightarrow 0\end{document} , the counts become increasingly dispersed relative to the Poisson distribution [135]. Admittedly, this model is limited as the count-valued nature of the molecule \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} X^{m^{\prime }}(s_k)\end{document} is not prioritized, appearing as a spatially varying predictor. Jointly modeling \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (X^{m}(s_k), X^{m^{\prime }}(s_k))\end{document} as bivariate NB (BNB) RVs is a possible approach that we do not pursue, as the existing definitions of the BNB distribution (outside of copula-based constructions) [202–205] vary considerably, often leading to restrictive correlation structures and inefficient MCMC sampling. To clarify, all applications in the manuscript assume no covariates, i.e., we do not include \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} C(s_k)\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \alpha m\end{document} , for simplicity. The models in Equations (1) and (2) are over-parameterized and do not incorporate spatial dependency between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{0}(s_k)\end{document} ’s and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s_k)\end{document} ’s, which we discuss next. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} G = (V, E)\end{document} denote the MST network between the locations constructed using the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} L^2\end{document} distance for a pair \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (s{k_1}, s_{k_2})\end{document} , where V and E are the sets of vertices and edges, respectively. Given a connected, weighted graph, an MST is an acyclic subgraph that connects all vertices and minimizes the sum of the weights of the included edges. Because of this property, MST is routinely used to develop transportation and telecommunication networks [206]. For a regular grid, MST is not unique, as the inter-point distances are not distinct. However, a unique random MST can be curated by simply adding small random values to the distances [160, 207]. Our simulation studies indicate that alternative adjacency graphs, such as k-nearest-neighbor (kNN) networks, can yield comparable performance and, in some settings, provide improved results.

Spatial modeling

Spatial fused lasso

In a recent study [157] of the temperature-salinity relationship in the Atlantic Ocean, Li et al. [157] consider Equation (1) and elegantly promote spatial homogeneity of the coefficients by considering fused lasso penalties [151, 208]: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} |\beta ^{mm^{\prime }}{0}(s{k_1}) - \beta ^{mm^{\prime }}{0}(s{k_2})| \approx 0\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} |\beta ^{mm^{\prime }}{1}(s{k_1}) - \beta ^{mm^{\prime }}{1}(s{k_2})| \approx 0\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (s_{k_1}, s_{k_2}) \in E\end{document} , in a frequentist setup. These constraints are intuitive, as it is reasonable to expect both the degree of co-expression, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s)\end{document} , and the effect of “unmeasured” factors, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{0}(s)\end{document} , to remain homogeneous across adjacent or connected locations. Extending this idea, a Bayesian fused lasso [153, 209] approach can be considered with Laplacian priors on the pair-wise differences of the coefficients. For a pair of locations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (s_{k_i^1}, s_{k_i^2})\end{document} connected by edge \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} i \in E\end{document} , let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \Delta \beta ^{(j)}{i} \equiv \beta ^{mm^{\prime }}{j}(s_{k_i^1}) - \beta ^{mm^{\prime }}{j}(s{k_i^2})\end{document} denote the difference, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} j\in \lbrace 0,1\rbrace\end{document} . The fused lasso prior can be imposed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \pi (\pmb {\beta }^{mm^{\prime }}_0|\ldots ) &\propto & \prod _{{i \in E}} \exp \left( -\frac{\lambda _0}{\sigma }\Delta \beta ^{(0)}_{i}\right), \\ &&\pmb {\beta }^{mm^{\prime }}_0 = \big(\beta ^{mm^{\prime }}_{0}(s_1), \ldots , \beta ^{mm^{\prime }}_{0}(s_n)\big)^\top , \\ \pi (\pmb {\beta }^{mm^{\prime }}_1|\ldots ) &\propto & \prod _{{i \in E}} \exp \left( -\frac{\lambda _1}{\sigma }\Delta \beta ^{(1)}_{i} \right), \\ &&\pmb {\beta }^{mm^{\prime }}_1 = \big(\beta ^{mm^{\prime }}_{1}(s_1), \ldots , \beta ^{mm^{\prime }}_{1}(s_n)\big)^\top , \end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sigma ^2\end{document} is the variance of the error term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \epsilon (s_k)\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lambda _0\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lambda _1\end{document} are regularization parameters that control the strength of fusion and are assumed to follow gamma priors. Note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sigma ^2\end{document} is only present in the Gaussian model (Equation 1) and could be omitted from the above exponents. We discuss the resemblance of the prior to the ICAR prior [146] and, more generally, the intrinsic GMRF (IGMRF) prior [150] in the Supplementary Material. Theoretically, using the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} L_1\end{document} distance seems appealing, as it has the potential to achieve better spatial smoothing by “exactly” fusing coefficient values at adjacent locations, unlike the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} L^2\end{document} distance implied by the ICAR prior. This is analogous to how lasso regression enforces sparsity in solutions, while ridge regression only shrinks effect sizes toward 0 [210]. For transparency, such a spatial fused lasso prior has already been proposed in the existing literature [194, 211].

Spatial fused horseshoe

In variable selection problems, failure of the Bayesian lasso or Laplacian prior to achieve exact sparsity, unlike the frequentist analog, has been reported while also underestimating larger effect sizes [212–214]. Consequently, the Bayesian fused lasso might struggle to promote spatial smoothness and preserve distinct local features simultaneously. For variable selection, the advantages of the horseshoe prior have been convincingly demonstrated to handle unknown sparsity and large outlying signals [154, 215, 216]. The horseshoe prior belongs to the class of global–local shrinkage priors [217], characterized by a “global” hyperparameter that controls overall shrinkage, while “local” hyperparameters control shrinkage per coefficient. Following recent developments on fused horseshoe priors [156, 218], we place horseshoe shrinkage on edgewise differences of the spatial coefficients. We specify

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \Delta \beta ^{(0)}_{i}\mid \Lambda ^{2}_{0i},\tau _{0}^{2},\sigma ^{2} &\sim & N \big (0,\ \Lambda ^{2}_{0i}\, \tau _{0}^{2}\, \sigma ^{2}\big ),\Lambda _{0i}\sim C^{+}(0,1),\\ &&\quad\tau _{0}\sim C^{+}(0,1),\\ \Delta \beta ^{(1)}_{i}\mid \Lambda ^{2}_{1i},\tau _{1}^{2},\sigma ^{2} &\sim & N \big (0,\ \Lambda ^{2}_{1i}\, \tau _{1}^{2}\, \sigma ^{2}\big ),\Lambda _{1i}\sim C^{+}(0,1),\\ &&\quad\tau _{1}\sim C^{+}(0,1), \end{eqnarray*}\end{document}

independently across edges \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} i=1,\dots ,p\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} |E| = p\end{document} (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} p = n - 1\end{document} for the MST). The error variance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sigma ^2 = 1\end{document} in the NB model. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} D\in \mathbb {R}^{p\times n}\end{document} be an oriented incidence matrix and define the weighted Laplacians \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} L(\Lambda _j)=D^\top !\mathrm{diag}(\Lambda _{j}^{-2}), D\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \Lambda _{j}^{-2} = {\lbrace \Lambda _{ji}^{-2}\rbrace _{i=1}^{p}}\end{document} being the vector of edgewise precisions. The above construction induces the following intrinsic GMRF prior:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \pi\big (\boldsymbol{\beta }^{mm^{\prime }}_j \mid \Lambda _j,\tau _j^{2},\sigma ^{2}\big ) &\propto & \big(\tau _j^{2}\sigma ^{2}\big)^{-\ell _{\mathrm{rank}}/2}\, \exp \left\lbrace -\frac{1}{2\, \tau _j^{2}\sigma ^{2}} \sum _{i=1}^{p}\frac{\big (\Delta \beta ^{(j)}_{i}\big )^{2}}{\Lambda ^{2}_{ji}}\right\rbrace \\ &=& \big(\tau _j^{2}\sigma ^{2}\big)^{-\ell _{\mathrm{rank}}/2}\, \exp \left\lbrace -\frac{1}{2\, \tau _j^{2}\sigma ^{2}}\, {\boldsymbol{\beta }^{mm^{\prime }}_j}^{\top } L(\Lambda _j)\boldsymbol{\beta }^{mm^{\prime }}_j\right\rbrace ,\\ \end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \ell _{\mathrm{rank}}=\mathrm{rank} (L(\Lambda _j))=n-C\end{document} , where C is the number of connected components of the graph ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \ell _{\mathrm{rank}}=n-1\end{document} for a connected graph). When all local scales are set to one, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \Lambda _{ji}\equiv 1\end{document} , the prior reduces to the standard ICAR prior (up to a scale factor) [97]. The deliberately omitted normalizing factor [150] (the generalized determinant of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} L(\Lambda _j)\end{document} ) depends on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \Lambda _j\end{document} but not on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \boldsymbol{\beta }^{mm^{\prime }}_j\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \tau _j^{2}\end{document} ; it therefore plays no direct role in the Gibbs updates for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \boldsymbol{\beta }^{mm^{\prime }}_j\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \tau _j^{2}\end{document} . In our implementation, we update the local scales \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \Lambda _{ji}\end{document} with conjugate per-edge steps, which are exact on trees and constitute a composite or pseudo-likelihood approximation on general graphs [219–221]. The fused horseshoe’s half-Cauchy global and local scales produce heavy tails (allowing large jumps across edges) and an infinitely tall spike at zero (aggressively shrinking small differences), thereby preserving spatial homogeneity while still accommodating sharp local variations. This connection between locally adaptive fusion and GMRFs was emphasized by Faulkner and Minin [222] in a longitudinal context, encouraging fusion between coefficients across time points. In the Supplementary Material, we provide the details of the Gibbs sampling steps, which include the Pólya–Gamma data augmentation strategy [201, 223] for the NB model (Equation 2).

One crucial aspect that deserves elucidation is the working assumption of independence between edge-wise differences. Specifically, for any 2 edges \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} i, i^{\prime }\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s{k_i^1}) - \beta ^{mm^{\prime }}{1}(s{k_{i}^2})\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta ^{mm^{\prime }}{1}(s{k_{i^{\prime }}^1}) - \beta ^{mm^{\prime }}{1}(s{k_{i^{\prime }}^2})\end{document} are assumed to be independent in Equation (4), conditional on the hyperparameters. To see how this assumption could be problematic for a general graph with cycles (i.e., not the MST), we briefly highlight 1 example from Rue and Held [150] provided in the context of IGMRF priors. Suppose there are only 3 locations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} A, B,\end{document} and C, all neighbors of each other. Letting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} e_1 = \beta ^{mm^{\prime }}{1}(A) - \beta ^{mm^{\prime }}{1}(B)\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} e_2 = \beta ^{mm^{\prime }}{1}(B) - \beta ^{mm^{\prime }}{1}(C)\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} e_3 = \beta ^{mm^{\prime }}{1}(C) - \beta ^{mm^{\prime }}{1}(A)\end{document} , Equation (4) proceeds to assume \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} e_1, e_2, e_3\end{document} are independent and normally distributed with non-identical parameters, yet there is a “hidden” linear constraint, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} e_1 + e_2 + e_3 = 0\end{document} , that contradicts independence. Analogously, using a highly connected or dense spatial neighborhood graph G introduces numerous hidden constraints corresponding to the cycles in G. Interestingly, as shown in Theorem 3 of the Supplementary Material, these constraints need not be enforced explicitly: the posterior sampling distribution of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \boldsymbol{\beta }^0\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \boldsymbol{\beta }^1\end{document} under the constraints is unchanged. However, penalizing too many edge-wise differences in the presence of implicit dependencies can lead to oversmoothing and the loss of salient local structure. This consideration naturally favors the MST, which is acyclic (hence no hidden constraints) and removes redundant relationships; moreover, a sparser G yields faster Cholesky factorizations of the precision matrix and thus improved computational efficiency. That said, as we show in the “Simulation design 3” section, a kNN graph with a small k can significantly outperform the MST in practice.

Hypothesis testing

We consider 2 types of hypothesis tests: (1) global test: to determine the significance of average association across the entire tissue domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ( H_0: \overline{\beta _1^{mm^{\prime }}} = \frac{1}{n}\sum _{k = 1}^n {\beta }_1^{mm^{\prime }}(s_k) = 0)\end{document} , based on the credible interval [224] of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overline{\beta _1^{mm^{\prime }}}\end{document} , and (2) local test: to determine the significance of location-level association \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} (H_0^k: {\beta }_1^{mm^{\prime }}(s_k) = 0)\end{document} , directly based on the credible intervals of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\beta }1^{mm^{\prime }}(s_k)\end{document} ’s. Additionally, in the genomic context, having a measure analogous to the frequentist P-value is often beneficial. To this end, we utilize a metric termed the probability of direction ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} p_d\end{document} ), which quantifies the probability (between 0.5 and 1) that a parameter has an effect in a specific direction, either positive or negative [225, 226]. Mathematically, it is defined as the proportion of the posterior distribution that shares the same sign as the median. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} p_d\end{document} resembles a 2-sided frequentist P-value as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} P{\mathrm{two-sided}} = 2(1 - p_d)\end{document} . It is implemented in the R package bayestestR [225]. For FDR control, we apply the Benjamini–Hochberg procedure using the p.adjust function in R.

Simulation design

We consider 2 different simulation designs as outlined below, for assessing the Type 1 error and power of the model proposed in Equation (2). The locations at which the variables are simulated are the same as the previously discussed cutaneous melanoma dataset ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} n = 293\end{document} ). We have observed that the results remain unaffected when a randomly generated set of locations or other real data-based sets of locations are used.

Simulation design 1

In the first design, we directly consider the model from Equation (2) to generate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ( {\bf X} ^m, {\bf X} ^{m^{\prime }})\end{document} based on 2 steps. First, we generate an NB-distributed RV, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} , using Gaussian copula [127], incorporating spatial dependency between the observations via a kernel covariance matrix H with an exponential kernel and varying lengthscale (l) parameters [114]. Then, based on the simulated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} , we generate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^m\end{document} following Equation (2) with a fixed slope \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta 1^{mm^{\prime }}(s_k) = \nu\end{document} . More specifically, for a fixed choice of l, failure probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \psi {m^{\prime }}\end{document} , dispersion parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r{m^{\prime }}\end{document} for variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} m^{\prime }\end{document} , and dispersion parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r{m}\end{document} for variable m, we consider the following steps:

Simulate a spatially autocorrelated normal RV of size n using a GP model:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \mathbf {Z}^{m^{\prime }} &\sim & MVN\left(\mathbf {0}, H\right), \\ H_{k_1k_2} &=& \exp \left(-\frac{||s_{k_1} - s_{k_2}||_1}{l}\right), ||.||_1 \text{denotes the} L^1 \mathrm{norm.} \end{eqnarray*}\end{document}

Transform to a vector of uniform RVs using the standard normal CDF ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \Phi\end{document} ):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} {\bf U} ^{m^{\prime }} = \Phi (\mathbf {Z}^{m^{\prime }}). \end{eqnarray*}\end{document}

Convert to a vector of NB RVs using the inverse CDF of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\mathrm{ NB}(\psi {m^{\prime }}, r{m^{\prime }})}\end{document} , denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} F^{-1}_{NB(\psi {m^{\prime }}, r{m^{\prime }})}\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} {\bf X} ^{m^{\prime }} = F_{\mathrm{ NB}(\psi _{m^{\prime }}, r_{m^{\prime }})}^{-1} ( {\bf U} ^{m^{\prime }}). \end{eqnarray*}\end{document}

Each element of the resulting vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} retains the marginal NB distribution, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\mathrm{ NB}(\psi {m^{\prime }}, r{m^{\prime }})}\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \psi {m^{\prime }}\end{document} is the failure probability and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r{m^{\prime }}\end{document} is the dispersion.Generate the link function to simulate variable m with a fixed slope of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _1^{mm^{\prime }}(s_k) = \nu\end{document} (2):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \mathrm{\boldsymbol \eta } _{m} = \mathrm{\boldsymbol \beta } _0^{mm^{\prime }} + \nu \log ( {\bf X} ^{m^{\prime }} + 1), \quad \mathrm{\boldsymbol \beta } _0^{mm^{\prime }} \sim MVN\left(\mathbf {0}, 0.5H\right). \end{eqnarray*}\end{document}

Convert the link vector to failure probabilities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \pmb {\psi }_m\end{document} and simulate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^m\end{document} from the NB distribution:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \pmb {\psi }_m = \frac{\exp ( \mathrm{\boldsymbol \eta } _{m})}{1+ \exp ( \mathrm{\boldsymbol \eta } _{m})}, {\bf X} ^m \sim \mathrm{ NB}(\pmb {\psi }_m, r_m), \end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_{m}\end{document} is a prefixed dispersion parameter. The kth element of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^m\end{document} follows \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{ NB}(\psi {mk}, r_m)\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \pmb {\psi }{m} = (\psi _{m1}, \ldots , \psi _{mn})^\top\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{\boldsymbol \eta } _{m} = (\eta _{m1}, \ldots , \eta _{mn})^\top\end{document} .

Three values of the lengthscale l are considered, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} l = 3.6, 7.2, 18\end{document} , with the corresponding structure of H displayed in Fig. 4A. The failure probability of variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} m^{\prime }\end{document} and dispersion parameters are kept fixed, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \psi {m^{\prime }} = 0.5\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_m = r{m^{\prime }} = 1\end{document} . The slope parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu\end{document} is varied between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lbrace -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75\rbrace\end{document} , with negative and positive values representing negative and positive association, respectively. Higher absolute value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu\end{document} dictates the strength of association, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu = 0\end{document} corresponds to the null model, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m}\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} are independent.

Simulation design 2

In this design, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ( {\bf X} ^m, {\bf X} ^{m^{\prime }})\end{document} are simulated jointly using a bivariate Gaussian copula and spatial dependency incorporated using a bivariate GP framework, where the joint covariance matrix has a Kronecker product structure, comprising a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} 2 \times 2\end{document} correlation matrix and the distance kernel covariance matrix H (Equation 9.11 from Banerjee et al. [97] and Equation 6 from the Supplementary Material). Specifically, we consider the following steps:

Simulate spatially cross-correlated normal RVs:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} (\mathbf {Z}^{m}, \mathbf {Z}^{m^{\prime }})^\top &\sim & MVN\left({\begin{bmatrix}\mathbf {0} \\ \mathbf {0} \end{bmatrix}}, \mathbf {\Sigma }= {\begin{bmatrix}1 & \nu \\ \nu & 1 \end{bmatrix}} \otimes H\right), \\ H_{k_1k_2} &=& \exp \left(-\frac{||s_{k_1} - s_{k_2}||_1}{l}\right). \end{eqnarray*}\end{document}

Transform to uniform RVs using the standard normal CDF ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \Phi\end{document} ):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} {\bf U} ^{m} = \Phi (\mathbf {Z}^{m}), \quad {\bf U} ^{m^{\prime }} = \Phi (\mathbf {Z}^{m^{\prime }}). \end{eqnarray*}\end{document}

Convert to NB RVs using the inverse CDFs:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} {\bf X} ^{m} = F^{-1}_{\mathrm{ NB}(\psi _{m}, r_{m})}( {\bf U} ^{m}), \quad {\bf X} ^{m^{\prime }} =F^{-1}_{\mathrm{ NB}(\psi _{m^{\prime }}, r_{m^{\prime }})} ( {\bf U} ^{m^{\prime }}). \end{eqnarray*}\end{document}

Each element of the resulting vectors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m}\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} retain the marginal NB distributions, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{ NB}(\psi {m}, r{m})\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{ NB}(\psi {m^{\prime }}, r{m^{\prime }})\end{document} , respectively, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \psi _{m}, \psi {m^{\prime }}\end{document} are failure probabilities and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_m, r{m^{\prime }}\end{document} are dispersions.

The lengthscale l is varied between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lbrace 0.6, 1.8, 3.6, 7.2\rbrace\end{document} . The failure probabilities and dispersion parameters are kept fixed, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \psi _{m} = \psi {m^{\prime }} = 0.5\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_m = r{m^{\prime }} = 1\end{document} . The parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu\end{document} is varied between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lbrace -0.75, -0.5, -0.25, 0,\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} 0.25, 0.5, 0.75\rbrace\end{document} , with similar implications on the direction and strength of association as before.

Simulation design 3

We consider the model from Equation (2) to generate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ( {\bf X} ^m, {\bf X} ^{m^{\prime }})\end{document} based on 2 steps. First, we generate an NB-distributed RV, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} , using Gaussian copula [127], incorporating spatial autocorrelation via H with a large lengthscale, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} l = 7.2\end{document} . Then, based on the simulated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} , we generate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^m\end{document} following Equation (2), this time with a spatially varying slope \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _1^{mm^{\prime }}(s_k)\end{document} . Specifically, we consider the following steps:

Simulate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} following steps 1, 2, and 3 from the “Simulation design 1” section.Simulate the slope \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{\boldsymbol \beta } 1^{mm^{\prime }}\end{document} in 2 ways: based on (a) linear partition boundary and (b) circular boundary. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} c_x=\mathrm{median}(s{k}^x)\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} c_y=\mathrm{median}(s_k^{y})\end{document} denote the median of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} xy\end{document} -coordinates, respectively.

• Linear partition boundary: Define 2 partitioning subsets of the spatial domain as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \mathcal {S}_1=\lbrace k:\ s_{k}^x > c_x,\ s_{k}^y\le c_y\rbrace ,\qquad \mathcal {S}_2 = \lbrace k:\ s_{k}^x\le c_x,\ s_{k}^y > c_y/r\rbrace . \end{eqnarray*}\end{document}

Draw \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sigma _k\sim \mathrm{Unif}(0.3,0.6)\end{document} . For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} k\in \mathcal {S}_1\end{document} , draw \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _{1k}^{mm^{\prime }}(s_k)\sim N(\nu ,\sigma _k^2)\end{document} ; for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} k\in \mathcal {S}_2\end{document} , draw \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _{1k}^{mm^{\prime }}(s_k)\sim N(-\nu ,\sigma _k^2)\end{document} ; for other k’s \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _{1k}^{mm^{\prime }}(s_k)=0\end{document} . Vary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r \in \lbrace 1, 2\rbrace\end{document} to create 2 different partitioning configurations.

• Circular boundary: Define the indicator for being inside a circle with radius \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r \in \lbrace 1, 2\rbrace\end{document} and centered at the median

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \mathbb {I}_k = \mathbf {1}\left\lbrace \frac{(x_k-c_x)^2+(y_k-c_y)^2}{r^2} \le 1 \right\rbrace . \end{eqnarray*}\end{document}

Set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sigma _k = \sigma _{\mathrm{in}}, \mathbb {I}_k + \sigma _{\mathrm{out}}, (1-\mathbb {I}_k)\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \sigma _{\mathrm{in}} = 0.3, \sigma _{\mathrm{out}}= 0.6\end{document} , and draw \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \varepsilon _k \sim N(0,\sigma _k^2)\end{document} . Finally, define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \beta ^{mm^{\prime }}_{1k}(s_k) \;=\; \nu \, (2\mathbb {I}_k-1) \;+\; \varepsilon _k . \end{eqnarray*}\end{document}

Thus, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathbb {E}[\beta ^{mm^{\prime }}_{1k}(s_k)\mid \mathbb {I}k=1]=\nu\end{document} (inside the circle) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathbb {E}[\beta ^{mm^{\prime }}{1k}(s_k)\mid \mathbb {I}_k=0]=-, \nu\end{document} (outside). Generate the link function to simulate variable m with the slope \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{\boldsymbol \beta } _1^{mm^{\prime }}\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \mathrm{\boldsymbol \eta } _{m} = \mathrm{\boldsymbol \beta } _0^{mm^{\prime }} + \mathrm{\boldsymbol \beta } _1^{mm^{\prime }} \log ( {\bf X} ^{m^{\prime }} + 1), \quad \mathrm{\boldsymbol \beta } _0^{mm^{\prime }} \sim MVN\left(\mathbf {0}, 0.5H\right), \end{eqnarray*}\end{document}

Convert the link vector to failure probabilities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \pmb {\psi }_m\end{document} and simulate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^m\end{document} from the NB distribution:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \pmb {\psi }_m = \frac{\exp ( \mathrm{\boldsymbol \eta } _{m})}{1+ \exp ( \mathrm{\boldsymbol \eta } _{m})}, {\bf X} ^m \sim \mathrm{ NB}(\pmb {\psi }_m, r_m), \end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_{m}\end{document} is a prefixed dispersion parameter. The kth element of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^m\end{document} follows \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{ NB}(\psi {mk}, r_m)\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \pmb {\psi }{m} = (\psi _{m1}, \ldots , \psi _{mn})^\top\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{\boldsymbol \eta } _{m} = (\eta _{m1}, \ldots , \eta _{mn})^\top\end{document} . Each element of the resulting vectors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m}\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\bf X} ^{m^{\prime }}\end{document} retain the marginal NB distributions, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{ NB}(\psi {m}, r{m})\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \mathrm{ NB}(\psi {m^{\prime }}, r{m^{\prime }})\end{document} , respectively, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \psi _{m}, \psi {m^{\prime }}\end{document} are failure probabilities and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_m, r{m^{\prime }}\end{document} are dispersions.

The dispersion parameters are kept fixed, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} r_m = r_{m^{\prime }} = 1\end{document} , and the failure probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \psi _{m^{\prime }} = 1\end{document} . The boundary parameter (or radius) r is varied between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lbrace 1, 2\rbrace\end{document} . The boundary (radius) parameter r varies over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lbrace 1,2\rbrace\end{document} , and the effect-size parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu\end{document} varies over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \lbrace 2,4\rbrace\end{document} . Results for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu =2\end{document} are provided in the Supplementary Material, while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \nu =4\end{document} results are shown in Figs 7 and 8.

SVC simulation based on circular boundary from the “Simulation design 3” section with effect size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu = 4$\end{document}. (A) A smaller circle with positive \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu$\end{document} values, radius \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r = 4$\end{document}, and (B) a larger circle with positive \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\nu$\end{document} values, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r = 6$\end{document}.

Graphical summary of the proposed approach.

Competing methods

We compare SpaceBF against 5 methods: (1) MERINGUE [98] with a Delaunay-triangulation graph; (2) SpatialDM [102] using a Gaussian kernel weight matrix (lengthscale 1.2, as in their original melanoma analysis); (3) Lee’s L [110] on an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \epsilon\end{document} -neighborhood network (implemented via the R package spdep [227]), where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \epsilon\end{document} is set to the maximum nearest-neighbor distance to ensure connectivity; (4) SpaGene [100]; and (5) PearsonCorr, the standard Pearson correlation. We do not include LIANA+ [104] or Voyager [117], as the former essentially wraps SpatialDM and the latter directly applies Lee’s L. We adopt an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \epsilon\end{document} -neighborhood graph rather than a fixed kNN graph to mitigate oversmoothing in Lee’s L while maintaining graph connectivity; however, this modification did not improve performance, as evidenced by our simulation studies. Table 1 summarizes the methods in terms of their assumptions and limitations. Note that, in the Supplementary Material, we derive the asymptotic mean and variance of the bivariate Moran’s I statistic and show that even under true independence, marginal spatial autocorrelation in each variable can inflate the estimated value. Moreover, this inflation is most pronounced when the autocorrelation patterns are aligned in the same direction. We also attempted to evaluate the performance of SpatialCorr [178] and Copulacci [105]. SpatialCorr was straightforward to use, but it proved highly sensitive to the choice of the lengthscale l in its innovative use of the spatial covariance matrix H. Copulacci was slightly difficult to use and will be benchmarked in a future study.

We further evaluate the performance of SpaceBF (the spatial horseshoe) on general spatial graphs, alongside standard spatial priors within our SVC framework (sections “Simulation design 3: comparison between spatial priors under SVC framework” and “Simulation design 3”). The configurations considered are summarized in Table 2, with brief notes on their limitations.

Table 2: Methods, graphs, priors, and limitations. HS can be placed on any spatial adjacency graph; denser graphs (e.g., Delaunay or large-k kNN) increase the risk of over-smoothing. HS-Del and HS-kNN (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k{=}3$\end{document}) probe denser and intermediate backbones than HS-MST. ICAR variants are included for comparison. The sdmTMB Matérn models differ by mesh density; a coarser mesh leads to a smoother spatial field.

Runtime comparison and convergence diagnostics

In most analyses, we ran \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} , 5{,}000,\end{document} MCMC iterations with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} , 2{,}500,\end{document} burn-in. We compared runtimes for our package SpaceBF across priors and spatial backbones (from sparser to denser). Figure 10 shows that HS and ICAR have comparable runtimes, scaling approximately linearly with n. Denser graphs (e.g., k-NN with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} k=9\end{document} ) are marginally slower. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} n=5{,}000\end{document} , SpaceBF completes in about 20 min on a Mac Pro (M3 Max). For substantially larger datasets, a practical alternative is to consider sdmTMB [162], which fits an NB SVC model via a Laplace-approximate maximum likelihood approach. It is extremely fast but can be less precise, may fail to converge, and often requires tuning the mesh density for interpretable results.

Run-time comparison of SpaceBF, with different priors: the horseshoe GMRF and ICAR, with varying spatial adjacency graphs.

For the convergence diagnostics, we computed the Geweke statistic [228] for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _1^{mm^{\prime }}(s_k)\end{document} , implemented in the R package coda [229], and investigated the trace plots of a few randomly chosen \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _1^{mm^{\prime }}(s_k)\end{document} ’s (see the Supplementary Material). When either the variable m or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} m^{\prime }\end{document} is highly sparse ( >75% zeroes), imposing additional normal priors on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _0^{mm^{\prime }}(s_k)\end{document} ’s and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _1^{mm^{\prime }}(s_k)\end{document} ’s with a moderate variance, such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} N(0, 10)\end{document} , drastically improves mixing and overall convergence performance.

Availability of source code and requirements

Project name: SpaceBFProject homepage: https://github.com/sealx017/SpaceBF/License: GPL-3.0Operating system: tested on macOS, WindowsProgramming language: RPackage management: GitHubHardware requirements: verified to run on laptops with 10 cores and 64 GB RAM

Additional files

Supplementary Figure S1: Cutaneous melanoma data analysis. (A) kNN adjacency network. (B) Clustering of spots based on centered and scaled estimates of slope surfaces of 72 statistically significant LR pairs. (C) The 4 main spatial patterns of the estimated surfaces. (D) Enrichment of LR interactions in 3 major cell types, with LR names arranged and color-coded according to their respective patterns. (E) The first 2 columns show the expression of 3 LR pairs. The third column displays the centered and scaled slope surfaces. In the fourth column, insignificant spot-level slope estimates are grayed.

Supplementary Figure S2: Cutaneous squamous cell carcinoma data analysis. (A) The 3 main spatial patterns of the estimated surfaces. (B) Bipartite graphs between Type I and Type 2 keratins based on their spatial pattern. (C) Study of coexpression between the Type 1 keratin KRT17 and 3 different Type 2 keratins. The insignificant spot-level slope estimates are grayed in the last column.

Supplementary Figure S3: Trace plot of 4 randomly chosen \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _{1}^{mm^{\prime }}(s_k)\end{document} ’s in the analysis of the LR pair: (IGF2, IGF1R) from the melanoma dataset.

Supplementary Figure S4: Trace plot of 4 randomly chosen \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \beta _{1}^{mm^{\prime }}(s_k)\end{document} ’s in the analysis of the LR pair: (SPP1, CD44) from the melanoma dataset.

Supplementary Material

giag006_Supplemental_File

giag006_Authors_Response_To_Reviewer_Comments_Original_Submission

giag006_GIGA-D-25-00259_Original_Submission

giag006_GIGA-D-25-00259_Revision_1

giag006_Reviewer_1_Report_Original_Submission Satwik Acharyya -- 7/27/2025

giag006_Reviewer_1_Report_Revision_1Satwik Acharyya -- 12/13/2025

giag006_Reviewer_2_Report_Original_SubmissionDaniel Domovic -- 8/22/2025

giag006_Reviewer_2_Report_Revision_1Daniel Domovic -- 1/3/2026