附录2  用于演示和模拟的R语言代码
Appendix 2   R scripts used for demonstration and simulation

##本演示计算的数据来自vegan包中的“mite”数据集(Borcard & Legendre, 1994)
library(vegan)#用于加载数据集以及计算PCoA
library(CommEcol)#用于计算Chao-jaccard指数
source("ESS.R")#导入ESS()函数(见Zou & Axmacher, 2020)

##1. 导入与查看“mite”数据集
  data(mite)#导入“mite”数据集
  row.names(mite) <- as.character(1:70)#将数据集的行名命名为1–70, 作为70个样方
  
  #查看数据集样方中的个体数
  summary(apply(mite,1,sum))#最小的样方拥有8个个体, 平均值为140个个体
  
  #查看数据集中各个物种多度并作图
  Sum <- colSums(mite)
  par(mar=c(5,5,2,1),mgp=c(2.2,1,0),family="STXihei")
  Abu <- data.frame(Sum)
  Abu$Species <- rownames(Abu)
  Abu$Rank <- rank(Abu$Sum)
  Abu <- Abu[order(Abu$Sum,decreasing=T),]
  barplot(Abu$Sum,
          ylim=c(0,2500),
          ylab="个体数\n Number of individuals",
          xlab="物种排序\n Ranked Species")
  
##2. 基于CNESS以及Chao-Jaccard指数, 使用PCoA(等同于经典的MDS)进行排序
  
  #计算CNESS指数并且进行PCoA排序
  Mat <- data.frame()#创建数据框
  for (m in c(1,8,140)) #CNESS的参数m选择1,8,140 
  {
      mi <- ESS(mite,m=m)#计算CNESS指数, 该指数自动移除数据集中样本量小于m的样方
      MDSm <-cmdscale(mi,  eig=TRUE)#PCoA排序
      Cm <- as.data.frame(MDSm$points)
      Cm$Var <- as.factor(paste("CNESS"," ","m=",m,sep=''))
      r.m <- which(apply(mite,1,sum)<m)#移除数据集中样本量小于m的样方, 以匹配ESS结果中的样方名称
      if (length(r.m)==0)
      {M.sub <- mite}else
      {M.sub <- mite[-r.m,]} 
      Cm$Plot <- row.names(M.sub)
      Mat <- rbind(Mat,Cm)
  }
  
  #计算Chao-Jaccard指数并且进行PCoA排序
  Ja <- dis.chao(mite,index="jaccard")
  M_Ja <- cmdscale(Ja, eig=TRUE)#PCoA排序
  Sdf <- as.data.frame(M_Ja$points)
  Sdf$Var <- "Chao-Jaccard"
  Sdf$Plot <- row.names(Sdf)
  colnames(Mat) <- colnames(Sdf)
  Dist <- rbind(Mat,Sdf)#合并两个数据框
 
 
##3. 基于两个不同指数得到的两个轴分别进行线性回归, 检验CNESS指数与Chao-Jaccard之间二维化之后所表达信息的关系
    Output <- data.frame()
    for (m in seq(from=1,to=140,by=2))
    {
        mi <- ESS(mite,m=m)
        MDSm <- cmdscale(mi,  eig=TRUE)
        Cm <- as.data.frame(MDSm$points)
        r.m <- which(apply(mite,1,sum)<m)#
        if (length(r.m)==0) 
          { S.sub <- Sdf} 
            else{S.sub <- Sdf[-r.m,]}#使用-r.m移除J.sub数据中的相应的行
        n=70-length(r.m)#剩余的样方
      CorJa1 <- cor.test(S.sub[,1],Cm[,1])
      CorJa2 <- cor.test(S.sub[,2],Cm[,2])
      CorJa <- data.frame(r=c(CorJa1$estimate,CorJa2$estimate),
                          p=c(CorJa1$p.value,CorJa2$p.value),
                          Var=c("MDS1","MSD2"),
                          Index="Chao-Jaccad")
      CorJa$m <- m
      CorJa$n <- n
      Output <- rbind(Output,CorJa)
    }

##4. 通过随机移除个体数从100依次递增至2,000, 分别计算新的样本的CNESS与Chao-Jaccard的平均值, 以及与完整数据集的差异
    LS <- list()#建立一个列表用于储存数据
    Tst <- as.matrix(mite)#将数据框改成矩阵, 有利于加快运算速度
    for (j in 1:20)
    {
      for (i in 1:100)#每次随机移除100个体
      {
        Int <- which(Tst>0) #移除只针对正数
        Sp <- sample(Int,1) #随机移除一个个体
        Tst[Sp] <- Tst[Sp]-1
      }
      LS[[j]] <- Tst
    }

    #针对每个新的数据, 重新计算CNESS以及Chao-Jaccard值
    Df  <- data.frame()
    for (k in 1:length(LS))
    {
    LS.k <- LS[[k]]
    m1 <- mean(ESS(LS.k,m=1))
    m8 <- mean(ESS(LS.k,m=8))
    m140 <- mean(ESS(LS.k,m=140))
    JD <- mean(dis.chao(LS.k,index="jaccard"))
    Out <- data.frame(Sample=k,m1,m8,m140,JD)
    Df <- rbind(Out,Df)
    }

    #计算与完整样本之间的差异
    Base <- subset(Df,Sample==1)
    Df$cm1 <- (Df$m1/Base$m1-1)*100
    Df$cm8 <- (Df$m8/Base$m8-1)*100
    Df$cm140 <- (Df$m140/Base$m140-1)*100
    Df$cJD <- (Df$JD/Base$JD-1)*100
    Df$Sample <- Df$Sample*100

    
参考文献
Borcard D, Legendre P (1994) Environmental control and spatial structure in ecological communities: An example using oribatid mites (Acari, Oribatei). Environmental and Ecological Statistics, 1, 37–61.
Zou Y, Axmacher JC (2020) The Chord-Normalized Expected Species Shared (CNESS)-distance represents a superior measure of species turnover patterns. Methods in Ecology and Evolution, 11, 273–280.
邹怡 (2021) 样本量不一致时的β多样性计算. 生物多样性, 29, 790–797. http://www.biodiversity-science.net/CN/10.17520/biods.2021011




16