如何在R中通过置换检验完成双矩阵比较 分析林木空间关联与功能性状关系

检验逻辑说明 #

我们会严格保留空间相互作用矩阵的正负关联结构（即值为3的正关联位置、值为1的负关联位置固定不变），仅将性状距离矩阵的数值分别在正关联单元格、负关联单元格两个集合内独立置换，生成无关联的零分布，再对比观测值与零分布的差异显著性。

前置校验 #

你需要先确认两份csv文件满足以下要求：

• 两个矩阵的行名、列名均为物种名，且顺序完全一致
• 矩阵为对称矩阵（物种两两配对的数值上下三角一致）
• 没有缺失值

完整R脚本 #

# 1. 读入数据，请替换引号内为你的文件实际路径
spatial_mat <- as.matrix(read.csv("空间矩阵文件路径.csv", row.names = 1, check.names = FALSE))
trait_mat <- as.matrix(read.csv("性状距离矩阵文件路径.csv", row.names = 1, check.names = FALSE))

# 2. 校验矩阵一致性，不一致会直接报错提示
if(!identical(dim(spatial_mat), dim(trait_mat)) || !identical(rownames(spatial_mat), rownames(trait_mat)) || !identical(colnames(spatial_mat), colnames(trait_mat))){
  stop("两个矩阵的维度、行名或列名不匹配，请检查输入文件")
}

# 3. 提取下三角索引避免重复计算成对数据（对称矩阵无需重复统计上下三角）
lower_tri_idx <- lower.tri(spatial_mat)
spatial_lower <- spatial_mat[lower_tri_idx]
trait_lower <- trait_mat[lower_tri_idx]

# 4. 分别提取正、负相互作用对应的性状值索引
pos_spatial_idx <- which(spatial_lower == 3)
neg_spatial_idx <- which(spatial_lower == 1)

# 5. 计算观测统计量：这里我们用正负组性状平均距离的差值作为统计量，你也可以按需换成相关系数等其他统计量
obs_pos_mean <- mean(trait_lower[pos_spatial_idx])
obs_neg_mean <- mean(trait_lower[neg_spatial_idx])
obs_diff <- obs_pos_mean - obs_neg_mean

# 6. 置换检验设置
n_perm <- 999 # 置换次数，可调整为9999获得更精确的p值，运行时间会对应增加
perm_diff <- numeric(n_perm) # 存储每次置换的统计量结果

# 7. 置换循环
for(i in 1:n_perm){
  # 正负组独立置换性状值，空间关联位置完全保留
  perm_trait <- trait_lower
  perm_trait[pos_spatial_idx] <- sample(perm_trait[pos_spatial_idx])
  perm_trait[neg_spatial_idx] <- sample(perm_trait[neg_spatial_idx])
  
  # 计算本次置换的统计量
  perm_pos_mean <- mean(perm_trait[pos_spatial_idx])
  perm_neg_mean <- mean(perm_trait[neg_spatial_idx])
  perm_diff[i] <- perm_pos_mean - perm_neg_mean
}

# 8. 计算显著性p值（双尾检验）
p_value <- (sum(abs(perm_diff) >= abs(obs_diff)) + 1) / (n_perm + 1)

# 9. 输出结果
cat("观测正关联组平均性状距离：", round(obs_pos_mean, 4), "\n")
cat("观测负关联组平均性状距离：", round(obs_neg_mean, 4), "\n")
cat("两组平均距离差值：", round(obs_diff, 4), "\n")
cat("置换检验p值：", round(p_value, 4), "\n")

# 可选：绘制置换分布直方图直观查看
hist(perm_diff, main = "置换检验零分布", xlab = "正负组平均性状距离差值", col = "lightblue")
abline(v = obs_diff, col = "red", lwd = 2, lty = 2)

结果解读 #

• 若p值 < 0.05，说明空间关联类型和功能性状距离的关联达到统计显著性，不是随机产生的
• 若观测差值为负，说明正关联的物种对平均性状距离更小，性状更相似；若为正则反之
• 直方图中红色虚线为观测差值，若落在置换分布的尾部则说明关联显著

注意事项 #

• 如果你的统计需求是计算空间矩阵和性状矩阵的整体相关而非组间均值差，可以将第5步的统计量替换为obs_cor <- cor(spatial_lower, trait_lower, method = "spearman")，对应置换步骤的统计量也同步替换即可
• 若矩阵存在缺失值，需要先做缺失值处理再运行脚本

内容的提问来源于stack exchange，提问作者Mike Harvey Salazar Villegas