在R包EcoSimR中设置na.rm=TRUE以解决全零行报错问题

问题场景 #

使用R包EcoSimR的ra4算法逐月逐站点计算Pianka生态位重叠指标，当数据集存在整行全为0的物种记录时，循环运行报错：

Error: Error in quantile.default(model$Sim, 0.05) :
若na.rm为FALSE，则不允许存在缺失值和NaN值

原代码如下：

#install.packages("EcoSimR")
library(EcoSimR)

set.seed(111)
month <- rep(c("J","J","J","F"), each = 4)
site <- rep(c("1","2","3","1"), each = 4)
species <- rep(c("A","B","C","D"), rep = 4)
q1 <- rtruncnorm(n=14, a=0, b=10, mean=0, sd=1) # 修正原代码多余的括号
q2 <- rtruncnorm(n=16, a=0, b=10, mean=1, sd=1)
q3 <- rtruncnorm(n=16, a=0, b=10, mean=0, sd=1)
q4 <- rtruncnorm(n=16, a=0, b=10, mean=1, sd=1)
q5 <- rtruncnorm(n=16, a=0, b=10, mean=1, sd=1)

df <- data.frame(month, site, species,q1,q2,q3,q4,q5)
df[1,c(4:8)] <- 0
df[6,c(4:8)] <- 0
df[15,c(4:8)] <- 0

get_eco_sim_result <- function(spd, algo= "ra4", metric = "pianka", nReps=500) {
  model = niche_null_model(speciesData = spd,
                           algo = algo,metric =metric, nReps = nReps, suppressProg = TRUE
  )
  return(list(
    Obs = model$Obs,
    Sim = mean(model$Sim, na.rm = TRUE),
    lower_1tailp = quantile(model$Sim,0.05),
    SES = (model$Obs - mean(model$Sim))/sd(model$Sim)
  ))
}

output <- do.call(
  rbind, lapply(split(df, list(month,site), drop=T), \(d) {
    data.frame(get_eco_sim_result(d[,-c(1,2,3)], nReps=5000))
  })
)

output 

报错原因 #

当物种数据行全为0时，niche_null_model运行模拟过程中会生成NaN值。原代码仅在计算Sim的均值时开启了na.rm=TRUE，但quantile()和计算SES时的mean()、sd()都未设置该参数，触发报错。

解决方案 #

1. 修正函数中的统计计算参数 #

在quantile()、mean(model$Sim)和sd(model$Sim)中都添加na.rm=TRUE，确保忽略NaN值：

修改后的get_eco_sim_result函数：

get_eco_sim_result <- function(spd, algo= "ra4", metric = "pianka", nReps=500) {
  model = niche_null_model(speciesData = spd,
                           algo = algo,metric =metric, nReps = nReps, suppressProg = TRUE
  )
  # 统一处理模拟值的均值和标准差，开启na.rm
  sim_mean <- mean(model$Sim, na.rm = TRUE)
  sim_sd <- sd(model$Sim, na.rm = TRUE)
  
  return(list(
    Obs = model$Obs,
    Sim = sim_mean,
    lower_1tailp = quantile(model$Sim, 0.05, na.rm = TRUE), # 添加na.rm=TRUE
    SES = (model$Obs - sim_mean)/sim_sd
  ))
}

2. 从根源避免NaN：提前过滤全0物种行 #

全0行的物种本身没有有效生态位数据，可在传入niche_null_model前过滤掉这类行，彻底避免NaN生成，同时增加物种数校验（至少2个物种才能计算生态位重叠）：

修改循环中的调用逻辑：

output <- do.call(
  rbind, lapply(split(df, list(month,site), drop=T), \(d) {
    sp_data <- d[,-c(1,2,3)]
    # 过滤所有列都为0的行
    sp_data_filtered <- sp_data[rowSums(sp_data) != 0, ]
    # 物种数不足时返回NA标记
    if(nrow(sp_data_filtered) >=2){
      data.frame(get_eco_sim_result(sp_data_filtered, nReps=5000))
    } else {
      data.frame(Obs=NA, Sim=NA, lower_1tailp=NA, SES=NA)
    }
  })
)

完整修正后代码 #

#install.packages("EcoSimR")
library(EcoSimR)

set.seed(111)
month <- rep(c("J","J","J","F"), each = 4)
site <- rep(c("1","2","3","1"), each = 4)
species <- rep(c("A","B","C","D"), rep = 4)
q1 <- rtruncnorm(n=14, a=0, b=10, mean=0, sd=1)
q2 <- rtruncnorm(n=16, a=0, b=10, mean=1, sd=1)
q3 <- rtruncnorm(n=16, a=0, b=10, mean=0, sd=1)
q4 <- rtruncnorm(n=16, a=0, b=10, mean=1, sd=1)
q5 <- rtruncnorm(n=16, a=0, b=10, mean=1, sd=1)

df <- data.frame(month, site, species,q1,q2,q3,q4,q5)
df[1,c(4:8)] <- 0
df[6,c(4:8)] <- 0
df[15,c(4:8)] <- 0

get_eco_sim_result <- function(spd, algo= "ra4", metric = "pianka", nReps=500) {
  model = niche_null_model(speciesData = spd,
                           algo = algo,metric =metric, nReps = nReps, suppressProg = TRUE
  )
  sim_mean <- mean(model$Sim, na.rm = TRUE)
  sim_sd <- sd(model$Sim, na.rm = TRUE)
  
  return(list(
    Obs = model$Obs,
    Sim = sim_mean,
    lower_1tailp = quantile(model$Sim, 0.05, na.rm = TRUE),
    SES = (model$Obs - sim_mean)/sim_sd
  ))
}

output <- do.call(
  rbind, lapply(split(df, list(month,site), drop=T), \(d) {
    sp_data <- d[,-c(1,2,3)]
    sp_data_filtered <- sp_data[rowSums(sp_data) != 0, ]
    if(nrow(sp_data_filtered) >=2){
      data.frame(get_eco_sim_result(sp_data_filtered, nReps=5000))
    } else {
      data.frame(Obs=NA, Sim=NA, lower_1tailp=NA, SES=NA)
    }
  })
)

output 

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