Skip to contents

Hilbert Space Gaussian Process (HSGP) Spatial Structure

Source: R/spatial.R
spatial_hsgp.Rd

Specify a Hilbert Space Gaussian Process approximation for spatial effects. HSGP approximates a GP using Laplacian eigenfunctions, providing O(NM^2) complexity with analytical gradients instead of O(Nk^2) with numerical gradients for NNGP.

This gives approximately 50x speedup over standard NNGP while maintaining high accuracy for smooth spatial fields.

Usage

spatial_hsgp(coords, m = 6, c = 1.5, shared = TRUE, scale_coords = TRUE)

# S3 method for class 'ratiod_hsgp'
print(x, ...)

Arguments

coords

A one-sided formula specifying coordinate columns (e.g., ~ lon + lat), or a character vector of length 2 with column names.

m

Number of basis functions per dimension. Total basis functions will be m^2. Default 6. Higher values give better approximation but slower computation. Recommended range: 5-15.

c

Boundary factor controlling domain extension beyond data range. Default 1.5. The domain is extended to \([-cL, cL]\) where \(L\) is half the data range. Larger values improve approximation at boundaries.

shared

Logical; if TRUE (default), spatial effect enters both numerator and denominator. Set to FALSE for process-specific spatial effects (triggers warning about potential confounding).

scale_coords

Logical; if TRUE (default), coordinates are scaled to unit variance before computing basis functions.

x

A ratiod_hsgp object

...

Ignored

Value

A ratiod_hsgp object

Details

HSGP approximates a GP as:

$$f(x) = \sum_{j=1}^{M^2} \phi_j(x) \sqrt{S(\lambda_j)} \beta_j$$

where:

  • \(\phi_j(x)\) are Laplacian eigenfunctions (products of sines)

  • \(S(\lambda_j)\) is the spectral density of the squared exponential kernel

  • \(\beta_j \sim N(0, 1)\) are basis coefficients

The hyperparameters are:

  • \(\sigma^2\): marginal variance (PC prior: P(sigma > 1) = 0.01)

  • \(\ell\): lengthscale (LogNormal(0, 1) prior)

Advantages over NNGP:

  • Analytical gradients enable ~50x speedup

  • Simple parameter interpretation

  • Works well for smooth spatial fields

Limitations:

  • Assumes squared exponential (smooth) covariance

  • Less accurate for rough fields (Matern with low nu)

  • Approximation quality depends on m and c choices

References

Riutort-Mayol, G., Bürkner, P. C., Andersen, M. R., Solin, A., & Vehtari, A. (2023). Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. Statistics and Computing, 33(1), 17.

See also

spatial_gp() for NNGP-based GP, spatial_car() for areal effects

Examples

# Create HSGP spatial structure
hsgp <- spatial_hsgp(~ lon + lat)
print(hsgp)
#> tulpaRatio Hilbert Space GP (HSGP) spatial specification
#> =====================================================
#> 
#> Coordinates: lon, lat 
#> Basis functions: 6 per dim ( 36 total )
#> Boundary factor: 1.5 
#> Shared: Yes (enters both processes) 

# \donttest{
# Generate synthetic spatial data
set.seed(42)
n <- 100
lon <- runif(n, 0, 10)
lat <- runif(n, 0, 10)
w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2)
x <- rnorm(n)
df <- data.frame(
  lon = lon,
  lat = lat,
  x = x,
  count = rpois(n, exp(3 + 0.3 * x + w)),
  effort = rgamma(n, shape = 4, rate = exp(-w))
)

# Fast spatial effect with HSGP (much faster than spatial_gp)
fit <- tratio(
  count | effort ~ x,
  data = df,
  family = ratiod_poisson_gamma(),
  spatial = spatial_hsgp(~ lon + lat, m = 6),
  mode = "hmc",
  control = list(iter = 500, warmup = 250, chains = 2, seed = 1)
)
summary(fit)
#> ratio model summary
#> ===================
#> 
#> Inference: Exact (Tier 1) via hmc
#> Family: poisson_gamma 
#> 
#> Fixed effects (numerator):
#>        parameter  mean    sd  2.5% 97.5%  rhat ess_bulk
#>  num:(Intercept) 2.900 0.193 2.453 3.239 1.072       36
#>            num:x 0.286 0.023 0.244 0.336 1.005      342
#> 
#> Fixed effects (denominator):
#>          parameter   mean    sd   2.5% 97.5%  rhat ess_bulk
#>  denom:(Intercept)  1.227 0.199  0.737 1.606 1.054       40
#>            denom:x -0.032 0.065 -0.167 0.085 1.020      147
#> 
#> Variance components:
#>  parameter  mean    sd  2.5% 97.5%  rhat ess_bulk
#>      shape 3.691 0.544 2.713 4.854 1.007      260
#> 
#> Diagnostics:
#>   Divergences: 0 
#>   Avg. acceptance: 0.919 
#>   Chains: 2 
#>   Iterations: 500 
#>   Warmup: 250 
#>   R-hat > 1.01 for 3 parameter(s), largest 1.072 (num:(Intercept)).
#>   Bulk ESS < 200 for 3 parameter(s), smallest 36 (num:(Intercept)).
#>   Run longer chains or reparameterise before using these estimates.
# }