Spatially-Varying-Coefficient Regression with Spike-and-Slab Group Lasso

2026 Joint Statistical Meetings, Boston, MA

Dr. Cheng-Han Yu
Department of Mathematical and Statistical Sciences
Marquette University

2026-08-03

Why Spatially-Varying-Coefficients (SVC) Regression?

• A single global coefficient may mask or distort spatially heterogeneous relationships.

• SVC models let regression effects change with location: \beta_j(s), not a single \beta_j.

Air pollution (Hamm et al. 2015)

Forestry (Babcock et al. 2015)

Hydrology (Roksvåg, Steinsland, and Engeland 2022)

• Spatial variation is scientifically meaningful only when it represents a real departure from the global association.

Spatially varying coefficient models are not only a theoretical idea. Smooth coefficient surfaces already appear in real geostatistical applications. For example, Hamm and coauthors mapped spatially varying calibration slopes for PM10 across Europe. Babcock and coauthors showed smoothly varying LiDAR effects in forest biomass modeling. Roksvag and coauthors used Gaussian random field coefficient surfaces in runoff modeling. These examples motivate our focus on learning beta surfaces directly, not only fitting responses.

Not Every Association Is Spatially Varying

Statistical question: which predictors truly need a spatially varying coefficient surface?

• Not every predictor needs a full spatial coefficient surface.

\small \beta_1(s) = 1

\small \beta_2(s) = 1 + u_2(s)

• A nonzero coefficient surface \{\beta(s) \ne 0, s \in \mathcal{S} \} does NOT necessarily imply meaningful spatial heterogeneity.

The Blind Spot of Whole-Surface Selection

• Traditional (whole-surface) variable selection detects whether a predictor has any effect.

• Our proposed Global-Deviation Spike-and-Slab Group Lasso SVC (GD-SSGL-SVC) detects whether the effect varies over space.

  • \beta(s) = \theta + u(s): global average association + spatial deviation

The Blind Spot of Whole-Surface Selection

• Our proposed Global-Deviation Spike-and-Slab Group Lasso SVC (GD-SSGL-SVC) detects whether the effect varies over space.
  • \beta(s) = \theta + u(s): global average association + spatial deviation
• Whole-surface selection cannot distinguish global-only and spatially-varying.

Gap in Existing Methods

• Existing methods usually either suppress spatial heterogeneity completely, or select the entire coefficient surface.

• They do not directly quantify whether a predictor needs a spatial deviation beyond its global association.

Method	\theta explicit?	u(s) explicit?	Selects u(s)?
GD-SSGL	Yes	Yes	Yes
Global Bayesian lasso	Yes	No	No
Full (Bayesian) SVC	Possibly	Yes	No
Whole-surface selection	No separate	No separate	No
Generalized Additive Models (GAM)	Possibly	Yes	No

GD-SSGL-SVC Regression Model

y(s_i) = \sum_{j=1}^{p} x_{j}(s_i)\,\beta_j(s_i) + \varepsilon(s_i), \qquad \varepsilon(s_i) \stackrel{iid}{\sim} N(0,\sigma^2), \qquad i = 1, \dots, n

\boxed{\beta_j(s_i) = \theta_j + u_j(s_i)}

Global component: \theta_j

• Spatially constant

• Average association with y

• Bayesian lasso (BL) shrinkage

Spatial deviation component: u_j(s_i)

• Departure from the average effect

• Centered over observed locations, i.e., \sum_{i}u_j(s_i) = 0

• Spike-and-slab group lasso (SSGL) selection

over the study domain

Deviation from the Average

• Centered deviation u_j(s) oscillates around the global effect \theta_j.

Speed-Up: Low-Rank Basis Expansion

• Traditional Gaussian process SVC dimension grows like O(pn^3)

• Approximate each spatial deviation surface with H basis functions, where H \ll n with complexity O(p(nH^2 + H^3))

u_j(s) = \sum_{h=1}^H B^h(s)\alpha_{jh}, \quad j = 1, \dots, p

Dimension changes from pn location-level coefficients u_j(s_i) to pH grouped basis coefficients \alpha_{jh}

The previous slide identified the methodological gap. This slide presents the proposed solution. Instead of estimating one coefficient at every location, SSGL-SVC first represents each coefficient surface using a lower-dimensional basis expansion. The model then applies spike-and-slab group lasso shrinkage to the basis coefficients associated with each predictor. This means selection is made at the predictor-surface level. If a predictor has meaningful spatially varying effects, its surface is retained. If the predictor is globally weak, the entire surface is shrunk toward zero. This gives us a model that is spatially flexible, easier to interpret than an unregularized SVC model, and more scalable because selection happens after reducing the raw spatial dimension.

