Apply Restricted Spatial Regression to mitigate spatial confounding. RSR orthogonalizes the spatial effect to the covariate space, preventing the spatial random effect from absorbing covariate information.
This is important when covariates are spatially smooth (e.g., climate variables, elevation) because the spatial random effect can "steal" variance from these covariates, leading to biased coefficient estimates.
Usage
spatial_rsr(spatial, restrict_to)
# S3 method for class 'ratiod_rsr'
print(x, ...)Arguments
- spatial
A spatial specification (
spatial_gp,spatial_car, etc.)- restrict_to
Formula specifying which covariates to orthogonalize against (e.g.,
~ depth + temp). The spatial effect will be constrained to be orthogonal to the column space of these covariates.- x
A ratiod_rsr object
- ...
Passed to underlying print method
Details
The RSR approach (Reich et al., 2006; Hodges & Reich, 2010) modifies the spatial random effect to be orthogonal to the fixed effect design matrix:
$$w_{RSR} = (I - P_X) w$$
where \(P_X = X(X'X)^{-1}X'\) is the projection matrix onto the column space of X.
When to use RSR:
Covariates are spatially smooth (environmental gradients)
Interested in causal interpretation of covariate effects
Coefficients appear attenuated toward zero
When NOT to use RSR:
Covariates are spatially uncorrelated
Spatial effect is the primary quantity of interest
Prediction is the main goal (not causal inference)
References
Reich, B. J., Hodges, J. S., & Zadnik, V. (2006). Effects of residual smoothing on the posterior of the fixed effects in disease-mapping models. Biometrics, 62(4), 1197-1206.
Hodges, J. S., & Reich, B. J. (2010). Adding spatially-correlated errors can mess up the fixed effect you love. The American Statistician, 64(4), 325-334.
Examples
# Create RSR spatial structure
rsr <- spatial_rsr(
spatial_gp(~ lon + lat),
restrict_to = ~ depth + temp
)
print(rsr)
#> tulpaRatio Gaussian Process spatial specification
#> =============================================
#>
#> Coordinates: lon, lat
#> Covariance: exponential
#> Neighbors (NNGP): 15
#> Shared: Yes (enters both processes)
#>
#> Restricted Spatial Regression (RSR):
#> Orthogonal to: ~depth + temp
#> (Spatial effect constrained to be orthogonal to covariate space)
# \donttest{
# Generate synthetic spatial data with confounding
set.seed(404)
n <- 40
lon <- runif(n, 0, 10)
lat <- runif(n, 0, 10)
# Smooth field w shared by count and effort
w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2)
# Make depth and temp spatially correlated
depth <- lon/5 + rnorm(n, 0, 0.5) # Correlated with lon
temp <- lat/5 + rnorm(n, 0, 0.5) # Correlated with lat
df <- data.frame(
lon = lon,
lat = lat,
depth = depth,
temp = temp,
count = rpois(n, exp(2.5 + 0.3 * depth - 0.2 * temp + w)),
effort = rgamma(n, shape = 4, rate = exp(-w))
)
# Standard GP (may have spatial confounding)
fit1 <- 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)
)
# RSR to protect depth and temp coefficients
fit2 <- tratio(
count | effort ~ depth + temp,
data = df,
family = ratiod_poisson_gamma(),
spatial = spatial_rsr(
spatial_gp(~ lon + lat),
restrict_to = ~ depth + temp
),
mode = "hmc",
control = list(iter = 400, chains = 1, seed = 1)
)
# Compare coefficient estimates
summary(fit1) # May be attenuated
#> 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.358 0.528 1.241 3.373 1.954 1
#> num:depth 0.316 0.096 0.129 0.491 1.057 26
#> num:temp -0.107 0.144 -0.365 0.138 1.026 21
#>
#> Fixed effects (denominator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> denom:(Intercept) 1.497 0.536 0.359 2.484 1.962 1
#> denom:depth -0.018 0.115 -0.265 0.177 1.098 19
#> denom:temp -0.005 0.145 -0.282 0.242 1.019 24
#>
#> Variance components:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> shape 7.970 2.215 4.569 12.653 1.069 25
#> sigma2_gp 0.843 0.496 0.260 2.149 0.998 66
#>
#> Other parameters:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> gp_w[10] 0.575 0.650 -0.275 1.682 1.700 2
#> gp_w[11] -0.119 0.594 -1.011 0.980 1.661 2
#> gp_w[12] 0.074 0.573 -1.020 1.058 1.804 2
#> gp_w[13] -0.068 0.631 -0.882 1.104 1.805 2
#> gp_w[14] -0.708 0.655 -1.613 0.508 1.605 2
#> gp_w[15] -0.723 0.659 -1.681 0.415 1.855 2
#> gp_w[16] -0.609 0.677 -1.559 0.516 1.486 2
#> gp_w[17] -0.040 0.649 -0.971 1.207 1.770 2
#> gp_w[18] 0.161 0.653 -0.717 1.288 1.539 2
#> gp_w[19] -0.965 0.656 -1.995 0.221 1.673 2
#> gp_w[1] -0.735 0.677 -1.769 0.419 1.566 2
#> gp_w[20] -0.267 0.658 -1.220 0.796 1.492 2
#> gp_w[21] 0.377 0.698 -0.590 1.474 1.525 2
#> gp_w[22] -0.927 0.679 -2.050 0.238 1.732 2
#> gp_w[23] 0.903 0.560 -0.126 2.052 1.869 2
#> gp_w[24] 0.843 0.574 -0.066 1.825 1.751 2
#> gp_w[25] 0.606 0.657 -0.235 1.759 1.581 2
#> gp_w[26] -0.044 0.571 -0.879 0.990 1.768 2
#> gp_w[27] -0.390 0.635 -1.275 0.781 1.805 2
#> gp_w[28] 0.699 0.558 -0.331 1.955 1.678 2
#> gp_w[29] -0.301 0.599 -1.298 0.676 1.791 2
#> gp_w[2] 0.199 0.702 -0.759 1.463 1.538 2
#> gp_w[30] 0.434 0.555 -0.627 1.431 1.853 2
#> gp_w[31] 0.150 0.641 -0.743 1.275 1.599 2
#> gp_w[32] -0.852 0.648 -1.732 0.270 1.786 2
#> gp_w[33] -0.513 0.636 -1.447 0.679 1.778 2
#> gp_w[34] -0.671 0.647 -1.615 0.548 1.711 2
#> gp_w[35] -0.778 0.649 -1.752 0.380 1.588 2
#> gp_w[36] 0.427 0.558 -0.650 1.472 1.697 2
#> gp_w[37] -0.670 0.659 -1.586 0.426 1.605 2
#> gp_w[38] -0.954 0.666 -1.999 0.200 1.643 2
#> gp_w[39] 1.110 0.574 -0.044 2.211 1.927 1
#> gp_w[3] 0.422 0.578 -0.470 1.455 1.698 2
#> gp_w[40] -0.782 0.664 -1.593 0.454 1.605 2
#> gp_w[4] 0.597 0.563 -0.461 1.535 1.804 2
#> gp_w[5] 1.056 0.545 0.184 2.102 1.809 2
#> gp_w[6] 1.145 0.559 0.256 2.181 1.833 2
#> gp_w[7] 0.253 0.634 -0.574 1.374 1.720 2
#> gp_w[8] -0.076 0.560 -1.118 1.042 1.865 2
#> gp_w[9] 0.148 0.623 -0.761 1.216 1.587 2
#> phi_gp 4.865 2.953 1.627 12.071 1.000 69
#>
#> Diagnostics:
#> Divergences: 0
#> Avg. acceptance: 0.917
#> Chains: 1
#> Iterations: 400
#> Warmup: 200
#> R-hat > 1.01 for 47 parameter(s), largest 1.962 (denom:(Intercept)).
#> Bulk ESS < 100 for 49 parameter(s), smallest 1 (denom:(Intercept)).
#> Run longer chains or reparameterise before using these estimates.
summary(fit2) # Protected from spatial confounding
#> 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.091 0.092 2.903 3.245 1.008 90
#> num:depth -0.079 0.067 -0.218 0.048 1.118 77
#> num:temp -0.384 0.048 -0.466 -0.284 1.010 118
#>
#> Fixed effects (denominator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> denom:(Intercept) 2.232 0.136 1.901 2.522 1.054 64
#> denom:depth -0.419 0.084 -0.572 -0.253 1.035 67
#> denom:temp -0.291 0.075 -0.435 -0.153 1.026 79
#>
#> Variance components:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> shape 7.568 1.665 4.808 11.045 1.001 74
#> sigma2_gp 0.919 0.878 0.323 2.537 1.057 33
#>
#> Other parameters:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> gp_w[10] 0.704 0.440 -0.093 1.461 1.151 5
#> gp_w[11] 0.064 0.471 -0.942 0.820 1.141 6
#> gp_w[12] 0.300 0.461 -0.594 1.029 1.194 5
#> gp_w[13] 0.068 0.452 -0.881 0.838 1.235 7
#> gp_w[14] -0.580 0.489 -1.532 0.198 1.138 6
#> gp_w[15] -0.551 0.470 -1.480 0.266 1.144 8
#> gp_w[16] -0.526 0.467 -1.342 0.336 1.176 5
#> gp_w[17] 0.114 0.464 -0.736 0.866 1.151 5
#> gp_w[18] 0.327 0.451 -0.579 1.118 1.137 8
#> gp_w[19] -0.823 0.491 -1.865 -0.026 1.128 6
#> gp_w[1] -0.595 0.456 -1.475 0.207 1.153 5
#> gp_w[20] -0.127 0.470 -1.037 0.693 1.171 5
#> gp_w[21] 0.539 0.447 -0.311 1.365 1.177 6
#> gp_w[22] -0.734 0.501 -1.804 0.140 1.121 8
#> gp_w[23] 1.116 0.447 0.136 1.799 1.139 8
#> gp_w[24] 1.054 0.448 0.080 1.691 1.151 5
#> gp_w[25] 0.726 0.454 -0.184 1.552 1.158 7
#> gp_w[26] 0.147 0.452 -0.809 0.918 1.181 4
#> gp_w[27] -0.227 0.471 -1.031 0.561 1.167 5
#> gp_w[28] 0.961 0.477 -0.096 1.648 1.129 6
#> gp_w[29] -0.106 0.450 -1.041 0.641 1.171 4
#> gp_w[2] 0.322 0.452 -0.453 1.164 1.138 7
#> gp_w[30] 0.648 0.460 -0.306 1.305 1.174 8
#> gp_w[31] 0.296 0.449 -0.630 1.087 1.137 6
#> gp_w[32] -0.677 0.493 -1.592 0.155 1.129 6
#> gp_w[33] -0.352 0.446 -1.240 0.435 1.279 7
#> gp_w[34] -0.519 0.519 -1.605 0.431 1.163 5
#> gp_w[35] -0.661 0.464 -1.537 0.134 1.123 6
#> gp_w[36] 0.602 0.475 -0.326 1.346 1.208 9
#> gp_w[37] -0.551 0.463 -1.498 0.224 1.132 6
#> gp_w[38] -0.793 0.481 -1.721 0.040 1.093 8
#> gp_w[39] 1.344 0.471 0.329 2.007 1.186 8
#> gp_w[3] 0.611 0.433 -0.341 1.294 1.144 6
#> gp_w[40] -0.642 0.517 -1.679 0.243 1.155 8
#> gp_w[4] 0.781 0.441 -0.134 1.516 1.128 7
#> gp_w[5] 1.290 0.455 0.425 1.968 1.160 6
#> gp_w[6] 1.357 0.461 0.349 2.021 1.203 5
#> gp_w[7] 0.374 0.440 -0.493 1.175 1.147 6
#> gp_w[8] 0.195 0.492 -0.854 0.979 1.165 5
#> gp_w[9] 0.261 0.461 -0.696 1.057 1.147 6
#> phi_gp 6.158 6.228 1.619 20.720 0.998 47
#>
#> Diagnostics:
#> Divergences: 0
#> Avg. acceptance: 0.859
#> Chains: 1
#> Iterations: 400
#> Warmup: 200
#> R-hat > 1.01 for 45 parameter(s), largest 1.279 (gp_w[33]).
#> Bulk ESS < 100 for 48 parameter(s), smallest 4 (gp_w[26]).
#> Run longer chains or reparameterise before using these estimates.
# }