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,
  shared = TRUE,
  scale_coords = TRUE,
  parameterization = c("centered", "noncentered", "collapsed")
)

# S3 method for class 'ratiod_gp'
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.

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.

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.

x

A ratiod_gp object

...

Ignored

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 
#> Shared: Yes (enters both processes) 

# \donttest{
# Generate synthetic spatial data with a smooth field w shared by
# count and effort
set.seed(789)
n <- 40
lon <- runif(n, 0, 10)
lat <- runif(n, 0, 10)
w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2)
depth <- rnorm(n)
temp <- rnorm(n)
df <- data.frame(
  lon = lon,
  lat = lat,
  depth = depth,
  temp = temp,
  count = rpois(n, exp(3 + 0.3 * depth - 0.2 * temp + w)),
  effort = rgamma(n, shape = 4, rate = exp(-w))
)

# 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 = 400, chains = 1, seed = 1, adapt_delta = 0.95)
)
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.377 0.401  1.844  3.331 1.201        9
#>        num:depth  0.333 0.039  0.268  0.399 1.030       30
#>         num:temp -0.248 0.034 -0.305 -0.174 1.001       41
#> 
#> Fixed effects (denominator):
#>          parameter  mean    sd   2.5% 97.5%  rhat ess_bulk
#>  denom:(Intercept) 0.947 0.404  0.375 1.836 1.194        6
#>        denom:depth 0.040 0.109 -0.150 0.301 1.359       33
#>         denom:temp 0.011 0.072 -0.120 0.147 0.995       47
#> 
#> Variance components:
#>  parameter  mean    sd  2.5% 97.5%  rhat ess_bulk
#>      shape 3.727 0.788 2.600 5.615 0.995       31
#>  sigma2_gp 0.764 0.529 0.201 2.248 1.022       27
#> 
#> Other parameters:
#>  parameter   mean    sd   2.5%  97.5%  rhat ess_bulk
#>   gp_w[10]  0.648 0.427 -0.312  1.315 1.151        9
#>   gp_w[11]  1.224 0.412  0.268  1.793 1.154        5
#>   gp_w[12]  1.368 0.434  0.456  2.003 1.167        5
#>   gp_w[13]  0.794 0.401 -0.140  1.376 1.117        7
#>   gp_w[14]  1.052 0.446 -0.020  1.717 1.211        7
#>   gp_w[15]  0.919 0.428  0.012  1.572 1.104        7
#>   gp_w[16]  0.420 0.449 -0.515  1.052 1.195        4
#>   gp_w[17]  0.404 0.432 -0.623  1.111 1.149        5
#>   gp_w[18]  0.158 0.430 -0.821  0.806 1.110        7
#>   gp_w[19] -0.041 0.431 -0.835  0.675 1.213        4
#>    gp_w[1]  0.403 0.424 -0.490  1.101 1.194        4
#>   gp_w[20]  0.962 0.405  0.127  1.538 1.160        6
#>   gp_w[21]  1.055 0.421 -0.037  1.660 1.135        5
#>   gp_w[22]  0.651 0.411 -0.386  1.235 1.138        6
#>   gp_w[23]  1.443 0.411  0.564  2.098 1.164        5
#>   gp_w[24]  0.633 0.397 -0.346  1.285 1.187       10
#>   gp_w[25]  0.658 0.419 -0.315  1.323 1.134        8
#>   gp_w[26]  0.505 0.430 -0.477  1.171 1.137        6
#>   gp_w[27]  0.718 0.396 -0.183  1.273 1.127        7
#>   gp_w[28]  0.523 0.429 -0.466  1.271 1.122        8
#>   gp_w[29]  0.088 0.428 -0.804  0.737 1.228        4
#>    gp_w[2]  0.436 0.441 -0.608  1.166 1.205        4
#>   gp_w[30]  1.042 0.410  0.046  1.620 1.195        9
#>   gp_w[31]  1.604 0.406  0.715  2.270 1.172        8
#>   gp_w[32]  0.692 0.411 -0.175  1.293 1.166        5
#>   gp_w[33]  0.678 0.437 -0.229  1.363 1.142        6
#>   gp_w[34]  1.222 0.417  0.220  1.765 1.148        6
#>   gp_w[35]  0.778 0.412 -0.139  1.449 1.116        7
#>   gp_w[36]  1.007 0.425  0.060  1.670 1.154        9
#>   gp_w[37]  0.706 0.413 -0.139  1.305 1.160        5
#>   gp_w[38]  0.700 0.395 -0.180  1.336 1.119        8
#>   gp_w[39]  0.439 0.426 -0.464  1.133 1.168        5
#>    gp_w[3]  0.329 0.428 -0.540  1.021 1.149        6
#>   gp_w[40]  0.254 0.442 -0.649  1.040 1.143        6
#>    gp_w[4]  0.846 0.412 -0.101  1.424 1.187       10
#>    gp_w[5]  1.209 0.415  0.282  1.843 1.132        6
#>    gp_w[6]  1.113 0.447  0.157  1.758 1.170        5
#>    gp_w[7]  0.337 0.429 -0.566  1.040 1.114        7
#>    gp_w[8]  0.973 0.445  0.012  1.654 1.132        6
#>    gp_w[9]  1.738 0.428  0.821  2.335 1.162        6
#>     phi_gp  8.171 5.508  1.823 27.023 1.051       28
#> 
#> Diagnostics:
#>   Divergences: 0 
#>   Avg. acceptance: 0.968 
#>   Chains: 1 
#>   Iterations: 400 
#>   Warmup: 200 
#>   R-hat > 1.01 for 46 parameter(s), largest 1.359 (denom:depth).
#>   Bulk ESS < 100 for 49 parameter(s), smallest 4 (gp_w[19]).
#>   Run longer chains or reparameterise before using these estimates.

# 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 = 400, chains = 1, seed = 1, adapt_delta = 0.95)
)
# }