Space-Varying Regression Model

Definition

A Space-Varying Regression Model (SVRM) is a Bayesian regression model in which each spatial unit $i$ has its own coefficient vector $\beta_i$ , replacing the global $\beta$ of a standard linear model. Spatial coherence is enforced by a pairwise difference prior — a multivariate Markov random field (MRF) prior that penalises differences between the coefficient vectors of neighbouring units. The result is a fully Bayesian spatial model in which coefficient heterogeneity is data-driven and regularised by the neighbourhood graph.

Key Ideas

Global regression $\mu_i = x_i'\beta$ becomes unit-specific $\mu_i = x_i'\beta_i$ ; the full parameter $\beta = \text{vec}(B)$ where $B$ has $n$ rows $\beta_1',\ldots,\beta_n'$ .
The pairwise difference prior penalises deviations between neighbours but imposes no level constraint — it is improper (the graph Laplacian $W$ has all-zero row sums). Proper posteriors follow from proper hyperpriors on $\kappa$ and $\eta$ .
The model nests the standard fixed-effect regression (taking the prior very tight, $\eta \to \infty$ ) and the unpooled unit-by-unit OLS (taking the prior very diffuse, $\eta \to 0$ ).
The neighbourhood weight matrix and prior precision $\eta$ together determine how strongly coefficients are smoothed across space.

How It Works

Pairwise Difference Prior

Neighbours $i \sim j$ are defined by shared borders; $w_{ij} = 1$ if $i \sim j$ , 0 otherwise, and $w_{i+} = \sum_j w_{ij}$ is the degree. Define the graph Laplacian $W$ with $W_{ii} = w_{i+}$ and $W_{ij} = -w_{ij}$ . The prior is:

$f(\beta|\kappa,\eta) \propto \kappa^{nr/2}|\eta|^{n/2}\exp\!\left\{-\frac{\kappa}{2}\,\beta'(W \otimes \eta)\beta\right\}$

where $\eta$ is an $r \times r$ precision hyperparameter. The quadratic form equals $\sum_{i \sim j}(\beta_i - \beta_j)'\eta(\beta_i - \beta_j)$ . The observational precision $\kappa \sim \text{Gamma}(\nu_\kappa/2, \nu_\kappa S_\kappa/2)$ and $\eta \sim \text{Wishart}(\nu/2, \nu S/2)$ .

Full Conditional of $\beta_i$ (Scheme A — Single-Site Gibbs)

$(\beta_i|\beta_{\partial i}, \kappa, \eta, y) \sim N\!\bigl(a_i,\; \kappa^{-1}R_i\bigr)$ $R_i = (x_ix_i' + w_{i+}\eta)^{-1}, \qquad a_i = R_i\!\left(x_i y_i + w_{i+}\eta\bar\beta_{\partial i}\right)$

where $\bar\beta_{\partial i} = w_{i+}^{-1}\sum_{j\sim i}w_{ij}\beta_j$ is the weighted neighbour mean. Each draw requires only an $r\times r$ matrix inversion. Slow mixing because neighbouring $\beta_i$ 's are strongly correlated.

Block Update (Scheme B)

Jointly draw $(\beta, \kappa)$ from the Normal-Gamma marginal:

$(\beta, \kappa | \eta, y) \sim \text{NG}(a, R, \hat\nu_\kappa/2, \hat S_\kappa/2)$

where $R = (X'X + W \otimes \eta)^{-1}$ , $a = RX'y$ , $\hat\nu_\kappa = \nu_\kappa + n$ . The $nr \times nr$ precision matrix $X'X + W \otimes \eta$ is sparse; Rue (2001) sparse Cholesky makes this feasible for $n \sim 10^2$ – $10^3$ .

Analytical Marginalisation (Scheme C)

Integrate out $(\beta, \kappa)$ analytically to obtain the marginal $\pi(\eta|y)$ , then draw $\eta$ via Metropolis-Hastings or sampling-importance-resampling (SIR). This is the preferred scheme when $r$ is small: it achieves the lowest chain autocorrelation by eliminating the high-dimensional $\beta$ from the sampler entirely.

MCMC Efficiency Ranking

C $>$ B $>$ A in effective sample size per unit time. Scheme A is approximately 20× faster per iteration but has much lower acceptance per effective sample. The gap between schemes is smaller than in state-space models because the observation likelihood is highly informative about each $\beta_i$ — single-site updates are less hampered by prior correlation when data overwhelm prior.

Proper Prior Variant

Add a pinning term to break impropriety:

$f(\beta|\kappa,\eta,\delta) \propto \exp\!\left\{-\frac{\kappa}{2}\left[\sum_{i\sim j}(\beta_i-\beta_j)'\eta(\beta_i-\beta_j) + \delta\sum_i(\beta_i - b_i)'\eta(\beta_i-b_i)\right]\right\}$

The full conditional becomes $(\beta_i|\ldots) \sim N\!\bigl(w_i\bar\beta_{\partial i} + (1-w_i)b_i,\; (1/(\kappa(w_{i+}+\delta)))\eta^{-1}\bigr)$ where $w_i = w_{i+}/(w_{i+}+\delta)$ .

Why It Matters

Provides a single coherent Bayesian framework for spatially heterogeneous regression; the strength of spatial smoothing is estimated from data rather than fixed.
Extends the univariate Conditional Autoregressive Model (Besag et al. 1991) to multivariate coefficient vectors — relevant for multivariate outcomes and multi-covariate regressions.
Amazon land use application (Gamerman-Moreira-Rue 2003): spatially varying model achieves $D = 0.016$ vs. $D = 0.035$ for the fixed-coefficient model (Gelfand-Ghosh 1998 criterion), confirming meaningful coefficient heterogeneity; clear North-South gradient in cropland-to-pasture conversion effects.
The analytical-marginalisation scheme (C) generalises to non-Gaussian likelihoods as a Laplace approximation framework — a precursor to the INLA approach (Rue-Martino-Chopin 2009).

Open Questions

Optimal choice of neighbourhood matrix: first-order contiguity ignores distance decay and higher-order proximity; distance-based or economic-flow weights may be more appropriate.
Performance for large $n$ (thousands of spatial units): the $nr \times nr$ block in Scheme B becomes prohibitive; approximate methods (INLA, variational Bayes) may be needed.
Extension to non-Gaussian responses: Schemes A and B need Metropolis acceptance steps; Scheme C loses analytical tractability.