Centering the Spatial Basis

We center the spatial bases B^h(s) to make \theta_j and u_j(s) interpretable and identifiable

u_j(s) = \sum_{h=1}^H B_c^h(s)\alpha_{jh}, \qquad \dfrac{1}{n}\sum_{i=1}^nu_j(s_i) = 0, \qquad \dfrac{1}{n}\sum_{i=1}^n\beta_j(s_i) = \theta_j

• \theta_j becomes the average effect over the training domain.
• u_j(s) contains only spatial departures from that average.
• Constant variation cannot be arbitrarily exchanged between \theta_j and u_j(s).

After centering, each basis column has mean zero over the training locations; therefore, any fitted spatial deviation also has training-location mean zero. \theta_j is the average coefficient over the observed training locations, and u_j(s) is the centered spatial deviation from that average.

GD-SSGL-SVC: Prior Setting

• Spike-and-Slab Group Lasso on basis coefficients: If \boldsymbol \alpha_j \approx 0, u_j(s) \approx 0 for all s.

p(\boldsymbol \alpha_j\mid \boldsymbol{\gamma}) = \prod_{j=1}^{p}\Big[(1-\gamma_j)\,\Psi(\boldsymbol \alpha_j\mid \lambda_0) + \gamma_j\,\Psi(\boldsymbol \alpha_j\mid \lambda_1)\Big], \quad \lambda_0 > \lambda_1 > 0

\Psi(\boldsymbol \alpha_j\mid \lambda) \propto \exp\big(-\lambda \lVert \boldsymbol \alpha_j \rVert_2\big)

\eta_j \mid \tau_j,\sigma^2 \sim N(\bm{0},\sigma^2\tau_j\hbox{\bf I}_H), \qquad \tau_j\mid \gamma_j \sim Ga\!\left(\frac{H+1}{2},\frac{\lambda_{\gamma_j}^2}{2}\right),\qquad \lambda_{\gamma_j}= \begin{cases} \lambda_0,&\gamma_j=0,\\ \lambda_1,&\gamma_j=1. \end{cases}

• The posterior inclusion probability (PIP) P(\gamma_j=1\mid\mathbf y) summarizes how strongly do the data support a spatial deviation for predictor j, not any association with the response.

• Other priors \theta_j \sim BL(\omega_j, \lambda, a_{\lambda}, b_{\lambda}), \quad \sigma^2\sim IG(a_\sigma,b_\sigma), \quad \gamma_j\mid \pi \sim \mathrm{Bernoulli}(\pi), \quad \pi \sim \mathrm{Beta}(a,b)

• Treat each predictor-specific coefficient surface as one group, so selection is performed at the predictor level, not basis by basis.
• If \gamma_j = 0, the spike component with large \lambda_0 strongly shrinks the whole surface \beta_j(\cdot) toward zero.
• If \gamma_j = 1, the slab component with smaller \lambda_1 allows a non-negligible and spatially varying effect surface.
• The shared inclusion probability \theta provides multiplicity control and yields posterior inclusion probabilities as natural summaries of selection uncertainty.
• This prior combines smooth coefficient estimation, group sparsity, and explicit model uncertainty quantification in one Bayesian framework. The prior acts on whole predictor blocks. The spike component shrinks an entire coefficient surface toward zero, while the slab component keeps a predictor active and allows spatial variation. This is the core reason the method can do predictor-level selection and still estimate smooth coefficient maps.

Posterior Computation

Do not derive the full conditionals. The audience needs to know the computation is posterior sampling, and that the basis representation is what makes the grouped selection problem computationally manageable.

• Gibbs style posterior sampling with structured block updates
• Update coefficient blocks, shrinkage parameters, inclusion indicators, and variance terms
• Basis reduction keeps the computation manageable
• Designed for practical spatial data analysis I will not go into full conditional details here, but the main computational strategy is a Gibbs style sampler with block updates. The basis representation is important not only for modeling but also for computation, because it reduces the effective dimensionality of the spatial problem.

Simulation Design

Predictors	Structure	Scientific role
X_1,X_2	global-only: \theta \ne 0, \,u(s) = 0	Avoid inventing spatial variation?
X_3,X_4	spatial-only: \theta = 0, \,u(s) \ne 0	Detect pure spatial heterogeneity?
X_5,X_6	global+spatial: \theta \ne 0, \,u(s) \ne 0	Recover both components?
X_7,\ldots,X_{10}	null: \theta = 0, \,u(s) = 0	Suppress unnecessary effects?

