如何迭代创建站点专属矩阵/向量，批量计算Moran's I并生成结果矩阵？

针对你12个站点、2400个分类群的批量计算需求，不要手动创建大量命名对象，用tidyverse的分组嵌套+映射操作可以高效完成，同时避免环境变量混乱。以下是具体实现步骤：

1. 加载依赖包 #

library(tidyverse)
library(ape) # 用于Moran.I函数

2. 预处理：生成每个站点的逆距离权重矩阵 #

先从原始数据中提取每个站点的坐标，批量生成逆距离矩阵，存储在嵌套数据框中（无需单独命名每个矩阵）：

# 提取每个站点的坐标，过滤样本数≥2的站点（避免距离计算报错）
site_coords <- data %>%
  select(Site, x_coord, y_coord) %>% # 注意：你原代码里的y_cood应为y_coord，需修正拼写
  group_by(Site) %>%
  filter(n() >= 2) %>%
  nest() %>%
  # 批量生成逆距离权重矩阵
  mutate(
    dist_inv = map(data, ~{
      dist_mat <- as.matrix(dist(select(.x, x_coord, y_coord)))
      dist_inv_mat <- 1 / dist_mat
      diag(dist_inv_mat) <- 0 # 对角线设为0（自身无距离权重）
      dist_inv_mat
    })
  ) %>%
  select(Site, dist_inv) # 保留站点和对应的逆距离矩阵

# 提取站点与处理组的对应关系（每个站点对应唯一Treatment）
site_treatment_map <- data %>%
  select(Site, Treatment) %>%
  distinct()

3. 转换残差数据为长格式 #

将宽格式的残差数据转为长格式，方便按站点+分类群分组处理：

resid_long <- Resid_data %>%
  pivot_longer(
    cols = starts_with("Abundance_Taxa"), # 匹配所有分类群残差列
    names_to = "Taxa",
    values_to = "Residual"
  )

4. 批量计算Moran's I并整理结果 #

结合逆距离矩阵和残差数据，按站点+分类群分组计算，最后转成你需要的宽格式结果矩阵：

# 批量计算Moran's I
moran_raw <- resid_long %>%
  left_join(site_coords, by = "Site") %>%
  left_join(site_treatment_map, by = "Site") %>%
  group_by(Site, Treatment, Taxa) %>%
  summarise(
    Moran_I = Moran.I(Residual, dist_inv[[1]])$observed, # 提取Moran's I值
    P_value = Moran.I(Residual, dist_inv[[1]])$p.value, # 提取P值
    .groups = "drop"
  )

# 转换为宽格式结果矩阵（每个分类群对应一列Moran_I和P_value）
final_result <- moran_raw %>%
  pivot_wider(
    id_cols = c(Site, Treatment),
    names_from = Taxa,
    values_from = c(Moran_I, P_value),
    names_glue = "{.value}_{str_remove(Taxa, 'Abundance_')}" # 重命名列，去掉Abundance_前缀
  )

关键说明 #

• 全程无需创建SiteI1_data这类命名对象，用嵌套数据框和映射操作批量处理，代码更简洁且易维护。
• 自动过滤样本数不足的站点，避免计算报错。
• 最终结果final_result包含Site、Treatment列，以及每个分类群对应的Moran_I_TaxaX和P_value_TaxaX列，完全符合你的需求。

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