Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Choosing Between Spatial Effect Models

Spatial Transcriptomics · Tutorial 640 of 700

Choosing Between Spatial Effect Models

Learn how domain geometry, neighborhood structure, computation, and boundary behavior guide the choice between Gaussian-process and graph-based spatial effects.

Advanced 11 min read

What You'll Learn

  • Distinguish coordinate-based Gaussian-process covariance from graph-based conditional dependence.
  • Compare the assumptions made by stationary distance kernels and conditional autoregressive models.
  • Calculate covariance and conditional expectations for a shared hypothetical tissue pattern.
  • Estimate the memory and factorization costs of dense Gaussian-process and sparse graph models.
  • Assess what each model implies at tissue edges, holes, and disconnected regions.
  • Recognize diagnostics and limitations that inform model choice.

Why the Spatial Domain Matters

A spatial random effect can account for residual similarity between nearby tissue locations, but “nearby” depends on how the tissue is represented. A Gaussian process (GP) defines dependence through coordinates and a covariance function. A graph-based effect defines local dependence through specified edges between spots. These choices can produce different spatial patterns, even when both models are applied to the same expression data.

This tutorial compares the two approaches using one hypothetical tissue section: a roughly grid-shaped set of spots with a localized expression hotspot and a narrow tear that separates two regions. A small neighborhood around the hotspot will make the calculations transparent; a larger version of the same section will illustrate computation and boundary behavior. The spatial effect is latent, so the pattern below is a hypothetical simulated effect, not a real study result.

Two Models for the Same Latent Effect

Let \(Y_i\) denote the count for a gene at spot \(i\), let \(x_i\) contain measured covariates, and let \(u_i\) be a latent spatial effect. A typical count-regression formulation is

$$ g\{\mathbb{E}(Y_i\mid u_i)\} = \log(L_i)+x_i^\mathsf{T}\beta+u_i, $$

where \(L_i\) is a library-size offset and \(g\) is a chosen link, often the log link for a count mean. The choice between a GP and a graph model concerns the joint distribution of \(u=(u_1,\ldots,u_n)^\mathsf{T}\), not the count likelihood itself. In either case, the spatial term does not by itself establish that location causes expression changes.

Definition: A Gaussian-process spatial effect assigns the latent effects a multivariate normal distribution whose covariance is determined by spot coordinates: \(u\sim N(0,K)\), with \(K_{ij}=k(s_i,s_j)\). A graph-based spatial effect specifies a graph of spots and a Gaussian distribution through a sparse precision matrix \(Q\), where \(Q\) is the inverse covariance matrix when the model is proper. The GP is built from a spatial domain and covariance kernel; the graph model is built from a domain and its declared neighborhood relations.

For a stationary, isotropic exponential GP, a common kernel is

$$ K_{ij} = \sigma_s^2\exp\left(-\frac{\|s_i-s_j\|}{\ell}\right), $$

where \(s_i\) is the coordinate of spot \(i\), \(\sigma_s^2\) is the marginal variance, and \(\ell>0\) is a distance scale. “Stationary” here means that covariance depends on the separation between locations rather than their absolute position; “isotropic” means it depends on distance, not direction. This kernel makes covariance decrease smoothly with Euclidean distance. It does not automatically encode tissue barriers, anisotropy, or changes in the spatial process across regions.

For a graph model, let \(W\) be a symmetric adjacency matrix with \(w_{ij}=1\) when spots \(i\) and \(j\) are connected and \(0\) otherwise, and let \(D\) be diagonal with \(d_i=\sum_jw_{ij}\). One proper conditional autoregressive (CAR) specification is

$$ u\sim N(0,Q^{-1}), \qquad Q=\tau(D-\rho W), $$

where \(\tau>0\) controls precision and \(0\leq\rho<1\) controls dependence. For a connected, undirected graph with no isolated vertices, this range gives a positive-definite precision matrix. More generally, propriety requires every eigenvalue of \(D^{-1/2}WD^{-1/2}\), multiplied by \(\rho\), to be less than 1. Isolated vertices need a separately specified effect rather than this CAR conditional.

The graph model has the conditional distribution

$$ u_i\mid u_{-i} \sim N\left( \frac{\rho}{d_i}\sum_jw_{ij}u_j,\; \frac{1}{\tau d_i} \right). $$

Thus each spot’s conditional mean is a scaled average of its graph neighbors. Non-neighbors are conditionally independent given the remaining effects, though they need not be marginally independent if paths connect them through the graph.

What the Models Assume

