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.
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.
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.
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
df <- data.frame(
lon = runif(n, 0, 10),
lat = runif(n, 0, 10),
x = rnorm(n),
count = rpois(n, 25),
effort = rgamma(n, shape = 4, rate = 1)
)
# 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)
)
#> Inference: Exact (Tier 1)
#> Backend: hmc
#> Fitting ratio model...
#> Family: poisson_gamma
#> Observations: 100
#> Iterations: 500 (warmup: 250)
#> Running NUTS sampler...
#> Parameters: 43
#> Iterations: 500 (warmup: 250)
#> Chains: 2 (cores: 2)
#> Spatial: hsgp (100 units)
#> Warning: 46 divergent transition(s) after warmup. Increase max_treedepth or reparameterize.
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
#> beta_num[1] 3.188 0.018 3.152 3.224 1.005 445
#> beta_num[2] -0.006 0.023 -0.047 0.038 1.006 329
#>
#> Fixed effects (denominator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> beta_denom[1] 1.317 0.056 1.211 1.434 1.002 271
#> beta_denom[2] 0.031 0.066 -0.086 0.164 1.021 241
#>
#> Diagnostics:
#> Divergences: 46
#> Avg. acceptance: 0.871
#> Warning: Divergent transitions detected.
#> Consider: reparameterization or increasing adapt_delta.
#> Chains: 2
#> Iterations: 500
#> Warmup: 250
# }