Gaussian Process spatial structure

Specify a Gaussian Process (GP) spatial random effect using coordinate-based distances. Uses Nearest Neighbor Gaussian Process (NNGP) approximation for computational efficiency with large datasets (scales to millions of observations).

This provides a continuous spatial effect that captures smooth spatial variation, unlike CAR/BYM2 which require discrete areal units.

Usage

spatial_gp(
  coords,
  cov = c("exponential", "matern", "gaussian", "spherical"),
  nu = 1.5,
  nn = 15,
  solver = c("auto", "cholesky", "cg", "pcg", "gpu"),
  cg_tol = 1e-06,
  cg_maxiter = 100,
  shared = TRUE,
  scale_coords = TRUE,
  parameterization = c("centered", "noncentered", "collapsed")
)

Arguments

coords

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

cov

Covariance function: "exponential" (default), "matern", "gaussian", or "spherical".

nu

Smoothness parameter for Matern covariance. Common values:

0.5: Equivalent to exponential (rough, once mean-square differentiable)
1.5: Once differentiable (moderate smoothness)
2.5: Twice differentiable (smooth) Ignored for non-Matern covariance functions.

nn

Number of nearest neighbors for NNGP approximation. Default 15. Larger values give better approximation but slower computation.

solver

Linear algebra solver for GP computations:

"auto" (default): Automatically select based on problem size. Uses Cholesky for N < 2000, CG for larger problems. Uses GPU if available and N > 5000.
"cholesky": Direct Cholesky decomposition. Exact but O(N*k^3). Best for smaller datasets or when high precision is critical.
"cg": Conjugate Gradient iterative solver. O(Nk^2iter). Better for large N (> 2000) with good preconditioning.
"pcg": Preconditioned CG with diagonal preconditioner. Faster convergence than CG for ill-conditioned systems.
"gpu": GPU-accelerated batched Cholesky using CUDA (if available). Requires tulpaRatio to be compiled with GPU support. Falls back to PCG if GPU is unavailable. Best for large datasets (N > 5000) with CUDA-capable GPU. Use gpu_available() to check support.

cg_tol

Convergence tolerance for CG/PCG solvers. Default 1e-6. Smaller values give more accurate solutions but slower convergence.

cg_maxiter

Maximum iterations for CG/PCG. Default 100.

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 distances.

parameterization

Parameterization of the latent field. One of "centered" (default), "noncentered" (Matt trick, better for diffuse prior regimes), or "collapsed" (marginalize spatial variance for improved mixing). Centered is robust for informative priors.

Value

A ratiod_gp object

Details

The GP spatial model adds a spatially-correlated random effect to the linear predictor:

$$\eta(s) = X\beta + w(s)$$

where $w(s)$ follows a Gaussian process: $$w(s) \sim GP(0, \sigma^2 C(d; \phi))$$

The correlation function $C(d; \phi)$ depends on distance $d$ and range parameter $\phi$:

Exponential: $C(d) = \exp(-d/\phi)$
Matern: $C(d) = \frac{2^{1-\nu}}{\Gamma(\nu)} (\sqrt{2\nu} d/\phi)^\nu K_\nu(\sqrt{2\nu} d/\phi)$
Gaussian: $C(d) = \exp(-(d/\phi)^2)$
Spherical: $C(d) = 1 - 1.5(d/\phi) + 0.5(d/\phi)^3$ for $d < \phi$

NNGP approximation: For computational tractability, we use the Nearest Neighbor Gaussian Process (Datta et al., 2016), which conditions each location on its k nearest neighbors. This reduces complexity from O(n^3) to O(n*k^2), enabling models with millions of observations.

References

Datta, A., Banerjee, S., Finley, A. O., & Gelfand, A. E. (2016). Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 111(514), 800-812.

Examples

# Create GP spatial structure
gp <- spatial_gp(~ lon + lat)
print(gp)
#> tulpaRatio Gaussian Process spatial specification
#> =============================================
#> 
#> Coordinates: lon, lat 
#> Covariance: exponential 
#> Neighbors (NNGP): 15 
#> Solver: auto (Cholesky<2k, PCG<5k, GPU/CG for larger) 
#> Shared: Yes (enters both processes) 

# \donttest{
# Generate synthetic spatial data
set.seed(789)
n <- 50
df <- data.frame(
  lon = runif(n, 0, 10),
  lat = runif(n, 0, 10),
  depth = rnorm(n),
  temp = rnorm(n),
  count = rpois(n, 25),
  effort = rgamma(n, shape = 4, rate = 1)
)

# Continuous spatial effect with exponential covariance
fit <- tratio(
  count | effort ~ depth + temp,
  data = df,
  family = ratiod_poisson_gamma(),
  spatial = spatial_gp(~ lon + lat),
  mode = "hmc",
  control = list(iter = 200, warmup = 100, chains = 1)
)
#> Inference: Exact (Tier 1)
#>   Backend: hmc
#> Fitting ratio model...
#>   Family: poisson_gamma
#>   Observations: 50
#>   Iterations: 200 (warmup: 100)
#>   Auto-capping max_treedepth to 8 for GP spatial model (use max_treedepth = 10 to override)
#> Running NUTS sampler...
#>   Parameters: 59
#>   Iterations: 200 (warmup: 100)
#>   Chains: 1 (cores: 1)
#>   Spatial: gp (50 units)
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]  71.629 0.003  71.624  71.636 1.690        2
#>  beta_num[2]  28.246 0.001  28.245  28.247 1.125        8
#>  beta_num[3] -44.863 0.003 -44.869 -44.857 1.468        2
#> 
#> Fixed effects (denominator):
#>      parameter    mean    sd    2.5%   97.5%  rhat ess_bulk
#>  beta_denom[1] -45.804 0.005 -45.813 -45.798 2.090        1
#>  beta_denom[2]  44.232 0.007  44.222  44.246 1.664        2
#>  beta_denom[3] -13.883 0.001 -13.885 -13.881 1.580        2
#> 
#> Diagnostics:
#>   Divergences: 0 
#>   Avg. acceptance: 1 
#>   Chains: 1 
#>   Iterations: 200 
#>   Warmup: 100 

# Smoother spatial field with Matern covariance
fit2 <- tratio(
  count | effort ~ depth + temp,
  data = df,
  family = ratiod_poisson_gamma(),
  spatial = spatial_gp(~ lon + lat, cov = "matern", nu = 1.5),
  mode = "hmc",
  control = list(iter = 200, warmup = 100, chains = 1)
)
#> Inference: Exact (Tier 1)
#>   Backend: hmc
#> Fitting ratio model...
#>   Family: poisson_gamma
#>   Observations: 50
#>   Iterations: 200 (warmup: 100)
#>   Auto-capping max_treedepth to 8 for GP spatial model (use max_treedepth = 10 to override)
#> Running NUTS sampler...
#>   Parameters: 59
#>   Iterations: 200 (warmup: 100)
#>   Chains: 1 (cores: 1)
#>   Spatial: gp (50 units)
#> Warning: 8 divergent transition(s) after warmup. Increase max_treedepth or reparameterize.
# }