在R语言中如何检验两类坐标数据框的统计相似性？

针对经纬度坐标这类空间点数据，以下是几种实用的统计检验方法，附R代码实现：

1. 核密度估计+Kolmogorov-Smirnov检验 #

通过核密度估计拟合两个点集的空间分布，再用KS检验比较密度分布是否存在显著差异。注意先将经纬度转换为平面投影（如UTM），避免球面距离干扰。

library(sf)
library(ks)

# 转换为sf对象并投影到UTM
A_sf <- st_as_sf(A, coords = c("lon", "lat"), crs = 4326)
B_sf <- st_as_sf(B, coords = c("lon", "lat"), crs = 4326)
utm_crs <- st_transform(A_sf, crs = st_crs(A_sf)$utm)$crs
A_utm <- st_transform(A_sf, crs = utm_crs)
B_utm <- st_transform(B_sf, crs = utm_crs)

# 提取平面坐标
A_coords <- st_coordinates(A_utm)
B_coords <- st_coordinates(B_utm)

# 拟合核密度
kd_A <- kde(A_coords)
kd_B <- kde(B_coords)

# 生成公共网格并计算密度值
grid <- expand.grid(
  x = seq(min(c(A_coords[,1], B_coords[,1])), max(c(A_coords[,1], B_coords[,1])), length.out = 100),
  y = seq(min(c(A_coords[,2], B_coords[,2])), max(c(A_coords[,2], B_coords[,2])), length.out = 100)
)
dens_A <- predict(kd_A, x = grid)
dens_B <- predict(kd_B, x = grid)

# KS检验密度分布差异
ks.test(dens_A, dens_B)

若检验结果的p值>0.05，无法拒绝“两个点集密度分布相似”的原假设。

2. 空间点格局的蒙特卡洛检验 #

基于最近邻距离的G函数，通过蒙特卡洛模拟比较两个点集的空间聚集/分散格局是否一致。

library(spatstat)

# 创建spatstat所需的ppp对象
win <- owin(
  xrange = range(c(A_coords[,1], B_coords[,1])),
  yrange = range(c(A_coords[,2], B_coords[,2]))
)
A_ppp <- ppp(A_coords[,1], A_coords[,2], window = win)
B_ppp <- ppp(B_coords[,1], B_coords[,2], window = win)

# 计算G函数（最近邻距离累积分布）
G_A <- Gest(A_ppp)
G_B <- Gest(B_ppp)

# 蒙特卡洛模拟检验
n_sim <- 999
diff_stats <- numeric(n_sim)
combined_coords <- rbind(A_coords, B_coords)
n_A <- nrow(A_coords)

for(i in 1:n_sim){
  idx <- sample(nrow(combined_coords), n_A)
  sim_A <- ppp(combined_coords[idx,1], combined_coords[idx,2], window = win)
  sim_B <- ppp(combined_coords[-idx,1], combined_coords[-idx,2], window = win)
  sim_GA <- Gest(sim_A)
  sim_GB <- Gest(sim_B)
  diff_stats[i] <- max(abs(sim_GA$rs - sim_GB$rs))
}

# 计算实际差异与p值
actual_diff <- max(abs(G_A$rs - G_B$rs))
p_val <- (sum(diff_stats >= actual_diff) + 1)/(n_sim + 1)
cat("蒙特卡洛检验p值：", p_val, "\n")

p值>0.05说明两个点集的空间格局无显著差异。

3. 多元方差分析（MANOVA） #

将经度、纬度作为两个连续变量，直接检验两组坐标的均值向量是否存在显著差异，适用于比较分布中心的相似性。

# 合并数据并添加分组标签
combined_df <- rbind(
  cbind(A, group = "A"),
  cbind(B, group = "B")
)

# 执行MANOVA
manova_result <- manova(cbind(lon, lat) ~ group, data = combined_df)
summary(manova_result)

MANOVA的p值>0.05表示两组的坐标中心无显著差异，但该方法不考虑分布的离散程度或空间格局。

4. 空间置换检验 #

通过随机置换分组标签，检验两组的空间统计量（如平均中心距离）是否由随机因素导致，可灵活适配不同的统计指标。

# 计算实际平均中心距离
mean_center_A <- colMeans(A_coords)
mean_center_B <- colMeans(B_coords)
actual_dist <- sqrt(sum((mean_center_A - mean_center_B)^2))

# 置换模拟
n_sim <- 999
sim_dists <- numeric(n_sim)
group_labels <- combined_df$group

for(i in 1:n_sim){
  permuted_labels <- sample(group_labels)
  perm_A <- combined_df[permuted_labels == "A", c("lon", "lat")]
  perm_B <- combined_df[permuted_labels == "B", c("lon", "lat")]
  perm_center_A <- colMeans(perm_A)
  perm_center_B <- colMeans(perm_B)
  sim_dists[i] <- sqrt(sum((perm_center_A - perm_center_B)^2))
}

# 计算双侧检验p值
p_val <- (sum(sim_dists >= actual_dist) + sum(sim_dists <= actual_dist) + 1)/(n_sim + 1)
cat("置换检验p值：", p_val, "\n")

p值>0.05说明两组的平均中心差异不显著，分布中心相似。

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