Specify a multi-scale spatial random effect that decomposes spatial variation into local (fine-scale) and regional (broad-scale) components. Each scale has its own range and variance parameters, and both are expanded on one Hilbert-space basis (Solin and Sarkka 2020).
Arguments
- coords
A one-sided formula specifying coordinate columns (e.g.,
~ lon + lat), or a character vector of length 2 with column names.- scales
Character vector specifying scale names. Default:
c("local", "regional").- range_local
Prior range for local scale as
c(lower, upper)in coordinate units. Default:c(0.01, 1)(after scaling).- range_regional
Prior range for regional scale as
c(lower, upper). Default:c(1, 10)(after scaling).- m
Basis functions per dimension, giving
m^2coefficients per scale. Default 10.- c_boundary
Boundary factor: the domain is extended by this multiple before the basis is built. Default 1.5.
Logical; if TRUE (default), spatial effects enter both numerator and denominator.
- scale_coords
Logical; if TRUE (default), coordinates are scaled to unit variance before computing distances.
- sampler
Warmup mass-matrix adaptation for the two GP blocks:
"auto"(default) leaves the sampler's own adaptation in place,"lbfgs"runs an L-BFGS pass during warmup to initialize it.- x
A ratiod_multiscale object
- ...
Ignored
Details
The multi-scale model decomposes spatial variation additively:
$$\eta(s) = X\beta + w_{local}(s) + w_{regional}(s)$$
where each component follows an independent Gaussian process with a squared-exponential covariance: $$w_{local}(s) \sim GP(0, \sigma^2_{local} C(\ell_{local}))$$ $$w_{regional}(s) \sim GP(0, \sigma^2_{regional} C(\ell_{regional}))$$
The range priors centre each scale's lengthscale prior.
Each field is the expansion of m^2 basis functions of the domain with
standard normal coefficients scaled by the covariance's spectral density,
so its cost does not grow with the number of locations.
Identifiability: the data inform the sum of the two fields, and the priors on the variances and ranges decide how it is split between them. Non-overlapping range priors give each scale its own band of the spectrum. A field with one value per location for each scale (a nearest-neighbour GP per scale) leaves that split along a ridge in as many dimensions as there are locations, which NUTS does not traverse, so the basis is the only representation offered.
References
Solin, A. and Sarkka, S. (2020). Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing, 30(2), 419-446.
See also
spatial_gp() for single-scale GP, temporal_multiscale() for
multi-scale temporal effects
Examples
# Create multi-scale spatial structure
ms <- spatial_multiscale(
~ lon + lat,
range_local = c(0.1, 0.5),
range_regional = c(1, 5)
)
print(ms)
#> tulpaRatio Multi-Scale spatial specification
#> ========================================
#>
#> Coordinates: lon, lat
#> Scales: local + regional
#>
#> Local scale:
#> Range prior: [0.1, 0.5]
#>
#> Regional scale:
#> Range prior: [1, 5]
#>
#> Approximation: HSGP, 10 basis functions per dim ( 100 per scale ), boundary factor 1.5
#> Covariance: squared exponential
#> Shared: Yes (enters both processes)
# \donttest{
# A regional trend plus local variation, shared by count and effort
set.seed(101)
n <- 60
lon <- runif(n, 0, 10)
lat <- runif(n, 0, 10)
w <- 0.5 * sin(lon / 3) + 0.5 * cos(lat / 3) + 0.3 * sin(2 * lon) * cos(2 * lat)
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 = 5, rate = exp(-w))
)
# Multi-scale spatial with local and regional components
fit <- tratio(
count | effort ~ depth + temp,
data = df,
family = ratiod_poisson_gamma(),
spatial = spatial_multiscale(
~ lon + lat,
range_local = c(0.1, 0.5),
range_regional = c(1, 5),
m = 6
),
control = list(iter = 400, chains = 1, 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) 3.260 0.183 2.999 3.750 1.034 19
#> num:depth 0.362 0.027 0.314 0.417 1.008 121
#> num:temp -0.139 0.031 -0.199 -0.068 1.012 84
#>
#> Fixed effects (denominator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> denom:(Intercept) 1.884 0.185 1.591 2.443 1.043 20
#> denom:depth 0.070 0.077 -0.094 0.240 0.999 84
#> denom:temp 0.130 0.055 0.035 0.235 0.997 105
#>
#> Variance components:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> shape 4.600 0.844 3.047 6.335 1.005 67
#> sigma2_hsgp_local 0.016 0.025 0.000 0.075 1.034 21
#> sigma2_hsgp_regional 0.625 0.452 0.179 1.856 1.026 70
#>
#> Other parameters:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> hsgp_beta_local[10] -0.081 0.990 -1.750 2.014 1.100 70
#> hsgp_beta_local[11] -0.219 0.801 -1.738 1.280 1.008 116
#> hsgp_beta_local[12] 0.210 0.907 -1.494 1.914 0.996 68
#> hsgp_beta_local[13] -0.385 0.829 -2.007 1.089 0.997 169
#> hsgp_beta_local[14] 0.426 0.986 -1.542 2.637 1.022 105
#> hsgp_beta_local[15] -0.295 0.911 -1.844 1.234 0.999 124
#> hsgp_beta_local[16] -0.100 1.073 -1.686 2.213 1.035 32
#> hsgp_beta_local[17] 0.056 1.009 -2.120 1.996 1.023 76
#> hsgp_beta_local[18] -0.332 0.956 -1.795 1.629 1.004 105
#> hsgp_beta_local[19] 0.252 0.957 -1.906 2.132 1.007 89
#> hsgp_beta_local[1] 0.019 1.075 -1.894 2.174 1.010 124
#> hsgp_beta_local[20] 0.233 0.795 -1.398 1.866 1.079 38
#> hsgp_beta_local[21] 0.022 0.950 -1.680 1.583 1.029 91
#> hsgp_beta_local[22] 0.013 0.988 -1.806 1.825 1.015 105
#> hsgp_beta_local[23] -0.073 1.058 -2.187 1.850 1.007 87
#> hsgp_beta_local[24] 0.024 0.874 -1.952 1.461 1.002 100
#> hsgp_beta_local[25] -0.165 0.856 -1.915 1.387 1.012 56
#> hsgp_beta_local[26] 0.677 1.058 -1.147 2.534 1.031 130
#> hsgp_beta_local[27] -0.111 0.955 -1.997 1.685 1.016 56
#> hsgp_beta_local[28] 0.160 0.868 -1.618 1.687 1.039 103
#> hsgp_beta_local[29] 0.490 0.847 -1.189 1.991 1.000 113
#> hsgp_beta_local[2] 0.383 1.051 -1.420 2.427 0.995 125
#> hsgp_beta_local[30] 0.139 0.930 -1.586 2.031 1.026 108
#> hsgp_beta_local[31] -0.040 0.941 -2.062 1.652 1.047 117
#> hsgp_beta_local[32] -0.002 0.936 -1.991 1.760 1.005 98
#> hsgp_beta_local[33] -0.085 1.044 -2.073 1.803 0.997 82
#> hsgp_beta_local[34] 0.117 0.883 -1.571 1.740 1.014 144
#> hsgp_beta_local[35] 0.193 0.917 -1.593 1.908 1.055 68
#> hsgp_beta_local[36] 0.053 0.873 -1.874 1.599 1.081 116
#> hsgp_beta_local[3] 0.152 1.202 -2.234 2.603 1.010 110
#> hsgp_beta_local[4] 0.434 0.833 -1.272 2.313 1.007 58
#> hsgp_beta_local[5] 0.515 0.896 -1.484 2.196 1.009 80
#> hsgp_beta_local[6] 0.316 1.014 -1.408 2.419 1.021 108
#> hsgp_beta_local[7] -0.149 0.901 -1.736 1.489 1.079 83
#> hsgp_beta_local[8] -0.003 0.926 -1.478 1.822 1.054 90
#> hsgp_beta_local[9] -0.358 0.879 -1.935 1.453 1.003 91
#> hsgp_beta_regional[10] -0.101 0.964 -2.042 1.625 1.013 111
#> hsgp_beta_regional[11] -0.192 0.987 -2.191 1.691 1.000 153
#> hsgp_beta_regional[12] 0.050 1.009 -1.873 1.682 1.006 121
#> hsgp_beta_regional[13] -1.112 0.584 -2.216 0.002 1.065 36
#> hsgp_beta_regional[14] 0.906 0.575 -0.298 2.032 1.014 83
#> hsgp_beta_regional[15] -0.435 0.815 -1.836 1.154 0.999 110
#> hsgp_beta_regional[16] -0.127 1.026 -2.044 2.021 1.035 144
#> hsgp_beta_regional[17] -0.118 1.033 -2.401 1.585 1.002 140
#> hsgp_beta_regional[18] -0.080 1.014 -2.036 2.038 1.002 102
#> hsgp_beta_regional[19] 0.297 0.767 -1.200 1.710 1.039 54
#> hsgp_beta_regional[1] 0.241 0.748 -1.662 1.512 1.023 21
#> hsgp_beta_regional[20] 0.026 0.874 -2.066 1.626 1.022 64
#> hsgp_beta_regional[21] -0.213 0.969 -2.229 1.680 0.997 100
#> hsgp_beta_regional[22] -0.071 1.064 -1.951 2.081 1.046 66
#> hsgp_beta_regional[23] -0.065 1.097 -2.073 2.068 1.014 128
#> hsgp_beta_regional[24] 0.110 0.924 -1.573 2.007 1.007 102
#> hsgp_beta_regional[25] -0.383 0.936 -2.310 1.295 1.021 84
#> hsgp_beta_regional[26] 0.584 1.134 -1.589 2.875 1.010 46
#> hsgp_beta_regional[27] 0.139 1.115 -2.199 2.273 0.998 74
#> hsgp_beta_regional[28] 0.098 0.973 -1.726 1.844 1.017 52
#> hsgp_beta_regional[29] 0.080 1.205 -2.894 2.460 1.017 71
#> hsgp_beta_regional[2] 2.342 0.540 1.466 3.364 1.081 80
#> hsgp_beta_regional[30] -0.027 0.974 -1.885 1.969 1.018 112
#> hsgp_beta_regional[31] -0.001 0.862 -1.689 1.797 1.006 62
#> hsgp_beta_regional[32] 0.027 0.744 -1.229 1.483 0.998 122
#> hsgp_beta_regional[33] -0.130 1.097 -2.008 2.022 1.031 99
#> hsgp_beta_regional[34] 0.008 0.934 -1.696 1.633 1.031 100
#> hsgp_beta_regional[35] 0.025 0.970 -1.674 1.898 1.002 144
#> hsgp_beta_regional[36] -0.073 0.978 -1.843 1.786 1.018 90
#> hsgp_beta_regional[3] -0.118 0.689 -1.635 1.127 1.030 31
#> hsgp_beta_regional[4] 0.580 0.796 -1.014 1.976 1.004 69
#> hsgp_beta_regional[5] 0.387 1.144 -1.676 2.332 1.006 36
#> hsgp_beta_regional[6] 0.050 0.894 -1.758 1.806 1.009 96
#> hsgp_beta_regional[7] 0.059 0.248 -0.282 0.681 1.016 75
#> hsgp_beta_regional[8] -0.178 0.359 -0.878 0.470 1.001 76
#> hsgp_beta_regional[9] -0.900 0.666 -2.429 0.184 1.015 67
#> lengthscale_hsgp_local 0.316 0.274 0.090 1.083 1.023 51
#> lengthscale_hsgp_regional 1.142 0.213 0.796 1.609 1.005 24
#>
#> Diagnostics:
#> Divergences: 0
#> Avg. acceptance: 0.954
#> Chains: 1
#> Iterations: 400
#> Warmup: 200
#> R-hat > 1.01 for 47 parameter(s), largest 1.100 (hsgp_beta_local[10]).
#> Bulk ESS < 100 for 52 parameter(s), smallest 19 (num:(Intercept)).
#> Run longer chains or reparameterise before using these estimates.
# }