空间依赖数据组检验：自然保护区与周边缓冲区干扰率差异分析

嘿，针对你这个100个森林保护区+4层非重叠缓冲区的干扰率对比需求，我整理了一套从方法选择到结果解读的完整思路，帮你高效验证零假设、定位高干扰区域：

一、先把核心目标理清楚 #

我们要做的是：检验保护区本体和100m/500m/1000m/2000m四个缓冲区这5组区域的年度干扰率是否存在显著差异，核心搞清楚哪块区域干扰率最高。零假设是「所有区域的干扰率完全相同」，而且你已经搞定了Global Forest Watch数据的标准化（这点太重要了，避免了面积大小干扰结果），省去了不少预处理麻烦。

二、选对统计方法是关键 #

因为是多组独立样本的均值/分布对比，优先推荐这些方法：

• 单因素方差分析（One-way ANOVA）：先整体检验5组间干扰率是否存在显著差异。如果ANOVA的p值<0.05，说明至少有一组和其他组不一样，但没法直接知道是哪几组。
• 事后多重比较：ANOVA显著后，就得用这个定位具体差异：
  • 要是你核心关注「缓冲区 vs 保护区」的对比，选Dunnett's检验，它专门针对“多处理组vs单一对照组”的场景，能减少误判概率；
  • 要是还想搞清楚缓冲区之间的差异，选Tukey's HSD检验，能覆盖所有组的两两对比。
• 如果数据不符合ANOVA的正态/方差齐性要求，就用非参数替代：
  • 整体检验用Kruskal-Wallis H检验；
  • 事后检验用Dunn's检验搭配Bonferroni校正。

三、数据准备要踩准细节 #

虽然你说数据已经标准化，但还是要确认这几点：

• 每个保护区的5个区域都有完整数据，没有缺失值；要是有缺失，要么删掉对应保护区样本，要么用同类型保护区的均值填充（别瞎填就行）；
• 检查数据分布：用直方图、Q-Q图看正态性，用Levene/Bartlett检验看方差齐性，这直接决定你用参数还是非参数方法；
• 再确认一遍缓冲区是非重叠的（你已经提过，这点必须保证，不然数据会重复计算，结果全错）。

四、拿R举个实操例子 #

假设你的数据框叫disturbance_data，列包括reserve_id（保护区ID）、core（保护区干扰率）、buf100到buf2000四个缓冲区的干扰率：

1. 先把宽格式转成长格式（统计方法都认长格式）：

library(tidyr)
long_data <- pivot_longer(disturbance_data, 
                          cols = c(core, buf100, buf500, buf1000, buf2000),
                          names_to = "region",
                          values_to = "disturbance_rate")

2. 跑单因素ANOVA：

anova_result <- aov(disturbance_rate ~ region, data = long_data)
summary(anova_result)

看输出里的p值，小于0.05就说明组间有显著差异，拒绝零假设。

3. 用Dunnett's检验做事后对比（以保护区为对照组）：

library(emmeans)
emmeans_result <- emmeans(anova_result, pairwise ~ region, adjust = "dunnett", ref = "core")
print(emmeans_result)

结果里会明确显示每个缓冲区和保护区的差异是否显著，以及均值高低。

4. 要是用非参数方法：

# Kruskal-Wallis整体检验
kruskal_result <- kruskal.test(disturbance_rate ~ region, data = long_data)
print(kruskal_result)

# Dunn's事后检验
library(FSA)
dunn_result <- dunnTest(disturbance_rate ~ region, data = long_data, method = "bonferroni")
print(dunn_result)

五、结果解读抓重点 #

• 先看整体检验的p值：显著=干扰率存在组间差异；
• 盯着事后检验的结果：找p值显著的组，对比均值就能知道哪个缓冲区干扰率显著高于（或低于）保护区，以及缓冲区之间的排名；
• 画个箱线图更直观：

library(ggplot2)
ggplot(long_data, aes(x = region, y = disturbance_rate)) +
  geom_boxplot(fill = "#4CAF50") +
  labs(title = "标准化干扰率对比：保护区 vs 缓冲区", 
       x = "区域类型", 
       y = "年度标准化干扰率") +
  theme_minimal()

六、几个容易踩坑的点 #

• 别搞混重复测量和独立样本：这里每个保护区的5个区域是独立空间单元，不是同一单元的重复测量，所以不用重复测量ANOVA；
• 标准化方法要靠谱：确认是「干扰面积/区域总面积」这种合理的标准化方式，不然结果毫无意义；
• 要是保护区有空间聚类（比如同一省份的保护区干扰模式类似），可以加个混合效应模型（用lme4包），把区域作为随机效应，控制空间自相关的影响。

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