QuestionGaussian processGraph-based CAR effect
What defines proximity?Coordinates and a covariance kernel, commonly Euclidean distance.Edges in \(W\), chosen from a neighborhood graph or tissue structure.
Dependence assumptionCovariance follows the selected kernel; a stationary isotropic kernel treats equal distances alike.Conditional dependence is local to graph neighbors; the precision matrix is sparse.
Domain flexibilityNatural for continuous coordinates and irregularly spaced samples, subject to the kernel assumptions.Natural when adjacency or barriers are well represented by a graph, including non-grid layouts.
Typical computationDense covariance operations can require quadratic memory and cubic factorization work.Sparse precision operations can be much cheaper, but factorization cost depends on graph structure and fill-in.
Boundary behaviorNearby points remain correlated across a boundary unless the kernel or domain is modified.Dependence follows the declared edges; missing cross-boundary edges prevent direct conditional dependence.

Worked Example: A Shared Hotspot Neighborhood

Worked Example: A Shared Hotspot Neighborhood

Consider a hypothetical \(3\times3\) square-grid neighborhood from the tissue section. Its simulated latent effects are

$$ \begin{matrix} 0 & 0.2 & 0\\ 0.2 & 1.0 & 0.2\\ 0 & 0.2 & 0 \end{matrix} $$

The central spot has effect \(1.0\), and its four edge-sharing neighbors each have effect \(0.2\). Suppose the GP uses the exponential kernel with \(\ell=1\) grid unit. The GP prior correlation between the center and each edge-sharing neighbor is

$$ \operatorname{Corr}(u_{\mathrm{center}},u_{\mathrm{edge}}) = \exp(-1/1) = 0.3679. $$

A corner spot is \(\sqrt{2}\) grid units from the center, so its prior correlation with the center is

$$ \exp(-\sqrt{2}/1)=0.2431. $$

These are prior correlations implied by the coordinates and kernel, not correlations calculated from this nine-value pattern. The GP distinguishes distance continuously: the corner is less correlated with the center than an edge-sharing spot.

Now use the four-edge neighborhood graph and a proper CAR model with \(\rho=0.8\) and \(\tau=2\). The center has degree \(d_i=4\). Conditional on its four neighbors, its CAR mean and variance are

$$ \mathbb{E}(u_{\mathrm{center}}\mid u_{-\mathrm{center}}) = \frac{0.8}{4}(0.2+0.2+0.2+0.2) = 0.16, \qquad \operatorname{Var}(u_{\mathrm{center}}\mid u_{-\mathrm{center}}) = \frac{1}{2\cdot4} = 0.125. $$

The simulated center value \(1.0\) is far above that conditional mean. This does not make the model invalid: the likelihood can support a high center effect. It does show that a sharply localized peak is not the same as a smooth field, and that the assumed spatial prior affects how surprising the peak is. In this graph, diagonal corners are not direct neighbors of the center; in the GP they still have nonzero covariance with it.

Why the CAR Conditional Formula Holds

The CAR conditional mean and variance follow directly from the precision matrix. For a multivariate normal distribution with precision \(Q\), the conditional distribution of component \(u_i\), given all other components, has variance \(1/Q_{ii}\) and mean \(-\sum_{j\ne i}Q_{ij}u_j/Q_{ii}\). Here \(Q_{ii}=\tau d_i\) and \(Q_{ij}=-\tau\rho w_{ij}\) for \(j\ne i\). Substitution gives

$$ \operatorname{Var}(u_i\mid u_{-i}) = \frac{1}{\tau d_i}, \qquad \mathbb{E}(u_i\mid u_{-i}) = -\frac{\sum_{j\ne i}(-\tau\rho w_{ij})u_j}{\tau d_i} = \frac{\rho}{d_i}\sum_jw_{ij}u_j. $$

This also explains the graph model’s local conditional structure: if \(w_{ij}=0\), the corresponding off-diagonal precision entry is zero, so \(u_i\) and \(u_j\) are conditionally independent given the other effects.

Worked Example: Computational Scale

Worked Example: Computational Scale

Suppose the full hypothetical section contains \(n=5{,}000\) spots. A dense GP covariance matrix has \(n^2=25{,}000{,}000\) entries. Storing each entry as an 8-byte number requires

$$ 25{,}000{,}000\times8 = 200{,}000{,}000\text{ bytes} \approx200\text{ MB}. $$

This is only the storage for one matrix, before temporary arrays or other model quantities. A dense Cholesky factorization has leading work of approximately \(n^3/3\), which here is

$$ \frac{5{,}000^3}{3} \approx 4.17\times10^{10}\text{ floating-point operations}. $$

For comparison, place the same number of spots on a \(100\times50\) grid and connect horizontal and vertical neighbors. The number of edges is \(99\cdot50+100\cdot49=9{,}850\). The symmetric CAR precision has \(5{,}000+2(9{,}850)=24{,}700\) nonzero entries, including the diagonal. At 8 bytes per value, those values alone use \(197{,}600\) bytes, about \(0.198\) MB. Sparse index storage and factorization workspace add to this figure.

The sparse representation is a substantial advantage, not a guarantee of linear-time fitting. Sparse Cholesky factorization can create fill-in, and its cost depends on graph connectivity and ordering. Conversely, GP approximations, low-rank representations, and structured kernels can reduce dense computation, but they change or approximate the original covariance model.

