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.
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
#> Solver: auto (Cholesky<2k, PCG<5k, GPU/CG for larger)
#> 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 <- 50
lon <- runif(n, 0, 10)
lat <- runif(n, 0, 10)
# Make depth and temp spatially correlated
df <- data.frame(
lon = lon,
lat = lat,
depth = lon/5 + rnorm(n, 0, 0.5), # Correlated with lon
temp = lat/5 + rnorm(n, 0, 0.5), # Correlated with lat
count = rpois(n, 25),
effort = rgamma(n, shape = 4, rate = 1)
)
# Standard GP (may have spatial confounding)
fit1 <- tratio(
count | effort ~ depth + temp,
data = df,
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: negbin_negbin
#> 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: 60
#> Iterations: 200 (warmup: 100)
#> Chains: 1 (cores: 1)
#> Spatial: gp (50 units)
# RSR to protect depth and temp coefficients
fit2 <- tratio(
count | effort ~ depth + temp,
data = df,
spatial = spatial_rsr(
spatial_gp(~ lon + lat),
restrict_to = ~ depth + temp
),
mode = "hmc",
control = list(iter = 200, warmup = 100, chains = 1)
)
#> Inference: Exact (Tier 1)
#> Backend: hmc
#> Fitting ratio model...
#> Family: negbin_negbin
#> 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: 60
#> Iterations: 200 (warmup: 100)
#> Chains: 1 (cores: 1)
#> Spatial: gp (50 units)
# Compare coefficient estimates
summary(fit1) # May be attenuated
#> ratio model summary
#> ===================
#>
#> Inference: Exact (Tier 1) via hmc
#> Family: negbin_negbin
#>
#> Fixed effects (numerator):
#> parameter mean sd 2.5% 97.5%
#> beta_num[1] 755.970 0.005 755.960 755.977
#> beta_num[2] -787.023 0.007 -787.043 -787.014
#> beta_num[3] 4496.873 0.018 4496.849 4496.909
#>
#> Fixed effects (denominator):
#> parameter mean sd 2.5% 97.5%
#> beta_denom[1] 1454.784 0.009 1454.767 1454.798
#> beta_denom[2] -562.495 0.008 -562.510 -562.479
#> beta_denom[3] 951.196 0.007 951.181 951.206
#>
#> Variance components:
#> parameter mean sd 2.5% 97.5%
#> phi_num 8.339098e+253 Inf 7.290979e+253 9.824072e+253
#> phi_denom 0.000000e+00 0 0.000000e+00 0.000000e+00
#> sigma2_gp 1.786672e+248 Inf 1.255865e+248 2.349596e+248
#>
#> Diagnostics:
#> Divergences: 0
#> Avg. acceptance: 1
#> Chains: 1
#> Iterations: 200
#> Warmup: 100
summary(fit2) # Protected from spatial confounding
#> ratio model summary
#> ===================
#>
#> Inference: Exact (Tier 1) via hmc
#> Family: negbin_negbin
#>
#> Fixed effects (numerator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> beta_num[1] 3.008 0.016 2.960 3.027 1.136 8
#> beta_num[2] 0.021 0.010 0.005 0.046 1.056 10
#> beta_num[3] 0.085 0.028 0.042 0.133 1.814 2
#>
#> Fixed effects (denominator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> beta_denom[1] 1.276 0.026 1.243 1.332 1.947 2
#> beta_denom[2] -0.086 0.020 -0.129 -0.041 1.381 3
#> beta_denom[3] -0.242 0.029 -0.286 -0.168 1.217 4
#>
#> Diagnostics:
#> Divergences: 0
#> Avg. acceptance: 0.832
#> Chains: 1
#> Iterations: 200
#> Warmup: 100
# }