含缺失值的网格数据计算Geary's C/Moran's I空间自相关指数

针对12×3网格数据中存在缺失值、spdep包geary()函数直接返回NA的问题，最快的解决方式是基于Geary's C的原始公式手动实现，仅处理有效观测单元及其空间邻接关系，步骤如下：

1. 数据预处理与有效单元筛选 #

先将网格数据转为一维向量，标记并提取非缺失的观测单元：

library(spdep)

# 原始网格数据
loc <- structure(c(0.653, 1.361, 1.392, 0.991, 1.228, 3.882, 1.46, 1.056, 
1.636, 1.535, 1.578, 2.106, NA, 0.801, NA, 1.51, 1.514, 1.846, 
0.791, 1.567, 0.902, 2.011, 1.683, 1.541, 3.437, 1.095, 1.469, 
1.897, 1.866, 1.498, 1.156, 1.342, 1.809, 1.56, 1.341, 1.005), dim = c(12L, 
3L), dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8", 
"9", "10", "11", "12"), c("1", "2", "3")))

# 转为一维向量，筛选非缺失单元
x_vec <- as.numeric(t(loc))
valid_idx <- which(!is.na(x_vec))
x_valid <- x_vec[valid_idx]
n_valid <- length(x_valid)

2. 过滤空间邻接关系 #

从原始queen邻接关系中，仅保留有效单元之间的邻接对：

# 生成原始网格邻接关系
adj <- cell2nb(nrow = nrow(loc), ncol = ncol(loc), type = "queen")

# 过滤邻接：仅保留有效单元的邻接，且邻接对象也为有效单元
adj_valid <- lapply(adj[valid_idx], function(neighbors) {
  intersect(neighbors, valid_idx)
})

# 转为空间权重矩阵（二进制权重）
ww_valid <- nb2listw(adj_valid, style = "B", zero.policy = TRUE)
S0_valid <- Szero(ww_valid)  # 权重总和

3. 手动计算Geary's C #

基于Geary's C的核心公式计算：
$$C = \frac{n-1}{2S_0} \times \frac{\sum_{i}\sum_{j} w_{ij}(x_i - x_j)^2}{\sum_{i}(x_i - \bar{x})^2}$$

# 计算分子：所有有效邻接对的权重加权平方差之和
numerator <- 0
for (i in seq_along(adj_valid)) {
  neighbors <- adj_valid[[i]]
  if (length(neighbors) > 0) {
    # 匹配邻接单元在有效数组中的位置
    neighbor_pos <- match(neighbors, valid_idx)
    numerator <- numerator + sum(ww_valid$weights[[i]] * (x_valid[i] - x_valid[neighbor_pos])^2)
  }
}

# 计算分母：有效观测值的离均差平方和
x_bar <- mean(x_valid)
denominator <- sum((x_valid - x_bar)^2)

# 计算最终Geary's C值
geary_c <- (n_valid - 1) / (2 * S0_valid) * numerator / denominator
print(geary_c)

方法优势 #

• 直接基于公式实现，跳过缺失值干扰，仅计算有效单元的空间关联；
• 针对12×3的小规模网格，循环计算速度极快，无需额外优化；
• 若需计算Moran's I，可通过替换对应公式快速修改代码。

内容的提问来源于stack exchange，提问作者Philopolis