Worked Example: A Tear and a Tissue Boundary

Worked Example: A Tear and a Tissue Boundary

In the same hypothetical section, suppose the nearest spots on opposite sides of a narrow tear are \(0.2\) coordinate units apart. With the GP kernel above, their prior correlation is

$$ \exp(-0.2/1)=0.8187. $$

A standard Euclidean GP therefore links the two sides strongly despite the tear. This may be appropriate if the tear is only a missing observation in an otherwise continuous process. It may be inappropriate if the tear separates distinct tissue compartments or blocks spatial continuity.

Now suppose the graph is intentionally constructed with no edges across the tear and with the two sides forming separate connected components. Under the proper CAR model, \(Q\) is block diagonal across those components, so their effects are independent under the prior and have zero cross-component covariance. The zero dependence is a consequence of the graph construction, not a universal property of CAR models. If the spots remain connected by another path, their effects can still be marginally dependent even without a direct edge.

This comparison is about a modeling decision: the GP uses geometric distance unless its covariance is adapted, while the graph uses declared connectivity. Neither representation discovers the correct tissue boundary automatically.

Diagnostics, Failure Modes, and Choice

A model can fit expression counts while describing spatial dependence poorly. Inspect the fitted spatial effects on the tissue coordinates and compare their pattern with the domain and the scientific question. Look for abrupt jumps where the model assumes smoothness, smoothing across a tear or compartment boundary, isolated spots with extreme effects, and edge behavior unlike the interior. Residual spatial structure can also indicate that the chosen spatial effect has not captured relevant dependence, although residual patterns can have other causes, including omitted covariates or unmodeled composition.

A GP with an isotropic stationary kernel can be a good choice when continuous coordinates and smooth distance-decaying dependence are plausible. It is less suitable without modification when barriers, directional dependence, or abrupt regional changes are central. A graph CAR effect can be a good choice when a scientifically defensible adjacency structure is available and local conditional smoothing is appropriate. It can be misleading if the graph connects biologically separated spots, omits meaningful links, or leaves isolated nodes untreated.

Boundary behavior deserves explicit attention. A standard GP does not know that a tissue edge is an edge: it assigns covariance to any pair of spots according to the kernel and their distance. This can yield sensible interpolation near an edge, but can also imply unwanted dependence across holes or narrow gaps. A graph model has no direct dependence across a missing edge, but its behavior at a boundary depends on how the graph and degrees are defined. Fewer neighbors change the CAR conditional variance \(1/(\tau d_i)\), so edge spots can have different conditional uncertainty from interior spots under this parameterization.

Neither model is universally superior. If the scientific domain is naturally continuous and covariance by distance is credible, the GP provides a direct, flexible covariance description. If a sparse neighborhood structure or barrier-aware connectivity is central, a graph model can encode that structure efficiently. For large datasets, computational feasibility may favor sparse graph methods, while structured or approximate GP methods may be appropriate when their approximation assumptions are acceptable. The choice should be documented in terms of domain, covariance or graph construction, computational constraints, and the spatial patterns each model can represent.

Common Mistake

Wrong reasoning: “The two spots are physically close, so a GP and a graph model must assign them similar dependence.” This overlooks what each model uses. A distance-kernel GP assigns covariance based on its kernel and coordinates; a graph model uses the specified edges, and a pair without a direct edge can still be marginally associated through graph paths. Correction: inspect the actual kernel and the actual graph, including holes, boundaries, and disconnected components, before interpreting the fitted spatial effect.

Key takeaway: Choose the spatial effect that represents the relevant domain. A GP encodes covariance through coordinates and a kernel; a graph CAR model encodes local conditional dependence through adjacency. Their computational costs and boundary implications follow from those constructions.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

  1. Under an exponential GP kernel with \(\ell=2\), what is the prior correlation between spots 2 coordinate units apart?
  2. A CAR spot has degree 5, \(\rho=0.6\), and neighbor effects \(0.1, 0.2, 0.2, 0.3,\) and \(0.2\). What is its conditional mean?
  3. Why does a standard Euclidean GP potentially link spots across a narrow tissue tear?
  4. Does the absence of a direct edge between two spots imply marginal independence in a connected CAR graph? Explain.
  5. Give one reason a sparse CAR model may still require substantial computation.

Solutions. 1. The correlation is \(\exp(-2/2)=\exp(-1)=0.3679\). 2. The neighbor effects sum to \(1.0\), so the conditional mean is \((0.6/5)(1.0)=0.12\). 3. The kernel uses Euclidean distance and does not encode the tear unless it is modified, so nearby coordinates can retain high covariance. 4. No. A path through other nodes can induce marginal dependence, even though the pair is conditionally independent given the remaining effects when there is no direct edge. 5. Sparse factorization can produce fill-in, and its cost depends on graph structure and ordering.

Clinical Biostats