考虑空间自相关的两幅栅格数据相关性系数校正问题
好问题!确实,Dutilleul检验(也就是你用的SpatialPack里的方法)核心是通过调整有效样本量修正相关性检验的显著性(p值),但它不会改变皮尔逊相关系数的点估计——毕竟这个系数的计算还是基于所有原始栅格像素对的协方差。如果想要校正相关系数的估计值,你得把空间自相关的影响直接纳入到关联强度的建模中,下面给你几个实用的方案(基于R环境,和你用的SpatialPack生态兼容):
方案1:空间回归框架下的调整系数
空间回归模型(比如空间滞后模型SLM、空间误差模型SEM)会直接把空间依赖性作为模型的一部分,拟合后得到的自变量系数可以看作是控制了空间自相关后的关联强度估计。步骤如下:
把栅格数据转换为包含坐标的数据集:
library(raster) # 假设你的两个栅格是r1和r2 df <- as.data.frame(stack(r1, r2), xy = TRUE) colnames(df) <- c("x", "y", "var1", "var2")构建空间权重矩阵(定义像素的空间邻接关系):
library(spdep) # 基于距离的邻接(比如距离小于5个单位的像素为邻域) coords <- df[, c("x", "y")] nb <- dnearneigh(coords, d1 = 0, d2 = 5) w <- nb2listw(nb, style = "W") # 行标准化权重拟合空间滞后模型并提取调整后的系数:
# 拟合SLM模型,var2为因变量,var1为自变量 slm_model <- lagsarlm(var2 ~ var1, data = df, listw = w) summary(slm_model) # 提取var1的系数作为校正后的关联强度 corrected_rho <- coef(slm_model)["var1"]注意:这里的系数和皮尔逊相关系数不是直接等价的,但它反映了两个变量在去除空间自相关影响后的线性关联强度。
方案2:基于变异函数的协方差校正
变异函数可以刻画空间数据的自相关结构,我们可以用交叉变异函数的参数来调整相关系数,本质是只考虑空间依赖部分的协方差,而非原始的全局协方差。步骤如下:
计算并拟合交叉变异函数:
library(gstat) # 构建gstat对象,指定两个变量的交叉变异 vgm_cross <- variogram(var1 ~ var2, data = df, locations = ~x+y) # 拟合球状变异函数模型(也可以用指数、高斯模型) fit_vgm <- fit.variogram(vgm_cross, vgm("Sph"))基于变异函数参数计算校正后的相关系数:
皮尔逊相关系数是rho = cov(var1, var2)/(sd(var1)*sd(var2)),而校正后的系数可以只考虑空间结构部分的协方差(排除块金效应,也就是随机噪声部分):# 提取两个变量的自变异函数参数(单独拟合) vgm_var1 <- variogram(var1 ~ 1, data = df, locations = ~x+y) fit_var1 <- fit.variogram(vgm_var1, vgm("Sph")) vgm_var2 <- variogram(var2 ~ 1, data = df, locations = ~x+y) fit_var2 <- fit.variogram(vgm_var2, vgm("Sph")) # 计算空间结构部分的协方差和标准差 cov_spatial <- fit_vgm$psill[fit_vgm$model == "Sph"] sd_spatial_var1 <- sqrt(fit_var1$psill[fit_var1$model == "Sph"]) sd_spatial_var2 <- sqrt(fit_var2$psill[fit_var2$model == "Sph"]) # 校正后的相关系数 corrected_rho <- cov_spatial / (sd_spatial_var1 * sd_spatial_var2)这个方法的前提是数据满足空间平稳性假设。
方案3:空间Bootstrap的非参数校正
如果你的数据不满足空间平稳性,非参数的空间Bootstrap是个不错的选择——它通过空间块重采样来模拟空间自相关下的样本分布,进而得到校正后的相关系数估计:
library(spatialEco) # 空间块Bootstrap,计算相关系数的分布 boot_result <- spatial.bootstrap( x = df, formula = var2 ~ var1, statistic = function(data) cor(data$var1, data$var2), n = 1000, # 重采样次数 type = "block", # 空间块重采样 size = 5 # 块的大小 ) # 取Bootstrap样本的中位数或均值作为校正后的相关系数 corrected_rho <- median(boot_result$t) # 同时可以得到校正后的置信区间 boot_ci <- quantile(boot_result$t, c(0.025, 0.975))
关键注意事项
- 这些方法的核心都是分离空间自相关和变量间的真实关联,不同方法的假设不同,你需要根据你的栅格数据的空间特征(比如是否平稳、自相关的范围)选择合适的方法。
- 你可以把校正后的相关系数和Dutilleul检验结合使用:先用上述方法得到校正后的点估计,再用SpatialPack的
modified.ttest()计算校正后的显著性p值,这样既得到了准确的关联强度,也得到了可靠的显著性检验结果。
内容的提问来源于stack exchange,提问作者tsutsume

