使用vegan包处理存在-缺失数据时tswap空模型运行失败求助

问题描述 #

针对带样地层级结构的存在-缺失数据开展beta多样性分区的空模型分析，选用推荐的tswap方法时无法运行，仅r00等方法可正常执行，已知错误源于未通过合理性检查，但疑惑tswap本应适配存在-缺失矩阵。测试代码如下：

library(parallel)
library(vegan)

# Setting up data ----
# community
com <- matrix(rbinom(120, 1, 0.5),nrow=12)
# plot hirachy 
plt_hir <- data.frame(plot=c(1:nrow(com)),
                      treatment = c(rep("reg1_t1",nrow(com)/4),rep("reg1_t2",nrow(com)/4),rep("reg2_t1",nrow(com)/4),rep("reg2_t2",nrow(com)/4)),
                      region=c(rep("reg1",nrow(com)/2),rep("reg2",nrow(com)/2)),
                      country=rep("country",nrow(com)))
# to numeric
plt_hir <- data.frame(apply(plt_hir, 2, function(x) as.numeric(factor(x))))


# Multipart ----

# matrix
y = com
# level matrix
x = plt_hir
# number of null model computations
nsimul = 100
# null model to be computed
method = "tswap"
# get and set number of cores for parallel computation (do not use all of them otherwise computer will be slowed)
cors <- detectCores()-1


# NULL MODEL ----

# define null model
nm <- nullmodel(com, method)

## socket type cluster, works on all platforms
cl <- makeCluster(cors)
clusterEvalQ(cl, library(vegan))
clusterExport(cl, c("simulate", "x", "nm"))
set.seed(3)
smlist <- parLapply(cl, 1:nsimul, function(i)  simulate(nm, nsim = 1, thin = 1))
stopCluster(cl)
smlist <- smbind(smlist, MARGIN=3)

问题排查与解决 #

• 核心原因：tswap对矩阵的行/列物种数有最低要求
  tswap算法依赖物种在样地间的交换操作，要求每个样地（行）至少包含2个物种，同时每个物种（列）至少出现在2个样地中。你的模拟数据用rbinom(120,1,0.5)生成，大概率会出现仅含0/1个物种的样地，或仅出现1次的物种，触发vegan内置的合理性检查，导致运行失败。

• 验证方式
  执行以下代码检查矩阵的行/列物种数分布：

  # 查看每个样地的物种数量
  rowSums(com)
  # 查看每个物种的出现次数
  colSums(com)
  

  若输出结果中存在0或1的数值，即为触发错误的直接原因。

• 解决办法

  1. 过滤不符合要求的行/列：移除不满足条件的样地和物种，确保剩余矩阵符合tswap的运行要求：

     # 保留至少含2个物种的样地
     com_filtered <- com[rowSums(com) >= 2, ]
     # 保留至少出现在2个样地的物种
     com_filtered <- com_filtered[, colSums(com_filtered) >= 2]
     # 同步更新样地层级数据，保证行数匹配
     plt_hir_filtered <- plt_hir[rowSums(com) >= 2, ]
     

  2. 调整模拟数据生成逻辑：若为测试场景，可循环生成符合要求的矩阵：

     repeat {
       com <- matrix(rbinom(120, 1, 0.6), nrow=12)
       # 确保所有样地至少2个物种，所有物种至少出现在2个样地
       if (all(rowSums(com)>=2) && all(colSums(com)>=2)) break
     }
     

  3. 并行代码优化：原代码中clusterExport无需导出simulate（vegan加载后已存在），可简化为：

     clusterExport(cl, c("nm"))
     

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