Use this as a clean design table. Do not spend time on all simulation details. The design should signal that the comparisons target both estimation and selection.

Identify Spatial Heterogeneity

The proposed GD-SSGL correctly distinguished global-only effects from genuinely spatially varying effects.

• X_1,X_2：global only，PIP near 0
• X_3,\ldots,X_6：spatial deviations，PIP near 1
• X_7,\ldots,X_{10}：null，PIP near 0

• WS-SSGL: No separations on \beta, \beta_j(s) = \sum_{h=1}^H B^h(s)\alpha_{jh}
• Gaussian SVC: No SSGL on basis coefficient, \alpha_{jh} \sim \text{Normal}

GAM had the lowest point-estimation error overall. Among the Bayesian spatial methods, proposed SSGL performed best.

Simulation Results: Surface Recovery

All methods well recover the main spatial patterns of active predictors.

\theta = 1

\beta(s) = \theta + u(s) surface with PIP \approx 1

Selection alone is not enough in this problem. We also need good recovery of the spatial coefficient surfaces. This slide shows that the proposed method not only identifies relevant predictors, but also reconstructs their spatial patterns well, while shrinking irrelevant effects toward zero.

Simulation Results: Variable Selection

GD-SSGL-SVC separates signal predictors from noise predictors more sharply than other methods.

\beta(s) = 0 surface with PIP \approx 0

Here the main message is that the proposed method distinguishes signal predictors from noise predictors more effectively across the simulation settings we considered. I would emphasize the overall separation and stability rather than any single number on this slide.

Real Data Application - MODIS

• MODIS (Moderate Resolution Imaging Spectroradiometer) by NASA.

• Remote sensing instruments on NASA’s Terra and Aqua satellites.

Real Data Application - Data Summary

• Response: Enhanced Vegetation Index (EVI) measuring vegetation greenness (Huete et al. 2002) on Western US.

• 10 predictors:

  • Spectral Reflectance: Red, NIR, Blue, MIR
  • Ecosystem: GPP (amount of carbon through photosynthesis), LE (heat flux)
  • Satellite Observation Geometry: View/Sun zenith angle, Azimuth angle
  • Land Surface: LC Type4

Which predictors have spatially varying associations with EVI?

• Effects may differ across locations, so a global coefficient can be misleading
• The proposed method selects a parsimonious set of predictors and reveals their spatial patterns This real data example illustrates the practical value of the method. In this application, the scientific question is not only which predictors matter, but also whether their effects differ across space. The SSGL framework gives a parsimonious model and a more interpretable spatial story than a single coefficient model.

MODIS: Global and Variation Surface and Selection

• GD-SSGL-SVC leads to a larger coefficient size for the variables having effects on EVI

• Red Reflectance: Chlorophyll absorbs red light.

• \text{PIP} \approx 1 shows the effect of Red Reflectance on EVI is spatially dependent.

MODIS: Global and Variation Surface and Selection

• GD-SSGL-SVC leads to a larger coefficient size for the variables having effects on EVI

• Red Reflectance: Chlorophyll absorbs red light.

• \text{PIP} \approx 1 shows the effect of Red Reflectance on EVI is spatially dependent.

MODIS: Null Surface and Selection

• GD-SSGL-SVC shrinks coefficients more toward to 0 when the spatial deviation is negligible.

• Given all other variables, satellite observation geometry effect is negligible

  • \text{PIP} \approx 0, having negligible spatial variation
  • posterior mean of \theta \approx 0

Take Home Messages

The novel GD-SSGL-SVC

• Select deviations, not complete surfaces: The spatial-deviation PIP quantifies evidence that the association must vary over space beyond its global component.

• Suppress unnecessary spatial complexity: GD-SSGL avoids fitting spatial variation for global-only and null predictors.

• Strong Bayesian spatial performance: GD-SSGL improves prediction and coefficient-surface recovery over competing methods with narrower uncertainty intervals.

• Work in process and future work includes
  • spatiotemporal modeling: AR noises, \beta(s, t), \gamma_j(t), etc
  • non-Gaussian responses: binary, counts, zero-inflated counts, etc
  • irregular or constrained space: smooth within region, abrupt across boundaries, within-surface local sparsity selection
  • scalable online learning: variational, EM + Laplace, NNGP / Vecchia / SPDE / sketching

To conclude, our SSGL framework addresses an important gap in spatial modeling by combining flexible coefficient surfaces, groupwise selection, and Bayesian uncertainty quantification. The method works well in both simulation and real data settings, and we see it as a practical tool for modern spatial data analysis.